From e11a6a70d0ccf829ccaead4e367ffffd642c85bc Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Sun, 13 Feb 2022 21:29:45 +0100 Subject: [PATCH 1/6] Make eigen(::Symmetric) work --- src/dual.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dual.jl b/src/dual.jl index 31a7ce99..30ae95a4 100644 --- a/src/dual.jl +++ b/src/dual.jl @@ -682,7 +682,7 @@ end function LinearAlgebra.eigen(A::Symmetric{<:Dual{Tg,T,N}}) where {Tg,T<:Real,N} λ = eigvals(A) - _,Q = eigen(SymTridiagonal(value.(parent(A).dv),value.(parent(A).ev))) + _,Q = eigen(Symmetric(value.(parent(A)))) parts = ntuple(j -> Q*_lyap_div!(Q' * getindex.(partials.(A), j) * Q - Diagonal(getindex.(partials.(λ), j)), value.(λ)), N) Eigen(λ,Dual{Tg}.(Q, tuple.(parts...))) end From becc2607a6c0278c6d35ebc426090fada95df6fa Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Sun, 13 Feb 2022 21:30:22 +0100 Subject: [PATCH 2/6] Add tests for eigen --- test/JacobianTest.jl | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl index 9ea3b0a7..647690f5 100644 --- a/test/JacobianTest.jl +++ b/test/JacobianTest.jl @@ -235,6 +235,24 @@ end @testset "eigen" begin @test ForwardDiff.jacobian(x -> eigvals(SymTridiagonal(x, x[1:end-1])), [1.,2.]) ≈ [(1 - 3/sqrt(5))/2 (1 - 1/sqrt(5))/2 ; (1 + 3/sqrt(5))/2 (1 + 1/sqrt(5))/2] @test ForwardDiff.jacobian(x -> eigvals(Symmetric(x*x')), [1.,2.]) ≈ [0 0; 2 4] + + function numdiff(f, x, Δ=1.e-5); + f0 = f(x); + dx = zero(x); + res = zeros(length(f0), length(x)); + for j=1:length(x) + fill!(dx, 0) + dx[j] = Δ + res[:,j] = (f(x+dx)-f(x-dx))/(2Δ) + end + return res + end + + x0 = [1.0, 2.0]; + ev1(x) = eigen(Symmetric(x*x')).vectors[:,1] + @test ForwardDiff.jacobian(ev1, x0) ≈ numdiff(ev1, x0) + ev2(x) = eigen(SymTridiagonal(x, x[1:end-1])).vectors[:,1] + @test ForwardDiff.jacobian(ev2, x0) ≈ numdiff(ev2, x0) end @testset "type stability" begin From 750f1fe8f77205f0ec7894e025cf8e6f76053e24 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Mon, 9 May 2022 21:13:44 +0200 Subject: [PATCH 3/6] Add support for autodiff of static symmetric tensors --- src/dual.jl | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/src/dual.jl b/src/dual.jl index 30ae95a4..f2c33bf7 100644 --- a/src/dual.jl +++ b/src/dual.jl @@ -658,6 +658,12 @@ function LinearAlgebra.eigvals(A::Symmetric{<:Dual{Tg,T,N}}) where {Tg,T<:Real,N Dual{Tg}.(λ, tuple.(parts...)) end +function LinearAlgebra.eigvals(A::Symmetric{<:Dual{Tg,T,N}, <:StaticArrays.StaticMatrix}) where {Tg,T<:Real,N} + λ,Q = eigen(Symmetric(value.(parent(A)))) + parts = ntuple(j -> diag(Q' * getindex.(partials.(A), j) * Q), N) + Dual{Tg}.(λ, tuple.(parts...)) +end + function LinearAlgebra.eigvals(A::Hermitian{<:Complex{<:Dual{Tg,T,N}}}) where {Tg,T<:Real,N} λ,Q = eigen(Hermitian(value.(real.(parent(A))) .+ im .* value.(imag.(parent(A))))) parts = ntuple(j -> diag(real.(Q' * (getindex.(partials.(real.(A)) .+ im .* partials.(imag.(A)), j)) * Q)), N) @@ -687,6 +693,13 @@ function LinearAlgebra.eigen(A::Symmetric{<:Dual{Tg,T,N}}) where {Tg,T<:Real,N} Eigen(λ,Dual{Tg}.(Q, tuple.(parts...))) end +function LinearAlgebra.eigen(A::Symmetric{<:Dual{Tg,T,N}, <:StaticArrays.StaticMatrix}) where {Tg,T<:Real,N} + λ = eigvals(A) + _,Q = eigen(Symmetric(value.(parent(A)))) + parts = ntuple(j -> Q*_lyap_div!(Q' * getindex.(partials.(A), j) * Q - Diagonal(getindex.(partials.(λ), j)), value.(λ)), N) + Eigen(λ,Dual{Tg}.(Q, tuple.(parts...))) +end + function LinearAlgebra.eigen(A::SymTridiagonal{<:Dual{Tg,T,N}}) where {Tg,T<:Real,N} λ = eigvals(A) _,Q = eigen(SymTridiagonal(value.(parent(A)))) From 1e8eaafd4ec994c99d3fe7e6c6e1bf8c05f52e8f Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Mon, 9 May 2022 21:20:22 +0200 Subject: [PATCH 4/6] Add tests for sym static eigen diff --- test/JacobianTest.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl index 647690f5..b590f124 100644 --- a/test/JacobianTest.jl +++ b/test/JacobianTest.jl @@ -253,6 +253,8 @@ end @test ForwardDiff.jacobian(ev1, x0) ≈ numdiff(ev1, x0) ev2(x) = eigen(SymTridiagonal(x, x[1:end-1])).vectors[:,1] @test ForwardDiff.jacobian(ev2, x0) ≈ numdiff(ev2, x0) + x0_static = SVector{2}(x0) + @test ForwardDiff.jacobian(ev1, x0_static) ≈ numdif(ev1, x0) end @testset "type stability" begin From 85c79a00933bf4d38d289889f3282947e3af76cd Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Mon, 9 May 2022 21:28:51 +0200 Subject: [PATCH 5/6] Fix typo in test --- test/JacobianTest.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl index b590f124..75f0285e 100644 --- a/test/JacobianTest.jl +++ b/test/JacobianTest.jl @@ -254,7 +254,7 @@ end ev2(x) = eigen(SymTridiagonal(x, x[1:end-1])).vectors[:,1] @test ForwardDiff.jacobian(ev2, x0) ≈ numdiff(ev2, x0) x0_static = SVector{2}(x0) - @test ForwardDiff.jacobian(ev1, x0_static) ≈ numdif(ev1, x0) + @test ForwardDiff.jacobian(ev1, x0_static) ≈ numdiff(ev1, x0) end @testset "type stability" begin From 40e25a2a4db37897b1d8659b12f21b446588d4cf Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Thu, 12 May 2022 17:05:59 +0200 Subject: [PATCH 6/6] Use Calculus.jl for test eigendiff --- test/JacobianTest.jl | 20 ++++---------------- 1 file changed, 4 insertions(+), 16 deletions(-) diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl index 75f0285e..69f9409c 100644 --- a/test/JacobianTest.jl +++ b/test/JacobianTest.jl @@ -235,26 +235,14 @@ end @testset "eigen" begin @test ForwardDiff.jacobian(x -> eigvals(SymTridiagonal(x, x[1:end-1])), [1.,2.]) ≈ [(1 - 3/sqrt(5))/2 (1 - 1/sqrt(5))/2 ; (1 + 3/sqrt(5))/2 (1 + 1/sqrt(5))/2] @test ForwardDiff.jacobian(x -> eigvals(Symmetric(x*x')), [1.,2.]) ≈ [0 0; 2 4] - - function numdiff(f, x, Δ=1.e-5); - f0 = f(x); - dx = zero(x); - res = zeros(length(f0), length(x)); - for j=1:length(x) - fill!(dx, 0) - dx[j] = Δ - res[:,j] = (f(x+dx)-f(x-dx))/(2Δ) - end - return res - end - + x0 = [1.0, 2.0]; ev1(x) = eigen(Symmetric(x*x')).vectors[:,1] - @test ForwardDiff.jacobian(ev1, x0) ≈ numdiff(ev1, x0) + @test ForwardDiff.jacobian(ev1, x0) ≈ Calculus.finite_difference_jacobian(ev1, x0) ev2(x) = eigen(SymTridiagonal(x, x[1:end-1])).vectors[:,1] - @test ForwardDiff.jacobian(ev2, x0) ≈ numdiff(ev2, x0) + @test ForwardDiff.jacobian(ev2, x0) ≈ Calculus.finite_difference_jacobian(ev2, x0) x0_static = SVector{2}(x0) - @test ForwardDiff.jacobian(ev1, x0_static) ≈ numdiff(ev1, x0) + @test ForwardDiff.jacobian(ev1, x0_static) ≈ Calculus.finite_difference_jacobian(ev1, x0) end @testset "type stability" begin