algorithm: estimate causal credit in vectorizer space

5 files changed, 403 insertions, 7 deletions

diff --git a/README.md b/README.md
index 4d4d3df..61f42de 100644
--- a/README.md
+++ b/README.md
@@ -110,6 +110,13 @@ 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.
 
+The subsequent V3 mechanism estimates the required A/G matrix statistics
+directly by perturbing the vectorizer parameter subspace. It remains
+forward-only and uses two causal loss queries per direction. Its exact and
+BatchNorm JVP identities, delta rule, and matched synthetic variance comparison
+are covered by `experiments/conv_local_smoke.py`; no task-level claim is made
+until a new frozen gate is completed.
+
 `experiments/finalize_accept.sh` additionally requires strict imports of the
 frozen BurstCCN and Dual Propagation author-code runs. Both records now pass;
 the accept-bar mechanical audit is green. The standard ResNet branch stopped
diff --git a/THEORY.md b/THEORY.md
index 8d20659..676d28e 100644
--- a/THEORY.md
+++ b/THEORY.md
@@ -132,6 +132,62 @@ 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.
 
+### Perturbing the vectorizer parameter subspace
+
+The channel-moment estimator still samples one independent base/gate
+coefficient for every example and channel and only afterward regresses those
+noisy coefficients on the output context.  The quantity the local delta rule
+actually needs is much smaller.  For
+
+```text
+p_ics = (A c_i)_c + t_ics (G c_i)_c
+```
+
+define the normalized prediction objective
+
+```text
+J(A,G) = (1 / 2B) sum_i mean_s ||p_i - q_i||^2,
+q_i = -d ell_i / dh_i.
+```
+
+Its causal target terms are matrices, not per-example coefficient fields:
+
+```text
+U_A = (1/B) sum_i mean_s(q_ics) c_i^T,
+U_G = (1/B) sum_i mean_s(q_ics t_ics) c_i^T.
+```
+
+Draw independent Rademacher matrices `Xi_A,Xi_G` with the same shapes as
+`A,G`, independently across layers, and intervene along
+
+```text
+d_ics = ((Xi_A c_i)_c + t_ics (Xi_G c_i)_c) / sqrt(2).
+```
+
+Let `D` be the centered derivative of the summed minibatch loss under all
+layer interventions. Since `D=-sum_ics q_ics d_ics`, independence gives
+
+```text
+Uhat_A = -sqrt(2) D Xi_A / (B S),     E[Uhat_A] = U_A,
+Uhat_G = -sqrt(2) D Xi_G / (B S),     E[Uhat_G] = U_G.
+```
+
+Cross-layer and A/G cross terms vanish in expectation. The prediction terms
+`(1/B) sum mean(p)c^T` and `(1/B) sum mean(pt)c^T` are locally available
+exactly, so subtracting them from `Uhat_A,Uhat_G` gives an unbiased stochastic
+gradient of `J`. BatchNorm changes neither identity: `D` is explicitly the
+summed minibatch derivative and therefore includes its true cross-example
+Jacobian.
+
+This estimator perturbs only the shared mapping that will survive
+amortization. It uses the same two loss queries per direction as the earlier
+estimators. It does not uniformly dominate them for every context covariance
+or batch size, but at the frozen batch-128 synthetic audit its one-direction
+MSE is `0.03381x` the coefficient-then-regress estimator's MSE. The executable
+mean has cosine `0.96928` and norm ratio `1.04888` to the exact A/G target.
+Those mechanics do not establish task utility; a separately frozen gate is
+still required.
+
 ## 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 1ea0cea..8070d9f 100644
--- a/experiments/conv_local_smoke.py
+++ b/experiments/conv_local_smoke.py
@@ -9,7 +9,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,
                        channel_subspace_apical_calibration, conv_local_step,
-                       simultaneous_conv_node_perturbation)
+                       simultaneous_conv_node_perturbation,
+                       vectorizer_subspace_apical_calibration)
 
 
 def architecture_checks():
@@ -391,6 +392,7 @@ def channel_subspace_estimator_check():
                        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,
@@ -399,6 +401,202 @@ def channel_subspace_estimator_check():
     }
 
 
