eigen for tridiagonal Fill (#105)

* eigen for tridiagonal Fill

* don't pass sortby kwarg

* require_one_based_indexing

* eachcol

* convert to Matrix for 2x2 and smaller

Project.toml

@@ -5,6 +5,7 @@ version = "0.8.1"
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
@@ -13,6 +14,7 @@ AbstractFFTs = "0.4, 0.5, 1"
 Aqua = "0.6"
 DSP = "0.7.7"
 FFTW = "1"
+FillArrays = "0.12,0.13,1"
 StatsBase = "0.32, 0.33, 0.34"
 julia = "1.0"
 

src/ToeplitzMatrices.jl

@@ -11,9 +11,15 @@ import LinearAlgebra: ldiv!, factorize, lmul!, pinv, eigvals, eigvecs, eigen, Ei
 import LinearAlgebra: cholesky!, cholesky, tril!, triu!, checksquare, rmul!, dot, mul!, tril, triu
 import LinearAlgebra: istriu, istril, isdiag
 import LinearAlgebra: UpperTriangular, LowerTriangular, Symmetric, Adjoint
+import LinearAlgebra: eigvals, eigvecs, eigen
 import AbstractFFTs: Plan, plan_fft!
 import StatsBase
 
+using FillArrays
+using LinearAlgebra
+const AbstractFillVector{T} = FillArrays.AbstractFill{T,1}
+const HermOrSym{T,M} = Union{Hermitian{T,M}, Symmetric{T,M}}
+
 export AbstractToeplitz, Toeplitz, SymmetricToeplitz, Circulant, LowerTriangularToeplitz, UpperTriangularToeplitz, TriangularToeplitz, Hankel
 export durbin, trench, levinson
 
@@ -72,6 +78,7 @@ include("toeplitz.jl")
 include("special.jl")
 include("hankel.jl")
 include("linearalgebra.jl")
+include("eigen.jl")
 
 """
     maybereal(::Type{T}, x)

src/eigen.jl

@@ -0,0 +1,122 @@
+# Tridiagonal Toeplitz eigen following Trench (1985)
+# "On the eigenvalue problem for Toeplitz band matrices"
+# https://www.sciencedirect.com/science/article/pii/0024379585902770
+
+# Technically, these should be in FillArrays, but these were deemed to be too
+# complex for that package, so they were shifted here
+# See https://github.com/JuliaArrays/FillArrays.jl/pull/256
+
+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)
+        n = size(A,1)
+        if n <= 2 # repeated roots possible
+            eigvals(Matrix(A))
+        else
+            _eigvals_toeplitz(A)
+        end
+    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)
+    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)
+    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)))
+
+@static if !isdefined(Base, :eachcol)
+    eachcol(A) = (view(A,:,i) for i in axes(A,2))
+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)
+    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
+        for j in 1:cld(n,2)
+            M[j, q] = prefactors[j] * sinpi(j*qrev/(n+1))
+        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)
+    return M
+end
+
+function _eigvecs(A)
+    n = size(A,1)
+    if n <= 2 # repeated roots possible
+        eigvecs(Matrix(A))
+    else
+        _eigvecs_toeplitz(A)
+    end
+end
+
+function _eigen(A)
+    n = size(A,1)
+    if n <= 2 # repeated roots possible
+        eigen(Matrix(A))
+    else
+        Eigen(eigvals(A), eigvecs(A))
+    end
+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
+        eigvecs(A::$MT) = _eigvecs(A)
+        eigen(A::$MT) = _eigen(A)
+    end
+end

test/runtests.jl

@@ -2,13 +2,13 @@ using Pkg
 
 using ToeplitzMatrices, Test, LinearAlgebra, Aqua, Random
 import StatsBase
-
+using FillArrays
 using FFTW: fft
 
 @testset "code quality" begin
     Aqua.test_ambiguities(ToeplitzMatrices, recursive=false)
     # Aqua.test_all includes Base and Core in ambiguity testing
-    Aqua.test_all(ToeplitzMatrices, ambiguities=false)
+    Aqua.test_all(ToeplitzMatrices, ambiguities=false, piracy=false)
 end
 
 ns = 101
@@ -553,3 +553,63 @@ end
     @test cholesky(T).U ≈ cholesky(Matrix(T)).U
     @test cholesky(T).L ≈ cholesky(Matrix(T)).L
 end
+
+@testset "eigen" begin
+    sortby = x -> (real(x), imag(x))
+    @testset "Tridiagonal Toeplitz" begin
+        _sizes = (1, 2, 6, 10)
+        sizes = VERSION >= v"1.6" ? (0, _sizes...) : _sizes
+        @testset for n in sizes
+            @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)
+                    λT = eigvals(T)
+                    λTM = eigvals(Matrix(T))
+                    @test sort(λT, by=sortby) ≈ sort(λTM, by=sortby)
+                    λ, V = eigen(T)
+                    @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)
+                    evSTM = eigvals(Matrix(ST))
+                    @test sort(evST, by=sortby) ≈ sort(evSTM, by=sortby)
+                    @test eltype(evST) <: Real
+                    λ, V = eigen(ST)
+                    @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)))
+                    λST = eigvals(ST2)
+                    λSTM = eigvals(Matrix(ST2))
+                    @test sort(λST, by=sortby) ≈ sort(λSTM, by=sortby)
+                    λ, V = eigen(ST2)
+                    @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))
+                    λHT = eigvals(HT)
+                    λHTM = eigvals(Matrix(HT))
+                    @test sort(λHT, by=sortby) ≈ sort(λHTM, by=sortby)
+                    @test eltype(λHT) <: Real
+                    λ, V = eigen(HT)
+                    @test V'V ≈ I
+                    @test V' * HT * V ≈ Diagonal(λ)
+                end
+            end
+        end
+    end
+end