From a4aa531c0e692851478565c9faf7e96ace9a3794 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Wed, 3 Dec 2025 14:52:15 +0100
Subject: [PATCH 01/15] differentiate general eigen and eigvals

---
 src/dual.jl | 32 ++++++++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)

diff --git a/src/dual.jl b/src/dual.jl
index f7aeb1d7..84c93d36 100644
--- a/src/dual.jl
+++ b/src/dual.jl
@@ -850,6 +850,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    #
 #---------------------------------------------------#

From a5c4b0c489608d5331be64c36ef7b89ee5f65f34 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Wed, 11 Feb 2026 18:09:54 +0100
Subject: [PATCH 02/15] add some basic tests

---
 test/JacobianTest.jl | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl
index 9cc5024c..f3b01227 100644
--- a/test/JacobianTest.jl
+++ b/test/JacobianTest.jl
@@ -259,6 +259,19 @@ 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)
 end
 
 @testset "type stability" begin

From 55fa69c640921cde8752e1ecbf572e4423dab036 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Wed, 11 Feb 2026 18:17:15 +0100
Subject: [PATCH 03/15] also add test for eigenvector

---
 test/JacobianTest.jl | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl
index f3b01227..83698069 100644
--- a/test/JacobianTest.jl
+++ b/test/JacobianTest.jl
@@ -272,6 +272,13 @@ end
     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

From 77f70a1a38264bfbf05f984388da2e1d19f3f553 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Mon, 10 Aug 2026 22:17:31 +0200
Subject: [PATCH 04/15] Compute general eigen derivatives at the value level

`eigen(A::StridedMatrix{<:Dual})` built the derivative from
`A_values_eig.vectors \ A * A_values_eig.vectors`, which mixes the
eigenvectors of `value.(A)` with `A` itself. `ForwardDiff.hessian` seeds
`Dual{T,Dual{T,V,N},N}` with the same tag on both levels, so the scalar
`*` in the generic matmul treats the two as one level and applies the
product rule at the outer level. First derivatives were unaffected;
second and higher ones were silently wrong.

Compute the derivatives entirely from `value.(A)` and assemble the
`Dual`s only at the end, the way the `Symmetric` and `SymTridiagonal`
methods already do, so the two levels never meet.

Along the way:

- `eigvals` no longer computes the eigenvector derivatives and throws
  them away; it is a separate method rather than `eigen(A).values`.
- reuse `_lyap_div!!` instead of building `F` by hand, which drops an
  n^2 allocation and a row-major traversal of a column-major array.
- share one `lu` factorization across the partials.
- rename `make_eigen_dual` to `_make_eigen_dual`, carry the tag as a
  type parameter, and give the section its `#---#` header line.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 src/dual.jl          | 68 +++++++++++++++++++++++++-------------------
 test/HessianTest.jl  | 33 +++++++++++++++++++++
 test/JacobianTest.jl | 18 ++++++++++++
 3 files changed, 89 insertions(+), 30 deletions(-)

diff --git a/src/dual.jl b/src/dual.jl
index ac266f2c..6c166b75 100644
--- a/src/dual.jl
+++ b/src/dual.jl
@@ -844,38 +844,46 @@ 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)
+# General eigvals and eigen #
+#---------------------------#
+
+# Assemble a value and its `N` partials into a `Dual` of type `D = Dual{Tg}`. Eigenvalues
+# and eigenvectors of a real matrix can be complex, in which case real and imaginary part
+# each become a `Dual` of their own.
+_make_eigen_dual(D::Type, val::Real, parts::NTuple{N,Real}) where {N} = D(val, parts...)
+_make_eigen_dual(D::Type, val::Complex, parts::NTuple{N,Number}) where {N} =
+    Complex(D(real(val), real.(parts)...), D(imag(val), imag.(parts)...))
+
+# The derivatives are computed entirely from the values of `A`, i.e. one `Dual` level
+# below the entries of `A`, and are only assembled into `Dual`s at the very end. Mixing
+# the two levels in a single expression (e.g. multiplying the eigenvectors of `value.(A)`
+# by `A` itself) would make the same-tag product rule apply at the outer level and hence
+# give wrong results for nested `Dual`s, as in `ForwardDiff.hessian`.
+#
+# The formulas are the ones of https://people.maths.ox.ac.uk/gilesm/files/NA-08-01.pdf:
+# with `A = U * Diagonal(λ) * inv(U)` and `M = inv(U) * Ȧ * U`, we have `λ̇ = diag(M)`
+# and `U̇ = U * (F .* M)`, where `F[i,j] = inv(λ[j] - λ[i])` off the diagonal and zero on it.
+
+LinearAlgebra.eigvals(A::StridedMatrix{Dual{Tg,T,N}}) where {Tg,T<:Real,N} = _eigvals_general(A)
+function _eigvals_general(A::StridedMatrix{Dual{Tg,T,N}}) where {Tg,T<:Real,N}
+    λ, U = eigen(value.(A))
+    luU = lu(U)
+    parts = ntuple(j -> diag(luU \ (getindex.(partials.(A), j) * U)), N)
+    return map((val, p) -> _make_eigen_dual(Dual{Tg}, val, p), λ, tuple.(parts...))
+end
+
+LinearAlgebra.eigen(A::StridedMatrix{Dual{Tg,T,N}}) where {Tg,T<:Real,N} = _eigen_general(A)
+function _eigen_general(A::StridedMatrix{Dual{Tg,T,N}}) where {Tg,T<:Real,N}
+    λ, U = eigen(value.(A))
+    luU = lu(U)
+    M = ntuple(j -> luU \ (getindex.(partials.(A), j) * U), N)
+    λ_parts = map(diag, M)
+    U_parts = ntuple(j -> U * _lyap_div!!(M[j] - Diagonal(λ_parts[j]), λ), N)
+    λ_dual = map((val, p) -> _make_eigen_dual(Dual{Tg}, val, p), λ, tuple.(λ_parts...))
+    U_dual = map((val, p) -> _make_eigen_dual(Dual{Tg}, val, p), U, tuple.(U_parts...))
+    return Eigen(λ_dual, U_dual)
 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/HessianTest.jl b/test/HessianTest.jl
index 8be72ee5..27995586 100644
--- a/test/HessianTest.jl
+++ b/test/HessianTest.jl
@@ -163,6 +163,39 @@ end
     @test ForwardDiff.hessian(x->dot(x,H,x), zeros(3)) ≈ [2 6 10; 6 10 14; 10 14 18]
 end
 
+@testset "nested duals in general eigen" begin
+    # The eigenvalue derivatives have to be computed from `value.(A)` alone; mixing that
+    # level with `A` itself applies the product rule at the wrong level for nested `Dual`s
+    B(w) = [3.0+w[1] 1.0+w[2]; 0.4+w[2] 2.0-2*w[1]]
+    w = [0.11, -0.07]
+    # sum(eigvals(B(w))) == tr(B(w)) == 5 - w[1] is linear in `w`
+    @test ForwardDiff.hessian(w -> sum(eigvals(B(w))), w) ≈ zeros(2, 2) atol=1e-12
+    # sum(eigvals(B(w)) .^ 2) == tr(B(w)^2) is quadratic in `w`
+    @test ForwardDiff.hessian(w -> sum(eigvals(B(w)) .^ 2), w) ≈ [10 0; 0 4]
+
+    # complex eigenvalues: λ = w[1] ± im*(1 + w[2]), i.e. sum(abs2, λ) == 2*(w[1]^2 + (1 + w[2])^2)
+    C(w) = [w[1] -1.0-w[2]; 1.0+w[2] w[1]]
+    @test ForwardDiff.hessian(w -> sum(abs2, eigvals(C(w))), [0.3, 0.2]) ≈ [4 0; 0 4]
+
+    # https://github.com/JuliaDiff/ForwardDiff.jl/issues/111
+    S(w) = [w[1]^2 w[1]*w[2]*w[3]; w[1]*w[2]*w[3] w[2]^2]
+    g(w) = sum(log, eigvals(S(w)))
+    gsym(w) = sum(log, eigvals(Symmetric(S(w))))
+    w111 = [0.9, 1.4, 0.3]
+    @test ForwardDiff.hessian(g, w111) ≈ ForwardDiff.hessian(gsym, w111)
+
+    # eigenvectors: renormalising removes the difference between the eigenvector
+    # convention of `eigen` and the one the derivatives are derived for
+    A0 = [2.0 1.0 0.5; 0.5 3.0 1.5; 0.25 0.75 4.0]
+    function v1(x)
+        v = eigen(reshape(x, 3, 3)).vectors[:, 1]
+        return sum(abs2, v / norm(v) .- [1.0, 0.5, -0.2])
+    end
+    x0 = vec(A0)
+    @test ForwardDiff.hessian(v1, x0) ≈ ForwardDiff.jacobian(x -> ForwardDiff.gradient(v1, x), x0)
+    @test ForwardDiff.hessian(v1, x0) ≈ Calculus.finite_difference_jacobian(x -> ForwardDiff.gradient(v1, x), x0) atol=1e-5
+end
+
 #https://github.com/JuliaDiff/ForwardDiff.jl/issues/720
 @testset "allocation-free hessian with StaticArrays" begin
     function hessian_allocs()
diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl
index f35a169f..e72c0042 100644
--- a/test/JacobianTest.jl
+++ b/test/JacobianTest.jl
@@ -279,6 +279,24 @@ end
     end
     x3 = [2.0, 1.0, 0.5, 3.0]
     @test ForwardDiff.jacobian(h, x3) ≈ Calculus.finite_difference_jacobian(h, x3)
+
+    # larger than 2x2, non-symmetric with real eigenvalues
+    f3(x) = eigvals(reshape(x, 3, 3))
+    x4 = vec([2.0 1.0 0.5; 0.5 3.0 1.5; 0.25 0.75 4.0])
+    @test ForwardDiff.jacobian(f3, x4) ≈ Calculus.finite_difference_jacobian(f3, x4)
+    h3(x) = begin
+        v = eigen(reshape(x, 3, 3)).vectors[:,2]
+        v = v / norm(v)
+    end
+    @test ForwardDiff.jacobian(h3, x4) ≈ Calculus.finite_difference_jacobian(h3, x4)
+
+    # eltypes of the general path
+    A_dual = Dual{TestTag}.([1.0 2.0; 3.0 4.0], [1.0 0.0; 0.0 0.0])
+    @test eigvals(A_dual) isa Vector{Dual{TestTag,Float64,1}}
+    @test eigen(A_dual).vectors isa Matrix{Dual{TestTag,Float64,1}}
+    A_dual_complex = Dual{TestTag}.([0.0 -1.0; 1.0 0.0], [1.0 0.0; 0.0 0.0])
+    @test eigvals(A_dual_complex) isa Vector{Complex{Dual{TestTag,Float64,1}}}
+    @test eigen(A_dual_complex).vectors isa Matrix{Complex{Dual{TestTag,Float64,1}}}
 end
 
 @testset "type stability" begin

From 6dd5491a8796c0eca5f0ea602c461bdb6061bc60 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Mon, 10 Aug 2026 22:41:30 +0200
Subject: [PATCH 05/15] Match `eigen`'s eigenvector normalization in the
 general case
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

`U̇ = U * (F .* M)` leaves the diagonal of `F` free; setting it to zero
normalizes the eigenvectors by `diag(inv(U) * U̇) == 0`. LAPACK's `geevx`
instead returns eigenvectors of unit 2-norm whose largest entry is real.
For a non-normal matrix those are different conventions, so `.vectors`
carried LAPACK values with partials belonging to a different
parameterization -- e.g. `dot(v, v̇)` came out as 0.5 rather than 0.

Pick the free diagonal such that the derivatives satisfy the two
constraints LAPACK's convention imposes, differentiated at `t = 0`.
This mirrors `_eigen_norm_phase_fwd!` in ChainRules, so the two systems
agree on `eigen`.

The eigenvector tests no longer renormalize; `v / norm(v)` was a no-op on
the values but removed exactly the scalar the two conventions differ by,
which is why they passed before.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 src/dual.jl          | 43 ++++++++++++++++++++++++++++++++++++++++++-
 test/HessianTest.jl  |  5 ++---
 test/JacobianTest.jl | 38 ++++++++++++++++++++++++++++++--------
 3 files changed, 74 insertions(+), 12 deletions(-)

diff --git a/src/dual.jl b/src/dual.jl
index 6c166b75..33397faf 100644
--- a/src/dual.jl
+++ b/src/dual.jl
@@ -864,6 +864,47 @@ _make_eigen_dual(D::Type, val::Complex, parts::NTuple{N,Number}) where {N} =
 # with `A = U * Diagonal(λ) * inv(U)` and `M = inv(U) * Ȧ * U`, we have `λ̇ = diag(M)`
 # and `U̇ = U * (F .* M)`, where `F[i,j] = inv(λ[j] - λ[i])` off the diagonal and zero on it.
 
+# Index of the entry of largest magnitude among the real entries of `v`. LAPACK normalizes
+# the eigenvectors it returns to unit 2-norm with their largest entry real, so this is the
+# entry that carries the phase convention. Modelled after `_findrealmaxabs2` in ChainRules.
+function _findrealmaxabs2(v)
+    imax = firstindex(v)
+    amax = abs2(zero(eltype(v)))
+    for i in eachindex(v)
+        vi = v[i]
+        isreal(vi) || continue
+        a = abs2(vi)
+        a < amax && continue
+        amax, imax = a, i
+    end
+    return imax
+end
+
+# The diagonal of `F` is a free gauge parameter: `U̇ + U * Diagonal(ċ)` solves the same
+# differentiated eigenvalue equation for any `ċ`. Setting it to zero, as `U̇ = U * (F .* M)`
+# does, normalizes the eigenvectors by `diag(inv(U) * U̇) == 0`, which for a non-normal `A`
+# is not the convention `eigen` itself returns. Choose `ċ` instead such that the derivatives
+# belong to LAPACK's convention, i.e. differentiate its two constraints at `t = 0`:
+#
+#   `u' * u == 1`             ⇒  `real(ċ) = -real(u' * u̇)`
+#   `imag(u[k]) == 0`         ⇒  `imag(ċ) = -imag(u̇[k]) / real(u[k])`,  `k = _findrealmaxabs2(u)`
+#
+# This mirrors `_eigen_norm_phase_fwd!` in ChainRules, so that both agree on `eigen`.
+function _eigen_norm_phase!(U̇, U)
+    for i in axes(U, 2)
+        u, u̇ = view(U, :, i), view(U̇, :, i)
+        ċ_norm = -real(dot(u, u̇))
+        if eltype(U) <: Real
+            ċ = ċ_norm
+        else
+            k = _findrealmaxabs2(u)
+            ċ = complex(ċ_norm, -imag(u̇[k]) / real(u[k]))
+        end
+        u̇ .+= u .* ċ
+    end
+    return U̇
+end
+
 LinearAlgebra.eigvals(A::StridedMatrix{Dual{Tg,T,N}}) where {Tg,T<:Real,N} = _eigvals_general(A)
 function _eigvals_general(A::StridedMatrix{Dual{Tg,T,N}}) where {Tg,T<:Real,N}
     λ, U = eigen(value.(A))
@@ -878,7 +919,7 @@ function _eigen_general(A::StridedMatrix{Dual{Tg,T,N}}) where {Tg,T<:Real,N}
     luU = lu(U)
     M = ntuple(j -> luU \ (getindex.(partials.(A), j) * U), N)
     λ_parts = map(diag, M)
-    U_parts = ntuple(j -> U * _lyap_div!!(M[j] - Diagonal(λ_parts[j]), λ), N)
+    U_parts = ntuple(j -> _eigen_norm_phase!(U * _lyap_div!!(M[j] - Diagonal(λ_parts[j]), λ), U), N)
     λ_dual = map((val, p) -> _make_eigen_dual(Dual{Tg}, val, p), λ, tuple.(λ_parts...))
     U_dual = map((val, p) -> _make_eigen_dual(Dual{Tg}, val, p), U, tuple.(U_parts...))
     return Eigen(λ_dual, U_dual)
diff --git a/test/HessianTest.jl b/test/HessianTest.jl
index 27995586..2261451d 100644
--- a/test/HessianTest.jl
+++ b/test/HessianTest.jl
@@ -184,12 +184,11 @@ end
     w111 = [0.9, 1.4, 0.3]
     @test ForwardDiff.hessian(g, w111) ≈ ForwardDiff.hessian(gsym, w111)
 
-    # eigenvectors: renormalising removes the difference between the eigenvector
-    # convention of `eigen` and the one the derivatives are derived for
+    # eigenvectors
     A0 = [2.0 1.0 0.5; 0.5 3.0 1.5; 0.25 0.75 4.0]
     function v1(x)
         v = eigen(reshape(x, 3, 3)).vectors[:, 1]
-        return sum(abs2, v / norm(v) .- [1.0, 0.5, -0.2])
+        return sum(abs2, v .- [1.0, 0.5, -0.2])
     end
     x0 = vec(A0)
     @test ForwardDiff.hessian(v1, x0) ≈ ForwardDiff.jacobian(x -> ForwardDiff.gradient(v1, x), x0)
diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl
index e72c0042..c1dc83b2 100644
--- a/test/JacobianTest.jl
+++ b/test/JacobianTest.jl
@@ -273,23 +273,45 @@ 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
+    # eigenvectors, deliberately without renormalizing: the derivatives have to belong to
+    # the normalization `eigen` itself returns, i.e. unit 2-norm with largest entry real
+    h(x) = vec(eigen(reshape(x, 2, 2)).vectors)
     x3 = [2.0, 1.0, 0.5, 3.0]
     @test ForwardDiff.jacobian(h, x3) ≈ Calculus.finite_difference_jacobian(h, x3)
 
