baseline: add convolutional hierarchical feedback alignment

6 files changed, 277 insertions, 8 deletions

diff --git a/BASELINES.md b/BASELINES.md
index 9c3a320..de7dfe8 100644
--- a/BASELINES.md
+++ b/BASELINES.md
@@ -15,6 +15,11 @@ weak in-house reimplementations from defining the state of the art.
 The audited implementations and completed results are in `RESULTS.md`:
 
 - BP, FA, and DFA share the SDIL forward architecture and data order;
+- convolutional hierarchical FA (`hfa`) mirrors the ResNet residual DAG with
+  independent random 3x3 feedback tensors, exact parameter-free shortcut
+  adjoints, and local ReLU/BatchNorm Jacobians; it never copies or reads a
+  downstream forward convolution and is the mandatory baseline for any learned
+  hierarchical extension;
 - direct node perturbation uses a causal target on every hidden update and no
   learned vectorizer;
 - the no-traffic, `P=0` SDIL backbone is learned direct feedback by node
diff --git a/README.md b/README.md
index c06e9b1..64abf25 100644
--- a/README.md
+++ b/README.md
@@ -115,6 +115,12 @@ 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 same smoke suite also audits a convolutional hierarchical-FA baseline. Its
+random feedback tensors follow the real residual DAG and use local activation
+Jacobians, but remain independent of forward weights. An audit-only symmetric
+copy reproduces exact hidden gradients to `5.3e-16` relative error and the BP
+parameter update to `2.98e-8`; actual HFA never performs that copy.
+
 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
diff --git a/THEORY.md b/THEORY.md
index 676d28e..93ef88c 100644
--- a/THEORY.md
+++ b/THEORY.md
@@ -356,6 +356,7 @@ work `S_l=sum_{j>l} F_j`.
 |:--|:--|:--|:--|:--|:--|
 | BP | forward + reverse | exact transposed Jacobians | ordinary supervised loss | about one reverse pass | saved activations / graph |
 | FA | forward + reverse-like serial feedback | fixed random feedback | ordinary loss | `O(F)` feedback work | activations + feedback states |
+| convolutional hierarchical FA | forward + residual-DAG feedback | independent random 3x3 feedback tensors; exact parameter-free shortcut graph | ordinary loss | one feedback convolution per non-stem forward convolution (`0.988x` forward MACs at ResNet-20, plus local correlations) | hidden feedback maps + forward caches |
 | DFA | forward + direct feedback | fixed random output maps | ordinary loss | `O(sum_l d_l d_out)` | activations + feedback states |
 | learned NP, simultaneous | ordinary forward; one duplicate clean plus `2K` noisy full forwards every `e` batches | learned direct maps; no weight transport | `2KB/e` per ordinary batch | `(1+2K)F/e` | up to `2KB` noisy examples in the vectorized implementation |
 | learned NP, layerwise | ordinary forward; clean baseline plus antithetic suffix replays | learned direct maps | `(1+2KH)B/e` | `(F + 2K sum_l S_l)/e`, quadratic in depth for balanced layers | serial suffix replay, no `2K` batch expansion |
diff --git a/experiments/conv_local_smoke.py b/experiments/conv_local_smoke.py
index 8070d9f..9207e04 100644
--- a/experiments/conv_local_smoke.py
+++ b/experiments/conv_local_smoke.py
@@ -7,8 +7,10 @@ import torch
 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,
+from sdil.conv import (CIFARHierarchicalFAResNet, CIFARLocalResNet,
+                       CIFARSDILResNet, ConvSDILConfig,
+                       channel_subspace_apical_calibration,
+                       conv_hierarchical_step, conv_local_step,
                        simultaneous_conv_node_perturbation,
                        vectorizer_subspace_apical_calibration)
 
@@ -597,6 +599,80 @@ def vectorizer_subspace_estimator_check():
     }
 
 
