Cone Bi-Rayleigh Stability Regularizer

Integrate this into the training step of a tanh RNN, state-space model, or implicit residual layer with hidden update \(h_{t+1}=F_\theta(h_t,x_t)\). The exact integration point is the local transition Jacobian \(A_\theta=\partial F_\theta/\partial h\), accessed through Jacobian-vector products rather than materialized. Initialize unconstrained probe parameters \(a,b\), set \(u=\operatorname{softplus}(a)+\delta\) and \(v=\operatorname{softplus}(b)+\delta\), and normalize both by their Euclidean norms. Compute \(Au\) using a forward-mode JVP and \(A^\top v\) using a reverse-mode VJP. Evaluate \(n=\langle u-Au,v\rangle\), \(d=\langle u,v\rangle+\varepsilon\), and \(\lambda=n/d\). Perform 3-10 alternating probe updates minimizing the right and left residuals \(r_R=\|u-Au-\lambda u\|_2^2\) and \(r_L=\|v-A^\top v-\lambda v\|_2^2\), then detach the probes before the outer model update. Add \(\rho\operatorname{softplus}(\mu-\mu_{\max})^2\) to the task loss, with \(\mu=1-\lambda\). Reject probes when \(d<10\varepsilon\). The quotient and residual equations come from the paper; softplus, thresholds, probe refresh rate, and penalty weight are empirical. First test on adding-task and copy-memory benchmarks using a 128-unit tanh RNN. Compare vanilla training, spectral-norm regularization, and this method at matched JVP/VJP cost. Track long-horizon loss, exploding-gradient frequency, residual norms, and wall-clock overhead. Success is lower long-horizon error and fewer unstable runs at equal parameter count; failure is no residual convergence or no improvement over spectral normalization.