+    # complex eigenvectors
+    hc(x) = begin
+        V = eigen(reshape(x, 2, 2)).vectors
+        vcat(real(vec(V)), imag(vec(V)))
+    end
+    x3c = vec([0.3 -1.2; 1.7 0.5])
+    @test ForwardDiff.jacobian(hc, x3c) ≈ Calculus.finite_difference_jacobian(hc, x3c)
+
     # larger than 2x2, non-symmetric with real eigenvalues
     f3(x) = eigvals(reshape(x, 3, 3))
     x4 = vec([2.0 1.0 0.5; 0.5 3.0 1.5; 0.25 0.75 4.0])
     @test ForwardDiff.jacobian(f3, x4) ≈ Calculus.finite_difference_jacobian(f3, x4)
-    h3(x) = begin
-        v = eigen(reshape(x, 3, 3)).vectors[:,2]
-        v = v / norm(v)
-    end
+    h3(x) = vec(eigen(reshape(x, 3, 3)).vectors)
     @test ForwardDiff.jacobian(h3, x4) ≈ Calculus.finite_difference_jacobian(h3, x4)
 
+    # 3x3 with one real and one complex conjugate pair of eigenvalues
+    g3(x) = begin
+        vals = eigvals(reshape(x, 3, 3))
+        vcat(real(vals), imag(vals))
+    end
+    hc3(x) = begin
+        V = eigen(reshape(x, 3, 3)).vectors
+        vcat(real(vec(V)), imag(vec(V)))
+    end
+    x4c = vec([0.5 -1.3 0.2; 1.1 0.4 -0.6; 0.3 0.7 2.0])
+    @test ForwardDiff.jacobian(g3, x4c) ≈ Calculus.finite_difference_jacobian(g3, x4c)
+    @test ForwardDiff.jacobian(hc3, x4c) ≈ Calculus.finite_difference_jacobian(hc3, x4c)
+
+    # the eigenvector derivatives used to belong to the normalization
+    # `diag(inv(U) * U̇) == 0` instead, which differs for a non-normal matrix
+    A_gauge = Dual{TestTag}.([1.0 1.0; 0.0 2.0], [1.0 0.0; 0.0 0.0])
+    @test ForwardDiff.partials.(eigen(A_gauge).vectors, 1) ≈ [0.0 1/(2*sqrt(2)); 0.0 -1/(2*sqrt(2))]
+
     # eltypes of the general path
     A_dual = Dual{TestTag}.([1.0 2.0; 3.0 4.0], [1.0 0.0; 0.0 0.0])
     @test eigvals(A_dual) isa Vector{Dual{TestTag,Float64,1}}

From 1666fb5f360dc2533239e16066e425c8328c2faa Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Mon, 10 Aug 2026 23:10:13 +0200
Subject: [PATCH 06/15] Forward keyword arguments in the general `eigen` and
 `eigvals`

`eigen(A; sortby=nothing)`, `eigvals(A; sortby=nothing)` and
`eigen(A; permute=false)` threw a `MethodError` for a `Dual` matrix,
while the internal decomposition of `value.(A)` silently applied the
default `sortby`.

Forward `kwargs...` to that decomposition. The derivatives are assembled
in whatever order it returns, so they follow the requested ordering, and
for nested `Dual`s every level is decomposed with the same keywords.
`eigvecs` picks this up through its generic fallback.

`eigen!` and `eigvals!` remain undefined.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 src/dual.jl          | 19 +++++++++++++------
 test/HessianTest.jl  |  4 ++++
 test/JacobianTest.jl | 15 +++++++++++++++
 3 files changed, 32 insertions(+), 6 deletions(-)

diff --git a/src/dual.jl b/src/dual.jl
index 33397faf..2109b8b8 100644
--- a/src/dual.jl
+++ b/src/dual.jl
@@ -905,17 +905,24 @@ function _eigen_norm_phase!(U̇, U)
     return U̇
 end
 
-LinearAlgebra.eigvals(A::StridedMatrix{Dual{Tg,T,N}}) where {Tg,T<:Real,N} = _eigvals_general(A)
-function _eigvals_general(A::StridedMatrix{Dual{Tg,T,N}}) where {Tg,T<:Real,N}
-    λ, U = eigen(value.(A))
+# `permute`, `scale` and `sortby` are forwarded to the underlying decomposition of
+# `value.(A)`; the derivatives are assembled in whatever order it returns, and for nested
+# `Dual`s every level is decomposed with the same keyword arguments
+function LinearAlgebra.eigvals(A::StridedMatrix{Dual{Tg,T,N}}; kwargs...) where {Tg,T<:Real,N}
+    return _eigvals_general(A; kwargs...)
+end
+function _eigvals_general(A::StridedMatrix{Dual{Tg,T,N}}; kwargs...) where {Tg,T<:Real,N}
+    λ, U = eigen(value.(A); kwargs...)
     luU = lu(U)
     parts = ntuple(j -> diag(luU \ (getindex.(partials.(A), j) * U)), N)
     return map((val, p) -> _make_eigen_dual(Dual{Tg}, val, p), λ, tuple.(parts...))
 end
 
-LinearAlgebra.eigen(A::StridedMatrix{Dual{Tg,T,N}}) where {Tg,T<:Real,N} = _eigen_general(A)
-function _eigen_general(A::StridedMatrix{Dual{Tg,T,N}}) where {Tg,T<:Real,N}
-    λ, U = eigen(value.(A))
+function LinearAlgebra.eigen(A::StridedMatrix{Dual{Tg,T,N}}; kwargs...) where {Tg,T<:Real,N}
+    return _eigen_general(A; kwargs...)
+end
+function _eigen_general(A::StridedMatrix{Dual{Tg,T,N}}; kwargs...) where {Tg,T<:Real,N}
+    λ, U = eigen(value.(A); kwargs...)
     luU = lu(U)
     M = ntuple(j -> luU \ (getindex.(partials.(A), j) * U), N)
     λ_parts = map(diag, M)
diff --git a/test/HessianTest.jl b/test/HessianTest.jl
index 2261451d..7e21c7a3 100644
--- a/test/HessianTest.jl
+++ b/test/HessianTest.jl
@@ -173,6 +173,10 @@ end
     # sum(eigvals(B(w)) .^ 2) == tr(B(w)^2) is quadratic in `w`
     @test ForwardDiff.hessian(w -> sum(eigvals(B(w)) .^ 2), w) ≈ [10 0; 0 4]
 
+    # keyword arguments reach every level of the nesting
+    @test ForwardDiff.hessian(w -> sum(eigvals(B(w); sortby = nothing)), w) ≈ zeros(2, 2) atol=1e-12
+    @test ForwardDiff.hessian(w -> sum(eigvals(B(w); permute = false, scale = false) .^ 2), w) ≈ [10 0; 0 4]
+
     # complex eigenvalues: λ = w[1] ± im*(1 + w[2]), i.e. sum(abs2, λ) == 2*(w[1]^2 + (1 + w[2])^2)
     C(w) = [w[1] -1.0-w[2]; 1.0+w[2] w[1]]
     @test ForwardDiff.hessian(w -> sum(abs2, eigvals(C(w))), [0.3, 0.2]) ≈ [4 0; 0 4]
diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl
index c1dc83b2..95537554 100644
--- a/test/JacobianTest.jl
+++ b/test/JacobianTest.jl
@@ -312,6 +312,21 @@ end
     A_gauge = Dual{TestTag}.([1.0 1.0; 0.0 2.0], [1.0 0.0; 0.0 0.0])
     @test ForwardDiff.partials.(eigen(A_gauge).vectors, 1) ≈ [0.0 1/(2*sqrt(2)); 0.0 -1/(2*sqrt(2))]
 
+    # keyword arguments are forwarded to the decomposition of the values
+    A_kw = reshape(x4, 3, 3)
+    A_kw_dual = Dual{TestTag}.(A_kw, Matrix(1.0I, 3, 3))
+    for kwargs in ((), (sortby = nothing,), (permute = false,), (scale = false,),
+                   (permute = false, scale = false), (sortby = λ -> -real(λ),))
+        @test ForwardDiff.value.(eigvals(A_kw_dual; kwargs...)) ≈ eigvals(A_kw; kwargs...)
+        @test eigvals(A_kw_dual; kwargs...) ≈ eigen(A_kw_dual; kwargs...).values
+        f_kw(x) = eigvals(reshape(x, 3, 3); kwargs...)
+        @test ForwardDiff.jacobian(f_kw, x4) ≈ Calculus.finite_difference_jacobian(f_kw, x4)
+    end
+    # and the derivatives follow the ordering they produce: `A_kw` has a real spectrum, so
+    # sorting by `-real` reverses the order `LinearAlgebra.eigsortby` gives
+    @test ForwardDiff.jacobian(x -> eigvals(reshape(x, 3, 3); sortby = λ -> -real(λ)), x4) ≈
+          ForwardDiff.jacobian(x -> eigvals(reshape(x, 3, 3)), x4)[3:-1:1, :]
+
     # eltypes of the general path
     A_dual = Dual{TestTag}.([1.0 2.0; 3.0 4.0], [1.0 0.0; 0.0 0.0])
     @test eigvals(A_dual) isa Vector{Dual{TestTag,Float64,1}}

From cee2467469ccf2547b789dedcc3b6b37a0f16447 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Mon, 10 Aug 2026 23:45:11 +0200
Subject: [PATCH 07/15] Error on repeated eigenvalues instead of returning
 `NaN`
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

