Improve 2x2 eigen #694

tfgast · 2019-11-27T23:04:28Z

Using the fact that the eigenvectors must be orthonormal, v21 and v22, can be computed from v21 and v11. In addition to being more performant, this is actually more robust in some cases, for instance with [1.0 -9.465257061381386e-20; -9.465257061381386e-20 1.0], the previous code would return duplicated eigenvectors.

coveralls · 2019-11-27T23:47:38Z

Coverage decreased (-0.04%) to 81.857% when pulling 54afe61 on tfgast:master into 103e9d4 on JuliaArrays:master.

c42f · 2019-11-27T23:56:47Z

Looks great to me, thanks!

Could you add a test or tests which cover the cases which the previous code got wrong? Is this characterized by the case where the off diagonal elements are less than eps(eltype(A))?

Would this fix also allow us to remove the special case for where A is detected to be diagonal? Just looking at the original code, the diagonal test seems unlikely to be numerically reliable... and ultimately the root cause of this bug.

tfgast · 2019-12-12T00:22:49Z

I've improved the handling in the case of underflow as well. I wasn't able to remove the special case for diagonal matrices, as we still need to avoid dividing by zero. I've also added some more tests.

c42f

Ok, no worries I'm happy to merge this more or less as is. Though I'd like it if you could clarify which cases you are targeting with the underflow tests (so we know what they are meant to test in the future).

c42f · 2019-12-12T01:44:44Z

test/eigen.jl

-        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


Suggested change

# issue #523

# issues #523, #694

c42f · 2019-12-12T02:04:15Z

test/eigen.jl

+        @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


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...

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?

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.

c42f

This seems so close, I'd love to get it in the next release but there's one problems still to resolve - the use of zip in the tests seems confusing and might be wrong.

Also, how does this change as a whole impact performance?

Anyway it's great that you've managed to delete half the code from the original implementation. That's what I like to see!

c42f · 2020-04-17T09:41:30Z

test/eigen.jl

+    smallest_normal = floatmin(zero)
+    largest_subnormal = prevfloat(smallest_normal)
+    epsilon = eps(1.0)
+    one_p_epsilon = 1.0 + epsilon


Quibble: Better written nextfloat(1.0) — it's the same in this case, but if you'd written 1.0 - epsilon that's not equal to prevfloat(1.0) which is a bit of a gotcha given that the boundary in floating point discretization density lies on the powers of two.

c42f · 2020-04-17T09:43:29Z

test/eigen.jl

+    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)


I don't think this zip does what you want because i,j,k are always the same? Though I'm not sure what you want!

Perhaps you meant to use for (i, j, k) in Iterators.product(degenerate,degenerate,degenerate)?

You're totally right, that's what I meant

andyferris · 2020-05-11T00:25:58Z

Where is this one at?

Also, how does this change as a whole impact performance?

Do we have an answer to this one?

thchr · 2022-01-06T16:09:29Z

I tested the performance impact of this just now: as far as I can see, there is a speed-up for Hermitian{ComplexF64, ...} matrices but a slow-down for Hermitian{Float64, ...} matrices (not completely clear why to me):

Hermitian{ComplexF64, ...} matrices

Master:

