From 01271c90bde572a454680aade0453bd5221e2b27 Mon Sep 17 00:00:00 2001 From: Theodore Gast Date: Fri, 13 Sep 2019 15:17:10 -0700 Subject: [PATCH 1/5] Improve 2x2 eigen --- src/eigen.jl | 18 ++++-------------- 1 file changed, 4 insertions(+), 14 deletions(-) diff --git a/src/eigen.jl b/src/eigen.jl index 7332459a..30386dc6 100644 --- a/src/eigen.jl +++ b/src/eigen.jl @@ -170,13 +170,8 @@ end 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 @@ -194,13 +189,8 @@ end 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 From e6d9f0158d7c7568d64183a7f94722fea989882f Mon Sep 17 00:00:00 2001 From: Theodore Gast Date: Wed, 11 Dec 2019 16:18:08 -0800 Subject: [PATCH 2/5] improve handling of underflow for 2x2 eigen --- src/eigen.jl | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/src/eigen.jl b/src/eigen.jl index 30386dc6..53c11083 100644 --- a/src/eigen.jl +++ b/src/eigen.jl @@ -157,16 +157,17 @@ 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 @@ -176,16 +177,17 @@ end 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 @@ -196,7 +198,8 @@ end 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 From 54afe616d04536d4606c03e71ca403e17b16a3b7 Mon Sep 17 00:00:00 2001 From: Theodore Gast Date: Wed, 11 Dec 2019 16:18:31 -0800 Subject: [PATCH 3/5] add more degenerate tests for 2x2 eigen --- test/eigen.jl | 42 +++++++++++++++++++++++++++--------------- 1 file changed, 27 insertions(+), 15 deletions(-) diff --git a/test/eigen.jl b/test/eigen.jl index fc8c98e8..1ad9bd50 100644 --- a/test/eigen.jl +++ b/test/eigen.jl @@ -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 + @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 + 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]) From bb4a8f395f0f80f8681843e7efe04f470d3ecd36 Mon Sep 17 00:00:00 2001 From: Theodore Gast Date: Thu, 19 Mar 2020 18:00:05 -0700 Subject: [PATCH 4/5] improve 2x2 eigen robustness --- src/eigen.jl | 66 +++++++++++++++++++++------------------------------ test/eigen.jl | 25 ++++++++----------- 2 files changed, 37 insertions(+), 54 deletions(-) diff --git a/src/eigen.jl b/src/eigen.jl index 53c11083..97bb1135 100644 --- a/src/eigen.jl +++ b/src/eigen.jl @@ -156,52 +156,40 @@ end a = A.data TA = eltype(A) - @inbounds if A.uplo == 'U' - 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] - 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(abs2(v11) + abs2a3) # always > 0 + A11 = real(a[1]) + A22 = real(a[4]) + if A.uplo == 'U' + @inbounds A21 = a[3]' + else + @inbounds A21 = a[2] + end + @inbounds if !iszero(A21) # A is not diagonal + t_half = (A11 + A22) / 2 + d_half = (A11 - A22) / 2 + + tmp = hypot(d_half, A21) + vals = SVector(t_half - tmp, t_half + tmp) + + v11 = d_half - tmp + n1 = hypot(v11, A21) # always > 0 + if n1 < floatmin(T) # n1 subnormal + scale = inv(floatmin(T)) + n1 *= scale + v11 = (v11 * scale) / n1 + v12 = (A21 * scale) / n1 + else v11 = v11 / n1 - v12 = a[3]' / n1 - - vecs = @SMatrix [ v11 -v12' ; - v12 v11' ] - - return Eigen(vals, vecs) + v12 = A21 / n1 end - else # A.uplo == 'L' - 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] - 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(abs2(v11) + abs2a2) # always > 0 - v11 = v11 / n1 - v12 = a[2] / n1 - vecs = @SMatrix [ v11 -v12' ; - v12 v11' ] + vecs = @SMatrix [ v11 -v12' ; + v12 v11' ] - return Eigen(vals,vecs) - end + return Eigen(vals, vecs) end - # A must be diagonal (or computation of abs2 underflowed) if we reached this point; + # A must be diagonal 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 vals = SVector(A11, A22) vecs = @SMatrix [convert(TA, 1) convert(TA, 0); diff --git a/test/eigen.jl b/test/eigen.jl index 1ad9bd50..0970afd8 100644 --- a/test/eigen.jl +++ b/test/eigen.jl @@ -77,21 +77,16 @@ using StaticArrays, Test, LinearAlgebra @test eigvals(Hermitian(m_c)) ≈ sort(m_d) end - # issue #523 - @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 - 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)) + # issue #523, #694 + zero = 0.0 + smallest_non_zero = nextfloat(zero) + smallest_normal = floatmin(zero) + largest_subnormal = prevfloat(smallest_normal) + epsilon = eps(1.0) + one_p_epsilon = 1.0 + epsilon + degenerate = (zero, -1, 1, smallest_non_zero, smallest_normal, largest_subnormal, epsilon, one_p_epsilon, -one_p_epsilon) + @testset "2×2 degenerate cases" for (i, j, k) in zip(degenerate,degenerate,degenerate), uplo in (:U, :L) + A = SMatrix{2,2,Float64}((i, k, k, j)) E = eigen(Symmetric(A, uplo)) @test eigvecs(E) * SDiagonal(eigvals(E)) * eigvecs(E)' ≈ A end From d8a91db4ec88804d956ddcffd9c0b26a6f3751de Mon Sep 17 00:00:00 2001 From: Theodore Gast Date: Wed, 6 May 2020 16:12:08 -0700 Subject: [PATCH 5/5] Improve eigen 2x2 test --- test/eigen.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/eigen.jl b/test/eigen.jl index 0970afd8..86b53767 100644 --- a/test/eigen.jl +++ b/test/eigen.jl @@ -83,9 +83,9 @@ using StaticArrays, Test, LinearAlgebra smallest_normal = floatmin(zero) largest_subnormal = prevfloat(smallest_normal) epsilon = eps(1.0) - one_p_epsilon = 1.0 + epsilon + one_p_epsilon = nextfloat(1.0) degenerate = (zero, -1, 1, smallest_non_zero, smallest_normal, largest_subnormal, epsilon, one_p_epsilon, -one_p_epsilon) - @testset "2×2 degenerate cases" for (i, j, k) in zip(degenerate,degenerate,degenerate), uplo in (:U, :L) + @testset "2×2 degenerate cases" for (i, j, k) in Iterators.product(degenerate,degenerate,degenerate), uplo in (:U, :L) A = SMatrix{2,2,Float64}((i, k, k, j)) E = eigen(Symmetric(A, uplo)) @test eigvecs(E) * SDiagonal(eigvals(E)) * eigvecs(E)' ≈ A