ForwardDiff.jl
ForwardDiff.jl copied to clipboard
JuliaDiff
Fix implementations of `eigen` and `eigvals`
The PR adds missing definitions of eigvals and eigen for Hermitian{<:Dual} and of eigen for Hermitian{<:Complex{<:Dual}}.
Additionally, the PR adds more tests of eigen and eigvals, unifies the existing implementations, and removes duplicate calculations in the definitions of eigen.
Fixes #756 without GenericLinearAlgebra.
For the example in #756, I get on [email protected] without GenericLinearAlgebra:
julia> test_hessian(Float64)
ERROR: MethodError: no method matching eigvals!(::Hermitian{ForwardDiff.Dual{…}, Matrix{…}}; sortby::Nothing)
The function `eigvals!` exists, but no method is defined for this combination of argument types.
...
julia> test_hessian(ComplexF64)
ERROR: MethodError: no method matching eigen!(::Hermitian{Complex{ForwardDiff.Dual{…}}, Matrix{Complex{…}}}; sortby::Nothing)
The function `eigen!` exists, but no method is defined for this combination of argument types.
...
And on [email protected] with import GenericLinearAlgebra: eigen:
julia> @time test_hessian(Float64)
0.000071 seconds (38 allocations: 33.750 KiB)
9×9 Matrix{Float64}:
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 4.44089e-16 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 4.44089e-16 0.0 0.0 2.66454e-15 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
julia> @time test_hessian(ComplexF64)
0.000556 seconds (1.97 k allocations: 1.387 MiB)
18×18 Matrix{Float64}:
0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 -6.93889e-18 0.0 0.0 0.0 0.0 -1.38778e-16 0.0
0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 4.33681e-19 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 1.38778e-17 0.0 0.0 0.0 0.0 2.77556e-17 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
On this PR I get without GenericLinearAlgebra:
julia> @time test_hessian(Float64)
0.000071 seconds (215 allocations: 55.547 KiB)
9×9 Matrix{Float64}:
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.77556e-17 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 2.77556e-17 0.0 0.0 0.0 -2.77556e-17
0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.77556e-17 0.0
julia> @time test_hessian(ComplexF64)
0.000491 seconds (1.00 k allocations: 477.500 KiB)
18×18 Matrix{Float64}:
0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 6.93889e-18 0.0 0.0 0.0 0.0 0.0 -1.38778e-17 0.0
0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 -6.93889e-18 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.77556e-17 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
Some benchmarks against master:
master
julia> using ForwardDiff, Chairmarks, LinearAlgebra
## Symmetric{<:Real}
julia> @b rand(10) ForwardDiff.gradient(x -> sum(eigvals(Symmetric(x * x'))), _)
11.312 μs (110 allocs: 49.297 KiB)
julia> @b rand(10) ForwardDiff.gradient(x -> sum(eigen(Symmetric(x * x')).values), _)
29.083 μs (273 allocs: 115.469 KiB)
## Hermitian{<:Complex}
julia> @b rand(10) ForwardDiff.gradient(x -> sum(eigvals(Hermitian(complex.(x * x', x' * x)))), _)
71.625 μs (4136 allocs: 542.594 KiB)
julia> @b rand(10) ForwardDiff.gradient(x -> sum(eigen(Hermitian(complex.(x * x', x' * x))).values), _)
ERROR: MethodError: no method matching eigen!(::Hermitian{Complex{ForwardDiff.Dual{…}}, Matrix{Complex{…}}}; sortby::Nothing)
...
## SymTridiagonal{<:Real}
julia> @b rand(10) ForwardDiff.gradient(x -> sum(eigvals(SymTridiagonal(x, x[begin:(end-1)]))), _)
14.500 μs (128 allocs: 37.625 KiB)
julia> @b rand(10) ForwardDiff.gradient(x -> sum(eigen(SymTridiagonal(x, x[begin:(end-1)])).values), _)
32.625 μs (314 allocs: 101.422 KiB)
## Hermitian{<:Real}
julia> @b rand(10) ForwardDiff.gradient(x -> sum(eigvals(Hermitian(x * x'))), _)
ERROR: MethodError: no method matching eigvals!(::Hermitian{ForwardDiff.Dual{ForwardDiff.Tag{…}, Float64, 10}, Matrix{ForwardDiff.Dual{…}}}; sortby::Nothing)
The function `eigvals!` exists, but no method is defined for this combination of argument types.
...
julia> @b rand(10) ForwardDiff.gradient(x -> sum(eigen(Hermitian(x * x')).values), _)
ERROR: MethodError: no method matching eigen!(::Hermitian{ForwardDiff.Dual{ForwardDiff.Tag{…}, Float64, 10}, Matrix{ForwardDiff.Dual{…}}}; sortby::Nothing)
The function `eigen!` exists, but no method is defined for this combination of argument types.
...
This PR
julia> using ForwardDiff, Chairmarks, LinearAlgebra
## Symmetric{<:Real}
julia> @b rand(10) ForwardDiff.gradient(x -> sum(eigvals(Symmetric(x * x'))), _)
15.834 μs (110 allocs: 49.297 KiB)
julia> @b rand(10) ForwardDiff.gradient(x -> sum(eigen(Symmetric(x * x')).values), _)
19.500 μs (133 allocs: 67.594 KiB)
## Hermitian{<:Complex}
julia> @b rand(10) ForwardDiff.gradient(x -> sum(eigvals(Hermitian(complex.(x * x', x' * x)))), _)
32.625 μs (136 allocs: 98.062 KiB)
julia> @b rand(10) ForwardDiff.gradient(x -> sum(eigen(Hermitian(complex.(x * x', x' * x))).values), _)
39.875 μs (159 allocs: 132.047 KiB)
## Symtridiagonal{<:Real}
julia> @b rand(10) ForwardDiff.gradient(x -> sum(eigvals(SymTridiagonal(x, x[begin:(end-1)]))), _)
10.938 μs (148 allocs: 31.062 KiB)
julia> @b rand(10) ForwardDiff.gradient(x -> sum(eigen(SymTridiagonal(x, x[begin:(end-1)])).values), _)
14.125 μs (171 allocs: 49.359 KiB)
## Hermitian{<:Real}
julia> @b rand(10) ForwardDiff.gradient(x -> sum(eigvals(Hermitian(x * x'))), _)
15.917 μs (110 allocs: 49.297 KiB)
julia> @b rand(10) ForwardDiff.gradient(x -> sum(eigen(Hermitian(x * x')).values), _)
19.250 μs (133 allocs: 67.594 KiB)
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 Hermitian, that is, if we change the function to
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)))
end
It doesn't matter to me, I never need a non-Hermitian matrix, I'm remarking just in case you wanted to know.
The PR fixes #780 as well: https://github.com/JuliaDiff/ForwardDiff.jl/issues/780#issuecomment-3373892904