Added Jacobian Stabilization as proposed by Eßl et al.

[Julpe]

Jacobian stabilization of the self-consistency loop

Adds a Jacobian-based stabilization of the ladder-DGA self-consistency cycle [1-3], behind the flag stabilization.use_jacobian_stabilization. The leading part of the spectrum of the Sigma map is estimated from the (iterate, proposal) pairs the mixing records anyway, which costs no extra evaluations of the map. The opt-in stabilization.use_exact_jacobian adds the exact Jacobian of the map: it certifies inside the loop the modes the estimate cannot, and replaces the estimate at the converged point.

What's new

• In-loop Jacobian tracker (dgamore/jacobian_stabilization.py):
  • Spectrum estimate: a secant Rayleigh-Ritz estimate of the leading eigenvalues lambda_Pi = 1 - lambda_J of the Sigma map, from the last pairs, refreshed every iteration. Only estimates whose residual passes a gate count as certified, and only certified ones drive decisions. Past the gate every value carries an error bound, the residual times the eigenvalue's condition number, so a non-normal map is not flipped on a backward error alone.
  • Reflection of unstable directions: on a direction certified with Re lambda_Pi < 0, the proposal residual is reflected before the mixing acts, from the iteration that certifies it on, which makes the physical fixed point attractive there [1, 2]. The projector onto the certified span is built in Schur form [3]. A reflection is released when a certified stable mode lies inside it or when the residual keeps growing, so a correct flip whose mode has contracted out of the estimate stays in place.
  • Predicted and pole crossings: a direction whose real part falls toward the stability boundary over consecutive matched estimates is flipped ahead of its crossing. One that jumps through a pole is flipped on the iteration of the jump [1, 2].
  • Damping bound: damped steps use the largest damping the certified spectrum allows [2], with the safety factor at the vertex of the stability parabola [3] and a floor of 0.01. Anderson and Pulay steps keep the configured mixing; while a reflection is installed, an accelerated step that points against the flipped damped step on the reflected subspace is replaced by the stabilized damped step, since a quasi-Newton step can converge to the unphysical fixed point [1, 2].
  • Per-direction damping: on damped steps, each certified direction takes the step its own eigenvalue allows [2, 3]. Blocks coupled too strongly to be damped apart share the smallest damping of their group.
  • Convergence and minimum iterations: the step residual of a damped step is read at the configured mixing, so the verdict means the same distance from the fixed point with and without the flag. The loop does not declare convergence before the minimum iteration count of the certified spectrum has passed [3], and warns when a mode inside its undecidable band asks for more.
  • Carry-over along a temperature ladder: the certified spectrum of the previous rung is re-gridded to the new frequencies and installed at the start (on another momentum grid or symmetry reduction, its eigenvalues and damping bound carry without vectors). Modes that the two preceding rungs extrapolate across the boundary are flipped preemptively [2, 3], and a mode that approached the boundary and receded, the cusp of [1], is reported. The saved spectrum keeps the flips the reflection in force still holds.
  • Symmetry sector: the tracker keeps its vectors on the weighted symmetry sector of the self-energy (irreducible momenta and positive frequencies, weighted so every norm is the one of the whole window), so its memory scales with the irreducible wedge: about 7 GB instead of 547 GB on rank 0 at 50^3 with 4 bands and 150 frequencies.
  • One-shot runs: the flag is disabled with a warning for max_iter: 1 as well as for the one-shot lambda correction, since a single iteration has no history.
• Exact Jacobian (dgamore/sigma_jacobian.py, flag stabilization.use_exact_jacobian): analytic Jacobian-vector products of the loop's map, with the pipeline's own functions run on linearized objects, chemical-potential and occupation feedback included (exact for V = 0; the V^q Hartree-Fock offset of the shell is held, and the loop warns about it). They act on the self-energies the loop can reach (the little-group average of every irreducible momentum and the populated orbital pairs, so multi-orbital and auto-symmetry lattices are covered) in the tracker's own coordinates:
  • Exact checks inside the loop: a least-stable mode that two consecutive estimates match below Re lambda_Pi = 0.1 without certifying it, and the near-unstable modes a carried exact spectrum hands over, are certified by a short block Arnoldi (at most ten products) before the iteration's step. Certified pairs act at once, and a negative real part is flipped. A run makes at most five such checks besides the one at its start; a check that certifies its lead and finds every certified pair stable gives its slot back, at most five times, so the budget is kept for a mode that turns out unstable.
  • Converged-point spectrum: at the pure fixed point, ARPACK computes the six eigenpairs with the largest real part and the six with the largest modulus, which replace the secant estimate in jacobian.npz for the next rung [3]. Both searches start from the tracker's least-stable and stiffest directions and share the products of their first Arnoldi factorization (41 products fewer per solve). A warning flags an unstable exact mode the reflection does not hold, and a search whose leading pairs are all unstable.
  • Block products: the Bethe-Salpeter systems of a block of products share one factorization per channel (invert_and_sum_over_last_vn_v2 takes a list of right-hand sides). The checks run their products in blocks whose width the memory detection sizes, at most four.
• Monitors:
  • the first-frequency density pole ratio R: a pole ring of the density ladder at w_1 lies inside the Brillouin zone once R >= 1
  • static and w_1 compound checks of chi_phys
  • the Eliashberg even-sector degeneracy hint (agreement within 1 %; expected for an SU(2N)-symmetric interaction)
  • the current filling, logged every iteration
  • the map residual |S(Sigma) - Sigma| / |Sigma| beside the step residual, with a stalled-step warning when the loop converges on the step while the map residual stays above 10x epsilon without falling
• Releasing scaffolds:
  • reworked use_chi_phys_restriction: a floor on the static inverse, and finite-frequency values clipped to the static maximum
  • a per-iteration lambda correction with a warm-started lambda, searched on the branch that keeps chi_lambda positive at every frequency
  • both release at 10x epsilon, so the result is always the pure fixed point
• Output: per-iteration iterates and the core window of the raw proposals in Sigma_Iterates/; the tracker's spectrum, certified modes and eigenvalue traces in jacobian.npz (uncompressed, written through a temporary file so a killed run leaves the previous version), which holds the exact spectrum (key exact) when the exact Jacobian ran.
• Memory estimate:
  • the tracker's history triples and its resident Ritz and reflector arrays count in every rank-0 slot, sized on the symmetry sector and widened by the exact checks' sets when the exact Jacobian runs
  • its secant transient counts in the mixing step
  • a branch for the exact Jacobian, including its block width; the flag is disabled with a warning when it does not fit
  • the pairing-vertex transient uses the per-rank peak
  • the k-resolved occupations are summed from Sigma in momentum chunks, without a full-box Green's function; whether they keep an imaginary part is decided on the whole dispersion, so the result does not depend on the rank count
• Garbage collection: deferred_collection is re-entrant per thread, so the products of the exact Jacobian, which wrap the pipeline's own deferred loops, no longer sweep the heap on every release.
• Docs: a new cooldown section on Jacobian tracking and the exact Jacobian; the configuration and output references are updated.

References

[1] H. Eßl, M. Reitner, E. Kozik, A. Toschi, Phys. Rev. Lett. 137, 016502 (2026).
[2] H. Eßl, S. Rohshap, M. Gievers, M. Wallerberger, A. Toschi, A. Kauch, Stabilizing the parquet problem, arXiv:2606.04936.
[3] M. Gievers, H. Eßl, S. Rohshap, A. Kauch, A. Toschi, Instabilities in self-consistent diagrammatic approaches and how to cure them, arXiv:2609.11405.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.

📢 Thoughts on this report? Let us know!