1 files changed, 53 insertions, 0 deletions
diff --git a/experiments/verify_theory.py b/experiments/verify_theory.py
index 7ee87bc..ba07487 100644
--- a/experiments/verify_theory.py
+++ b/experiments/verify_theory.py
@@ -108,6 +108,58 @@ def check_predictor_timescale():
             assert abs(empirical - theoretical) < 1e-12
 
 
+def check_conditional_projection():
+    """Finite-sample L2 projection realizes the innovation identities exactly."""
+    rng = np.random.default_rng(351)
+    examples = 4096
+    h = rng.uniform(-1.5, 1.5, size=(examples, 2))
+    basis = np.column_stack((
+        np.ones(examples), h[:, 0], h[:, 1], np.square(h[:, 0]),
+        np.square(h[:, 1]), h[:, 0] * h[:, 1], np.sin(h[:, 0]),
+        np.cos(h[:, 1])))
+    coefficients = rng.normal(size=(basis.shape[1], 3))
+    conditional_mean = basis @ coefficients
+    innovation = rng.normal(scale=0.4, size=conditional_mean.shape)
+    # Project the finite-sample noise off every soma-measurable basis vector.
+    innovation -= basis @ np.linalg.lstsq(basis, innovation, rcond=None)[0]
+    nuisance = conditional_mean + innovation
+
+    linear_basis = basis[:, :3]
+    restricted = linear_basis @ np.linalg.lstsq(
+        linear_basis, nuisance, rcond=None)[0]
+    lhs = np.square(nuisance - restricted).mean()
+    irreducible = np.square(innovation).mean()
+    approximation = np.square(conditional_mean - restricted).mean()
+    orthogonality = np.abs(basis.T @ innovation / examples).max()
+
+    teaching = rng.normal(size=conditional_mean.shape)
+    raw = teaching + nuisance
+    residual = teaching + innovation
+    alpha = (np.linalg.norm(residual, axis=1)
+             / np.linalg.norm(raw, axis=1).clip(min=1e-12))
+    matched = alpha[:, None] * raw
+
+    def row_cosine(left, right):
+        numerator = np.sum(left * right, axis=1)
+        denominator = (np.linalg.norm(left, axis=1)
+                       * np.linalg.norm(right, axis=1)).clip(min=1e-12)
+        return numerator / denominator
+
+    direction_difference = np.abs(
+        row_cosine(raw, teaching) - row_cosine(matched, teaching)).max()
+    norm_difference = np.abs(
+        np.linalg.norm(matched, axis=1) - np.linalg.norm(residual, axis=1)).max()
+    print("\nCONDITIONAL INNOVATION PROJECTION")
+    print(f"orthogonality={orthogonality:.3e} pythagorean_error="
+          f"{abs(lhs - irreducible - approximation):.3e}")
+    print(f"norm_match_error={norm_difference:.3e} "
+          f"direction_change={direction_difference:.3e}")
+    assert orthogonality < 2e-14
+    assert abs(lhs - irreducible - approximation) < 2e-14
+    assert norm_difference < 2e-14
+    assert direction_difference < 2e-14
+
+
 def squared_cosine(x, y):
     return float((x.ravel() @ y.ravel()) ** 2
                  / ((x.ravel() @ x.ravel()) * (y.ravel() @ y.ravel())))
@@ -148,6 +200,7 @@ def main():
     check_simultaneous_variance()
     check_sigma_bias()
     check_descent_threshold()
+    check_conditional_projection()
     check_predictor_timescale()
     check_innovation_identification()
     print("\nALL THEORY CHECKS PASSED")