+def hierarchical_feedback_checks():
+    """The residual feedback graph becomes exact only under an audit copy."""
+    torch.manual_seed(91)
+    exact = CIFARHierarchicalFAResNet(
+        depth=8, base_width=2, seed=92, dtype=torch.float64,
+        normalization="batchnorm", residual_scale=1.0)
+    exact.Q = [value.clone() for value in exact.W]
+    exact.R_out.copy_(-exact.W_out.t())
+    x = torch.randn(3, 3, 32, 32, dtype=torch.float64)
+    y = torch.tensor([1, 5, 8])
+    parameters = (exact.W + exact.gamma + exact.beta
+                  + [exact.W_out, exact.b_out])
+    for parameter in parameters:
+        parameter.requires_grad_(True)
+    forward = exact.forward(x, return_cache=True, training=True,
+                            update_stats=False)
+    gradients = torch.autograd.grad(
+        F.cross_entropy(forward["logits"], y), forward["hiddens"])
+    output_error = (torch.softmax(forward["logits"].detach(), dim=1)
+                    - F.one_hot(y, 10).to(torch.float64))
+    teaching = exact.hierarchical_teaching(output_error, forward)
+    relative = [float((signal + x.shape[0] * gradient).abs().max()
+                      / gradient.abs().max().clamp_min(1e-30)
+                      / x.shape[0])
+                for signal, gradient in zip(teaching, gradients)]
+    assert max(relative) < 2e-12
+    for parameter in parameters:
+        parameter.requires_grad_(False)
+
+    # With independent feedback the forward model is bitwise unchanged, while
+    # changing only the feedback seed changes Q/R. This is the actual baseline.
+    left = CIFARHierarchicalFAResNet(
+        depth=8, base_width=2, seed=93, feedback_seed=1001)
+    right = CIFARHierarchicalFAResNet(
+        depth=8, base_width=2, seed=93, feedback_seed=1002)
+    assert all(torch.equal(a, b) for a, b in zip(
+        left.W + [left.W_out], right.W + [right.W_out]))
+    assert any(not torch.equal(a, b) for a, b in zip(
+        left.Q + [left.R_out], right.Q + [right.R_out]))
+
+    # The local update must match exact BP when feedback is explicitly copied
+    # in this audit-only comparator. Actual HFA never performs this copy.
+    local = CIFARHierarchicalFAResNet(
+        depth=8, base_width=2, seed=94, normalization="batchnorm",
+        residual_scale=1.0)
+    bp = CIFARLocalResNet(
+        depth=8, base_width=2, seed=94, normalization="batchnorm",
+        residual_scale=1.0)
+    local.Q = [value.clone() for value in local.W]
+    local.R_out.copy_(-local.W_out.t())
+    xf = torch.randn(4, 3, 32, 32)
+    yf = torch.tensor([0, 2, 5, 9])
+    eta = 0.013
+    conv_hierarchical_step(
+        local, xf, yf, ConvSDILConfig(
+            eta=eta, eta_output=eta, eta_A=0.0, momentum=0.0,
+            weight_decay=0.0, learn_A=False))
+    bp.bp_step(xf, yf, eta, momentum=0.0, weight_decay=0.0)
+    parameter_error = max(float((a - b).abs().max()) for a, b in zip(
+        local.W + local.gamma + local.beta + [local.W_out, local.b_out],
+        bp.W + bp.gamma + bp.beta + [bp.W_out, bp.b_out]))
+    running_error = max(float((a - b).abs().max()) for a, b in zip(
+        local.running_mean + local.running_var,
+        bp.running_mean + bp.running_var))
+    assert parameter_error < 2e-7
+    assert running_error == 0.0
+    return {
+        "hierarchical_symmetric_hidden_relative_error": max(relative),
+        "hierarchical_symmetric_update_absolute_error": parameter_error,
+        "hierarchical_feedback_to_forward_mac_ratio": (
+            local.apical_macs_per_example / local.forward_macs_per_example),
+    }
+
+
 def apical_learning_checks():
     torch.manual_seed(11)
     net = CIFARSDILResNet(depth=8, base_width=2, seed=6)
