"""#2 — leading-abscissa control for the ept non-conservative operator (ports the aep-dynamics 'control the LEADING spectral signal, surgically' finding to C512). Why not jacreg: jacreg penalizes ||J_nc||_F^2 (Hutchinson) — the WHOLE Jacobian norm. That is blunt: it over-constrains directions that never threaten stability, and when the controller ramps it high it HIJACKS the task gradient (the known jr-hijack failure). The aep leading-vs-lagging result says the right knob is the leading SPECTRAL ABSCISSA, not the norm. What we control: the numerical abscissa omega(J_nc) = lambda_max( (J_nc + J_nc^T)/2 ) = the 2-norm log-norm mu_2(J_nc) = one-sided Lipschitz constant. It upper-bounds the spectral abscissa, governs the transient growth ||e^{Jt}||, and its crossing past (1+c) IS the free-phase Hopf (J_F = J_nc - (1+c)I). Power iteration on the SYMMETRIC PART -> matvec-only (jvp+vjp of nc_force, the same primitives jacreg and the AEP correction already call), so it scales and is analog-compatible (no eigendecomposition). One-sided ReLU penalty = a LEADING signal: acts only as the abscissa nears the margin, so unlike jacreg it does not over-contract / hijack when the operator is already safe. === VERDICT (eig_recheck.py, 2026-07-03): the omega premise above is REFUTED for this operator. === Gold-standard Lanczos on Sym(J_nc) at z*: omega = +5..+13 on EVERY operator — including the stablest warm source s2000 (omega=+12.8 while its TRUE spectral abscissa, eig_probe FD-JVP, is ~-0.02). The non-normality gap omega - alpha ~ 10: the numerical abscissa is hopelessly conservative here and even ANTI-correlates with actual stability (rand-init +5.0 -> s2000 +12.8). An omega-margin penalty would fire permanently and over-constrain worse than jacreg. KEEP num_abscissa/eig_penalty as DIAGNOSTICS ONLY. (The pre-fix plain power iteration also returned garbage mixtures on the indefinite Sym — the old 'warm -10.14 vs scratch +1.11' fingerprint story is retracted.) The CONTROL is lead_rho/spec_penalty below: the TRUE leading eigenvalue of the forward map M = I + eps*J_F (a real 2-D Rayleigh-Ritz block captures the complex Hopf pair that symmetric power iteration cannot see). rho = |lam_lead(M)| is the exact discrete stability criterion (rho > 1 = the relaxation spirals out — covers both the continuous Hopf and the Euler/Neimark-Sacker route), matches eig_probe's measurement, and is the aep-dynamics 'spectral' steering ported to C512 — applied SOFTLY (small coefficient, one-sided), per the aep finding that hard projection over-constrains. """ import torch from torch.autograd.functional import jvp, vjp def num_abscissa(blk, zs, cache, iters=3): """SHIFTED power iteration for the TOP (algebraic-max) eigenpair of Sym(J_nc)=(J_nc+J_nc^T)/2 at zs. Returns (v_detached, lambda_max_float) — lambda_max(Sym) = the numerical abscissa. Why shifted: plain power iteration converges to the largest-|lambda| eigenpair. Sym(J_nc) is indefinite and its NEGATIVE end can dominate in magnitude — plain PI then silently returns lambda_min, which reads 'deeply stable' exactly when lambda_max is over the margin, disarming the penalty. Phase 1 estimates the dominant magnitude s; phase 2 iterates (Sym + s·I), whose top eigenvalue lambda_max + s is non-negative and dominant, so the Rayleigh quotient converges to lambda_max. Both phases warm-start from the cache, so accuracy accrues across training steps.""" z = zs.detach() def smv(u): # Sym(J_nc) u return 0.5 * (jvp(blk.nc_force, z, u)[1] + vjp(blk.nc_force, z, u)[1]) def unit(key): u = cache.get(key) if u is None or u.shape != z.shape or u.dtype != z.dtype or u.device != z.device: u = torch.randn_like(z) return u / (u.norm() + 1e-12) with torch.no_grad(): u, s = unit('u'), 0.0 for _ in range(iters): # phase 1: dominant |lambda| -> the shift Su = smv(u); s = float(Su.norm()); u = Su / (s + 1e-12) s = 1.1 * s + 1e-3 # safe overestimate (need s >= |lambda_min|) v = unit('v') for _ in range(iters): # phase 2: top eigenpair of Sym + s·I w = smv(v) + s * v v = w / (w.norm() + 1e-12) lam = float((v * jvp(blk.nc_force, z, v)[1]).sum() / (v * v).sum()) # v^T J v = v^T Sym(J) v cache['u'] = u; cache['v'] = v return v, lam def lead_rho(blk, zs, eps, c, cache, iters=6): """TRUE leading eigenvalue of the forward map M = I + eps*J_F, J_F = J_nc - (1+c)I. Warm-started 2-D subspace (Rayleigh-Ritz) iteration: a real 2-D block captures the complex Hopf pair that plain power iteration cannot represent. Matvec-only (jvp of nc_force). Returns (V=[v0,v1] detached orthonormal, rho, mu_re): rho = |lam_lead(M)| (discrete stability: rho > 1 = relaxation spirals out), mu = (lam-1)/eps the continuous eigenvalue (Hopf: Re mu > 0).""" z = zs.detach() k = 1.0 - eps * (1.0 + c) def mv(u): # M u = u + eps*(J_nc - (1+c)I) u return k * u + eps * jvp(blk.nc_force, z, u)[1] V = cache.get('V') if V is None or V[0].shape != z.shape or V[0].dtype != z.dtype or V[0].device != z.device: V = [torch.randn_like(z), torch.randn_like(z)] with torch.no_grad(): for _ in range(iters): W0, W1 = mv(V[0]), mv(V[1]) W0 = W0 / (W0.norm() + 1e-12) # Gram-Schmidt: orthonormal 2-D block W1 = W1 - (W1 * W0).sum() * W0 W1 = W1 / (W1.norm() + 1e-12) V = [W0, W1] M0, M1 = mv(V[0]), mv(V[1]) h00, h01 = float((V[0] * M0).sum()), float((V[0] * M1).sum()) # H = V^T M V (Rayleigh-Ritz) h10, h11 = float((V[1] * M0).sum()), float((V[1] * M1).sum()) tr, det = h00 + h11, h00 * h11 - h01 * h10 disc = tr * tr / 4.0 - det if disc >= 0: # two real Ritz values: take the larger |.| l1, l2 = tr / 2.0 + disc ** 0.5, tr / 2.0 - disc ** 0.5 lam_re, rho = (l1, abs(l1)) if abs(l1) >= abs(l2) else (l2, abs(l2)) else: # complex pair a+-bi: |lam|=sqrt(det), Re=tr/2 rho, lam_re = det ** 0.5, tr / 2.0 cache['V'] = V return V, rho, (lam_re - 1.0) / eps def spec_penalty(blk, zs, eps, c, specreg, rho_target, cache, iters=6): """Grads of the one-sided TRUE-leading-eigenvalue penalty R = specreg * relu(rho - rho_target)^2 with rho = |lam_lead(I + eps*J_F)| via Rayleigh-Ritz on the converged (detached) 2-D subspace. This is the aep-dynamics 'spectral' steering at C512, applied softly (one-sided, small specreg). Returns ({id(p): grad}, rho, mu_re) — empty grads when rho below target.""" V, rho0, mu_re = lead_rho(blk, zs, eps, c, cache, iters) if rho0 <= rho_target: # contraction margin held: signal off return {}, rho0, mu_re z = zs.detach() k = 1.0 - eps * (1.0 + c) with torch.enable_grad(): # H(theta), V fixed (envelope/Rayleigh-Ritz) M0 = k * V[0] + eps * jvp(blk.nc_force, z, V[0], create_graph=True)[1] M1 = k * V[1] + eps * jvp(blk.nc_force, z, V[1], create_graph=True)[1] h00, h01 = (V[0] * M0).sum(), (V[0] * M1).sum() h10, h11 = (V[1] * M0).sum(), (V[1] * M1).sum() tr, det = h00 + h11, h00 * h11 - h01 * h10 disc = tr * tr / 4.0 - det rho = torch.where(disc >= 0, torch.maximum((tr / 2 + disc.clamp(min=0).sqrt()).abs(), (tr / 2 - disc.clamp(min=0).sqrt()).abs()), det.clamp(min=1e-12).sqrt()) R = specreg * torch.relu(rho - rho_target) ** 2 gr = torch.autograd.grad(R, blk.block, allow_unused=True) return {id(p): g for p, g in zip(blk.block, gr) if g is not None}, rho0, mu_re def eig_penalty(blk, zs, eigreg, margin, cache, iters=3): """[DIAGNOSTIC-ONLY — omega control refuted, see module docstring; use spec_penalty instead.] Grads of the one-sided leading-abscissa penalty R = eigreg * relu(omega(J_nc) - margin)^2. Returns ({id(p): grad}, omega) — omega logged as the controller signal. Empty grads when below margin.""" v, lam0 = num_abscissa(blk, zs, cache, iters) if lam0 <= margin: # below the stability margin: leading signal off return {}, lam0 z = zs.detach() with torch.enable_grad(): Jv = jvp(blk.nc_force, z, v, create_graph=True)[1] # differentiable in theta (nc_force params) lam = (v * Jv).sum() / (v * v).sum() # numerical abscissa, v fixed R = eigreg * torch.relu(lam - margin) ** 2 gr = torch.autograd.grad(R, blk.block, allow_unused=True) return {id(p): g for p, g in zip(blk.block, gr) if g is not None}, lam0