theory: expose residual coupling instability

1 files changed, 39 insertions, 0 deletions

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")