theory: expose residual coupling instability

2 files changed, 95 insertions, 0 deletions

diff --git a/THEORY.md b/THEORY.md
index 146621f..8299c75 100644
--- a/THEORY.md
+++ b/THEORY.md
@@ -499,6 +499,61 @@ The relevant timescale is therefore the product of learning rate and number of
 neutral updates. Too few neutral updates leave traffic; task-period updates
 converge to the wrong coefficient regardless of speed.
 
+### Incomplete residualization can create multiplicative local instability
+
+Residual power alone is not a stability certificate. Consider a linearized
+local population `h=W x` with input covariance `C=E[x x^T]`. If the true
+soma-proportional traffic coefficient is `M` and the neutral predictor is `P`,
+write the remaining coupling as `D=M-P`. The innovation used for the local
+update is
+
+```text
+r = s + D h = s + D W x.
+```
+
+Ignoring nonlinear gates for this local calculation, the expected forward
+update with decay is
+
+```text
+E[Delta W] = eta (S + D W C - lambda W),
+S = E[s x^T].
+```
+
+Thus predictor error does not enter as additive noise. It multiplies the
+current forward weights through `D W C`; a positive mode is a Hebbian feedback
+loop. In vector form its homogeneous operator is
+
+```text
+vec(W_(t+1))
+  = [I + eta (C^T kron D - lambda I)] vec(W_t).
+```
+
+The implemented momentum dynamics make the same point sharply. Along a scalar
+joint eigenmode let `k=d c-lambda`, where `d` and `c>=0` are eigenvalues of
+`D` and `C`. With momentum `mu`, the homogeneous state follows
+
+```text
+[w_(t+1)]   [1 + eta k, eta mu] [w_t]
+[m_(t+1)] = [k,          mu   ] [m_t].
+```
+
+The determinant is `mu`. Evaluating its characteristic polynomial at one
+gives `p(1)=-eta k`; therefore every `k>0` produces a real eigenvalue above
+one. Weight decay stabilizes the mode only if it makes the effective `k`
+nonpositive, with the usual additional discrete-time step-size constraints.
+Consequently a small aggregate residual RMS can coexist with an unstable
+worst-direction coefficient. Deep state dependence can amplify the same mode
+across layers before a slow neutral predictor catches it.
+
+This calculation explains the frozen MT-1 failure boundary without changing
+its result. The training-only audit finds the first active nonfinite tensors
+at the stem forward weight and momentum: step 23 for raw, 69 for norm-matched
+raw, and 74 for innovation. Residualization and norm matching delay the
+positive-feedback mode but do not stabilize it at the frozen four-to-one
+traffic setting. A scalable version needs an operator-level dissipativity or
+gain-control mechanism; the pretraining residual-RMS threshold is
+insufficient by itself.
+
 ### Intermittent feedback tracking can hide behind a final cosine
 
 An idealized mirror event every `k` task updates uses
@@ -576,4 +631,5 @@ The verifier checks the exact Rademacher MSE across depth, width, and K; the
 quadratic finite-sigma bias law; the smooth-descent step threshold; the
 conditional-projection Pythagorean identity and norm-matching direction
 invariance; neutral predictor convergence across `eta_P` and update counts;
+the multiplicative residual-coupling operator and its unstable positive mode;
 and the task-fit teaching-signal absorption predicted above.
diff --git a/experiments/verify_theory.py b/experiments/verify_theory.py
index cc9b6c3..ad0e9b9 100644
--- a/experiments/verify_theory.py
+++ b/experiments/verify_theory.py
@@ -196,6 +196,44 @@ def check_innovation_identification():
     assert abs(retained_fraction - (1.0 - predictable_fraction)) < 0.01
 
 
+def check_residual_coupling_instability():
+    rng = np.random.default_rng(551)
+    rows, columns = 4, 3
+    residual_coupling = rng.normal(size=(rows, rows))
+    weight = rng.normal(size=(rows, columns))
+    input_covariance = rng.normal(size=(columns, columns))
+    input_covariance = input_covariance @ input_covariance.T / columns
+    direct = residual_coupling @ weight @ input_covariance
+    operator = np.kron(input_covariance.T, residual_coupling)
+    vectorized = (operator @ weight.ravel(order="F")).reshape(
+        weight.shape, order="F")
+    identity_error = float(np.abs(direct - vectorized).max())
+
+    eta = 0.1
+    momentum = 0.9
+
+    def radius(k):
+        transition = np.asarray((
+            (1.0 + eta * k, eta * momentum),
+            (k, momentum)))
+        return float(np.abs(np.linalg.eigvals(transition)).max())
+
+    positive_k = 0.02
+    negative_k = -0.20
+    positive_radius = radius(positive_k)
+    negative_radius = radius(negative_k)
+    polynomial_at_one = -eta * positive_k
+    print("\nMULTIPLICATIVE RESIDUAL COUPLING")
+    print(f"vectorization_error={identity_error:.3e} "
+          f"rho(k={positive_k:+.3f})={positive_radius:.9f} "
+          f"rho(k={negative_k:+.3f})={negative_radius:.9f} "
+          f"p_positive(1)={polynomial_at_one:+.3e}")
+    assert identity_error < 2e-15
+    assert polynomial_at_one < 0.0
+    assert positive_radius > 1.0
+    assert negative_radius < 1.0
+
+
 def check_intermittent_feedback_tracking():
     eta_m = 0.1
     cadence = 16
@@ -264,6 +302,7 @@ def main():
     check_conditional_projection()
     check_predictor_timescale()
     check_innovation_identification()
+    check_residual_coupling_instability()
     check_intermittent_feedback_tracking()
     check_kolen_pollack_difference_dynamics()
     print("\nALL THEORY CHECKS PASSED")