+def vectorizer_subspace_estimator_check():
+    """Direct A/G perturbations are unbiased and lower variance at batch 128."""
+    torch.manual_seed(81)
+    batch = 128
+    output_dim = 5
+    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]
+    output_signal = torch.randn(batch, output_dim, dtype=torch.float64)
+    exact_pairs = []
+    for hidden, target in zip(hiddens, negative_gradients):
+        exact_pairs.extend([
+            target.mean(dim=(2, 3)).t() @ output_signal / batch,
+            (target * torch.tanh(hidden)).mean(dim=(2, 3)).t()
+            @ output_signal / batch,
+        ])
+    exact = torch.cat([value.flatten() for value in exact_pairs])
+    coefficient_sum = torch.zeros_like(exact)
+    vectorizer_sum = torch.zeros_like(exact)
+    coefficient_mse = 0.0
+    vectorizer_mse = 0.0
+    directions = 2048
+    generator = torch.Generator(device="cpu").manual_seed(401)
+    inverse_sqrt_two = 1.0 / (2.0 ** 0.5)
+    for _ in range(directions):
+        random_coefficients = []
+        hidden_directions = []
+        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)
+            random_coefficients.append((base, gate))
+            hidden_directions.append((
+                base[:, :, None, None]
+                + torch.tanh(hidden) * gate[:, :, None, None])
+                * inverse_sqrt_two)
+        directional = -sum((target * direction).sum()
+                           for target, direction in zip(
+                               negative_gradients, hidden_directions))
+        coefficient_values = []
+        for hidden, (base, gate) in zip(hiddens, random_coefficients):
+            spatial = hidden.shape[2] * hidden.shape[3]
+            base_target = -(2.0 ** 0.5) * directional * base / spatial
+            gate_target = -(2.0 ** 0.5) * directional * gate / spatial
+            coefficient_values.extend([
+                base_target.t() @ output_signal / batch,
+                gate_target.t() @ output_signal / batch,
+            ])
+        coefficient_sample = torch.cat(
+            [value.flatten() for value in coefficient_values])
+
+        random_matrices = []
+        hidden_directions = []
+        for hidden in hiddens:
+            shape = (hidden.shape[1], output_dim)
+            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)
+            base_field = output_signal @ base.t()
+            gate_field = output_signal @ gate.t()
+            random_matrices.append((base, gate))
+            hidden_directions.append((
+                base_field[:, :, None, None]
+                + torch.tanh(hidden) * gate_field[:, :, None, None])
+                * inverse_sqrt_two)
+        directional = -sum((target * direction).sum()
+                           for target, direction in zip(
+                               negative_gradients, hidden_directions))
+        vectorizer_values = []
+        for hidden, (base, gate) in zip(hiddens, random_matrices):
+            spatial = hidden.shape[2] * hidden.shape[3]
+            scale = -(2.0 ** 0.5) * directional / (batch * spatial)
+            vectorizer_values.extend([scale * base, scale * gate])
+        vectorizer_sample = torch.cat(
+            [value.flatten() for value in vectorizer_values])
+        coefficient_sum.add_(coefficient_sample)
+        vectorizer_sum.add_(vectorizer_sample)
+        coefficient_mse += float((coefficient_sample - exact).square().mean())
+        vectorizer_mse += float((vectorizer_sample - exact).square().mean())
+    coefficient_mean = coefficient_sum / directions
+    vectorizer_mean = vectorizer_sum / directions
+    vectorizer_cosine = float(F.cosine_similarity(
+        vectorizer_mean, exact, dim=0))
+    vectorizer_norm_ratio = float(vectorizer_mean.norm() / exact.norm())
+    variance_ratio = vectorizer_mse / coefficient_mse
+    assert vectorizer_cosine > 0.95
+    assert 0.90 < vectorizer_norm_ratio < 1.10
+    assert variance_ratio < 0.25
+
+    # Match the executable forward-only derivative and its exact A/G update.
+    torch.manual_seed(3)
+    net = CIFARSDILResNet(
+        depth=8, base_width=2, seed=82, dtype=torch.float64,
+        vectorizer_mode="channel_gated")
+    x = torch.randn(2, 3, 32, 32, dtype=torch.float64)
+    y = torch.tensor([1, 6])
+    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()
+    F.cross_entropy(clean["logits"], y).backward()
+    output_error = (torch.softmax(clean["logits"].detach(), dim=1)
+                    - F.one_hot(y, 10).to(torch.float64))
+    _, diagnostics = vectorizer_subspace_apical_calibration(
+        net, x, y, clean, output_error, sigma=1e-6,
+        n_directions=1, eta=0.0,
+        generator=torch.Generator(device="cpu").manual_seed(503),
+        return_diagnostics=True)
+    finite_difference = diagnostics["directional_derivatives"][0][
+        "scaled_directional"]
+    exact_directional = sum(
+        (x.shape[0] * hidden.grad * direction).sum()
+        for hidden, direction in zip(
+            clean["hiddens"], diagnostics["directions"][0]["hidden"]))
+    jvp_relative = float((finite_difference - exact_directional).abs()
+                         / exact_directional.abs().clamp_min(1e-12))
+    assert jvp_relative < 2e-6
+    for parameter in parameters:
+        parameter.requires_grad_(False)
+
+    eta = 0.0017
+    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_error @ before_a[index].t()
+        gate_coefficient = output_error @ before_g[index].t()
+        gate = torch.tanh(hidden.detach())
+        mean = gate.mean(dim=(2, 3))
+        second = gate.square().mean(dim=(2, 3))
+        base_prediction = (base + mean * gate_coefficient).t() @ output_error / 2
+        gate_prediction = (
+            mean * base + second * gate_coefficient).t() @ output_error / 2
+        expected_a.append(before_a[index] + eta * (
+            diagnostics["target_base"][index] - base_prediction))
+        expected_g.append(before_g[index] + eta * (
+            diagnostics["target_gate"][index] - gate_prediction))
+    vectorizer_subspace_apical_calibration(
+        net, x, y, clean, output_error, sigma=1e-6,
+        n_directions=1, eta=eta,
+        generator=torch.Generator(device="cpu").manual_seed(503))
+    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
+
+    batchnorm = CIFARSDILResNet(
+        depth=8, base_width=2, seed=83, dtype=torch.float64,
+        normalization="batchnorm", residual_scale=1.0,
+        vectorizer_mode="channel_gated")
+    xb = torch.randn(3, 3, 32, 32, dtype=torch.float64)
+    yb = torch.tensor([0, 4, 9])
+    parameters = (batchnorm.W + batchnorm.gamma + batchnorm.beta
+                  + [batchnorm.W_out, batchnorm.b_out])
+    for parameter in parameters:
+        parameter.requires_grad_(True)
+    clean_b = batchnorm.forward(xb, training=True, update_stats=False)
+    for hidden in clean_b["hiddens"]:
+        hidden.retain_grad()
+    F.cross_entropy(clean_b["logits"], yb).backward()
+    output_b = (torch.softmax(clean_b["logits"].detach(), dim=1)
+                - F.one_hot(yb, 10).to(torch.float64))
+    _, diagnostics_b = vectorizer_subspace_apical_calibration(
+        batchnorm, xb, yb, clean_b, output_b, sigma=1e-6,
+        n_directions=1, eta=0.0,
+        generator=torch.Generator(device="cpu").manual_seed(509),
+        return_diagnostics=True)
+    finite_b = diagnostics_b["directional_derivatives"][0][
+        "scaled_directional"]
+    exact_b = sum(
+        (xb.shape[0] * hidden.grad * direction).sum()
+        for hidden, direction in zip(
+            clean_b["hiddens"], diagnostics_b["directions"][0]["hidden"]))
+    batchnorm_relative = float(
+        (finite_b - exact_b).abs() / exact_b.abs().clamp_min(1e-12))
+    assert batchnorm_relative < 2e-6
+    for parameter in parameters:
+        parameter.requires_grad_(False)
+    return {
+        "vectorizer_subspace_mean_cosine": vectorizer_cosine,
+        "vectorizer_subspace_mean_norm_ratio": vectorizer_norm_ratio,
+        "vectorizer_vs_coefficient_mse_ratio": variance_ratio,
+        "vectorizer_subspace_jvp_relative_error": jvp_relative,
+        "vectorizer_subspace_batchnorm_jvp_relative_error": batchnorm_relative,
+        "vectorizer_subspace_delta_rule_absolute_error": update_error,
+    }
+
+
 def apical_learning_checks():
     torch.manual_seed(11)
     net = CIFARSDILResNet(depth=8, base_width=2, seed=6)
