Fix implementations of eigen and eigvals
#757
Conversation
|
Amazing. I can confirm it fixes not only the MWE I posted, but also the real code I extracted it from. Note that there is still a regression if the matrix is not wrapped in function barrier(point)
d2 = div(d, 2)
M = mat(point, T)
M[1:d2, d2+1:end] .= 0
M[d2+1:end, 1:d2] .= 0
return sum(real(eigvals(M)))
endIt doesn't matter to me, I never need a non-Hermitian matrix, I'm remarking just in case you wanted to know. |
| # analytic solutions | ||
| @test ForwardDiff.jacobian(eigvals_symreal, x) ≈ [0 0; 2 4] | ||
| @test ForwardDiff.jacobian(eigvals_hermreal, x) ≈ [0 0; 2 4] | ||
| if !(x isa StaticArrays.StaticArray) | ||
| @test ForwardDiff.jacobian(eigvals_symtridiag, x) ≈ [(1 - 3/sqrt(5))/2 (1 - 1/sqrt(5))/2 ; (1 + 3/sqrt(5))/2 (1 + 1/sqrt(5))/2] | ||
| end |
There was a problem hiding this comment.
I haven't looked very closely, but are there more tests somewhere that check the values?
In what's changed in this PR, I see one more line with Calculus.finite_difference_jacobian(ev, float.(x0)) where I think x0 = [1.0, 2.0].
There was a problem hiding this comment.
Yes, there are more in the lines below. Much more exhaustive than on master.
There was a problem hiding this comment.
Sorry that wasn't clear. I see the many tests of different paths against each other, which is good.
But what I meant is tests of any ForwardDiff path against something completely independent -- a known answer, or finite differences. To catch things like making the same logic error in writing both _eigvals_hermitian and _eigen_hermitian here.
There was a problem hiding this comment.
In the tests below, Jacobians of every test function with ForwardDiff are checked against finite differencing based Jacobians. The latter doesn't involve any ForwardDiff paths.
There was a problem hiding this comment.
Sorry to be dense but where are you looking? Besides line 293, as mentioned:
I see one more line with Calculus.finite_difference_jacobian(ev, float.(x0)) where I think x0 = [1.0, 2.0].
There was a problem hiding this comment.
It is line 293. Calculus doesn't support all vector types, hence the ForwardDiff results with x (the possibly static version of x0: lines 259 and 262) is compared with the finite differencing result using the regular array x0.
There was a problem hiding this comment.
Ok.
So the matrices being tested are these. Does it seem OK to test e.g. no negative eigenvalue for Symmetric{Real}, no zero eigenvalues for Hermitian{Complex}? Haven't thought hard but it just seemed a small sample.
julia> x
2-element Vector{Float64}:
1.0
2.0
julia> Symmetric(x*x') # also wrapped in Hermitian
2×2 Symmetric{Float64, Matrix{Float64}}:
1.0 2.0
2.0 4.0
julia> eigvals(ans)
2-element Vector{Float64}:
0.0
5.0
julia> Hermitian(complex.(x*x', x'*x))
2×2 Hermitian{ComplexF64, Matrix{ComplexF64}}:
1.0+0.0im 2.0+5.0im
2.0-5.0im 4.0+0.0im
julia> eigvals(ans)
2-element Vector{Float64}:
-3.0901699437494754
8.090169943749475There was a problem hiding this comment.
(That said, when I tried various tests locally, I did not manage to break this.)
The PR adds missing definitions of
eigvalsandeigenforHermitian{<:Dual}and ofeigenforHermitian{<:Complex{<:Dual}}.Additionally, the PR adds more tests of
eigenandeigvals, unifies the existing implementations, and removes duplicate calculations in the definitions ofeigen.Fixes #756 without GenericLinearAlgebra.
For the example in #756, I get on ForwardDiff@0.10 without GenericLinearAlgebra:
And on ForwardDiff@0.10 with
import GenericLinearAlgebra: eigen:On this PR I get without GenericLinearAlgebra:
Some benchmarks against master:
master
This PR