diff --git a/src/dual.jl b/src/dual.jl index 7552dd5f..d8b7698c 100644 --- a/src/dual.jl +++ b/src/dual.jl @@ -843,6 +843,38 @@ function LinearAlgebra.eigen(A::SymTridiagonal{<:Dual{Tg,T,N}}) where {Tg,T<:Rea Eigen(λ,Dual{Tg}.(Q, tuple.(parts...))) end +# General eigvals # +function make_eigen_dual(val::Real, partial) + Dual{tagtype(partial)}(val, partial.partials) +end + +function make_eigen_dual(val::Complex, partial::Complex) + Complex(Dual{tagtype(real(partial))}(real(val), real(partial).partials), + Dual{tagtype(imag(partial))}(imag(val), imag(partial).partials)) +end + +function LinearAlgebra.eigen(A::StridedMatrix{<:Dual}) + A_values = map(d -> d.value, A) + A_values_eig = eigen(A_values) + UinvAU = A_values_eig.vectors \ A * A_values_eig.vectors + vals_diff = diag(UinvAU) + F = similar(A_values, eltype(A_values_eig.values)) + for i in axes(A_values, 1), j in axes(A_values, 2) + if i == j + F[i, j] = 0 + else + F[i, j] = inv(A_values_eig.values[j] - A_values_eig.values[i]) + end + end + vectors_diff = A_values_eig.vectors * (F .* UinvAU) + for i in eachindex(vectors_diff) + vectors_diff[i] = make_eigen_dual(A_values_eig.vectors[i], vectors_diff[i]) + end + Eigen(vals_diff, vectors_diff) +end + +LinearAlgebra.eigvals(A::StridedMatrix{<:Dual}) = eigen(A).values + # Functions in SpecialFunctions which return tuples # # Their derivatives are not defined in DiffRules # #---------------------------------------------------# diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl index 9cc5024c..83698069 100644 --- a/test/JacobianTest.jl +++ b/test/JacobianTest.jl @@ -259,6 +259,26 @@ end x0_mvector = MVector{2}(x0) @test ForwardDiff.jacobian(ev1, x0_mvector) isa MMatrix{2, 2} @test ForwardDiff.jacobian(ev1, x0_mvector) ≈ Calculus.finite_difference_jacobian(ev1, x0) + + # real eigenvalues + f(x) = eigvals(reshape(x, 2, 2)) + x1 = [1.0, 2.0, 3.0, 4.0] + @test ForwardDiff.jacobian(f, x1) ≈ Calculus.finite_difference_jacobian(f, x1) + + # complex eigenvalues + g(x) = begin + vals = eigvals(reshape(x, 2, 2)) + vcat(real(vals), imag(vals)) + end + x2 = [0.0, -1.0, 1.0, 0.0] + @test ForwardDiff.jacobian(g, x2) ≈ Calculus.finite_difference_jacobian(g, x2) + + h(x) = begin + v = eigen(reshape(x, 2, 2)).vectors[:,1] + v = v / norm(v) + end + x3 = [2.0, 1.0, 0.5, 3.0] + @test ForwardDiff.jacobian(h, x3) ≈ Calculus.finite_difference_jacobian(h, x3) end @testset "type stability" begin