theory: connect local bias removal to CLLN scaling

1 files changed, 66 insertions, 0 deletions

diff --git a/experiments/verify_theory.py b/experiments/verify_theory.py
index c00a091..c047193 100644
--- a/experiments/verify_theory.py
+++ b/experiments/verify_theory.py
@@ -160,6 +160,71 @@ def check_conditional_projection():
     assert direction_difference < 2e-14
 
 
+def check_state_dependent_component_bias():
+    """Local conditional subtraction beats a constant and scales by edges."""
+    rng = np.random.default_rng(371)
+    levels = np.asarray((-1.5, -0.5, 0.5, 1.5))
+    repeats = 2048
+    state = np.repeat(levels, repeats)
+    predictable = 0.7 * state + 0.35 * (np.square(state) - 1.25)
+    unpredictable = rng.normal(scale=0.3, size=state.size)
+    # Make the finite-sample remainder orthogonal to every state cell, so the
+    # conditional-expectation identity is checked to floating-point accuracy.
+    for level in levels:
+        group = state == level
+        unpredictable[group] -= unpredictable[group].mean()
+    bias = predictable + unpredictable
+    constant = np.full_like(bias, bias.mean())
+    conditional = np.empty_like(bias)
+    for level in levels:
+        group = state == level
+        conditional[group] = bias[group].mean()
+
+    static_mse = np.square(bias - constant).mean()
+    conditional_mse = np.square(bias - conditional).mean()
+    predictable_mse = np.square(conditional - constant).mean()
+    identity_error = abs(static_mse - conditional_mse - predictable_mse)
+
+    base_residual = bias - constant
+    edge_counts = np.asarray((32, 128, 512, 2048))
+    aggregate_power = np.asarray([
+        np.square(np.resize(base_residual, edges)).sum()
+        for edges in edge_counts
+    ])
+    # Use exact periodic tiling for the scaling identity.
+    periodic = base_residual[:32]
+    periodic_power = np.asarray([
+        np.square(np.tile(periodic, edges // periodic.size)).sum()
+        for edges in edge_counts
+    ])
+    normalized_power = periodic_power / edge_counts
+    scaling_error = float(np.ptp(normalized_power))
+
+    hessian = np.diag((0.4, 0.9, 1.7, 3.2))
+    delta = np.asarray((0.2, -0.1, 0.3, 0.15))
+    displacement = np.linalg.solve(hessian, delta)
+    gradient_at_fixed_point = hessian @ displacement
+    excess = 0.5 * displacement @ hessian @ displacement
+    predicted_excess = 0.5 * delta @ np.linalg.solve(hessian, delta)
+    lower = delta @ delta / (2.0 * np.linalg.eigvalsh(hessian).max())
+    upper = delta @ delta / (2.0 * np.linalg.eigvalsh(hessian).min())
+
+    print("\nSTATE-DEPENDENT COMPONENT BIAS")
+    print(f"static_mse={static_mse:.6f} conditional_mse={conditional_mse:.6f} "
+          f"removed={predictable_mse:.6f} identity_error={identity_error:.3e}")
+    print(f"edge_power_per_component={normalized_power[0]:.6f} "
+          f"scaling_error={scaling_error:.3e}")
+    print(f"fixed_point_error={np.abs(gradient_at_fixed_point-delta).max():.3e} "
+          f"excess={excess:.6f} bounds=[{lower:.6f},{upper:.6f}]")
+    assert conditional_mse < static_mse
+    assert identity_error < 2e-15
+    assert aggregate_power.shape == edge_counts.shape
+    assert scaling_error < 2e-15
+    assert np.abs(gradient_at_fixed_point - delta).max() < 2e-15
+    assert abs(excess - predicted_excess) < 2e-15
+    assert lower <= excess <= upper
+
+
 def squared_cosine(x, y):
     return float((x.ravel() @ y.ravel()) ** 2
                  / ((x.ravel() @ x.ravel()) * (y.ravel() @ y.ravel())))
@@ -393,6 +458,7 @@ def main():
     check_sigma_bias()
     check_descent_threshold()
     check_conditional_projection()
+    check_state_dependent_component_bias()
     check_predictor_timescale()
     check_innovation_identification()
     check_residual_coupling_instability()