theory: connect local bias removal to CLLN scaling

2 files changed, 138 insertions, 0 deletions

diff --git a/THEORY.md b/THEORY.md
index 4aecc96..112c7eb 100644
--- a/THEORY.md
+++ b/THEORY.md
@@ -499,6 +499,78 @@ 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.
 
+### State-dependent component bias accumulates with system size
+
+The physical-network setting makes the nuisance-removal statement concrete.
+For component `e`, let `z_e` denote locally observable state, `g_e` the ideal
+credit signal, and let the measured local signal be
+
+```text
+a_e = g_e + b_e(z_e) + epsilon_e,
+```
+
+where `b_e(z_e)` is a repeatable component imperfection and `epsilon_e` is the
+remaining zero-mean variation. A one-time constant calibration subtracts
+`c_e=E[b_e]`. The best predictor from the local information `F_e` is
+`p_e=E[b_e | F_e]`. Conditional projection gives
+
+```text
+E[(b_e-c_e)^2] - E[(b_e-p_e)^2]
+  = E[(p_e-c_e)^2] >= 0.
+```
+
+Thus constant calibration is optimal only when the predictable imperfection is
+constant over the component's operating states. SDIL estimates `p_e` online
+from neutral observations and subtracts it locally; it does not require the
+simulator's imperfection parameters or a component-by-component oracle table.
+Its remaining error includes predictor approximation, tracking error, and the
+unpredictable part `epsilon_e`.
+
+For `E` components, the squared norm of the residual bias is exactly the sum
+of its component powers:
+
+```text
+E[||delta||^2] = sum_e E[delta_e^2].
+```
+
+If the average residual power per component remains bounded away from zero,
+the total residual power grows linearly in component count. Correlation between
+components can change particular projected modes, but it does not remove this
+sum-of-squares accumulation. Reducing the per-component residual therefore
+reduces the coefficient of size-dependent error growth; it need not make a
+large physical network perfect.
+
+The link from residual bias to optimization error can be stated exactly for a
+local quadratic model. Let
+
+```text
+F(k) = 1/2 (k-k_star)^T H (k-k_star),
+implemented update = -grad F(k) + delta.
+```
+
+For constant `delta` in the range of positive-semidefinite `H`, the displaced
+fixed point and its excess objective are
+
+```text
+k_infinity = k_star + H^+ delta,
+F(k_infinity)-F(k_star) = 1/2 delta^T H^+ delta.
+```
+
+On the positive-curvature subspace this implies
+
+```text
+||delta||^2/(2 lambda_max)
+ <= F(k_infinity)-F(k_star)
+ <= ||delta||^2/(2 lambda_min).
+```
+
+A component of `delta` in the nullspace of `H` causes drift rather than a
+displaced optimum. These statements explain why persistent bias can produce a
+size-dependent floor and why reducing residual power can reduce that floor.
+They do not by themselves determine classification error: the CLLN ladder
+tests the downstream effect using held-out logic tasks, while residual-bias
+measurements remain diagnostic.
+
 ### Incomplete residualization can create multiplicative local instability
 
 Residual power alone is not a stability certificate. Consider a linearized
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()