From b0c8b9cc9d2eac64f64bd9063ccfc5607867e87d Mon Sep 17 00:00:00 2001 From: YurenHao0426 Date: Wed, 22 Jul 2026 12:32:14 -0500 Subject: algorithm: calibrate feedback in representable subspace --- README.md | 6 +- THEORY.md | 56 +++++++++++++ experiments/conv_local_smoke.py | 149 ++++++++++++++++++++++++++++++++- experiments/conv_run.py | 5 ++ sdil/conv.py | 177 +++++++++++++++++++++++++++++++++++++--- 5 files changed, 381 insertions(+), 12 deletions(-) diff --git a/README.md b/README.md index a690477..d13f8f4 100644 --- a/README.md +++ b/README.md @@ -98,7 +98,11 @@ The convolutional infrastructure can be checked independently: These checks cover deterministic CIFAR splitting/augmentation, the standard CIFAR `6n+2` ResNet topology, exact local convolution and BatchNorm eligibilities, the BatchNorm-coupled perturbation objective, predictor fitting, -translation-shared feedback, and parameter/cost accounting. +translation-shared feedback, representable-subspace causal calibration, and +parameter/cost accounting. The latter is a transparent post-A3 development +path: it perturbs the channel-gated basis directly and reproduces the expected +full-field delta rule with far less irrelevant spatial variance. It does not +retroactively reopen or relabel the failed frozen Oral-A protocol. `experiments/finalize_accept.sh` additionally requires strict imports of the frozen BurstCCN and Dual Propagation author-code runs. Both records now pass; diff --git a/THEORY.md b/THEORY.md index da3e77a..8d20659 100644 --- a/THEORY.md +++ b/THEORY.md @@ -76,6 +76,62 @@ The important dependence is `O(sigma^2)`, with a dimension-dependent constant. Reducing K removes variance but not this bias; increasing width/depth without a sigma audit is therefore unsafe. +### Perturbing the representable convolutional subspace + +The failed ResNet-20 run exposed a separate projection problem. Its +translation-shared feedback field for example `i`, channel `c`, and spatial +site `s` is + +```text +p_ics = b_ic + t_ics v_ic, t_ics = tanh(h_ics). +``` + +The legacy estimator perturbs every `h_ics` independently, constructs a full +unit target, and only then averages its delta-rule update into the shared +`b/v` parameterization. At K1, almost all sampled energy can lie outside this +two-dimensional per-channel subspace. SDIL-v2 instead draws independent +Rademacher coefficients `z^b_ic,z^v_ic` and intervenes along + +```text +xi_ics = (z^b_ic + t_ics z^v_ic) / sqrt(2). +``` + +Let `q_ics=-dL/dh_ics`, let `S` be the number of spatial sites, and let `D` be +the centered directional loss derivative. Then + +```text +qhat^b_ic = -sqrt(2) D z^b_ic / S, +qhat^v_ic = -sqrt(2) D z^v_ic / S +``` + +are unbiased for the two negative-gradient moments + +```text +qbar^b_ic = mean_s q_ics, +qbar^v_ic = mean_s (q_ics t_ics). +``` + +All other examples, channels, and layers contribute zero-mean cross terms. +For BatchNorm, `D` is the derivative of the summed minibatch loss, so the same +identity targets each example's contribution including the true cross-example +Jacobian coupling. + +This does not mistake moments for expansion coefficients. The predicted +moments of `p=b+t v` are + +```text +pbar^b = b + mean(t) v, +pbar^v = mean(t) b + mean(t^2) v. +``` + +Consequently updating the two vectorizer pathways with +`qhat^b-pbar^b` and `qhat^v-pbar^v` is exactly the expected gradient of the +original full-field squared prediction error. It uses the same two causal +loss queries per direction as unit perturbation, but avoids estimating +unrepresentable spatial components. It remains a variance reduction, not a +guarantee of useful feedback or stable end-to-end optimization; those require +a separately frozen empirical gate. + ## 2. What is sufficient for a descent step Flatten all forward parameters into `theta`, let `p = grad L(theta)`, and let diff --git a/experiments/conv_local_smoke.py b/experiments/conv_local_smoke.py index 99ca3b6..1ea0cea 100644 --- a/experiments/conv_local_smoke.py +++ b/experiments/conv_local_smoke.py @@ -8,7 +8,8 @@ import torch.nn.functional as F sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) from sdil.conv import (CIFARLocalResNet, CIFARSDILResNet, ConvSDILConfig, - conv_local_step, simultaneous_conv_node_perturbation) + channel_subspace_apical_calibration, conv_local_step, + simultaneous_conv_node_perturbation) def architecture_checks(): @@ -266,6 +267,138 @@ def perturbation_estimator_check(): } +def channel_subspace_estimator_check(): + """The structured estimator targets representable base/gate moments.""" + torch.manual_seed(71) + batch = 3 + hiddens = [ + torch.randn(batch, 2, 4, 4, dtype=torch.float64), + torch.randn(batch, 3, 2, 2, dtype=torch.float64), + ] + negative_gradients = [torch.randn_like(value) for value in hiddens] + estimated_base = [torch.zeros( + batch, value.shape[1], dtype=value.dtype) for value in hiddens] + estimated_gate = [torch.zeros_like(value) for value in estimated_base] + generator = torch.Generator(device="cpu").manual_seed(211) + directions = 4096 + inverse_sqrt_two = 1.0 / (2.0 ** 0.5) + for _ in range(directions): + hidden_directions = [] + base_random = [] + gate_random = [] + for hidden in hiddens: + shape = (batch, hidden.shape[1]) + base = torch.empty(shape, dtype=hidden.dtype).bernoulli_( + 0.5, generator=generator).mul_(2).sub_(1) + gate = torch.empty_like(base).bernoulli_( + 0.5, generator=generator).mul_(2).sub_(1) + hidden_directions.append(( + base[:, :, None, None] + + torch.tanh(hidden) * gate[:, :, None, None]) + * inverse_sqrt_two) + base_random.append(base) + gate_random.append(gate) + # The exact loss derivative uses g=-negative_gradient and remains + # per-example without BatchNorm coupling. + directional = -sum( + (gradient * direction).flatten(1).sum(dim=1) + for gradient, direction in zip( + negative_gradients, hidden_directions)) + for index, (hidden, base, gate) in enumerate(zip( + hiddens, base_random, gate_random)): + spatial = hidden.shape[2] * hidden.shape[3] + scale = -(2.0 ** 0.5) / (spatial * directions) + estimated_base[index].add_(directional[:, None] * base, alpha=scale) + estimated_gate[index].add_(directional[:, None] * gate, alpha=scale) + exact_base = [value.mean(dim=(2, 3)) for value in negative_gradients] + exact_gate = [(value * torch.tanh(hidden)).mean(dim=(2, 3)) + for value, hidden in zip(negative_gradients, hiddens)] + estimated = torch.cat([ + value.flatten() for pair in zip(estimated_base, estimated_gate) + for value in pair]) + exact = torch.cat([ + value.flatten() for pair in zip(exact_base, exact_gate) + for value in pair]) + cosine = float(F.cosine_similarity(estimated, exact, dim=0)) + norm_ratio = float(estimated.norm() / exact.norm()) + assert cosine > 0.985 + assert 0.90 < norm_ratio < 1.10 + + # The executable antithetic implementation must match the same structured + # directional derivative, not an autograd surrogate. + # Use a fixed nondegenerate point. Width-one, zero-bias ReLU networks can + # contain structurally exact-zero preactivations, where central differences + # and PyTorch's chosen subgradient need not agree even as sigma -> 0. + torch.manual_seed(3) + net = CIFARSDILResNet( + depth=8, base_width=2, seed=72, dtype=torch.float64, + vectorizer_mode="channel_gated") + x = torch.randn(2, 3, 32, 32, dtype=torch.float64) + y = torch.tensor([2, 8]) + parameters = net.W + [net.W_out, net.b_out] + for parameter in parameters: + parameter.requires_grad_(True) + clean = net.forward(x, return_cache=True) + for hidden in clean["hiddens"]: + hidden.retain_grad() + loss = F.cross_entropy(clean["logits"], y) + loss.backward() + output_signal = (torch.softmax(clean["logits"].detach(), dim=1) + - F.one_hot(y, 10).to(torch.float64)) + _, diagnostics = channel_subspace_apical_calibration( + net, x, y, clean, output_signal, sigma=1e-6, + n_directions=1, eta=0.0, + generator=torch.Generator(device="cpu").manual_seed(307), + return_diagnostics=True) + hidden_direction = diagnostics["directions"][0]["hidden"] + finite_difference = diagnostics["directional_derivatives"][0][ + "scaled_directional"] + exact_directional = sum( + (x.shape[0] * hidden.grad * direction).flatten(1).sum(dim=1) + for hidden, direction in zip(clean["hiddens"], hidden_direction)) + relative = ((finite_difference - exact_directional).abs() + / exact_directional.abs().clamp_min(1e-12)) + assert float(relative.max()) < 2e-6 + for parameter in parameters: + parameter.requires_grad_(False) + + # The local A/G update must equal the full-field delta rule after replacing + # only its two target moments with the structured causal estimates. + eta = 0.0023 + before_a = [value.clone() for value in net.A] + before_g = [value.clone() for value in net.A_gate] + expected_a = [] + expected_g = [] + for index, hidden in enumerate(clean["hiddens"]): + base = output_signal @ before_a[index].t() + gate_coefficient = output_signal @ before_g[index].t() + gate = torch.tanh(hidden.detach()) + mean = gate.mean(dim=(2, 3)) + second = gate.square().mean(dim=(2, 3)) + base_error = diagnostics["target_base"][index] - ( + base + mean * gate_coefficient) + gate_error = diagnostics["target_gate"][index] - ( + mean * base + second * gate_coefficient) + expected_a.append(before_a[index] + eta * ( + base_error.t() @ output_signal / x.shape[0])) + expected_g.append(before_g[index] + eta * ( + gate_error.t() @ output_signal / x.shape[0])) + channel_subspace_apical_calibration( + net, x, y, clean, output_signal, sigma=1e-6, + n_directions=1, eta=eta, + generator=torch.Generator(device="cpu").manual_seed(307)) + update_error = max(float((actual - expected).abs().max()) + for actual, expected in zip( + net.A + net.A_gate, expected_a + expected_g)) + assert update_error < 1e-14 + return { + "channel_subspace_moment_cosine": cosine, + "channel_subspace_moment_norm_ratio": norm_ratio, + "channel_subspace_jvp_relative_error": float(relative.max()), + "channel_subspace_delta_rule_absolute_error": update_error, + } + + def apical_learning_checks(): torch.manual_seed(11) net = CIFARSDILResNet(depth=8, base_width=2, seed=6) @@ -342,6 +475,19 @@ def apical_learning_checks(): for before_weight, after_weight in zip(weights_before, net.W)) assert all(not parameter.requires_grad for parameter in net.W + [net.W_out, net.b_out]) + + gated_weights_before = [value.clone() for value in gated.A + gated.A_gate] + gated_result = conv_local_step( + gated, x[:2], y[:2], + ConvSDILConfig( + eta=1e-3, eta_A=1e-3, momentum=0.0, weight_decay=0.0, + pert_every=1, apical_calibration_mode="channel_subspace"), + step=0, generator=torch.Generator(device="cpu").manual_seed(17)) + assert gated_result["did_perturb"] + assert all(torch.isfinite(torch.tensor(value)) for value in + gated_result["calibration"].values()) + assert any(not torch.equal(before, after) for before, after in zip( + gated_weights_before, gated.A + gated.A_gate)) return {"apical_mse_ratio": after / before, "gated_apical_mse_ratio": gated_after / gated_before, "gated_vectorizer_parameter_reduction": ( @@ -356,6 +502,7 @@ def main(): report = exact_local_gradient_check() report.update(exact_batchnorm_local_gradient_check()) report.update(perturbation_estimator_check()) + report.update(channel_subspace_estimator_check()) report.update(apical_learning_checks()) print(report) print("ALL CONVOLUTIONAL LOCAL-ELIGIBILITY CHECKS PASSED") diff --git a/experiments/conv_run.py b/experiments/conv_run.py index 98c1f19..cb99545 100644 --- a/experiments/conv_run.py +++ b/experiments/conv_run.py @@ -111,6 +111,7 @@ def build(args): use_residual=bool(args.use_residual), nuisance_scale=args.nuisance_scale, pert_sigma=args.pert_sigma, pert_every=args.pert_every, pert_directions=args.pert_directions, + apical_calibration_mode=args.apical_calibration_mode, direct_node_perturbation=args.mode == "nodepert") config.validate() return net, config @@ -481,6 +482,10 @@ def parse_args(): parser.add_argument("--pert_sigma", type=float, default=0.01) parser.add_argument("--pert_every", type=int, default=4) parser.add_argument("--pert_directions", type=int, default=1) + parser.add_argument( + "--apical_calibration_mode", + choices=("unit_targets", "channel_subspace"), + default="unit_targets") parser.add_argument("--predictor_warmup_steps", type=int, default=0) parser.add_argument("--a_warmup_steps", type=int, default=0) parser.add_argument("--alignment_probe", type=int, default=0) diff --git a/sdil/conv.py b/sdil/conv.py index 279acac..7935027 100644 --- a/sdil/conv.py +++ b/sdil/conv.py @@ -640,6 +640,142 @@ def simultaneous_conv_node_perturbation(net, x, y, clean_forward, sigma=1e-2, return targets +@torch.no_grad() +def channel_subspace_apical_calibration( + net, x, y, clean_forward, output_signal, sigma=1e-2, + n_directions=1, eta=1e-3, generator=None, return_diagnostics=False): + """Calibrate channel-gated feedback in its representable causal subspace. + + The legacy estimator perturbs every spatial unit independently, estimates a + full hidden target, and only then averages that target into the shared + channel coefficients. With K=1, most of its variance lies outside the + vectorizer's representable subspace. Here each intervention is instead + ``(z_base + tanh(h) z_gate) / sqrt(2)`` with one Rademacher coefficient per + example and channel. If ``D`` is the antithetic loss derivative and ``S`` + is the number of spatial sites, ``-sqrt(2) D z/S`` is an unbiased estimate + of the corresponding negative-gradient moment. Cross-example and + cross-layer terms remain zero mean. Subtracting the predicted moments + gives exactly the expected local delta-rule update that full unit targets + would produce, without first estimating directions outside the feedback + model's representable subspace. + + This is still forward-only causal calibration: it consumes the same two + scalar loss queries per direction, never differentiates through the + network, and updates only the local A/A_gate tensors. + """ + if getattr(net, "vectorizer_mode", None) != "channel_gated": + raise ValueError("channel-subspace calibration requires channel_gated A") + if sigma <= 0 or n_directions < 1 or eta < 0: + raise ValueError("invalid channel-subspace calibration hyperparameters") + if len(clean_forward["hiddens"]) != net.n_hidden: + raise ValueError("clean forward does not match network hidden populations") + if generator is None: + generator = torch.Generator(device=x.device).manual_seed(0) + batch = x.shape[0] + target_base = [torch.zeros( + batch, hidden.shape[1], device=hidden.device, dtype=hidden.dtype) + for hidden in clean_forward["hiddens"]] + target_gate = [torch.zeros_like(value) for value in target_base] + diagnostic_directions = [] + diagnostic_derivatives = [] + inverse_sqrt_two = 1.0 / math.sqrt(2.0) + for _ in range(n_directions): + base_random = [] + gate_random = [] + directions = [] + for hidden in clean_forward["hiddens"]: + shape = (batch, hidden.shape[1]) + base = torch.empty( + shape, device=hidden.device, dtype=hidden.dtype).bernoulli_( + 0.5, generator=generator).mul_(2).sub_(1) + gate = torch.empty_like(base).bernoulli_( + 0.5, generator=generator).mul_(2).sub_(1) + direction = (base[:, :, None, None] + + torch.tanh(hidden) * gate[:, :, None, None]) + direction.mul_(inverse_sqrt_two) + base_random.append(base) + gate_random.append(gate) + directions.append(direction) + plus = F.cross_entropy(net.forward( + x, perturbations=[sigma * value for value in directions], + training=True, update_stats=False)["logits"], y, reduction="none") + minus = F.cross_entropy(net.forward( + x, perturbations=[-sigma * value for value in directions], + training=True, update_stats=False)["logits"], y, reduction="none") + if net.normalization == "batchnorm": + batch_directional = (plus.mean() - minus.mean()) / (2.0 * sigma) + directional = batch_directional.mul(batch).expand(batch) + else: + batch_directional = None + directional = (plus - minus) / (2.0 * sigma) + for index, (hidden, base, gate) in enumerate(zip( + clean_forward["hiddens"], base_random, gate_random)): + spatial = hidden.shape[2] * hidden.shape[3] + scale = -math.sqrt(2.0) / (spatial * n_directions) + target_base[index].add_(directional[:, None] * base, alpha=scale) + target_gate[index].add_(directional[:, None] * gate, alpha=scale) + if return_diagnostics: + diagnostic_directions.append({ + "hidden": directions, "base": base_random, "gate": gate_random}) + diagnostic_derivatives.append({ + "scaled_directional": directional, + "batch_mean_directional": batch_directional, + "coupling": ("batch_objective" if batch_directional is not None + else "per_example_objective"), + }) + + before_error = 0.0 + target_power = 0.0 + prediction_power = 0.0 + dot = 0.0 + update_power = 0.0 + coefficients = 0 + for index, (base_target, gate_target) in enumerate(zip( + target_base, target_gate)): + base_coefficient = output_signal @ net.A[index].t() + gate_coefficient = output_signal @ net.A_gate[index].t() + gate = torch.tanh(clean_forward["hiddens"][index]) + gate_mean = gate.mean(dim=(2, 3)) + gate_second_moment = gate.square().mean(dim=(2, 3)) + # For prediction b + tanh(h) g, these are its inner products with + # the two representable basis fields. The errors are therefore the + # exact stochastic gradients of full-field squared prediction error. + base_prediction = base_coefficient + gate_mean * gate_coefficient + gate_prediction = (gate_mean * base_coefficient + + gate_second_moment * gate_coefficient) + base_error = base_target - base_prediction + gate_error = gate_target - gate_prediction + base_update = base_error.t() @ output_signal / batch + gate_update = gate_error.t() @ output_signal / batch + net.A[index].add_(base_update, alpha=eta) + net.A_gate[index].add_(gate_update, alpha=eta) + for prediction, target, error in ( + (base_prediction, base_target, base_error), + (gate_prediction, gate_target, gate_error)): + before_error += float(error.square().sum()) + target_power += float(target.square().sum()) + prediction_power += float(prediction.square().sum()) + dot += float((prediction * target).sum()) + coefficients += target.numel() + update_power += float(base_update.square().sum() + gate_update.square().sum()) + denominator = math.sqrt(target_power * prediction_power) + calibration = { + "calibration_mse": before_error / coefficients, + "target_power": target_power / coefficients, + "prediction_target_cosine": dot / denominator if denominator else 0.0, + "parameter_update_rms": math.sqrt( + update_power / max(1, net.n_vectorizer_parameters)), + } + if return_diagnostics: + return calibration, { + "directions": diagnostic_directions, + "directional_derivatives": diagnostic_derivatives, + "target_base": target_base, + "target_gate": target_gate, + } + return calibration + + @dataclass class ConvSDILConfig: eta: float = 0.01 @@ -655,6 +791,7 @@ class ConvSDILConfig: pert_sigma: float = 1e-2 pert_every: int = 4 pert_directions: int = 1 + apical_calibration_mode: str = "unit_targets" direct_node_perturbation: bool = False def validate(self): @@ -664,8 +801,14 @@ class ConvSDILConfig: raise ValueError("apical learning rates must be nonnegative") if self.pert_every < 1 or self.pert_directions < 1: raise ValueError("perturbation cadence/directions must be positive") + if self.apical_calibration_mode not in ( + "unit_targets", "channel_subspace"): + raise ValueError("unknown apical calibration mode") if self.direct_node_perturbation and self.pert_every != 1: raise ValueError("direct node perturbation requires a target every step") + if (self.direct_node_perturbation + and self.apical_calibration_mode != "unit_targets"): + raise ValueError("direct node perturbation requires unit targets") def conv_local_step(net, x, y, config, step, generator=None): @@ -692,10 +835,18 @@ def conv_local_step(net, x, y, config, step, generator=None): did_perturb = ((config.learn_A or config.direct_node_perturbation) and step % config.pert_every == 0) targets = None + calibration = None if did_perturb: - targets = simultaneous_conv_node_perturbation( - net, x, y, forward, sigma=config.pert_sigma, - n_directions=config.pert_directions, generator=generator) + if config.apical_calibration_mode == "channel_subspace": + calibration = channel_subspace_apical_calibration( + net, x, y, forward, output_error, + sigma=config.pert_sigma, + n_directions=config.pert_directions, eta=config.eta_A, + generator=generator) + else: + targets = simultaneous_conv_node_perturbation( + net, x, y, forward, sigma=config.pert_sigma, + n_directions=config.pert_directions, generator=generator) weight_teaching = targets if config.direct_node_perturbation else teaching if weight_teaching is None: raise RuntimeError("direct perturbation target is unavailable") @@ -708,8 +859,8 @@ def conv_local_step(net, x, y, config, step, generator=None): momentum=config.momentum, weight_decay=config.weight_decay, gamma_directions=gamma_directions, beta_directions=beta_directions) - calibration = None - if did_perturb and config.learn_A: + if (did_perturb and config.learn_A + and config.apical_calibration_mode == "unit_targets"): calibration = net.calibrate_apical( output_error, forward["hiddens"], teaching, targets, config.eta_A) predictor_mse = None @@ -742,11 +893,17 @@ def conv_apical_calibration_step(net, x, y, config, generator=None): output_error, forward["hiddens"], config.nuisance_scale, config.use_residual) del raw, innovations - targets = simultaneous_conv_node_perturbation( - net, x, y, forward, sigma=config.pert_sigma, - n_directions=config.pert_directions, generator=generator) - calibration = net.calibrate_apical( - output_error, forward["hiddens"], teaching, targets, config.eta_A) + if config.apical_calibration_mode == "channel_subspace": + calibration = channel_subspace_apical_calibration( + net, x, y, forward, output_error, sigma=config.pert_sigma, + n_directions=config.pert_directions, eta=config.eta_A, + generator=generator) + else: + targets = simultaneous_conv_node_perturbation( + net, x, y, forward, sigma=config.pert_sigma, + n_directions=config.pert_directions, generator=generator) + calibration = net.calibrate_apical( + output_error, forward["hiddens"], teaching, targets, config.eta_A) return float(F.cross_entropy(logits, y)), calibration -- cgit v1.2.3