Tearing - BUGFIX! - Polish roots with a real 2x2 Jacobian so the anti-holomorphic SLAYER residual converges

src/Tearing/Dispersion/GrowthRateExtraction.jl

@@ -112,7 +112,10 @@ single-surface scans; `mc.surfaces[mc.ref_idx].tauk` for coupled scans).
     reports the raw marching-squares estimate, which is sensitive to grid and
     floating-point rounding.
   - `polish_maxit` -- max polish iterations per root (default 20). Each costs a
-    handful of `f` evaluations; failures fall back to the unpolished point.
+    handful of `f` evaluations (two for the real 2×2 Jacobian, one trial step,
+    at most three backtracks); failures fall back to the unpolished point. The
+    polish does not assume `f` is holomorphic, so it also converges on the SLAYER
+    residual, which is a function of conj(Q).
   - `validity_rtol` -- with `residual` supplied, a polished root is kept only if
     `|residual| < validity_rtol · median(|Δ|)` (default 1e-3). Intersections
     that don't polish to a true zero (contour near-misses, or surfaces whose
@@ -500,9 +503,12 @@ end
 
 # Bounded Gauss-Newton minimization of |f(Q)| within a disk of radius `R`
 # around `Q0`. Returns (Q_best, n_evals, improved). `improved == false` ⇒
-# Q0 is returned unchanged (hard fallback). Secant derivative + trust-region
-# projection + |f|-decrease backtracking keep it strictly "improve-or-no-op":
-# it cannot diverge, jump to another root, or be pulled into a pole.
+# Q0 is returned unchanged (hard fallback). Trust-region projection +
+# |f|-decrease backtracking keep it strictly "improve-or-no-op": it cannot
+# diverge, jump to another root, or be pulled into a pole.
+#
+# The local model is the real 2×2 forward-difference Jacobian of (Re f, Im f), exact for holomorphic
+# and anti-holomorphic f alike; SLAYER's residual is a function of conj(Q) (see Riccati.jl).
 function _polish_root(f, Q0::ComplexF64, R::Float64; maxit::Int=8,
     tol_step::Float64=1e-12, tol_f::Float64=1e-8)
     (R > 0) || return (Q0, 0, false)
