Skip to content

Commit

Permalink
Disable some eigen methods (#111)
Browse files Browse the repository at this point in the history 
* Disable some eigen methods

* Test non-extended eigen on Julia v1.9

* Test internal methods

* Add symmetric eigen test with negative dv

* Don't use sortby to compare

* separate out bands in Tridiagonal loops
  • Loading branch information
@jishnub
jishnub committed Aug 1, 2023
1 parent c9bea8f commit b40021e
Show file tree
Hide file tree
Showing 2 changed files with 107 additions and 67 deletions.
49 changes: 25 additions & 24 deletions src/eigen.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,31 +6,43 @@
# 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}})
# The methods aren't defined for general Tridiagonal or Hermitian, as the
# ordering of eigenvectors needs fixing
for MT in (:(SymTridiagonal{<:Union{Real,Complex}, <:AbstractFillVector}),
:(Symmetric{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)
_eigvals(A)
end
end
end

___eigvals_toeplitz(a, sqrtbc, n) = [a + 2 * sqrtbc * cospi(q/(n+1)) for q in n:-1:1]
for MT in (:(SymTridiagonal{<:Union{Real,Complex}, <:AbstractFillVector}),
:(Symmetric{T, <:Tridiagonal{T, <:AbstractFillVector{T}}} where {T<:Union{Real,Complex}}),
)

@eval begin
eigvecs(A::$MT) = _eigvecs(A)
eigen(A::$MT) = _eigen(A)
end
end


___eigvals(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)
__eigvals(::AbstractMatrix, a, b, c, n) =
___eigvals(a, √(complex(b*c)), n)
__eigvals(::Union{SymTridiagonal, Symmetric{<:Any, <:Tridiagonal}}, a, b, c, n) =
___eigvals(a, b, n)
__eigvals(::Hermitian{<:Any, <:Tridiagonal}, a, b, c, n) =
___eigvals(real(a), abs(b), n)

# triangular Toeplitz
function _eigvals_toeplitz(T)
function _eigvals(T)
require_one_based_indexing(T)
n = checksquare(T)
# extra care to handle 0x0 and 1x1 matrices
Expand All @@ -40,7 +52,7 @@ function _eigvals_toeplitz(T)
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)
vals = __eigvals(T, a, b, c, n)
return vals
end

Expand Down Expand Up @@ -152,14 +164,3 @@ function _eigen(A)
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
125 changes: 82 additions & 43 deletions test/runtests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -585,59 +585,98 @@ end
end

@testset "eigen" begin
sortby = x -> (real(x), imag(x))
@testset "Tridiagonal Toeplitz" begin
_sizes = (1, 2, 5, 6, 10, 15)
sizes = VERSION >= v"1.6" ? (0, _sizes...) : _sizes
sizes = (1, 2, 5, 6, 10, 15)
@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(λ)
if VERSION >= v"1.9"
@testset "Tridiagonal" begin
evm1r = Fill(2, max(0, n-1))
ev1r = Fill(3, max(0,n-1))
dvr = Fill(-4, n)
evm1c = Fill(2+3im, max(0, n-1))
dvc = Fill(-4+4im, n)
ev1c = Fill(3im, max(0,n-1))

for (dl, d, du) in (
(evm1r, dvr, ev1r),
(evm1r, dvr, -ev1r),
(evm1c, dvc, ev1c),
)

T = Tridiagonal(dl, d, du)
λT = eigvals(T)
λTM = eigvals(Matrix(T))
@test all(x -> any(y -> y ≈ x, λTM), λT)
λ, V = eigen(T)
@test T * V ≈ V * Diagonal(λ)

# Test that internal methods are correct,
# aside from the ordering of eigenvectors
λT2 = ToeplitzMatrices._eigvals(T)
@test all(x -> any(y -> y ≈ x, λTM), λT2)
V2 = ToeplitzMatrices._eigvecs(T)
for v in eachcol(V2)
w = T * v
@test any(λ -> w ≈ λ * v, λT2)
end
end
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(λ)
_dv = Fill(1, n)
_ev = Fill(3, max(0,n-1))
for dv in (_dv, -_dv), ev in (_ev, -_ev)
for ST in (SymTridiagonal(dv, ev), Symmetric(Tridiagonal(ev, dv, ev)))
λST = eigvals(ST)
λSTM = eigvals(Matrix(ST))
@test all(x -> any(y -> y ≈ x, λSTM), λST)
@test eltype(λST) <: Real
λ, V = eigen(ST)
@test V'V ≈ I
@test V' * ST * V ≈ Diagonal(λ)
end
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(λ)
_dv = Fill(-4+4im, n)
_ev = Fill(2+3im, max(0,n-1))
for dv in (_dv, -_dv, conj(_dv)), ev in (_ev, -_ev, conj(_ev))
for ST2 in (SymTridiagonal(dv, ev), Symmetric(Tridiagonal(ev, dv, ev)))
λST = eigvals(ST2)
λSTM = eigvals(Matrix(ST2))
@test all(x -> any(y -> y ≈ x, λSTM), λST)
λ, V = eigen(ST2)
@test ST2 * V ≈ V * Diagonal(λ)
end
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(λ)
if VERSION >= v"1.9"
@testset "Hermitian Tridiagonal" begin
_dvR = Fill(2, n)
_evR = Fill(3, max(0, n-1))
_dvc = complex(_dvR)
_evc = Fill(3-4im, max(0, n-1))
for (dv, ev) in ((_dvc, _evc), (_dvc, conj(_evc)),
(_dvR, _evR), (_dvR, -_evR))
HT = Hermitian(Tridiagonal(ev, dv, ev))
λHT = eigvals(HT)
λHTM = eigvals(Matrix(HT))
@test all(x -> any(y -> y ≈ x, λHTM), λHT)
@test eltype(λHT) <: Real
λ, V = eigen(HT)
@test V'V ≈ I
@test V' * HT * V ≈ Diagonal(λ)

# Test that internal methods are correct,
# aside from the ordering of eigenvectors
λHT2 = ToeplitzMatrices._eigvals(HT)
@test all(x -> any(y -> y ≈ x, λHTM), λHT2)
V2 = ToeplitzMatrices._eigvecs(HT)
for v in eachcol(V2)
w = HT * v
@test any(λ -> w ≈ λ * v, λHT2)
end
end
end
end
end
Expand Down

0 comments on commit b40021e

Please sign in to comment.