From f8a0cbee7fe599a75ab02f54021e45c1f19454cb Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Fri, 2 Oct 2026 08:57:36 +0200
Subject: [PATCH 1/2] Fix overflow in `eig_trunc_pullback!` and use
 `accelerative_smith_iteration!`

---
 src/pullbacks/eig.jl | 39 +++++++++++----------------------------
 1 file changed, 11 insertions(+), 28 deletions(-)

diff --git a/src/pullbacks/eig.jl b/src/pullbacks/eig.jl
index ad2fa2464..2e48d5f85 100755
--- a/src/pullbacks/eig.jl
+++ b/src/pullbacks/eig.jl
@@ -161,36 +161,19 @@ function eig_trunc_pullback!(
     Z = ViG * VᴴΔAV
 
     # add contribution from orthogonal complement
-    AP = mul!(complex.(A), V * Dmat, ViG', -1, 1)
-    X₀ = iszerotangent(ΔV₊) ? AP' * Z : mul!(ΔV₊, AP', Z, 1, 1)
+    # build the adjoint of AP directly, since that is what the series is summed with
+    APᴴ = mul!(complex.(A'), ViG, (V * Dmat)', -1, 1)
+    X₀ = iszerotangent(ΔV₊) ? APᴴ * Z : mul!(ΔV₊, APᴴ, Z, 1, 1)
     X₀ ./= D'
-    dabsmax = maximum(abs, D)
-    AP ./= dabsmax
-    D̄⁻¹ = dabsmax ./ conj.(D)
-    X₁ = rmul!(AP' * X₀, Diagonal(D̄⁻¹))
-    X₁ .+= X₀
-    Xₖ, Xₖ₊₁ = X₁, X₀
-    APₖ, APₖ₊₁ = AP * AP, AP
-    D̄⁻¹ₖ, D̄⁻¹ₖ₊₁ = D̄⁻¹ .^ 2, D̄⁻¹
-    for k in 1:maxiter
-        Xₖ₊₁ = rmul!(mul!(Xₖ₊₁, APₖ', Xₖ), Diagonal(D̄⁻¹ₖ))
-        if norm(Xₖ₊₁, Inf) < degeneracy_atol
-            break
-        end
-        Xₖ₊₁ .+= Xₖ
-        if k == maxiter
-            @warn "Sylvester iteration did not converge after $k iterations, final norm of X: $(norm(Xₖ₊₁, Inf)))"
-            break
-        end
-        D̄⁻¹ₖ₊₁ .= D̄⁻¹ₖ .^ 2
-        APₖ₊₁ = mul!(APₖ₊₁, APₖ, APₖ)
-        Xₖ, Xₖ₊₁ = Xₖ₊₁, Xₖ
-        APₖ, APₖ₊₁ = APₖ₊₁, APₖ
-        D̄⁻¹ₖ, D̄⁻¹ₖ₊₁ = D̄⁻¹ₖ₊₁, D̄⁻¹ₖ
-    end
-    Z .+= Xₖ
+    # Normalize by the smallest |eigenvalue|, which caps `max|D̄⁻¹|` at 1, so squaring can
+    # only shrink it.
+    dabsmin = minimum(abs, D)
+    APᴴ ./= dabsmin
+    D̄⁻¹ = dabsmin ./ conj.(D)
+    X = accelerative_smith_iteration!(X₀, similar(X₀), APᴴ, D̄⁻¹, degeneracy_atol, maxiter)
+    Z .+= X
     if eltype(ΔA) <: Real
-        ΔAc = mul!(AP, Z, V') # recycle AP
+        ΔAc = mul!(APᴴ, Z, V') # recycle APᴴ
         ΔA .+= real.(ΔAc)
     else
         ΔA = mul!(ΔA, Z, V', 1, 1)

From d0455b9130b6b02f95379d50949e48aedd09c5ba Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Fri, 2 Oct 2026 19:04:12 +0200
Subject: [PATCH 2/2] Remove superseded outer normalization

---
 src/pullbacks/eig.jl | 7 +------
 1 file changed, 1 insertion(+), 6 deletions(-)

diff --git a/src/pullbacks/eig.jl b/src/pullbacks/eig.jl
index 2e48d5f85..c865d88c0 100755
--- a/src/pullbacks/eig.jl
+++ b/src/pullbacks/eig.jl
@@ -165,12 +165,7 @@ function eig_trunc_pullback!(
     APᴴ = mul!(complex.(A'), ViG, (V * Dmat)', -1, 1)
     X₀ = iszerotangent(ΔV₊) ? APᴴ * Z : mul!(ΔV₊, APᴴ, Z, 1, 1)
     X₀ ./= D'
-    # Normalize by the smallest |eigenvalue|, which caps `max|D̄⁻¹|` at 1, so squaring can
-    # only shrink it.
-    dabsmin = minimum(abs, D)
-    APᴴ ./= dabsmin
-    D̄⁻¹ = dabsmin ./ conj.(D)
-    X = accelerative_smith_iteration!(X₀, similar(X₀), APᴴ, D̄⁻¹, degeneracy_atol, maxiter)
+    X = accelerative_smith_iteration!(X₀, similar(X₀), APᴴ, inv.(conj.(D)), degeneracy_atol, maxiter)
     Z .+= X
     if eltype(ΔA) <: Real
         ΔAc = mul!(APᴴ, Z, V') # recycle APᴴ