@@ -512,17 +518,21 @@ function _polish_root(f, Q0::ComplexF64, R::Float64; maxit::Int=8,
     a0 = abs(f0)
     a0 == 0 && return (Q0, nev, false)
 
-    # Second point for the secant slope: a small real offset scaled to R.
+    h = max(R * 1e-3, eps(Float64))                       # finite-difference offset, scaled to R
     Qk, fk = Q0, f0
-    Qp = Q0 + complex(max(R * 1e-3, eps(Float64)), 0.0)
-    fp = ComplexF64(f(Qp))
-    nev += 1
     Qbest, abest = Q0, a0
 
     for _ in 1:maxit
-        df = (Qp == Qk) ? ComplexF64(0) : (fp - fk) / (Qp - Qk)
-        (isfinite(df) && abs(df) > 0) || break
-        step = -fk / df                                   # zero of local linear model
+        fx = ComplexF64(f(Qk + h))
+        fy = ComplexF64(f(Qk + im * h))
+        nev += 2
+        (isfinite(fx) && isfinite(fy)) || break
+        j11, j12 = real(fx - fk) / h, real(fy - fk) / h
+        j21, j22 = imag(fx - fk) / h, imag(fy - fk) / h
+        jdet = j11 * j22 - j12 * j21
+        (isfinite(jdet) && jdet != 0) || break
+        # zero of the local linear model: J * [dRe, dIm] = -[Re f, Im f]
+        step = complex((-real(fk) * j22 + j12 * imag(fk)) / jdet, (-j11 * imag(fk) + j21 * real(fk)) / jdet)
         Qt = Qk + step
         if abs(Qt - Q0) > R                               # project into trust region
             Qt = Q0 + R * (Qt - Q0) / abs(Qt - Q0)
@@ -538,7 +548,6 @@ function _polish_root(f, Q0::ComplexF64, R::Float64; maxit::Int=8,
             bt += 1
         end
         if isfinite(ft) && abs(ft) < abest
-            Qp, fp = Qk, fk
             Qk, fk = Qt, ft
             Qbest, abest = Qt, abs(ft)
             (abs(step) < tol_step * (abs(Qk) + tol_step) || abest < tol_f * a0) && break

test/runtests_dispersion_polish.jl

@@ -14,7 +14,7 @@
         @test improved
         @test abs(Qp - r) < 1e-10
         @test abs(f(Qp)) < 1e-10
-        @test nev <= 3 * 8 + 2             # bounded eval budget
+        @test nev <= 1 + 6                 # one iteration: 2 differences, the trial step, ≤ 3 backtracks
     end
 
     @testset "_polish_root: nonlinear residual converges" begin
@@ -25,6 +25,43 @@
         @test abs(Qp - r) < 1e-7
     end
 
+    @testset "_polish_root: residual that is neither holomorphic nor anti-holomorphic converges" begin
+        r = 0.15 - 0.35im
+        f(Q) = (ComplexF64(Q) - r) + 0.5 * conj(ComplexF64(Q) - r)   # zero only at r
+        Qp, _, improved = D._polish_root(f, r + 0.01 + 0.012im, 0.1)
+        @test improved
+        @test abs(Qp - r) < 1e-10
+    end
+
+    # The SLAYER residual is a function of conj(Q) (solve_inner evaluates
+    # conj(Δ̂(ĝ))), so the polish must not assume a holomorphic f.
+    @testset "_polish_root: anti-holomorphic residual converges" begin
+        r = 0.42 + 0.27im
+        f(Q) = (2 - 1im) * conj(ComplexF64(Q) - r)
+        Qp, _, improved = D._polish_root(f, r + (0.013 - 0.009im), 0.1)
+        @test improved
+        @test abs(Qp - r) < 1e-10
+        @test abs(f(Qp)) < 1e-9
+    end
+
+    @testset "_polish_root: anti-holomorphic root offset along Im Q (mirrored complex step)" begin
+        # A complex-slope step on f = a·conj(Q - r) maps an offset i·y to 2i·y and
+        # every halving backtrack stays worse, so a complex-secant polish stalls here.
+        r = -0.75 - 0.009im
+        f(Q) = 1e3 * conj(ComplexF64(Q) - r)
+        Qp, _, improved = D._polish_root(f, r + 0.0017im, 0.01)
+        @test improved
+        @test abs(Qp - r) < 1e-10
+    end
+
+    @testset "_polish_root: steep nonlinear anti-holomorphic residual converges" begin
+        r = -0.75 - 0.009im
+        f(Q) = 1e3 * conj((ComplexF64(Q) - r) * (1 + 0.4 * (ComplexF64(Q) - r)))
+        Qp, _, improved = D._polish_root(f, r + (0.004 + 0.003im), 0.02)
+        @test improved
+        @test abs(Qp - r) < 1e-9
+    end
+
     @testset "_polish_root: no-op-or-better invariant (never worse)" begin
         # The polish must never return a point with larger |f| than the start.
         for r in (0.1 + 0.2im, -0.5 + 0.0im, 0.3 - 0.4im)
@@ -181,4 +218,17 @@
         @test !(:spurious in gr3.warning_flags)
         @test !isnan(real(gr3.Q_root))
     end
+
+    @testset "find_growth_rates: steep anti-holomorphic root passes the validity gate" begin
+        # A steep residual leaves the raw contour intersection far above the gate
+        # threshold, so the root is kept only if the polish reaches the true zero.
+        r = 0.42 + 0.27im
+        f(Q) = 1e3 * conj((ComplexF64(Q) - r) * (1 + 0.4 * (ComplexF64(Q) - r)))
+        scan = brute_force_scan(f, (0.0, 1.0), (0.0, 0.6); nre=9, nim=7, threaded=false)   # raw |f| above the gate
+        # the pole classifier's threshold is absolute: scale it with f, as SLAYER's adaptive one does
+        gr = find_growth_rates(scan, 1.0; residual=f, validity_rtol=1e-3, pole_threshold=1e3 * 10.0)
+        @test !(:spurious in gr.warning_flags)
+        @test !isnan(real(gr.Q_root))
+        @test abs(gr.Q_root - r) < 1e-8
+    end
 end