@@ -503,6 +701,7 @@ def main():
     report.update(exact_batchnorm_local_gradient_check())
     report.update(perturbation_estimator_check())
     report.update(channel_subspace_estimator_check())
+    report.update(vectorizer_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 a33005d..3c0de00 100644
--- a/experiments/conv_run.py
+++ b/experiments/conv_run.py
@@ -213,10 +213,11 @@ def run(args):
         "schema_version": 1,
         "protocol_family": "oral_a_cifar_local_resnet_development",
         "calibration_metric_space": (
-            None if config is None else
-            ("channel_basis_moments"
-             if config.apical_calibration_mode == "channel_subspace"
-             else "full_hidden_field")),
+            None if config is None else {
+                "unit_targets": "full_hidden_field",
+                "channel_subspace": "channel_basis_moments",
+                "vectorizer_subspace": "vectorizer_parameter_gradients",
+            }[config.apical_calibration_mode]),
         "args": vars(args),
         "provenance": provenance(),
         "split": split,
@@ -496,7 +497,7 @@ def parse_args():
     parser.add_argument("--pert_directions", type=int, default=1)
     parser.add_argument(
         "--apical_calibration_mode",
-        choices=("unit_targets", "channel_subspace"),
+        choices=("unit_targets", "channel_subspace", "vectorizer_subspace"),
         default="unit_targets")
     parser.add_argument("--predictor_warmup_steps", type=int, default=0)
     parser.add_argument("--a_warmup_steps", type=int, default=0)
