Eigen for AbstractFill and eigvals for Tridiagonal{<:Number,<:AbstractFillVector}

Project.toml

@@ -1,6 +1,6 @@
 name = "FillArrays"
 uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
-version = "1.1.1"
+version = "1.2.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/FillArrays.jl

@@ -10,7 +10,8 @@ import Base: size, getindex, setindex!, IndexStyle, checkbounds, convert,
 
 import LinearAlgebra: rank, svdvals!, tril, triu, tril!, triu!, diag, transpose, adjoint, fill!,
     dot, norm2, norm1, normInf, normMinusInf, normp, lmul!, rmul!, diagzero, AdjointAbsVec, TransposeAbsVec,
-    issymmetric, ishermitian, AdjOrTransAbsVec, checksquare, mul!
+    issymmetric, ishermitian, AdjOrTransAbsVec, checksquare, mul!, eigvals, eigen, eigvecs, Eigen,
+    HermOrSym
 
 
 import Base.Broadcast: broadcasted, DefaultArrayStyle, broadcast_shape

src/fillalgebra.jl

@@ -407,3 +407,194 @@ fillzero(::Type{Zeros{T,N,AXIS}}, n, m) where {T,N,AXIS} = Zeros{T,N,AXIS}((n, m
 fillzero(::Type{F}, n, m) where F = throw(ArgumentError("Cannot create a zero array of type $F"))
 
 diagzero(D::Diagonal{F}, i, j) where F<:AbstractFill = fillzero(F, axes(D.diag[i], 1), axes(D.diag[j], 2))
+
+_eigenind(λ0, n, sortby) = (isnothing(sortby) || sortby(λ0) <= sortby(zero(λ0))) ? 1 : n
+
+function eigvals(A::AbstractFillMatrix{<:Number}; sortby = nothing)
+    Base.require_one_based_indexing(A)
+    n = checksquare(A)
+    # only one non-trivial eigenvalue for a rank-1 matrix
+    λ0 = float(getindex_value(A)) * n
+    ind = _eigenind(λ0, n, sortby)
+    OneElement(λ0, ind, n)
+end
+
+function eigvecs(A::AbstractFillMatrix{<:Number}; sortby = nothing)
+    Base.require_one_based_indexing(A)
+    n = checksquare(A)
+    M = similar(A, float(eltype(A)))
+    n == 0 && return M
+    val = getindex_value(A)
+    ind = _eigenind(val, n, sortby)
+    # The non-trivial eigenvector is normalize(ones(n))
+    M[:, ind] .= 1/sqrt(n)
+    # eigenvectors corresponding to zero eigenvalues
+    for (i, j) in enumerate(axes(M,2)[(ind == 1) .+ (1:end-1)])
+        # The eigenvectors are v = normalize([ones(n-1); -(n-1)]), and sum(v) == 0
+        # The ordering is arbitrary,
+        # and we choose to order in terms of the number of non-zero elements
+        nrm = 1/sqrt(i*(i+1))
+        M[1:i, j] .= nrm
+        M[i+1, j] = -i * nrm
+        M[i+2:end, j] .= zero(eltype(M))
+    end
+    return M
+end
+
+function eigen(A::AbstractFillMatrix{<:Number}; sortby = nothing)
+    Eigen(eigvals(A; sortby), eigvecs(A; sortby))
+end
+
+# Tridiagonal Toeplitz eigen following Trench (1985)
+# "On the eigenvalue problem for Toeplitz band matrices"
+# https://www.sciencedirect.com/science/article/pii/0024379585902770
+
+for MT in (:(Tridiagonal{<:Union{Real, Complex}, <:AbstractFillVector}),
+            :(SymTridiagonal{<:Union{Real, Complex}, <:AbstractFillVector}),
+            :(HermOrSym{T, <:Tridiagonal{T, <:AbstractFillVector{T}}} where {T<:Union{Real, Complex}})
+            )
+    @eval function eigvals(A::$MT; sortby = nothing)
+        _eigvals_toeplitz(A; sortby)
+    end
+end
+
+___eigvals_toeplitz(a, sqrtbc, n) = [a + 2 * sqrtbc * cospi(q/(n+1)) for q in n:-1:1]
+
+__eigvals_toeplitz(::AbstractMatrix, a, b, c, n) =
+    ___eigvals_toeplitz(a, √(b*c), n)
+__eigvals_toeplitz(::Union{SymTridiagonal, Symmetric{<:Any, <:Tridiagonal}}, a, b, c, n) =
+    ___eigvals_toeplitz(a, b, n)
+__eigvals_toeplitz(::Hermitian{<:Any, <:Tridiagonal}, a, b, c, n) =
+    ___eigvals_toeplitz(real(a), abs(b), n)
+
+# triangular Toeplitz
+function _eigvals_toeplitz(T; sortby = nothing)
+    Base.require_one_based_indexing(T)
+    n = checksquare(T)
+    # extra care to handle 0x0 and 1x1 matrices
+    # diagonal
+    a = get(T, (1,1), zero(eltype(T)))
+    # subdiagonal
+    b = get(T, (2,1), zero(eltype(T)))
+    # superdiagonal
+    c = get(T, (1,2), zero(eltype(T)))
+    vals = __eigvals_toeplitz(T, a, b, c, n)
+    !isnothing(sortby) && sort!(vals, by = sortby)
+    return vals
+end
+
+_eigvec_prefactor(A, cm1, c1, m) = sqrt(complex(cm1/c1))^m
+_eigvec_prefactor(A::Union{SymTridiagonal, Symmetric{<:Any, <:Tridiagonal}}, cm1, c1, m) = oneunit(_eigvec_eltype(A))
+
+function _eigvec_prefactors(A, cm1, c1)
+    x = _eigvec_prefactor(A, cm1, c1, 1)
+    [x^(j-1) for j in axes(A,1)]
+end
+_eigvec_prefactors(A::Union{SymTridiagonal, Symmetric{<:Any, <:Tridiagonal}}, cm1, c1) =
+    Fill(_eigvec_prefactor(A, cm1, c1, 1), size(A,1))
+
+_eigvec_eltype(A::SymTridiagonal) = float(eltype(A))
+_eigvec_eltype(A) = complex(float(eltype(A)))
+
+# The functions swapcols! and permutecols!! are copied from Julia Base, which is under the MIT License
+# https://github.com/JuliaLang/julia/blob/master/LICENSE.md
+#=
+MIT License
+Copyright (c) 2009-2022: Jeff Bezanson, Stefan Karpinski, Viral B. Shah,
+and other contributors: https://github.com/JuliaLang/julia/contributors
+Permission is hereby granted, free of charge, to any person obtaining a
+copy of this software and associated documentation files (the "Software"),
+to deal in the Software without restriction, including without limitation the
+rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is furnished to do so,
+subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all copies
+or substantial portions of the Software.
+=#
+# swap columns i and j of a, in-place
+function swapcols!(a::AbstractMatrix, i, j)
+    i == j && return
+    cols = axes(a,2)
+    @boundscheck i in cols || throw(BoundsError(a, (:,i)))
+    @boundscheck j in cols || throw(BoundsError(a, (:,j)))
+    for k in axes(a,1)
+        @inbounds a[k,i],a[k,j] = a[k,j],a[k,i]
+    end
+end
+
+function permutecols!!(A::AbstractMatrix, p::AbstractVector{<:Integer})
+    Base.require_one_based_indexing(A, p)
+    count = 0
+    start = 0
+    while count < length(p)
+        ptr = start = findnext(!iszero, p, start+1)
+        next = p[start]
+        count += 1
+        while next != start
+            swapcols!(A, ptr, next)
+            p[ptr] = 0
+            ptr = next
+            next = p[next]
+            count += 1
+        end
+        p[ptr] = 0
+    end
+    A
+end
+
+_normalizecols!(M, T) = foreach(normalize!, eachcol(M))
+function _normalizecols!(M, T::Union{SymTridiagonal, Symmetric{<:Number, <:Tridiagonal}})
+    n = size(M,1)
+    invnrm = sqrt(2/(n+1))
+    M .*= invnrm
+    return M
+end
+
+function _eigvecs_toeplitz(T; sortby = nothing)
+    Base.require_one_based_indexing(T)
+    n = checksquare(T)
+    M = Matrix{_eigvec_eltype(T)}(undef, n, n)
+    n == 0 && return M
+    n == 1 && return fill!(M, oneunit(eltype(M)))
+    cm1 = T[2,1] # subdiagonal
+    c1 = T[1,2]  # superdiagonal
+    prefactors = _eigvec_prefactors(T, cm1, c1)
+    for q in axes(M,2)
+        qrev = n+1-q # match the default eigenvalue sorting
+        x = qrev/(n+1)
+        cs = cispi(x)
+        cs1 = copy(cs)
+        for j in 1:cld(n,2)
+            M[j, q] = prefactors[j] * imag(cs)
+            cs *= cs1
+        end
+        phase = iseven(n+q) ? 1 : -1
+        for j in cld(n,2)+1:n
+            M[j, q] = phase * prefactors[2j-n] * M[n+1-j,q]
+        end
+    end
+    _normalizecols!(M, T)
+    if !isnothing(sortby)
+        perm = sortperm(eigvals(T), by = sortby)
+        permutecols!!(M, perm)
+    end
+    return M
+end
+
+for MT in (:(Tridiagonal{<:Union{Real,Complex}, <:AbstractFillVector}),
+            :(SymTridiagonal{<:Union{Real,Complex}, <:AbstractFillVector}),
+            :(HermOrSym{T, <:Tridiagonal{T, <:AbstractFillVector{T}}} where {T<:Union{Real,Complex}}),
+            )
+
+    @eval begin
+        function eigvecs(A::$MT; sortby = nothing)
+            _eigvecs_toeplitz(A; sortby)
+        end
+        function eigen(T::$MT; sortby = nothing)
+            Eigen(eigvals(T; sortby), eigvecs(T; sortby))
+        end
+    end
+end
+
+

test/runtests.jl

@@ -1451,6 +1451,88 @@ end
     end
 end
 
+@testset "eigen" begin
+    sortby = x -> (real(x), imag(x))
+    @testset "AbstractFill" begin
+        @testset for val in (2.0, -2, 3+2im, 4 - 5im), n in (0, 1, 4)
+            F = Fill(val, n, n)
+            M = Matrix(F)
+            @test eigvals(F; sortby) ≈ eigvals(M; sortby)
+            λ, V = eigen(F; sortby)
+            @test λ == eigvals(F; sortby)
+            @test V'V ≈ I
+            @test V' * F * V ≈ Diagonal(λ)
+        end
+        @testset for MT in (Ones, Zeros), T in (Float64, Int, ComplexF64), n in (0, 1, 4)
+            F = MT{T}(n,n)
+            M = Matrix(F)
+            @test eigvals(F; sortby) ≈ eigvals(M; sortby)
+            λ, V = eigen(F; sortby)
+            @test λ == eigvals(F; sortby)
+            @test V'V ≈ I
+            @test V' * F * V ≈ Diagonal(λ)
+        end
+    end
+    @testset "Tridiagonal Toeplitz" begin
+        @testset for n in (0, 1, 2, 6)
+            @testset "Tridiagonal" begin
+                for (dl, d, du) in (
+                    (Fill(2, max(0, n-1)), Fill(-4, n), Fill(3, max(0,n-1))),
+                    (Fill(2+3im, max(0, n-1)), Fill(-4+4im, n), Fill(3im, max(0,n-1)))
+                    )
+                    T = Tridiagonal(dl, d, du)
+                    @test eigvals(T; sortby) ≈ eigvals(Matrix(T); sortby)
+                    λ, V = eigen(T)
+                    @test T * V ≈ V * Diagonal(λ)
+                    λ, V = eigen(T, sortby=abs)
+                    @test T * V ≈ V * Diagonal(λ)
+                end
+            end
+
+            @testset "SymTridiagonal/Symmetric" begin
+                dv = Fill(1, n)
+                ev = Fill(3, max(0,n-1))
+                for ST in (SymTridiagonal(dv, ev), Symmetric(Tridiagonal(ev, dv, ev)))
+                    evST = eigvals(ST; sortby)
+                    @test evST ≈ eigvals(Matrix(ST); sortby)
+                    @test eltype(evST) <: Real
+                    λ, V = eigen(ST)
+                    @test V'V ≈ I
+                    @test V' * ST * V ≈ Diagonal(λ)
+                    λ, V = eigen(ST, sortby=abs)
+                    @test V'V ≈ I
+                    @test V' * ST * V ≈ Diagonal(λ)
+                end
+                dv = Fill(-4+4im, n)
+                ev = Fill(2+3im, max(0,n-1))
+                for ST2 in (SymTridiagonal(dv, ev), Symmetric(Tridiagonal(ev, dv, ev)))
+                    @test eigvals(ST2; sortby) ≈ eigvals(Matrix(ST2); sortby)
+                    λ, V = eigen(ST2)
+                    @test ST2 * V ≈ V * Diagonal(λ)
+                    λ, V = eigen(ST2, sortby=abs)
+                    @test ST2 * V ≈ V * Diagonal(λ)
+                end
+            end
+
+            @testset "Hermitian Tridiagonal" begin
+                for (dv, ev) in ((Fill(2+0im, n), Fill(3-4im, max(0, n-1))),
+                                    (Fill(2, n), Fill(3, max(0, n-1))))
+                    HT = Hermitian(Tridiagonal(ev, dv, ev))
+                    evHT = eigvals(HT; sortby)
+                    @test evHT ≈ eigvals(Matrix(HT); sortby)
+                    @test eltype(evHT) <: Real
+                    λ, V = eigen(HT)
+                    @test V'V ≈ I
+                    @test V' * HT * V ≈ Diagonal(λ)
+                    λ, V = eigen(HT, sortby=abs)
+                    @test V'V ≈ I
+                    @test V' * HT * V ≈ Diagonal(λ)
+                end
+            end
+        end
+    end
+end
+
 @testset "count" begin
     @test count(Ones{Bool}(10)) == count(Fill(true,10)) == 10
     @test count(Zeros{Bool}(10)) == count(Fill(false,10)) == 0