Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Improve 2x2 eigen #694

Open
wants to merge 5 commits into
base: master
Choose a base branch
from
Open

Improve 2x2 eigen #694

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
01271c9
Improve 2x2 eigen
tfgast Sep 13, 2019
e6d9f01
improve handling of underflow for 2x2 eigen
tfgast Dec 12, 2019
54afe61
add more degenerate tests for 2x2 eigen
tfgast Dec 12, 2019
bb4a8f3
improve 2x2 eigen robustness
tfgast Mar 20, 2020
d8a91db
Improve eigen 2x2 test
tfgast May 6, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
35 changes: 14 additions & 21 deletions src/eigen.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -157,56 +157,49 @@ end
TA = eltype(A)

@inbounds if A.uplo == 'U'
if !iszero(a[3]) # A is not diagonal
abs2a3 = abs2(a[3])
if !iszero(abs2a3) # A is not diagonal
t_half = real(a[1] + a[4]) / 2
d = real(a[1] * a[4] - a[3]' * a[3]) # Should be real
d = real(a[1] * a[4] - abs2a3) # Should be real

tmp2 = t_half * t_half - d
tmp = tmp2 < 0 ? zero(tmp2) : sqrt(tmp2) # Numerically stable for identity matrices, etc.
vals = SVector(t_half - tmp, t_half + tmp)

v11 = vals[1] - a[4]
n1 = sqrt(v11' * v11 + a[3]' * a[3])
n1 = sqrt(abs2(v11) + abs2a3) # always > 0
v11 = v11 / n1
v12 = a[3]' / n1

v21 = vals[2] - a[4]
n2 = sqrt(v21' * v21 + a[3]' * a[3])
v21 = v21 / n2
v22 = a[3]' / n2

vecs = @SMatrix [ v11 v21 ;
v12 v22 ]
vecs = @SMatrix [ v11 -v12' ;
v12 v11' ]

return Eigen(vals, vecs)
end
else # A.uplo == 'L'
if !iszero(a[2]) # A is not diagonal
abs2a2 = abs2(a[2])
if !iszero(abs2a2) # A is not diagonal
t_half = real(a[1] + a[4]) / 2
d = real(a[1] * a[4] - a[2]' * a[2]) # Should be real
d = real(a[1] * a[4] - abs2a2) # Should be real

tmp2 = t_half * t_half - d
tmp = tmp2 < 0 ? zero(tmp2) : sqrt(tmp2) # Numerically stable for identity matrices, etc.
vals = SVector(t_half - tmp, t_half + tmp)

v11 = vals[1] - a[4]
n1 = sqrt(v11' * v11 + a[2]' * a[2])
n1 = sqrt(abs2(v11) + abs2a2) # always > 0
v11 = v11 / n1
v12 = a[2] / n1

v21 = vals[2] - a[4]
n2 = sqrt(v21' * v21 + a[2]' * a[2])
v21 = v21 / n2
v22 = a[2] / n2

vecs = @SMatrix [ v11 v21 ;
v12 v22 ]
vecs = @SMatrix [ v11 -v12' ;
v12 v11' ]

return Eigen(vals,vecs)
end
end

# A must be diagonal if we reached this point; treatment of uplo 'L' and 'U' is then identical
# A must be diagonal (or computation of abs2 underflowed) if we reached this point;
# treatment of uplo 'L' and 'U' is then identical
A11 = real(a[1])
A22 = real(a[4])
if A11 < A22
Expand Down
42 changes: 27 additions & 15 deletions test/eigen.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -75,24 +75,36 @@ using StaticArrays, Test, LinearAlgebra
@test vals::SVector ≈ sort(m_d)
@test eigvals(m_c) ≈ sort(m_d)
@test eigvals(Hermitian(m_c)) ≈ sort(m_d)
end

# issue #523
for (i, j) in ((1, 2), (2, 1)), uplo in (:U, :L)
A = SMatrix{2,2,Float64}((i, 0, 0, j))
E = eigen(Symmetric(A, uplo))
@test eigvecs(E) * SDiagonal(eigvals(E)) * eigvecs(E)' ≈ A
end

m1_a = randn(2,2)
m1_a = m1_a*m1_a'
m1 = SMatrix{2,2}(m1_a)
m2_a = randn(2,2)
m2_a = m2_a*m2_a'
m2 = SMatrix{2,2}(m2_a)
@test (@inferred_maybe_allow SVector{2,ComplexF64} eigvals(m1, m2)) ≈ eigvals(m1_a, m2_a)
@test (@inferred_maybe_allow SVector{2,ComplexF64} eigvals(Symmetric(m1), Symmetric(m2))) ≈ eigvals(Symmetric(m1_a), Symmetric(m2_a))
# issue #523
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
# issue #523
# issues #523, #694

@testset "2×2 degenerate cases" for (i, j) in ((1 , 1), (1, 2), (2, 1)), uplo in (:U, :L)
fmin = floatmin(Float64)
pfmin = prevfloat(fmin)
nfmin = nextfloat(fmin)
A = SMatrix{2,2,Float64}((i, 0, 0, j))
E = eigen(Symmetric(A, uplo))
@test eigvecs(E) * SDiagonal(eigvals(E)) * eigvecs(E)' ≈ A
A = SMatrix{2,2,Float64}((i, pfmin, pfmin, j))
E = eigen(Symmetric(A, uplo))
@test eigvecs(E) * SDiagonal(eigvals(E)) * eigvecs(E)' ≈ A
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Cool, so can you add some comment about exactly which cases we are straddling here by examining the boundary between pfmin vs nfmin? We seem to have have x^2 == 0 for both of those and i ± x == 1.0, i*j ± x == i*j etc...

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd love to merge this, but I'm unsure whether these tests are testing the right thing. Do you have some thoughts on the above query?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I looked into it again. You were right, they were testing the wrong thing. When I tested the right thing it wasn't working properly. I put some time in to make it work and simplify things, mostly by using the hypot function which already has a lot of work put into making it accurate.

A = SMatrix{2,2,Float64}((i, fmin, fmin, j))
E = eigen(Symmetric(A, uplo))
@test eigvecs(E) * SDiagonal(eigvals(E)) * eigvecs(E)' ≈ A
A = SMatrix{2,2,Float64}((i, nfmin, nfmin, j))
E = eigen(Symmetric(A, uplo))
@test eigvecs(E) * SDiagonal(eigvals(E)) * eigvecs(E)' ≈ A
end

m1_a = randn(2,2)
m1_a = m1_a*m1_a'
m1 = SMatrix{2,2}(m1_a)
m2_a = randn(2,2)
m2_a = m2_a*m2_a'
m2 = SMatrix{2,2}(m2_a)
@test (@inferred_maybe_allow SVector{2,ComplexF64} eigvals(m1, m2)) ≈ eigvals(m1_a, m2_a)
@test (@inferred_maybe_allow SVector{2,ComplexF64} eigvals(Symmetric(m1), Symmetric(m2))) ≈ eigvals(Symmetric(m1_a), Symmetric(m2_a))

@test_throws DimensionMismatch eigvals(SA[1 2 3; 4 5 6], SA[1 2 3; 4 5 5])
@test_throws DimensionMismatch eigvals(SA[1 2; 4 5], SA[1 2 3; 4 5 5; 3 4 5])

Expand Down