From 274044a5cee4ec58e151360a62c7f7fed89e8e0a Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Tue, 11 Jul 2023 19:13:08 +0530 Subject: [PATCH] `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 | 2 + src/ToeplitzMatrices.jl | 7 +++ src/eigen.jl | 122 ++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 64 ++++++++++++++++++++- 4 files changed, 193 insertions(+), 2 deletions(-) create mode 100644 src/eigen.jl diff --git a/Project.toml b/Project.toml index a51dc60..7a37216 100644 --- a/Project.toml +++ b/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" diff --git a/src/ToeplitzMatrices.jl b/src/ToeplitzMatrices.jl index d19ef27..e44287a 100644 --- a/src/ToeplitzMatrices.jl +++ b/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) diff --git a/src/eigen.jl b/src/eigen.jl new file mode 100644 index 0000000..14a36dd --- /dev/null +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index b87e9da..eb510c8 100644 --- a/test/runtests.jl +++ b/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