@@ -702,6 +778,7 @@ def main():
     report.update(perturbation_estimator_check())
     report.update(channel_subspace_estimator_check())
     report.update(vectorizer_subspace_estimator_check())
+    report.update(hierarchical_feedback_checks())
     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 3c0de00..162a8be 100644
--- a/experiments/conv_run.py
+++ b/experiments/conv_run.py
@@ -12,9 +12,11 @@ import time
 import torch
 
 sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
-from sdil.conv import (CIFARLocalResNet, CIFARSDILResNet, ConvSDILConfig,
+from sdil.conv import (CIFARHierarchicalFAResNet, CIFARLocalResNet,
+                       CIFARSDILResNet, ConvSDILConfig,
                        conv_alignment_report, conv_apical_calibration_step,
-                       conv_local_step, evaluate_conv)
+                       conv_hierarchical_alignment_report,
+                       conv_hierarchical_step, conv_local_step, evaluate_conv)
 from sdil.data import DATA_DIR, get_cifar_image_splits
 
 
@@ -101,6 +103,16 @@ def build(args):
         bn_momentum=args.bn_momentum, bn_eps=args.bn_eps)
     if args.mode == "bp":
         return CIFARLocalResNet(**common), None
+    if args.mode == "hfa":
+        net = CIFARHierarchicalFAResNet(
+            **common, feedback_seed=args.apical_seed,
+            feedback_scale=args.a_scale)
+        config = ConvSDILConfig(
+            eta=args.lr, eta_output=args.output_lr, eta_A=0.0,
+            momentum=args.momentum, weight_decay=args.weight_decay,
+            learn_A=False, learn_P=False)
+        config.validate()
+        return net, config
     net = CIFARSDILResNet(
         **common, a_scale=args.a_scale, apical_seed=args.apical_seed,
         vectorizer_mode=args.vectorizer_mode)
@@ -213,7 +225,7 @@ def run(args):
         "schema_version": 1,
         "protocol_family": "oral_a_cifar_local_resnet_development",
         "calibration_metric_space": (
-            None if config is None else {
+            None if config is None or args.mode == "hfa" else {
                 "unit_targets": "full_hidden_field",
                 "channel_subspace": "channel_basis_moments",
                 "vectorizer_subspace": "vectorizer_parameter_gradients",
@@ -238,6 +250,8 @@ def run(args):
             "vectorizer_mode": getattr(net, "vectorizer_mode", None),
             "fixed_traffic_coefficients": getattr(
                 net, "n_fixed_traffic_coefficients", 0),
+            "fixed_feedback_parameters": getattr(
+                net, "n_fixed_feedback_parameters", 0),
         },
         "epochs": [],
     }
@@ -328,6 +342,10 @@ def run(args):
                     x, y, lr, momentum=args.momentum,
                     weight_decay=args.weight_decay)
                 did_perturb = False
+            elif args.mode == "hfa":
+                result = conv_hierarchical_step(net, x, y, config)
+                loss = result["loss"]
+                did_perturb = False
             else:
                 result = conv_local_step(
                     net, x, y, config, step, generator=perturb_generator)
@@ -404,8 +422,11 @@ def run(args):
         probe = min(args.alignment_probe, train.x.shape[0])
         sync(args.device)
         diagnostic_start = time.time()
-        diagnostics = conv_alignment_report(
-            net, train.x[:probe], train.y[:probe], config)
+        diagnostics = (conv_hierarchical_alignment_report(
+            net, train.x[:probe], train.y[:probe])
+            if args.mode == "hfa" else
+            conv_alignment_report(
+                net, train.x[:probe], train.y[:probe], config))
         sync(args.device)
         diagnostics["wall_s"] = time.time() - diagnostic_start
         values = diagnostics["teaching_negative_gradient_cosine"]
@@ -449,7 +470,9 @@ def run(args):
 
 def parse_args():
     parser = argparse.ArgumentParser()
