theory: verify dynamic neutral projection stability

1 files changed, 49 insertions, 0 deletions

diff --git a/experiments/verify_theory.py b/experiments/verify_theory.py
index d4325e9..21a7c06 100644
--- a/experiments/verify_theory.py
+++ b/experiments/verify_theory.py
@@ -239,6 +239,54 @@ def check_residual_coupling_instability():
     assert too_negative_radius > 1.0
 
 
+def check_dynamic_neutral_projection():
+    """A paired neutral fit nulls affine coupling across covariance scales."""
+    rng = np.random.default_rng(571)
+    examples = 128
+    cells = 64
+    soma = rng.normal(size=(examples, cells))
+    coupling = rng.normal(scale=0.04, size=cells)
+    offset = rng.normal(scale=0.01, size=cells)
+    neutral_residual = soma * coupling + offset
+    centered_soma = soma - soma.mean(axis=0)
+    centered_residual = neutral_residual - neutral_residual.mean(axis=0)
+    variance = np.square(centered_soma).mean(axis=0)
+    correction = (centered_soma * centered_residual).mean(axis=0) / variance
+    remainder = centered_residual - correction * centered_soma
+    remainder_ratio = (np.linalg.norm(remainder)
+                       / np.linalg.norm(neutral_residual))
+    post_slope = ((centered_soma * remainder).mean(axis=0) / variance)
+    maximum_post_slope = float(np.abs(post_slope).max())
+
+    eta = 0.1
+    momentum = 0.9
+    decay = 1e-4
+
+    def radius(k):
+        transition = np.asarray((
+            (1.0 + eta * k, eta * momentum),
+            (k, momentum)))
+        return float(np.abs(np.linalg.eigvals(transition)).max())
+
+    # With zero residual coupling, covariance drops out and only decay remains.
+    null_radii = [radius(-decay) for _ in (1e-3, 1.0, 1e3, 1e6)]
+    # A fixed negative coefficient eventually crosses the lower Jury boundary
+    # as the nonnegative covariance eigenvalue grows.
+    fixed_coefficient = -0.03
+    high_covariance_k = fixed_coefficient * 1e6 - decay
+    high_covariance_radius = radius(high_covariance_k)
+    print("\nDYNAMIC NEUTRAL PROJECTION")
+    print(f"remainder_ratio={remainder_ratio:.3e} "
+          f"post_slope={maximum_post_slope:.3e} "
+          f"rho_null={max(null_radii):.9f} "
+          f"rho_fixed_high_cov={high_covariance_radius:.3f}")
+    assert remainder_ratio < 2e-15
+    assert maximum_post_slope < 2e-16
+    assert max(null_radii) < 1.0
+    assert high_covariance_k < -2.0 * (1.0 + momentum) / eta
+    assert high_covariance_radius > 1.0
+
+
 def check_intermittent_feedback_tracking():
     eta_m = 0.1
     cadence = 16
@@ -308,6 +356,7 @@ def main():
     check_predictor_timescale()
     check_innovation_identification()
     check_residual_coupling_instability()
+    check_dynamic_neutral_projection()
     check_intermittent_feedback_tracking()
     check_kolen_pollack_difference_dynamics()
     print("\nALL THEORY CHECKS PASSED")