Conditioned Numerical-Range Stability Regularizer
Implementation & benchmark of arXiv:2608.09866 — Sharp spectral constants for scaled $q$-numerical ranges
Source paper: Sharp spectral constants for scaled $q$-numerical ranges arXiv:2608.09866 ⓘ · analyzed Sep 1, 2026
AI-generated research hypothesis, automatically tested. Not peer-reviewed.
Idea description
Regularize a recurrent or state-space transition matrix using numerical ranges after bounded-condition-number similarity transforms, rather than only penalizing eigenvalues or the raw spectral norm. The resulting penalty targets nonnormal transient amplification and can certify bounds on powers or other polynomial functions of the transition matrix.
Formulas
Mathematical statement
For a matrix A in M_n(C), define the numerical range W(A)={<Ax,x>: ||x||=1}. The scaled q-numerical range is Ω_q(A)=q^{-1}W_q(A), where W_q(A)={⟨Ax,y⟩: ||x||=||y||=1, ⟨x,y⟩=q} for 0<|q|≤1. The paper proves the exact identity Ω_{η(γ)}(A)=⋃_{κ(S)≤γ}W(S^{-1}AS), with η(γ)=2/(γ+γ^{-1}) and κ(S)=||S||·||S^{-1}||. Thus a q-range samples the numerical ranges of every coordinate transform S whose condition number is at most γ. It also proves ||p(A)||≤C(q) max_{z∈Ω_q(A)}|p(z)| for every polynomial p, where C(q)=max{1, 2|q|/(1+√(1-|q|²))}. For p(z)=z^k, this gives ||A^k||≤C(q)(max_{z∈Ω_q(A)}|z|)^k. The adaptation estimates the maximum transformed numerical-range radius R_γ(A)=max_{κ(S)≤γ} max_{||x||=1}|⟨S^{-1}ASx,x⟩| through sampled invertible transforms, and penalizes R_γ(A) or its certified power-growth surrogate.
Implementation notes
(1) Integrate this into the recurrent/state-space transition matrix, not into the input/output projections. For an RNN use the hidden update h_{t+1}=\phi(Ah_t+Bu_t); for an SSM use the discretized state matrix A_d. Apply the regularizer to A or A_d once per optimizer step, with target radius ρ<1 for contractive dynamics or ρ≈1 for long-memory dynamics. The explicit theorem parameter is γ≥1; set q=η(γ)=2/(γ+γ^{-1}) and C(q) from the paper. (2) Use the following differentiable Monte Carlo estimator. Sample m random Gaussian matrices G_j, normalize their singular values so that S_j=U diag(s) V^* has condition number γ exactly or at most γ, and sample r random unit vectors x_{j,r}. Compute z_{j,r}=x_{j,r}^* S_j^{-1} A S_j x_{j,r}, then R_hat=max_{j,r}|z_{j,r}|, replacing max by τ log sum exp(|z|/τ) during training. Pseudocode: R=0; for j in 1..m: S=make_bounded_condition_transform(G_j,gamma); M=solve(S,A@S); for r in 1..r0: x=normalize(randn()); R=max_or_softmax(R, abs(x.H@M@x)); loss += beta*softplus(R-rho). (3) The identity Ω_{η(γ)}=union W(S^{-1}AS) and the polynomial constant C(q) come from the paper; the sampled R_hat, random transforms, and empirical radius are approximations. Estimate the actual transient envelope independently by computing ||A^k|| for k=1,...,K and compare it with C(q)R_hat^k. Use power iteration or exact SVD on small matrices to avoid confusing numerical-range shrinkage with true stability. (4) First experiment: compare vanilla spectral-radius clipping, raw spectral-norm regularization, and this penalty on a 1-layer tanh RNN or diagonal-plus-low-rank SSM trained on sequential MNIST or adding-task. Use hidden size 64, γ∈{2,5,10}, m=8 transforms, 16 probe vectors, and K=50. Measure gradient explosion frequency, maximum ||A^k||, validation accuracy, and training loss at equal FLOPs. Success is lower transient growth and fewer exploding gradients at matched accuracy, or improved long-sequence accuracy at matched one-step spectral radius. After validation, replace random transforms with adversarial gradient ascent on S to obtain a stronger worst-case regularizer.
Verification
This idea has not been verified yet.
Verification happens in two stages: Stage 1 — a mechanism check on a toy system confirms the claimed mathematical phenomenon reproduces; Stage 2 — a benchmark implements the idea on a real (small) neural network task and compares it against a tuned baseline over 8 paired seeds with a permutation test.
Artifacts
Artifacts unavailable.