-    parser.add_argument("--mode", choices=("bp", "dfa", "sdil", "nodepert"), required=True)
+    parser.add_argument(
+        "--mode", choices=("bp", "dfa", "hfa", "sdil", "nodepert"),
+        required=True)
     parser.add_argument("--out", required=True)
     parser.add_argument("--device", default="cpu")
     parser.add_argument("--data_dir", default=DATA_DIR)
diff --git a/sdil/conv.py b/sdil/conv.py
index 384584e..f0592a6 100644
--- a/sdil/conv.py
+++ b/sdil/conv.py
@@ -411,6 +411,163 @@ class CIFARLocalResNet:
         return float(loss.detach())
 
 
+class CIFARHierarchicalFAResNet(CIFARLocalResNet):
+    """Residual-DAG feedback alignment with independent convolutional weights.
+
+    Feedback follows the actual child edges of the forward residual graph and
+    uses locally available ReLU/BatchNorm Jacobians, but every convolutional
+    feedback tensor is initialized independently and never reads its forward
+    counterpart.  This is a baseline and an infrastructure step for learned
+    hierarchical dendritic feedback, not an SDIL novelty claim.
+    """
+
+    def __init__(self, *args, feedback_seed=None, feedback_scale=1.0, **kwargs):
+        model_seed = kwargs.get("seed", 0)
+        super().__init__(*args, **kwargs)
+        if feedback_scale <= 0:
+            raise ValueError("feedback_scale must be positive")
+        generator = torch.Generator(device="cpu").manual_seed(
+            model_seed + 30011 if feedback_seed is None else feedback_seed)
+        self.Q = []
+        for weight in self.W:
+            fan_in = weight.shape[1] * weight.shape[2] * weight.shape[3]
+            value = (torch.randn(weight.shape, generator=generator)
+                     * (feedback_scale * math.sqrt(2.0 / fan_in)))
+            self.Q.append(value.to(device=weight.device, dtype=weight.dtype))
+        channels = self.W_out.shape[1]
+        self.R_out = (torch.randn(
+            channels, self.n_classes, generator=generator)
+            * (feedback_scale / math.sqrt(channels))).to(
+                device=self.W_out.device, dtype=self.W_out.dtype)
+
+    @property
+    def n_fixed_feedback_parameters(self):
+        # Q[0] maps the stem to pixels and is not used for hidden credit.
+        return (sum(value.numel() for value in self.Q[1:])
+                + self.R_out.numel())
+
+    @property
+    def apical_macs_per_example(self):
+        conv = 0
+        for weight, spec in zip(self.Q[1:], self.layer_specs[1:]):
+            out_channels, in_channels, kh, kw = weight.shape
+            _, height, width = spec.hidden_shape
+            conv += out_channels * height * width * in_channels * kh * kw
+        return int(conv + self.R_out.numel())
+
+    @staticmethod
+    def _option_a_shortcut_transpose(value, in_channels, stride, output_shape):
+        """Adjoint of the parameter-free option-A shortcut."""
+        out_channels = value.shape[1]
+        if in_channels > out_channels:
+            raise ValueError("option-A transpose cannot recover reduced channels")
+        missing = out_channels - in_channels
+        before = missing // 2
+        selected = value[:, before:before + in_channels]
+        if stride == 1:
+            if tuple(selected.shape) != tuple(output_shape):
+                raise ValueError("shortcut transpose shape mismatch")
+            return selected
+        result = value.new_zeros(output_shape)
+        result[:, :, ::2, ::2] = selected
+        return result
+
+    @torch.no_grad()
+    def hierarchical_teaching(self, output_signal, forward):
+        """Propagate teaching fields through independent feedback convolutions."""
+        caches = forward.get("caches")
+        hiddens = forward.get("hiddens")
+        if caches is None or hiddens is None:
+            raise ValueError("hierarchical feedback requires cached hidden states")
+        if len(hiddens) != self.n_hidden:
+            raise ValueError("hierarchical hidden population mismatch")
+        teaching = [torch.zeros_like(value) for value in hiddens]
+        spatial = hiddens[-1].shape[2] * hiddens[-1].shape[3]
+        teaching[-1].copy_(
+            (output_signal @ self.R_out.t())[:, :, None, None] / spatial)
+        for block in reversed(self.blocks):
+            first = block["first"]
+            second = block["second"]
+            parent = first - 1
+
+            second_gate = caches[second]["gate"].to(teaching[second].dtype)
+            second_delta = teaching[second] * second_gate
+            branch_delta, _, _ = self._normalization_backward(
+                second, second_delta * self.residual_scale,
+                caches[second]["normalization"])
+            teaching[first].add_(F.conv_transpose2d(
+                branch_delta, self.Q[second], stride=1, padding=1))
+            teaching[parent].add_(self._option_a_shortcut_transpose(
+                second_delta, block["in_channels"], block["stride"],
+                hiddens[parent].shape))
+
+            first_gate = caches[first]["gate"].to(teaching[first].dtype)
+            first_delta = teaching[first] * first_gate
+            first_delta, _, _ = self._normalization_backward(
+                first, first_delta, caches[first]["normalization"])
+            teaching[parent].add_(F.conv_transpose2d(
+                first_delta, self.Q[first], stride=block["stride"], padding=1,
+                output_padding=block["stride"] - 1))
+        return teaching
+
+
+def conv_hierarchical_step(net, x, y, config):
+    """One hierarchical-FA update using no reverse-mode graph or weight transport."""
+    config.validate()
+    with torch.no_grad():
+        forward = net.forward(
+            x, return_cache=True, training=True, update_stats=True)
+        logits = forward["logits"]
+        loss = F.cross_entropy(logits, y)
+        output_error = (torch.softmax(logits, dim=1)
+                        - F.one_hot(y, net.n_classes).to(logits.dtype))
+        teaching = net.hierarchical_teaching(output_error, forward)
+        total_units = sum(value.numel() for value in teaching)
+        teaching_rms = math.sqrt(
+            sum(float(value.square().sum()) for value in teaching) / total_units)
+        (directions, gamma_directions, beta_directions,
+         output_weight, output_bias) = net.local_ascent_directions(
+            teaching, output_error, forward)
+        net.apply_ascent(
+            directions, output_weight, output_bias,
+            eta_hidden=config.eta, eta_output=config.eta_output,
+            momentum=config.momentum, weight_decay=config.weight_decay,
+            gamma_directions=gamma_directions,
+            beta_directions=beta_directions)
+        return {
+            "loss": float(loss), "did_perturb": False, "calibration": None,
+            "predictor_mse": None, "teaching_rms": teaching_rms,
+            "raw_apical_rms": teaching_rms, "innovation_rms": teaching_rms,
+        }
+
+
+def conv_hierarchical_alignment_report(net, x, y):
+    """Audit hierarchical teaching against exact hidden gradients."""
+    parameters = net.W + net.gamma + net.beta + [net.W_out, net.b_out]
+    for parameter in parameters:
+        parameter.requires_grad_(True)
+    forward = net.forward(x, return_cache=True, training=True, update_stats=False)
+    gradients = torch.autograd.grad(
+        F.cross_entropy(forward["logits"], y), forward["hiddens"])
+    batch = x.shape[0]
+    negative_gradients = [-batch * value.detach() for value in gradients]
+    with torch.no_grad():
+        output_error = (torch.softmax(forward["logits"], dim=1)
+                        - F.one_hot(y, net.n_classes).to(forward["logits"].dtype))
+        teaching = net.hierarchical_teaching(output_error, forward)
+        values = [float(F.cosine_similarity(
+            left.flatten(1), right.flatten(1), dim=1).mean())
+                  for left, right in zip(teaching, negative_gradients)]
+    for parameter in parameters:
+        parameter.requires_grad_(False)
+    return {
+        "normalization_state": "training_batch_stats_without_running_update",
+        "teaching_negative_gradient_cosine": values,
+        "raw_negative_gradient_cosine": values,
+        "innovation_negative_gradient_cosine": values,
+    }
+
+
 class CIFARSDILResNet(CIFARLocalResNet):
     """CIFAR local ResNet with per-unit apical vectorizers and predictors.