diff --git a/sdil/conv.py b/sdil/conv.py
index 7935027..384584e 100644
--- a/sdil/conv.py
+++ b/sdil/conv.py
@@ -776,6 +776,128 @@ def channel_subspace_apical_calibration(
     return calibration
 
 
+@torch.no_grad()
+def vectorizer_subspace_apical_calibration(
+        net, x, y, clean_forward, output_signal, sigma=1e-2,
+        n_directions=1, eta=1e-3, generator=None, return_diagnostics=False):
+    """Estimate the causal A/G delta rule directly in parameter space.
+
+    Channel-subspace calibration first estimates one causal coefficient target
+    per example/channel and then regresses those targets on ``output_signal``.
+    This estimator instead draws Rademacher matrices with the exact shapes of
+    A and A_gate.  Their induced hidden intervention already contains the
+    output context, so the antithetic scalar directly estimates the matrix
+    moments ``mean(q c^T)`` and ``mean(q tanh(h) c^T)`` required by the local
+    vectorizer delta rule.  It uses the same two loss queries per direction and
+    removes variance in coefficient directions that the shared A/G maps cannot
+    represent.
+    """
+    if getattr(net, "vectorizer_mode", None) != "channel_gated":
+        raise ValueError("vectorizer-subspace calibration requires channel_gated A")
+    if sigma <= 0 or n_directions < 1 or eta < 0:
+        raise ValueError("invalid vectorizer-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_like(value) for value in net.A]
+    target_gate = [torch.zeros_like(value) for value in net.A_gate]
+    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, base_map, gate_map in zip(
+                clean_forward["hiddens"], net.A, net.A_gate):
+            base = torch.empty_like(base_map).bernoulli_(
+                0.5, generator=generator).mul_(2).sub_(1)
+            gate = torch.empty_like(gate_map).bernoulli_(
+                0.5, generator=generator).mul_(2).sub_(1)
+            base_field = output_signal @ base.t()
+            gate_field = output_signal @ gate.t()
+            direction = (base_field[:, :, None, None]
+                         + torch.tanh(hidden)
+                         * gate_field[:, :, 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)
+        minus = F.cross_entropy(net.forward(
+            x, perturbations=[-sigma * value for value in directions],
+            training=True, update_stats=False)["logits"], y)
+        # A/G are shared across the minibatch, so their sufficient statistic is
+        # the derivative of the summed loss even when examples are uncoupled.
+        directional = (plus - minus) * batch / (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) / (
+                batch * spatial * n_directions)
+            target_base[index].add_(directional * base, alpha=scale)
+            target_gate[index].add_(directional * gate, alpha=scale)
+        if return_diagnostics:
+            diagnostic_directions.append({
+                "hidden": directions, "base": base_random,
+                "gate": gate_random})
+            diagnostic_derivatives.append({
+                "scaled_directional": directional,
+                "coupling": "summed_batch_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))
+        base_moment = base_coefficient + gate_mean * gate_coefficient
+        gate_moment = (gate_mean * base_coefficient
+                       + gate_second_moment * gate_coefficient)
+        base_prediction = base_moment.t() @ output_signal / batch
+        gate_prediction = gate_moment.t() @ output_signal / batch
+        base_error = base_target - base_prediction
+        gate_error = gate_target - gate_prediction
+        net.A[index].add_(base_error, alpha=eta)
+        net.A_gate[index].add_(gate_error, 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())
+            update_power += float(error.square().sum())
+            coefficients += target.numel()
+    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 / coefficients),
+    }
+    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
@@ -802,7 +924,7 @@ class ConvSDILConfig:
         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"):
+                "unit_targets", "channel_subspace", "vectorizer_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")
@@ -843,6 +965,12 @@ def conv_local_step(net, x, y, config, step, generator=None):
                     sigma=config.pert_sigma,
                     n_directions=config.pert_directions, eta=config.eta_A,
                     generator=generator)
+            elif config.apical_calibration_mode == "vectorizer_subspace":
+                calibration = vectorizer_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,
@@ -898,6 +1026,11 @@ def conv_apical_calibration_step(net, x, y, config, generator=None):
             net, x, y, forward, output_error, sigma=config.pert_sigma,
             n_directions=config.pert_directions, eta=config.eta_A,
             generator=generator)
+    elif config.apical_calibration_mode == "vectorizer_subspace":
+        calibration = vectorizer_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,