`_lyap_div!!` divides by `λ[j] - λ[i]`, so a repeated eigenvalue made
`.vectors` come out as `NaN` with no indication of what went wrong. The
eigenvectors of a repeated eigenvalue are not unique and hence not
differentiable, so throw an `ArgumentError` naming the pair instead.

The check compares the primal values rather than the `Dual`s: `==` and
`iszero` on a `Dual` take the partials into account, but the quotient
already breaks down when only the values coincide, which is what happens
one level up in a `hessian`.

This covers the `Symmetric` and `SymTridiagonal` methods as well, which
share the same weakness, and the `MMatrix` method in the StaticArrays
extension through `_lyap_div!`. `eigvals` never divides by the gaps and
keeps working.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 src/dual.jl          | 26 ++++++++++++++++++++++++++
 test/JacobianTest.jl | 17 +++++++++++++++++
 2 files changed, 43 insertions(+)

diff --git a/src/dual.jl b/src/dual.jl
index 2109b8b8..f48b0b75 100644
--- a/src/dual.jl
+++ b/src/dual.jl
@@ -805,9 +805,34 @@ function LinearAlgebra.eigvals(A::SymTridiagonal{<:Dual{Tg,T,N}}) where {Tg,T<:R
     Dual{Tg}.(λ, tuple.(parts...))
 end
 
+# Strip every `Dual` level, so that eigenvalues can be compared by their primal value alone
+_primalvalue(x::Real) = x
+_primalvalue(d::Dual) = _primalvalue(value(d))
+_primalvalue(z::Complex) = complex(_primalvalue(real(z)), _primalvalue(imag(z)))
+
+@noinline function _throw_repeated_eigvals(i, j, λ)
+    throw(ArgumentError(
+        "eigenvector derivatives are not defined for repeated eigenvalues, but " *
+        "λ[$i] == λ[$j] == $λ"
+    ))
+end
+
+# `_lyap_div!!` divides by `λ[j] - λ[i]`, which vanishes for repeated eigenvalues: there the
+# eigenvectors are not unique and hence not differentiable. Comparing the primal values
+# catches the cases in which the quotient would silently come out as `NaN`.
+function _check_distinct_eigvals(λ::AbstractVector)
+    for j in eachindex(λ), i in eachindex(λ)
+        i < j || continue
+        vi, vj = _primalvalue(λ[i]), _primalvalue(λ[j])
+        vi == vj && _throw_repeated_eigvals(i, j, vi)
+    end
+    return nothing
+end
+
 # A ./ (λ' .- λ) but with diag special cased
 # Default out-of-place method
 function _lyap_div!!(A::AbstractMatrix, λ::AbstractVector)
+    _check_distinct_eigvals(λ)
     return map(
         (a, b, idx) -> a / (idx[1] == idx[2] ? oneunit(b) : b),
         A,
@@ -818,6 +843,7 @@ end
 # For `Matrix` (and e.g. `StaticArrays.MMatrix`) we can use an in-place method
 _lyap_div!!(A::Matrix, λ::AbstractVector) = _lyap_div!(A, λ)
 function _lyap_div!(A::AbstractMatrix, λ::AbstractVector)
+    _check_distinct_eigvals(λ)
     for (j,μ) in enumerate(λ), (k,λ) in enumerate(λ)
         if k ≠ j
             A[k,j] /= μ - λ
diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl
index 95537554..49ffa38c 100644
--- a/test/JacobianTest.jl
+++ b/test/JacobianTest.jl
@@ -327,6 +327,23 @@ end
     @test ForwardDiff.jacobian(x -> eigvals(reshape(x, 3, 3); sortby = λ -> -real(λ)), x4) ≈
           ForwardDiff.jacobian(x -> eigvals(reshape(x, 3, 3)), x4)[3:-1:1, :]
 
+    # repeated eigenvalues: the eigenvectors are not unique, so their derivatives do not
+    # exist. This used to be a silent `NaN` in `.vectors`.
+    for A_rep in (Dual{TestTag}.([2.0 0.0; 0.0 2.0], [1.0 0.0; 0.0 0.0]),   # diagonalizable
+                  Dual{TestTag}.([2.0 1.0; 0.0 2.0], [1.0 0.0; 0.0 0.0]))   # defective
+        @test_throws ArgumentError eigen(A_rep)
+    end
+    @test_throws ArgumentError eigen(Symmetric(Dual{TestTag}.([2.0 0.0; 0.0 2.0], [1.0 0.0; 0.0 0.0])))
+    @test_throws ArgumentError eigen(SymTridiagonal(Dual{TestTag}.([2.0, 2.0], [1.0, 0.0]),
+                                                    Dual{TestTag}.([0.0], [0.0])))
+    # `eigvals` never divides by the eigenvalue gaps and is unaffected
+    @test eigvals(Dual{TestTag}.([2.0 0.0; 0.0 2.0], [1.0 0.0; 0.0 0.0])) ==
+          Dual{TestTag}.([2.0, 2.0], [1.0, 0.0])
+    # equal values with differing partials would divide by a `Dual` with a zero value
+    A_nested = Dual{TestTag}.(Dual{TestTag}.([2.0 0.0; 0.0 2.0], [1.0 0.0; 0.0 0.0]),
+                              Dual{TestTag}.([0.0 1.0; 1.0 0.0], [0.0 0.0; 0.0 0.0]))
+    @test_throws ArgumentError eigen(A_nested)
+
     # eltypes of the general path
     A_dual = Dual{TestTag}.([1.0 2.0; 3.0 4.0], [1.0 0.0; 0.0 0.0])
     @test eigvals(A_dual) isa Vector{Dual{TestTag,Float64,1}}

From 273a9652341606e3f20fc69e5311fd00ee545837 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Mon, 10 Aug 2026 23:45:33 +0200
Subject: [PATCH 08/15] Dispatch a Hermitian `Matrix` to the `Symmetric`
 methods

`LinearAlgebra.eigen!` and `eigvals!` check `issymmetric`/`ishermitian`
and route to the symmetric algorithm, but a plain `Matrix{<:Dual}` went
through the general path regardless -- which is the situation in #111.

Do the same check on the `Dual` matrix. `==` on a `Dual` compares the
partials as well, so `ishermitian` only accepts a matrix whose values
*and* partials are Hermitian; a Hermitian value with a non-Hermitian
perturbation keeps using the general path, since symmetrizing it would
change the derivative.

The eigenvalues are then real by construction rather than by whether
`geevx!` produced a nonzero imaginary part, and the derivative avoids the
LU factorization in favour of `Q'`.

The `Symmetric` methods take no keyword arguments. `permute` and `scale`
are balancing options the symmetric algorithm does not use and are
ignored the same way `eigen!` ignores them, but any `sortby` other than
the ascending order it returns anyway stays on the general path.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 src/dual.jl          | 16 ++++++++++++++++
 test/JacobianTest.jl | 21 +++++++++++++++++++++
 2 files changed, 37 insertions(+)

diff --git a/src/dual.jl b/src/dual.jl
index f48b0b75..33847655 100644
--- a/src/dual.jl
+++ b/src/dual.jl
@@ -931,10 +931,25 @@ function _eigen_norm_phase!(U̇, U)
     return U̇
 end
 
+# A matrix that is Hermitian in both its values and its partials is handled by the
+# `Symmetric` methods, which is what `LinearAlgebra.eigen!` and `eigvals!` do for the
+# values as well. Their eigenvalues are real, so the eltype no longer depends on whether
+# `geevx!` happens to produce a nonzero imaginary part.
+#
+# Those methods take no keyword arguments. `permute` and `scale` only control the balancing
+# that the symmetric algorithm does not use, so they can be ignored the way
+# `LinearAlgebra.eigen!` ignores them, but any `sortby` other than the ascending order that
+# is returned anyway has to go through the general path.
+function _use_symmetric(A; permute::Bool=true, scale::Bool=true,
+                        sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby)
+    return (sortby === nothing || sortby === LinearAlgebra.eigsortby) && ishermitian(A)
+end
+
 # `permute`, `scale` and `sortby` are forwarded to the underlying decomposition of
 # `value.(A)`; the derivatives are assembled in whatever order it returns, and for nested
 # `Dual`s every level is decomposed with the same keyword arguments
 function LinearAlgebra.eigvals(A::StridedMatrix{Dual{Tg,T,N}}; kwargs...) where {Tg,T<:Real,N}
+    _use_symmetric(A; kwargs...) && return _eigvals(Symmetric(A))
     return _eigvals_general(A; kwargs...)
 end
 function _eigvals_general(A::StridedMatrix{Dual{Tg,T,N}}; kwargs...) where {Tg,T<:Real,N}
@@ -945,6 +960,7 @@ function _eigvals_general(A::StridedMatrix{Dual{Tg,T,N}}; kwargs...) where {Tg,T
 end
 
 function LinearAlgebra.eigen(A::StridedMatrix{Dual{Tg,T,N}}; kwargs...) where {Tg,T<:Real,N}
+    _use_symmetric(A; kwargs...) && return _eigen(Symmetric(A))
     return _eigen_general(A; kwargs...)
 end
 function _eigen_general(A::StridedMatrix{Dual{Tg,T,N}}; kwargs...) where {Tg,T<:Real,N}
diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl
index 49ffa38c..b1d65be4 100644
--- a/test/JacobianTest.jl
+++ b/test/JacobianTest.jl
@@ -327,6 +327,27 @@ end
     @test ForwardDiff.jacobian(x -> eigvals(reshape(x, 3, 3); sortby = λ -> -real(λ)), x4) ≈
           ForwardDiff.jacobian(x -> eigvals(reshape(x, 3, 3)), x4)[3:-1:1, :]
 
+    # a plain `Matrix` that is Hermitian in its values *and* its partials goes through the
+    # `Symmetric` methods, so both agree exactly
+    S0 = [2.0 1.0 0.5; 1.0 3.0 1.5; 0.5 1.5 4.0]
+    Ṡ = [0.5 -1.0 0.25; -1.0 2.0 0.75; 0.25 0.75 -0.5]
+    A_herm = Dual{TestTag}.(S0, Ṡ)
+    @test ishermitian(A_herm)
+    @test eigvals(A_herm) == eigvals(Symmetric(A_herm))
+    @test eigen(A_herm).values == eigen(Symmetric(A_herm)).values
+    @test eigen(A_herm).vectors == eigen(Symmetric(A_herm)).vectors
+    @test eigvals(A_herm) isa Vector{Dual{TestTag,Float64,1}}
+    # but a `sortby` the `Symmetric` methods cannot honour has to stay on the general path
+    @test ForwardDiff.partials.(eigvals(A_herm; sortby = λ -> -real(λ)), 1) ≈
+          reverse(ForwardDiff.partials.(eigvals(A_herm), 1))
+
+    # Hermitian values with non-Hermitian partials must *not* take the `Symmetric` path,
+    # which would symmetrize the perturbation and give a different derivative
+    A_asym = Dual{TestTag}.([1.0 2.0; 2.0 1.0], [0.0 1.0; 0.0 0.0])
+    @test !ishermitian(A_asym)
+    @test ForwardDiff.partials.(eigvals(A_asym), 1) ≈ [-0.5, 0.5]
+    @test ForwardDiff.partials.(eigvals(Symmetric(A_asym)), 1) ≈ [-1.0, 1.0]
+
     # repeated eigenvalues: the eigenvectors are not unique, so their derivatives do not
     # exist. This used to be a silent `NaN` in `.vectors`.
     for A_rep in (Dual{TestTag}.([2.0 0.0; 0.0 2.0], [1.0 0.0; 0.0 0.0]),   # diagonalizable

From 6535c41c40ea93034fa6cb594f2274ae1fa1a9aa Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Tue, 11 Aug 2026 00:12:05 +0200
Subject: [PATCH 09/15] Consume the internal temporaries in the general `eigen`
 and `eigvals`
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

`value.(A)` and `Ȧ * U` are allocated by these methods and used nowhere
else, so the decomposition and the solve may overwrite them: `eigen!`
instead of `eigen` for the former, `ldiv!` instead of `\` for the latter.
Both `eigen` and `\` copy their argument, so each was costing an extra
n^2 matrix -- once for the decomposition and once per partial.

`eigen!` only exists for BLAS element types, which is the innermost level
of a nested `Dual`; `_eigen!!` falls back to the copying `eigen` above it
and for value types such as `Float16`.

Results are unchanged bit for bit. Measured at n = 40 with 8 partials,
allocations drop by 27.9% for `eigvals` (416 KiB to 300 KiB, 1.04 ms to
0.70 ms) and by 13.8% for `eigen` (841 KiB to 725 KiB).

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 src/dual.jl | 15 +++++++++++----
 1 file changed, 11 insertions(+), 4 deletions(-)

diff --git a/src/dual.jl b/src/dual.jl
index 33847655..af70e5fa 100644
--- a/src/dual.jl
+++ b/src/dual.jl
@@ -945,6 +945,12 @@ function _use_symmetric(A; permute::Bool=true, scale::Bool=true,
     return (sortby === nothing || sortby === LinearAlgebra.eigsortby) && ishermitian(A)
 end
 
+# `value.(A)` is a temporary of our own, so the decomposition may consume it. `eigen!` only
+# exists for BLAS element types, which is the innermost level of the nesting; above it the
+# recursion goes through the copying `eigen` below. `!!` as in `_lyap_div!!`: may mutate.
+_eigen!!(B::StridedMatrix{<:LinearAlgebra.BlasFloat}; kwargs...) = eigen!(B; kwargs...)
+_eigen!!(B::AbstractMatrix; kwargs...) = eigen(B; kwargs...)
+
 # `permute`, `scale` and `sortby` are forwarded to the underlying decomposition of
 # `value.(A)`; the derivatives are assembled in whatever order it returns, and for nested
 # `Dual`s every level is decomposed with the same keyword arguments
@@ -953,9 +959,10 @@ function LinearAlgebra.eigvals(A::StridedMatrix{Dual{Tg,T,N}}; kwargs...) where
     return _eigvals_general(A; kwargs...)
 end
 function _eigvals_general(A::StridedMatrix{Dual{Tg,T,N}}; kwargs...) where {Tg,T<:Real,N}
-    λ, U = eigen(value.(A); kwargs...)
+    λ, U = _eigen!!(value.(A); kwargs...)
     luU = lu(U)
-    parts = ntuple(j -> diag(luU \ (getindex.(partials.(A), j) * U)), N)
+    # `Ȧ * U` is a temporary as well, so the solve can overwrite it
+    parts = ntuple(j -> diag(ldiv!(luU, getindex.(partials.(A), j) * U)), N)
     return map((val, p) -> _make_eigen_dual(Dual{Tg}, val, p), λ, tuple.(parts...))
 end
 
@@ -964,9 +971,9 @@ function LinearAlgebra.eigen(A::StridedMatrix{Dual{Tg,T,N}}; kwargs...) where {T
     return _eigen_general(A; kwargs...)
 end
 function _eigen_general(A::StridedMatrix{Dual{Tg,T,N}}; kwargs...) where {Tg,T<:Real,N}
-    λ, U = eigen(value.(A); kwargs...)
+    λ, U = _eigen!!(value.(A); kwargs...)
     luU = lu(U)
-    M = ntuple(j -> luU \ (getindex.(partials.(A), j) * U), N)
+    M = ntuple(j -> ldiv!(luU, getindex.(partials.(A), j) * U), N)
     λ_parts = map(diag, M)
     U_parts = ntuple(j -> _eigen_norm_phase!(U * _lyap_div!!(M[j] - Diagonal(λ_parts[j]), λ), U), N)
     λ_dual = map((val, p) -> _make_eigen_dual(Dual{Tg}, val, p), λ, tuple.(λ_parts...))

From 2130a6dc8d53ef2a6a993a8044caa37104dab1bc Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Tue, 11 Aug 2026 10:34:59 +0200
Subject: [PATCH 10/15] Compare eigenvalues without stripping their primal
 values

Review feedback on the repeated-eigenvalue check:

- Do not extract primal values unconditionally. `_primalvalue` stripped
  every `Dual` level, including levels belonging to an independent,
  enclosing differentiation. It turns out no stripping is needed at all:
  when two eigenvalues are `!=` but share a primal, the decomposition of
  the values -- which every method computes first, one `Dual` level down
  -- has already seen them as exact duplicates and thrown. Each level of
  the nesting runs its own check, so `==` is enough.

- Make the error message lazy, as in `throw_cannot_dual`.

- The comment claimed `_lyap_div!!` did not handle a zero denominator. It
  does, but only on the diagonal, where the difference vanishes by
  construction; what is unhandled is an off-diagonal denominator that
  vanishes because two eigenvalues coincide. Say that, and say `Inf` or
  `NaN` rather than just `NaN`.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 src/dual.jl | 25 +++++++++++--------------
 1 file changed, 11 insertions(+), 14 deletions(-)

diff --git a/src/dual.jl b/src/dual.jl
index af70e5fa..51fea9ae 100644
--- a/src/dual.jl
+++ b/src/dual.jl
@@ -805,26 +805,23 @@ function LinearAlgebra.eigvals(A::SymTridiagonal{<:Dual{Tg,T,N}}) where {Tg,T<:R
     Dual{Tg}.(λ, tuple.(parts...))
 end
 
-# Strip every `Dual` level, so that eigenvalues can be compared by their primal value alone
-_primalvalue(x::Real) = x
-_primalvalue(d::Dual) = _primalvalue(value(d))
-_primalvalue(z::Complex) = complex(_primalvalue(real(z)), _primalvalue(imag(z)))
-
 @noinline function _throw_repeated_eigvals(i, j, λ)
-    throw(ArgumentError(
-        "eigenvector derivatives are not defined for repeated eigenvalues, but " *
-        "λ[$i] == λ[$j] == $λ"
-    ))
+    throw(ArgumentError(lazy"eigenvector derivatives are not defined for repeated eigenvalues, but λ[$i] == λ[$j] == $λ"))
 end
 
-# `_lyap_div!!` divides by `λ[j] - λ[i]`, which vanishes for repeated eigenvalues: there the
-# eigenvectors are not unique and hence not differentiable. Comparing the primal values
-# catches the cases in which the quotient would silently come out as `NaN`.
+# `_lyap_div!!` special cases only the diagonal, where `λ[j] - λ[i]` vanishes by
+# construction. Two equal eigenvalues make an off-diagonal denominator vanish as well, and
+# the eigenvector partials come out as `Inf` or `NaN`: the eigenvectors of a repeated
+# eigenvalue are not unique and hence not differentiable.
+#
+# The eigenvalues are compared as they are, without extracting primal values: for nested
+# `Dual`s two of them can be `!=` here and still divide to `Inf`, but then their primals
+# coincide, and the decomposition of the values -- which every method computes first, one
+# `Dual` level down -- has already seen them as exact duplicates and thrown.
 function _check_distinct_eigvals(λ::AbstractVector)
     for j in eachindex(λ), i in eachindex(λ)
         i < j || continue
-        vi, vj = _primalvalue(λ[i]), _primalvalue(λ[j])
-        vi == vj && _throw_repeated_eigvals(i, j, vi)
+        λ[i] == λ[j] && _throw_repeated_eigvals(i, j, λ[i])
     end
     return nothing
 end

From 7660612e3cf5b03a378c2271776ad9c0ba73e941 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Fri, 14 Aug 2026 22:16:33 +0200
Subject: [PATCH 11/15] Fix the eigenvector phase gauge at the wrong entry
 under nesting
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

`eigen` returns eigenvectors whose largest-magnitude entry is real, so
`imag(u[k]) == 0` holds identically along the curve and `imag(u̇[k])` is
exactly zero. It was left as the cancellation `q + fl(r * fl(-q / r))`
instead, which keeps a rounding residue of about `5.6e-17`.

At first order the residue is harmless. One differentiation level up it
sits in the partials, where `isreal` sees it -- `isreal` on a
`Complex{<:Dual}` is `iszero(imag(...))`, and `iszero` on a `Dual`
requires the partials to vanish too. `_findrealmaxabs2` then finds no
real entry in the column at all and silently returned `firstindex(v)`, so
the phase was fixed at the wrong entry and second and higher derivatives
through complex eigenvectors were wrong from `n = 4` on.

Assert the imaginary part instead of computing it, which is exact and so
fixes every order, and make the fallback in `_findrealmaxabs2` throw:
reaching it means the assumed normalization does not hold.

This is not inherited from ChainRules. `_eigen_norm_phase_fwd!` only ever
runs on `StridedMatrix{<:BlasFloat}` there, where `isreal` is exact.

Tested by asserting the two gauge constraints, `imag(u[k]) == 0` and
`u' * u == 1`, at nesting depths 1, 2 and 3 for n = 2 to 6, and by
comparing second derivatives of the eigenvectors against central
differences of the first derivatives. `eigen` does not pin the sign of
the real entry, so the columns are realigned before differencing; plain
central differences are off by O(1) and look like an AD bug. Before this
commit those tests fail 29 times, with the second derivative off by 55.6
at n = 4 and 7.0 at n = 5; n = 2 passes either way.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 src/dual.jl         | 23 +++++++++--
 test/HessianTest.jl | 99 +++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 118 insertions(+), 4 deletions(-)

diff --git a/src/dual.jl b/src/dual.jl
index 51fea9ae..4dc21bf5 100644
--- a/src/dual.jl
+++ b/src/dual.jl
@@ -887,19 +887,28 @@ _make_eigen_dual(D::Type, val::Complex, parts::NTuple{N,Number}) where {N} =
 # with `A = U * Diagonal(λ) * inv(U)` and `M = inv(U) * Ȧ * U`, we have `λ̇ = diag(M)`
 # and `U̇ = U * (F .* M)`, where `F[i,j] = inv(λ[j] - λ[i])` off the diagonal and zero on it.
 
+@noinline function _throw_no_real_eigvec_entry()
+    throw(ArgumentError("no entry of an eigenvector is real, so the entry that carries the phase convention cannot be identified; the derivatives assume the normalization of `eigen`, whose eigenvectors have a real entry of largest magnitude"))
+end
+
 # Index of the entry of largest magnitude among the real entries of `v`. LAPACK normalizes
 # the eigenvectors it returns to unit 2-norm with their largest entry real, so this is the
-# entry that carries the phase convention. Modelled after `_findrealmaxabs2` in ChainRules.
+# entry that carries the phase convention. Modelled after `_findrealmaxabs2` in ChainRules,
+# except that reaching the end without a real entry is an error rather than a silent
+# fallback: it means the assumed normalization does not hold.
 function _findrealmaxabs2(v)
     imax = firstindex(v)
     amax = abs2(zero(eltype(v)))
+    found = false
     for i in eachindex(v)
         vi = v[i]
         isreal(vi) || continue
         a = abs2(vi)
         a < amax && continue
         amax, imax = a, i
+        found = true
     end
+    found || _throw_no_real_eigvec_entry()
     return imax
 end
 
@@ -918,12 +927,18 @@ function _eigen_norm_phase!(U̇, U)
         u, u̇ = view(U, :, i), view(U̇, :, i)
         ċ_norm = -real(dot(u, u̇))
         if eltype(U) <: Real
-            ċ = ċ_norm
+            u̇ .+= u .* ċ_norm
         else
             k = _findrealmaxabs2(u)
-            ċ = complex(ċ_norm, -imag(u̇[k]) / real(u[k]))
+            u̇ .+= u .* complex(ċ_norm, -imag(u̇[k]) / real(u[k]))
+            # `imag(u[k]) == 0` holds identically along the curve, so `imag(u̇[k])` is exactly
+            # zero. Leaving it as the cancellation `q + fl(r * fl(-q / r))` keeps a rounding
+            # residue instead, and one differentiation level up that residue sits in the
+            # partials, where `isreal` sees it: `_findrealmaxabs2` then finds no real entry
+            # in the column and fixes the phase at the wrong one. Being exact rather than
+            # approximate, this also keeps third and higher derivatives right.
+            u̇[k] = complex(real(u̇[k]), zero(real(u̇[k])))
         end
-        u̇ .+= u .* ċ
     end
     return U̇
 end
diff --git a/test/HessianTest.jl b/test/HessianTest.jl
index 7e21c7a3..2e0166f5 100644
--- a/test/HessianTest.jl
+++ b/test/HessianTest.jl
@@ -199,6 +199,105 @@ end
     @test ForwardDiff.hessian(v1, x0) ≈ Calculus.finite_difference_jacobian(x -> ForwardDiff.gradient(v1, x), x0) atol=1e-5
 end
 
+# Helpers for the eigenvector phase gauge tests below
+struct GaugeTag end
+
+# A deterministic `n x n` matrix with a complex spectrum and well separated eigenvalues
+function complex_spectrum_matrix(n)
+    A = zeros(n, n)
+    for b in 1:(n ÷ 2)
+        i = 2b - 1
+        A[i, i] = A[i+1, i+1] = 0.5 + b / 4
+        A[i, i+1] = -1.0 - b / 8
+        A[i+1, i] = 1.0 + b / 8
+    end
+    isodd(n) && (A[n, n] = 2.0)
+    for i in 1:n, j in 1:n
+        A[i, j] += 0.15 * sinpi((i + 2j) / (n + 1))
+    end
+    return A
+end
+
+seed_matrix(n, i, j) = (S = zeros(n, n); S[i, j] = 1.0; S)
+
+# `A` lifted to a `Dual` of nesting depth `k` with vanishing partials
+lift_dual(A, k) = k == 0 ? A : Dual{GaugeTag}.(lift_dual(A, k - 1), lift_dual(zero(A), k - 1))
+# `A` seeded with one partial per level, innermost seed first, the way `hessian` nests them
+function nest_dual(A, seeds...)
+    M = A
+    for (k, seed) in enumerate(seeds)
+        M = Dual{GaugeTag}.(M, lift_dual(seed, k - 1))
+    end
+    return M
+end
+
+# The entry `eigen` made real, i.e. the one of largest magnitude
+phase_index(v) = argmax(j -> abs2(v[j]), eachindex(v))
+
+# All components of a (possibly nested) `Dual`: its value and every partial, recursively
+function dual_components!(out, x)
+    if x isa Dual
+        dual_components!(out, ForwardDiff.value(x))
+        for i in 1:ForwardDiff.npartials(x)
+            dual_components!(out, ForwardDiff.partials(x, i))
+        end
+    else
+        push!(out, x)
+    end
+    return out
+end
+dual_components(x) = dual_components!(Float64[], x)
+
+# The eigenvectors and their first derivative in direction `seed`, with the columns flipped
+# to match the signs of `Vref`
+function eigvecs_and_derivative(X, seed, Vref)
+    V = eigen(Dual{GaugeTag}.(X, seed)).vectors
+    V0 = map(z -> complex(ForwardDiff.value(real(z)), ForwardDiff.value(imag(z))), V)
+    V1 = map(z -> complex(ForwardDiff.partials(real(z), 1), ForwardDiff.partials(imag(z), 1)), V)
+    for i in axes(V0, 2)
+        k = phase_index(view(Vref, :, i))
+        if real(V0[k, i]) * real(Vref[k, i]) < 0
+            V0[:, i] .*= -1
+            V1[:, i] .*= -1
+        end
+    end
+    return V0, V1
+end
+
+@testset "eigenvector phase gauge under nesting, n = $n" for n in 2:6
+    A = complex_spectrum_matrix(n)
+    # otherwise the phase convention never comes up
+    @test !isreal(eigvals(A))
+    seeds = (seed_matrix(n, 1, 1), seed_matrix(n, 2, min(3, n)), seed_matrix(n, min(3, n), 1))
+
+    # `eigen` returns eigenvectors of unit 2-norm whose largest-magnitude entry is real, so
+    # `imag(u[k]) == 0` and `u' * u == 1` hold identically along the curve and every
+    # derivative of them has to vanish. A rounding residue in the first is invisible at
+    # first order, but one level up it sits in the partials, where `isreal` sees it, and
+    # `_findrealmaxabs2` then fixes the phase at the wrong entry.
+    @testset "nesting depth $depth" for depth in 1:3
+        V = eigen(nest_dual(A, seeds[1:depth]...)).vectors
+        for i in axes(V, 2)
+            v = view(V, :, i)
+            @test iszero(imag(V[phase_index(v), i]))
+            @test all(x -> abs(x) < 1e-12, dual_components(real(dot(v, v)) - 1))
+        end
+    end
+
+    # second derivatives of the eigenvectors, against central differences of the first
+    # derivatives. `eigen` does not pin the sign of the real entry, so the columns have to
+    # be realigned before differencing; plain central differences are off by `O(1)`.
+    Ea, Eb = seeds[1], seeds[2]
+    Vref = eigen(A).vectors
+    V = eigen(nest_dual(A, Ea, Eb)).vectors
+    ad = map(z -> complex(ForwardDiff.partials(ForwardDiff.partials(real(z), 1), 1),
+                          ForwardDiff.partials(ForwardDiff.partials(imag(z), 1), 1)), V)
+    h = 1e-5
+    _, Dp = eigvecs_and_derivative(A .+ h .* Eb, Ea, Vref)
+    _, Dm = eigvecs_and_derivative(A .- h .* Eb, Ea, Vref)
+    @test maximum(abs, ad .- (Dp .- Dm) ./ (2h)) < 1e-6
+end
+
 #https://github.com/JuliaDiff/ForwardDiff.jl/issues/720
 @testset "allocation-free hessian with StaticArrays" begin
     function hessian_allocs()

From ed1ad14105134e8110c88e3f4b18e23da7a533b3 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Fri, 14 Aug 2026 22:49:19 +0200
Subject: [PATCH 12/15] Test the general path for #111, and complex
 eigenvectors beyond 3x3

`S(w) = [w[1]^2 w[1]*w[2]*w[3]; w[1]*w[2]*w[3] w[2]^2]` is Hermitian in
its values *and* its partials, so `g` and `gsym` both take the
`Symmetric` shortcut and the two Hessians are bitwise identical. The test
pinned the Hermitian dispatch, not the general path that #111 needs.

Keep it, as the dispatch test it actually is, and add the same
log-determinant Hessian on a matrix that is not Hermitian.
`det(Bgen(w)) == w[1]^2 * (1 + w[2]^2 / 2)`, so there is a closed form to
compare against that never goes through `eigen`; the analytic Hessian is
matched to 3.6e-15, rather than to finite-difference accuracy.

The complex eigenvector tests were all 2x2 and 3x3, sizes at which the
entry carrying the phase convention can come out right by accident, so
add a 4x4 with two complex conjugate pairs. That matrix is also chosen so
that plain central differences are stable: `eigen` does not pin the sign
of the real entry, and on the 4x4 of the review the sign of column 3
flips under a 1e-6 perturbation, which makes the finite-difference
baseline wrong by 1e5 on a correct implementation.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 test/HessianTest.jl  | 19 +++++++++++++++++--
 test/JacobianTest.jl | 16 ++++++++++++++++
 2 files changed, 33 insertions(+), 2 deletions(-)

diff --git a/test/HessianTest.jl b/test/HessianTest.jl
index 2e0166f5..f2f0c408 100644
--- a/test/HessianTest.jl
+++ b/test/HessianTest.jl
@@ -181,12 +181,27 @@ end
     C(w) = [w[1] -1.0-w[2]; 1.0+w[2] w[1]]
     @test ForwardDiff.hessian(w -> sum(abs2, eigvals(C(w))), [0.3, 0.2]) ≈ [4 0; 0 4]
 
-    # https://github.com/JuliaDiff/ForwardDiff.jl/issues/111
+    # https://github.com/JuliaDiff/ForwardDiff.jl/issues/111. `S(w)` is Hermitian in its
+    # values *and* its partials, so both sides take the `Symmetric` shortcut: this pins the
+    # dispatch, not the general path, and the two Hessians are bitwise identical.
     S(w) = [w[1]^2 w[1]*w[2]*w[3]; w[1]*w[2]*w[3] w[2]^2]
     g(w) = sum(log, eigvals(S(w)))
     gsym(w) = sum(log, eigvals(Symmetric(S(w))))
     w111 = [0.9, 1.4, 0.3]
-    @test ForwardDiff.hessian(g, w111) ≈ ForwardDiff.hessian(gsym, w111)
+    @test ishermitian(S(w111))
+    @test ForwardDiff.hessian(g, w111) == ForwardDiff.hessian(gsym, w111)
+
+    # The same log-determinant Hessian on a matrix that is *not* Hermitian, so that the
+    # general path is what is under test. `det(Bgen(w)) == w[1]^2 * (1 + w[2]^2 / 2)`, which
+    # gives a closed-form reference that never goes through `eigen`.
+    Bgen(w) = [w[1]^2 w[1]*w[2]; 0.5*w[1]*w[2] w[2]^2+1]
+    hgen(w) = sum(log, eigvals(Bgen(w)))
+    wgen = [0.9, 1.4]
+    @test !ishermitian(Bgen(wgen))
+    @test isreal(eigvals(Bgen(wgen)))
+    @test ForwardDiff.hessian(hgen, wgen) ≈
+          [-2/wgen[1]^2 0; 0 (1 - wgen[2]^2/2)/(1 + wgen[2]^2/2)^2]
+    @test ForwardDiff.hessian(hgen, wgen) ≈ ForwardDiff.hessian(w -> log(det(Bgen(w))), wgen)
 
     # eigenvectors
     A0 = [2.0 1.0 0.5; 0.5 3.0 1.5; 0.25 0.75 4.0]
diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl
index b1d65be4..7fde9ea2 100644
--- a/test/JacobianTest.jl
+++ b/test/JacobianTest.jl
@@ -307,6 +307,22 @@ end
     @test ForwardDiff.jacobian(g3, x4c) ≈ Calculus.finite_difference_jacobian(g3, x4c)
     @test ForwardDiff.jacobian(hc3, x4c) ≈ Calculus.finite_difference_jacobian(hc3, x4c)
 
+    # 4x4 with two complex conjugate pairs. At 2x2 and 3x3 the entry that carries the phase
+    # convention can come out right by accident, so the complex eigenvector derivatives need
+    # to be pinned at a larger size as well.
+    g4(x) = begin
+        vals = eigvals(reshape(x, 4, 4))
+        vcat(real(vals), imag(vals))
+    end
+    hc4(x) = begin
+        V = eigen(reshape(x, 4, 4)).vectors
+        vcat(real(vec(V)), imag(vec(V)))
+    end
+    x5c = vec([1.0 -2.0 0.5 0.0; 2.0 1.0 0.0 0.5; 0.0 0.5 2.0 -1.0; 0.5 0.0 1.0 2.0])
+    @test !isreal(eigvals(reshape(x5c, 4, 4)))
+    @test ForwardDiff.jacobian(g4, x5c) ≈ Calculus.finite_difference_jacobian(g4, x5c)
+    @test ForwardDiff.jacobian(hc4, x5c) ≈ Calculus.finite_difference_jacobian(hc4, x5c)
+
     # the eigenvector derivatives used to belong to the normalization
     # `diag(inv(U) * U̇) == 0` instead, which differs for a non-normal matrix
     A_gauge = Dual{TestTag}.([1.0 1.0; 0.0 2.0], [1.0 0.0; 0.0 0.0])

From f30acb9d0d6c854e0cbeffa718fbab2d2a74f952 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Fri, 14 Aug 2026 23:16:18 +0200
Subject: [PATCH 13/15] Check for repeated eigenvalues once, in the `eigen`
 methods
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

`_check_distinct_eigvals` sat inside `_lyap_div!!`, which the `eigen`
methods call once per partial direction, so it ran `N` times per
decomposition instead of once. Hoist it into `_eigen`, `eigen` for
`SymTridiagonal` and `_eigen_general`, right after the decomposition,
which also leaves the `_lyap_div` kernel purely numeric.

In `_eigen_general` it now runs before `lu(U)`, so a repeated eigenvalue
cannot surface as a `SingularException` from factorizing a defective `U`.
That ordering is defensive rather than an observed fix: for every Jordan
block and repeated-diagonal case tried, `lu(U)` succeeds, because the
near-duplicate eigenvector columns LAPACK returns differ by about 1e-16
rather than being exactly equal.

Hoisting `value.(λ)` out of the `ntuple` in the two Hermitian methods
saves the `N - 1` redundant copies as well: 2800 bytes of 1060144 at
n = 40 with 8 partials.

All seven degenerate paths still throw the same `ArgumentError`: the
general one for a diagonalizable and for a defective matrix, `Symmetric`,
`SymTridiagonal`, a nested `Dual`, and both the `MMatrix` and `SMatrix`
paths of the StaticArrays extension.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 src/dual.jl | 17 ++++++++++++-----
 1 file changed, 12 insertions(+), 5 deletions(-)

diff --git a/src/dual.jl b/src/dual.jl
index 4dc21bf5..ae235b2b 100644
--- a/src/dual.jl
+++ b/src/dual.jl
@@ -812,7 +812,9 @@ end
 # `_lyap_div!!` special cases only the diagonal, where `λ[j] - λ[i]` vanishes by
 # construction. Two equal eigenvalues make an off-diagonal denominator vanish as well, and
 # the eigenvector partials come out as `Inf` or `NaN`: the eigenvectors of a repeated
-# eigenvalue are not unique and hence not differentiable.
+# eigenvalue are not unique and hence not differentiable. The `eigen` methods below call
+# this once, right after the decomposition and before anything can fail for another reason,
+# rather than from inside `_lyap_div!!`, which runs once per partial direction.
 #
 # The eigenvalues are compared as they are, without extracting primal values: for nested
 # `Dual`s two of them can be `!=` here and still divide to `Inf`, but then their primals
@@ -829,7 +831,6 @@ end
 # A ./ (λ' .- λ) but with diag special cased
 # Default out-of-place method
 function _lyap_div!!(A::AbstractMatrix, λ::AbstractVector)
-    _check_distinct_eigvals(λ)
     return map(
         (a, b, idx) -> a / (idx[1] == idx[2] ? oneunit(b) : b),
         A,
@@ -840,7 +841,6 @@ end
 # For `Matrix` (and e.g. `StaticArrays.MMatrix`) we can use an in-place method
 _lyap_div!!(A::Matrix, λ::AbstractVector) = _lyap_div!(A, λ)
 function _lyap_div!(A::AbstractMatrix, λ::AbstractVector)
-    _check_distinct_eigvals(λ)
     for (j,μ) in enumerate(λ), (k,λ) in enumerate(λ)
         if k ≠ j
             A[k,j] /= μ - λ
@@ -856,14 +856,18 @@ LinearAlgebra.eigen(A::Symmetric{<:Dual{Tg,T,N}}) where {Tg,T<:Real,N} = _eigen(
 function _eigen(A::Symmetric{<:Dual{Tg,T,N}}) 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)
+    λvals = value.(λ)
+    _check_distinct_eigvals(λvals)
+    parts = ntuple(j -> Q*_lyap_div!!(Q' * getindex.(partials.(A), j) * Q - Diagonal(getindex.(partials.(λ), j)), λvals), 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))))
-    parts = ntuple(j -> Q*_lyap_div!!(Q' * getindex.(partials.(A), j) * Q - Diagonal(getindex.(partials.(λ), j)), value.(λ)), N)
+    λvals = value.(λ)
+    _check_distinct_eigvals(λvals)
+    parts = ntuple(j -> Q*_lyap_div!!(Q' * getindex.(partials.(A), j) * Q - Diagonal(getindex.(partials.(λ), j)), λvals), N)
     Eigen(λ,Dual{Tg}.(Q, tuple.(parts...)))
 end
 
@@ -984,6 +988,9 @@ function LinearAlgebra.eigen(A::StridedMatrix{Dual{Tg,T,N}}; kwargs...) where {T
 end
 function _eigen_general(A::StridedMatrix{Dual{Tg,T,N}}; kwargs...) where {Tg,T<:Real,N}
     λ, U = _eigen!!(value.(A); kwargs...)
+    # before `lu`, so that a repeated eigenvalue is reported as such rather than surfacing
+    # as a `SingularException` from the factorization of a defective `U`
+    _check_distinct_eigvals(λ)
     luU = lu(U)
     M = ntuple(j -> ldiv!(luU, getindex.(partials.(A), j) * U), N)
     λ_parts = map(diag, M)

From 808abfb99d56f736d68a456ae0d72469f0dad6f1 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Fri, 11 Sep 2026 15:16:36 +0200
Subject: [PATCH 14/15] Say what `eigvals` actually returns for a repeated
 eigenvalue

The comment claimed `eigvals` "never divides by the eigenvalue gaps and
is unaffected". Not dividing is why it returns rather than throwing, but
the values are the diagonal of the compressed perturbation on the
degenerate eigenspace, not the derivatives of the eigenvalue branches.
They coincide only when the perturbation leaves the degeneracy unsplit,
which is the case in this test, and even there the sorted output is
one-sided, since `2` and `2 + t` swap order at `t = 0`.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 test/JacobianTest.jl | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl
index 744c9aac..0913501a 100644
--- a/test/JacobianTest.jl
+++ b/test/JacobianTest.jl
@@ -449,7 +449,11 @@ end
     @test_throws ArgumentError eigen(Symmetric(Dual{TestTag}.([2.0 0.0; 0.0 2.0], [1.0 0.0; 0.0 0.0])))
     @test_throws ArgumentError eigen(SymTridiagonal(Dual{TestTag}.([2.0, 2.0], [1.0, 0.0]),
                                                     Dual{TestTag}.([0.0], [0.0])))
-    # `eigvals` never divides by the eigenvalue gaps and is unaffected
+    # `eigvals` does not divide by the eigenvalue gaps, so it returns rather than throwing.
+    # What it returns is the diagonal of the compressed perturbation, which coincides with
+    # the derivative of the eigenvalue branches only when the perturbation leaves the
+    # degeneracy unsplit, as it does here -- and `2` and `2 + t` still swap order at
+    # `t = 0`, so even this result is one-sided. What is pinned is that `eigvals` still runs.
     @test eigvals(Dual{TestTag}.([2.0 0.0; 0.0 2.0], [1.0 0.0; 0.0 0.0])) ==
           Dual{TestTag}.([2.0, 2.0], [1.0, 0.0])
     # equal values with differing partials would divide by a `Dual` with a zero value

From 24f8600f2404cf4cf87d7bc9f02e9bc22bb19824 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Fri, 11 Sep 2026 15:37:13 +0200
Subject: [PATCH 15/15] Take the eigenvalue derivatives without the second
 matmul
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

`diag(inv(U) * Ȧ * U)` was formed as `inv(U) * (Ȧ * U)`, computing a full
`n^3` product to keep only its diagonal. Entry `i` of that diagonal is
row `i` of `inv(U) * Ȧ` against column `i` of `U`, so the second product
can be `n^2` multiplications instead, as in ChainRules' `eigvals` frule.
`_eigen_general` still needs the whole matrix and is unchanged.

Solving into one reusable `n^2` buffer rather than a fresh product per
direction removes the rest of the per-direction allocation.

                             before                 after
  n = 40,  N = 8,  real    646 us  296168 B     613 us  139728 B
  n = 40,  N = 8,  complex 909 us  752688 B     804 us  202456 B
  n = 150, N = 12, real     15.9 ms 5061096 B    13.7 ms 1121584 B
  n = 150, N = 12, complex  24.7 ms 14279120 B   20.5 ms 1904440 B

The eigenvalues are unchanged bit for bit, since they come straight from
the decomposition. The partials are not: summing `n` products in place of
a `gemm` reorders the additions, and they agree to a relative 1.5e-16 at
n = 3 and 1.0e-14 at n = 150.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 src/dual.jl | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

diff --git a/src/dual.jl b/src/dual.jl
index b066a944..f975b987 100644
--- a/src/dual.jl
+++ b/src/dual.jl
@@ -987,8 +987,16 @@ end
 function _eigvals_general(A::StridedMatrix{Dual{Tg,T,N}}; kwargs...) where {Tg,T<:Real,N}
     λ, U = _eigen!!(value.(A); kwargs...)
     luU = lu(U)
-    # `Ȧ * U` is a temporary as well, so the solve can overwrite it
-    parts = ntuple(j -> diag(ldiv!(luU, partials.(A, j) * U)), N)
+    # `diag(inv(U) * Ȧ * U)` without forming the second product: entry `i` is row `i` of
+    # `inv(U) * Ȧ` against column `i` of `U`, which is `n` multiplications rather than a
+    # matmul. Only the diagonal is wanted here, unlike in `_eigen_general`. One `n^2` buffer
+    # serves every direction, so the loop allocates only the result vectors.
+    B = similar(U)
+    parts = ntuple(N) do j
+        B .= partials.(A, j)
+        ldiv!(luU, B)
+        map((b, u) -> sum(prod, zip(b, u)), eachrow(B), eachcol(U))
+    end
     return _to_duals(Val(Tg), λ, parts)
 end