diff options
| -rw-r--r-- | THEORY.md | 72 | ||||
| -rw-r--r-- | experiments/verify_theory.py | 66 |
2 files changed, 138 insertions, 0 deletions
@@ -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() |