diff --git a/src/dual.jl b/src/dual.jl index 31a7ce99..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) @@ -682,7 +688,14 @@ 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 + +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 diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl index 9ea3b0a7..69f9409c 100644 --- a/test/JacobianTest.jl +++ b/test/JacobianTest.jl @@ -235,6 +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] + + x0 = [1.0, 2.0]; + ev1(x) = eigen(Symmetric(x*x')).vectors[:,1] + @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) ≈ Calculus.finite_difference_jacobian(ev2, x0) + x0_static = SVector{2}(x0) + @test ForwardDiff.jacobian(ev1, x0_static) ≈ Calculus.finite_difference_jacobian(ev1, x0) end @testset "type stability" begin