julia> @benchmark StaticArrays._eig(Size{(2,2)}(), A, true, true) setup=begin
           a = @SMatrix rand(ComplexF64,2,2)
           A = Hermitian(a + a')
       end
BenchmarkTools.Trial: 10000 samples with 940 evaluations.
 Range (min … max):  101.915 ns … 324.681 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     102.766 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   105.722 ns ±   9.921 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▅█▁          ▁▃▁  ▁▁                                          ▁
  ███▆▅▄▆▆▆▅▅▅▆███▇█████▇▇▇▆▅▅▅▅▅▅▄▅▅▅▄▅▄▄▅▄▅▄▄▄▅▅▄▄▄▅▅▅▄▄▄▄▃▄▃ █
  102 ns        Histogram: log(frequency) by time        149 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

PR: (speed-up)

julia> @benchmark StaticArrays._eig(Size{(2,2)}(), A, true, true) setup=begin
           a = @SMatrix rand(ComplexF64,2,2)
           A = Hermitian(a + a')
       end
BenchmarkTools.Trial: 10000 samples with 979 evaluations.
 Range (min … max):  64.147 ns … 324.821 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     66.394 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   67.999 ns ±   7.887 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▄▆▆█▂            ▁▂▁▁  ▁                                     ▁
  █████▆▄▅▃▆▅▆▆▆▄▆██████████▇▇▇▇▆▆▅▅▅▅▅▅▄▅▅▅▄▆▄▅▅▄▅▁▄▅▄▅▃▅▅▄▄▄ █
  64.1 ns       Histogram: log(frequency) by time       104 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

Hermitian{Float64, ...} matrices

Master:

julia> @benchmark StaticArrays._eig(Size{(2,2)}(), A, true, true) setup=begin
           a = @SMatrix rand(Float64,2,2)
           A = Hermitian(a + a')
       end
BenchmarkTools.Trial: 10000 samples with 999 evaluations.
 Range (min … max):  8.909 ns … 69.670 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     9.009 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   9.230 ns ±  1.770 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▄                                                         ▁
  ██▆▄▄▄▄▁▃▄▅▄▃▃▃▄▁▁▁▄▅▇▇▆▅▄▅▄▄▃▄▄▅▄▄▄▁▃▃▁▃▃▁▁▁▃▁▁▁▁▁▃▁▄▃▃▄▅ █
  8.91 ns      Histogram: log(frequency) by time     18.9 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

PR: (slow-down)

julia> @benchmark StaticArrays._eig(Size{(2,2)}(), A, true, true) setup=begin
           a = @SMatrix rand(Float64,2,2)
           A = Hermitian(a + a')
       end
BenchmarkTools.Trial: 10000 samples with 994 evaluations.
 Range (min … max):  29.879 ns … 204.930 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     31.489 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   32.438 ns ±   5.469 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▁▆▃ █▃                                                       ▁
  ███████▆▆▅▄▄▄▅▁▃▄▄▅▆▅▆▅▅▄▁▄▅▇█▇█▇▇▆▇▇▆▆▇▅▇▅▇▇▅▅▅▆▆▆█▇▇▆▅▆▆▅▆ █
  29.9 ns       Histogram: log(frequency) by time      51.6 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

c42f added linear-algebra numerical-robustness numerical robustness and accuracy labels Nov 27, 2019

c42f added the bugfix label Nov 27, 2019

tfgast added 3 commits December 11, 2019 16:24

Improve 2x2 eigen

01271c9

improve handling of underflow for 2x2 eigen

e6d9f01

add more degenerate tests for 2x2 eigen

54afe61

tfgast force-pushed the master branch from de286f6 to 54afe61 Compare December 12, 2019 00:27

c42f reviewed Dec 12, 2019

View reviewed changes

improve 2x2 eigen robustness

bb4a8f3

c42f added this to the v0.12.2 milestone Apr 14, 2020

c42f requested changes Apr 17, 2020

View reviewed changes

Improve eigen 2x2 test

d8a91db

fredrikekre mentioned this pull request Jan 6, 2022

Incorrect eigen vectors on symmetric tensors Ferrite-FEM/Tensors.jl#173

Closed

mschauer mentioned this pull request Jan 13, 2022

Fix for eigen not working properly for 2x2 Hermitian SMatrices #957

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Improve 2x2 eigen #694

Improve 2x2 eigen #694

tfgast commented Nov 27, 2019

coveralls commented Nov 27, 2019 •

edited

c42f commented Nov 27, 2019 •

edited

tfgast commented Dec 12, 2019

c42f left a comment

c42f Dec 12, 2019

c42f Dec 12, 2019

c42f Mar 17, 2020

tfgast Mar 20, 2020

c42f left a comment •

edited

c42f Apr 17, 2020

tfgast May 6, 2020

c42f Apr 17, 2020

tfgast May 6, 2020

andyferris commented May 11, 2020

thchr commented Jan 6, 2022 •

edited

Improve 2x2 eigen #694

Are you sure you want to change the base?

Improve 2x2 eigen #694

Conversation

tfgast commented Nov 27, 2019

coveralls commented Nov 27, 2019 • edited

c42f commented Nov 27, 2019 • edited

tfgast commented Dec 12, 2019

c42f left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

c42f left a comment • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

andyferris commented May 11, 2020

thchr commented Jan 6, 2022 • edited

Hermitian{ComplexF64, ...} matrices

Hermitian{Float64, ...} matrices

coveralls commented Nov 27, 2019 •

edited

c42f commented Nov 27, 2019 •

edited

c42f left a comment •

edited

thchr commented Jan 6, 2022 •

edited