diff options
Diffstat (limited to 'ep_run/eig_control.py')
| -rw-r--r-- | ep_run/eig_control.py | 121 |
1 files changed, 109 insertions, 12 deletions
diff --git a/ep_run/eig_control.py b/ep_run/eig_control.py index 6e3f598..0e62068 100644 --- a/ep_run/eig_control.py +++ b/ep_run/eig_control.py @@ -13,30 +13,127 @@ Power iteration on the SYMMETRIC PART -> matvec-only (jvp+vjp of nc_force, the s 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): - """Power-iterate Sym(J_nc)=(J_nc+J_nc^T)/2 at zs for the leading eigenpair. Returns (v_detached, lambda_float). - lambda = v^T Sym(J_nc) v = v^T J_nc v (Rayleigh quotient at the leading eigenvector) = numerical abscissa.""" + """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() - v = cache.get('v') - if v is None or v.shape != z.shape or v.dtype != z.dtype or v.device != z.device: - v = torch.randn_like(z) - v = v / (v.norm() + 1e-12) + + 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(): - for _ in range(iters): - Sv = 0.5 * (jvp(blk.nc_force, z, v)[1] + vjp(blk.nc_force, z, v)[1]) # Sym(J_nc) v - v = Sv / (Sv.norm() + 1e-12) - lam = float((v * jvp(blk.nc_force, z, v)[1]).sum() / (v * v).sum()) # v^T J_nc v - cache['v'] = v + 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): - """Grads of the one-sided leading-abscissa penalty R = eigreg * relu(omega(J_nc) - margin)^2. + """[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 |