Add contribution roadmap and phase transition plot

3 files changed, 589 insertions, 0 deletions

diff --git a/notes/18_contribution_roadmap.md b/notes/18_contribution_roadmap.md
new file mode 100644
index 0000000..179bd5d
--- /dev/null
+++ b/notes/18_contribution_roadmap.md
@@ -0,0 +1,182 @@
+# Contribution Roadmap
+
+Working title:
+
+`Distributional Capacity Bounds for Feedback Alignment in Multilayer Perceptrons`
+
+The current paper should use five contributions. The fifth contribution is not
+extra decoration: it separates structural capacity theory from finite-time
+loss-gap prediction.
+
+## 1. Capacity Formalization
+
+We define the capacity burden induced by feedback alignment.
+
+Layerwise matrix alignment:
+
+```text
+Q_l = cos^2(W_{l+1}^T, B_l)
+```
+
+For isotropic random feedback in dimension `D_l`:
+
+```text
+Q_l ~ Beta(1/2, (D_l - 1)/2)
+```
+
+Alignment threshold `q_l` gives a log-volume capacity cost:
+
+```text
+C_l(q_l) = -log P(Q_l >= q_l)
+```
+
+This formalizes what it means for feedback alignment to consume parameter
+direction volume.
+
+## 2. Scaling Law and Redundancy Exhaustion
+
+Across independent layers, raw feasible volume multiplies while log-capacity
+cost adds:
+
+```text
+p_all = product_l p_l
+C_all = sum_l C_l
+```
+
+For an equal-width MLP with width `n`, depth `L`, and fixed alignment threshold:
+
+```text
+C_all = Θ(L n^2)
+p_all = exp[-Θ(L n^2)]
+```
+
+Functional loss does not need to appear immediately. If total parameter
+dimension is `P`, task dimension is `d`, and alignment imposes `k` generic
+constraints, then hard functional rank loss is:
+
+```text
+Δd_hard = max(0, k - (P - d))
+```
+
+So the FA/BP gap should open when redundant directions are exhausted.
+
+This is the phase-transition contribution.
+
+## 3. Prior-Free Minimax Initialization Bound
+
+Let `b` be the normalized feedback direction in `D` dimensions and `a` be the
+unknown normalized target backward direction.
+
+For any feedback initialization distribution `μ`, define:
+
+```text
+M_μ = E[b b^T]
+trace(M_μ) = 1
+```
+
+Then:
+
+```text
+inf_a E_μ[(a^T b)^2] = λ_min(M_μ) <= 1/D
+```
+
+Therefore:
+
+```text
+sup_μ inf_a E_μ[(a^T b)^2] = 1/D
+```
+
+Isotropic random feedback reaches the bound. Without a prior over `W` or task
+directions, no initialization can beat isotropic random feedback in worst-case
+expected squared alignment.
+
+## 4. Tangent-Operator Loss/Gap Estimator
+
+This is the bridge from capacity regime to finite-time loss gap.
+
+For BP squared loss, standard tangent-kernel residual dynamics gives:
+
+```text
+r_{t+1}^{BP} ≈ (I - η K_t^{BP}/N) r_t
+K_t^{BP} = J_t J_t^T
+```
+
+For FA, the parameter update uses a surrogate backward Jacobian `J_tilde_t`,
+but output change is still measured by the true forward Jacobian `J_t`:
+
+```text
+r_{t+1}^{FA} ≈ (I - η K_t^{FA}/N) r_t
+K_t^{FA} = J_t J_tilde_t^T
+```
+
+`K_t^{FA}` is better called a tangent operator, not a PSD kernel.
+
+The local estimator freezes `K_0`. The finite-time estimator uses an early
+operator velocity:
+
+```text
+K_hat_t = K_0 + t (K_s - K_0) / s
+```
+
+Then:
+
+```text
+L_hat_T = ||r_hat_T||^2 / (2N)
+gap_hat_T = L_hat_T^FA - L_hat_T^BP
+```
+
+This contribution is conditional on early operator observations `(K_0, K_s)`.
+It is not an architecture-only theorem.
+
+## 5. Distributional Empirical Validation
+
+Experiments should validate three levels:
+
+1. static alignment distributions:
+   ```text
+   Q_l ~ Beta(1/2, (D_l - 1)/2)
+   ```
+2. capacity and redundancy transition:
+   FA/BP gap opens near the hard FA capacity margin crossing;
+3. tangent-operator trajectory distributions:
+   predicted gap distributions overlap empirical trajectory distributions.
+
+Current strong evidence:
+
+- 256-trajectory finite-time overlap plot for `d=2,w=64,T=50,s=20`;
+- stress grid over depth, width, and horizon;
+- systematic low-bias analysis showing residual error is kernel-path curvature,
+  not normalization.
+
+## Important Separation
+
+Do not claim:
+
+```text
+capacity alone predicts exact loss gap
+```
+
+Correct claim:
+
+```text
+capacity bounds identify the structural regime;
+tangent-operator dynamics quantify the finite-time gap within that regime.
+```
+
+This separation avoids overclaiming while still giving a coherent theory chain.
+
+## Next Priority
+
+The weakest current visual is the phase-transition contribution.
+
+The desired figure should show:
+
+```text
+capacity margin decreases -> redundant directions exhausted -> FA/BP gap opens
+and then grows
+```
+
+The plot should use hard FA capacity margin on the x-axis, place the zero
+margin as a vertical reference, and show BP/FA train loss or FA-BP train gap.
+The most direct target is a width sweep on a random-label task, because random
+labels make the task dimension controllable and force memorization capacity.
diff --git a/notes/19_phase_transition_visualization.md b/notes/19_phase_transition_visualization.md
new file mode 100644
index 0000000..a2672c2
--- /dev/null
+++ b/notes/19_phase_transition_visualization.md
@@ -0,0 +1,162 @@
+# Phase Transition Visualization
+
+This note records the first clean phase-transition visualization for the
+capacity-exhaustion contribution.
+
+## Goal
+
+We want a figure that shows:
+
+```text
+available FA capacity decreases
+-> redundant directions are exhausted
+-> FA/BP train gap opens
+-> gap grows while BP remains capable
+```
+
+The key x-axis is the hard FA capacity margin:
+
+```text
+M_FA = P - K_FA - N*out
+```
+
+where:
+
+- `P` is the parameter count;
+- `K_FA` is the hard feedback-alignment constraint count;
+- `N*out` is the random-label task dimension.
+
+The predicted transition is at:
+
+```text
+M_FA = 0
+```
+
+## Experiment
+
+This run uses SGD, not Adam.
+
+Configuration:
+
+```text
+task: random-label regression
+architecture: 16 -> width -> width -> 4
+widths: 8, 12, 16, 24, 32, 48, 64, 96
+train samples: 128
+optimizer: full-batch SGD
+learning rate: 0.01
+steps: 3000
+init seeds: 4
+feedback seeds per init: 8
+FA trajectories: 256
+```
+
+Output directory:
+
+```text
+outputs/phase_transition_sgd_256_lr001_T3000
+```
+
+Main figures:
+
+```text
+outputs/phase_transition_sgd_256_lr001_T3000/phase_transition_capacity_exhaustion.png
+outputs/phase_transition_sgd_256_lr001_T3000/phase_transition_gap_only.png
+```
+
+## Result
+
+Summary:
+
+| width | FA margin | BP train | FA train | FA-BP train gap |
+|---:|---:|---:|---:|---:|
+| 96 | 1026 | 0.000249 | 0.009436 | 0.009187 |
+| 64 | 514 | 0.003354 | 0.053192 | 0.049838 |
+| 48 | 258 | 0.016732 | 0.146680 | 0.129947 |
+| 32 | 2 | 0.089511 | 0.386777 | 0.297266 |
+| 24 | -126 | 0.201897 | 0.600796 | 0.398899 |
+| 16 | -254 | 0.547408 | 0.993219 | 0.445812 |
+| 12 | -318 | 0.872346 | 1.204686 | 0.332340 |
+| 8 | -382 | 1.260234 | 1.443789 | 0.183555 |
+
+The gap is small in the redundant region:
+
+```text
+M_FA > 0
+```
+
+It rises sharply around the hard FA capacity boundary:
+
+```text
+M_FA ≈ 0
+```
+
+and becomes largest in the FA-deficient but BP-capable region:
+
+```text
+M_FA < 0
+and
+M_BP = P - N*out > 0
+```
+
+At the smallest widths, BP itself also becomes under-capacity:
+
+```text
+M_BP < 0
+```
+
+In that both-deficient regime, the FA-BP gap decreases because both methods
+fail to fit the random labels.
+
+## Interpretation
+
+The visual story is not simply "smaller network means larger gap." The clean
+phase-transition story is:
+
+1. high redundancy:
+   ```text
+   M_FA >> 0
+   ```
+   BP and FA both have enough effective room, so the train gap is small.
+
+2. FA redundancy exhausted:
+   ```text
+   M_FA crosses 0
+   ```
+   FA starts paying the feedback-alignment burden in task-relevant directions,
+   so the gap opens.
+
+3. FA-deficient but BP-capable:
+   ```text
+   M_FA < 0, M_BP > 0
+   ```
+   the gap grows and peaks.
+
+4. both-deficient:
+   ```text
+   M_BP < 0
+   ```
+   BP also cannot memorize, so the difference between FA and BP no longer
+   grows monotonically.
+
+This is the right phase-transition framing for the paper.
+
+## Caveat
+
+The hard margin is intentionally simple. It uses:
+
+```text
+K_FA = (width*width - 1) + (width*out - 1)
+```
+
+for the two-hidden-layer MLP. It predicts the transition location well enough
+to provide a clean regime variable, but the exact peak is shifted because real
+FA alignment is soft, not a hard rank constraint.
+
+Do not claim that hard capacity margin alone predicts the exact final train
+gap. The correct claim is:
+
+```text
+hard FA margin predicts the regime where the gap opens;
+tangent/operator dynamics or empirical trajectories determine the gap size.
+```
diff --git a/scripts/plot_phase_transition.py b/scripts/plot_phase_transition.py
new file mode 100644
index 0000000..025e267
--- /dev/null
+++ b/scripts/plot_phase_transition.py
@@ -0,0 +1,245 @@
+#!/usr/bin/env python3
+"""Plot the capacity-exhaustion phase transition for BP/FA."""
+
+from __future__ import annotations
+
+import argparse
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "--run-dir",
+        type=Path,
+        default=Path("outputs/phase_transition_sgd_256_lr001_T3000"),
+    )
+    parser.add_argument(
+        "--outdir",
+        type=Path,
+        default=Path("outputs/phase_transition_sgd_256_lr001_T3000"),
+    )
+    return parser.parse_args()
+
+
+def summarize_runs(runs: pd.DataFrame, summary: pd.DataFrame) -> pd.DataFrame:
+    fa = runs[runs["run_type"] == "fa"].copy()
+    records = []
+    for margin, group in fa.groupby("fa_capacity_margin"):
+        gaps = group["train_gap_to_bp"].to_numpy(dtype=np.float64)
+        records.append(
+            {
+                "fa_capacity_margin": float(margin),
+                "gap_mean": float(gaps.mean()),
+                "gap_q05": float(np.quantile(gaps, 0.05)),
+                "gap_q25": float(np.quantile(gaps, 0.25)),
+                "gap_q75": float(np.quantile(gaps, 0.75)),
+                "gap_q95": float(np.quantile(gaps, 0.95)),
+            }
+        )
+    gap_summary = pd.DataFrame(records)
+    merged = summary.merge(gap_summary, on="fa_capacity_margin", how="left")
+    return merged.sort_values("fa_capacity_margin", ascending=False)
+
+
+def trajectory_id(frame: pd.DataFrame) -> pd.Series:
+    return frame["init_seed"].astype(str) + ":" + frame["feedback_seed"].astype(str)
+
+
+def add_regime_shading(ax: plt.Axes, summary: pd.DataFrame) -> None:
+    max_margin = float(summary["fa_capacity_margin"].max())
+    min_margin = float(summary["fa_capacity_margin"].min())
+    bp_bad_margin = float(
+        summary.loc[summary["bp_capacity_margin"] < 0, "fa_capacity_margin"].max()
+    )
+    ax.axvspan(max_margin + 80, 0, color="#e8f1fb", alpha=0.55, linewidth=0)
+    ax.axvspan(0, bp_bad_margin, color="#fff1d8", alpha=0.55, linewidth=0)
+    ax.axvspan(bp_bad_margin, min_margin - 80, color="#fde2e2", alpha=0.55, linewidth=0)
+    ax.axvline(0, color="black", linestyle="--", linewidth=1.1)
+    ax.axvline(bp_bad_margin, color="#8b1a1a", linestyle=":", linewidth=1.2)
+
+
+def plot_transition(run_dir: Path, outdir: Path) -> list[Path]:
+    runs = pd.read_csv(run_dir / "runs.csv")
+    summary = pd.read_csv(run_dir / "width_summary.csv")
+    merged = summarize_runs(runs, summary)
+    fa_runs = runs[runs["run_type"] == "fa"].copy()
+    fa_runs["_trajectory_id"] = trajectory_id(fa_runs)
+
+    outdir.mkdir(parents=True, exist_ok=True)
+    paths: list[Path] = []
+
+    fig, (ax_loss, ax_gap) = plt.subplots(
+        2,
+        1,
+        figsize=(11.2, 8.4),
+        dpi=180,
+        sharex=True,
+        height_ratios=[1.0, 1.12],
+    )
+
+    add_regime_shading(ax_loss, merged)
+    add_regime_shading(ax_gap, merged)
+
+    x = merged["fa_capacity_margin"].to_numpy(dtype=np.float64)
+    ax_loss.plot(
+        x,
+        merged["bp_train_mean"],
+        color="#174f91",
+        marker="o",
+        linewidth=2.1,
+        label="BP train loss",
+    )
+    ax_loss.plot(
+        x,
+        merged["fa_train_mean"],
+        color="#a84708",
+        marker="o",
+        linewidth=2.1,
+        label="FA train loss",
+    )
+    ax_loss.fill_between(
+        x,
+        np.maximum(merged["bp_train_mean"] - merged["bp_train_std"], 1e-8),
+        merged["bp_train_mean"] + merged["bp_train_std"],
+        color="#174f91",
+        alpha=0.13,
+        linewidth=0,
+    )
+    ax_loss.fill_between(
+        x,
+        np.maximum(merged["fa_train_mean"] - merged["fa_train_std"], 1e-8),
+        merged["fa_train_mean"] + merged["fa_train_std"],
+        color="#a84708",
+        alpha=0.13,
+        linewidth=0,
+    )
+    ax_loss.set_yscale("log")
+    ax_loss.set_ylabel("train MSE")
+    ax_loss.set_title("Capacity exhaustion opens the FA/BP train-gap")
+    ax_loss.legend(loc="upper left", frameon=True)
+    ax_loss.grid(alpha=0.18)
+
+    for _tid, group in fa_runs.groupby("_trajectory_id"):
+        group = group.sort_values("fa_capacity_margin", ascending=False)
+        ax_gap.plot(
+            group["fa_capacity_margin"],
+            group["train_gap_to_bp"],
+            color="#c65f16",
+            alpha=0.13,
+            linewidth=0.95,
+        )
+    ax_gap.fill_between(
+        x,
+        merged["gap_q05"],
+        merged["gap_q95"],
+        color="#c65f16",
+        alpha=0.13,
+        linewidth=0,
+        label="FA trajectories 5-95%",
+    )
+    ax_gap.fill_between(
+        x,
+        merged["gap_q25"],
+        merged["gap_q75"],
+        color="#c65f16",
+        alpha=0.24,
+        linewidth=0,
+        label="FA trajectories 25-75%",
+    )
+    ax_gap.plot(
+        x,
+        merged["train_gap_mean"],
+        color="#a84708",
+        marker="o",
+        linewidth=2.2,
+        label="mean gap",
+    )
+    ax_gap.axhline(0, color="black", linewidth=1.0)
+    ax_gap.set_ylabel("FA train MSE - BP train MSE")
+    ax_gap.set_xlabel("hard FA capacity margin P - K_FA - N*out (capacity decreases ->)")
+    ax_gap.grid(alpha=0.18)
+    ax_gap.legend(loc="upper left", frameon=True)
+
+    max_margin = float(merged["fa_capacity_margin"].max())
+    min_margin = float(merged["fa_capacity_margin"].min())
+    ax_gap.set_xlim(max_margin + 80, min_margin - 80)
+
+    ax_loss.text(
+        0.18,
+        0.92,
+        "redundant",
+        transform=ax_loss.transAxes,
+        ha="center",
+        va="center",
+        fontsize=10,
+        color="#174f91",
+    )
+    ax_loss.text(
+        0.57,
+        0.92,
+        "FA-deficient\nBP-capable",
+        transform=ax_loss.transAxes,
+        ha="center",
+        va="center",
+        fontsize=10,
+        color="#9a5b00",
+    )
+    ax_loss.text(
+        0.90,
+        0.92,
+        "both-deficient",
+        transform=ax_loss.transAxes,
+        ha="center",
+        va="center",
+        fontsize=10,
+        color="#8b1a1a",
+    )
+
+    fig.tight_layout()
+    path = outdir / "phase_transition_capacity_exhaustion.png"
+    fig.savefig(path)
+    plt.close(fig)
+    paths.append(path)
+
+    detail_fig, ax = plt.subplots(figsize=(8.4, 5.4), dpi=180)
+    add_regime_shading(ax, merged)
+    ax.errorbar(
+        x,
+        merged["train_gap_mean"],
+        yerr=merged["train_gap_std"],
+        color="#a84708",
+        marker="o",
+        linewidth=2.1,
+        capsize=3,
+        label="mean +/- sd",
+    )
+    ax.axhline(0, color="black", linewidth=1.0)
+    ax.set_xlim(max_margin + 80, min_margin - 80)
+    ax.set_xlabel("hard FA capacity margin P - K_FA - N*out (capacity decreases ->)")
+    ax.set_ylabel("FA train MSE - BP train MSE")
+    ax.set_title("FA/BP gap increases after redundant capacity is exhausted")
+    ax.legend(loc="upper left", frameon=True)
+    ax.grid(alpha=0.18)
+    detail_fig.tight_layout()
+    detail_path = outdir / "phase_transition_gap_only.png"
+    detail_fig.savefig(detail_path)
+    plt.close(detail_fig)
+    paths.append(detail_path)
+
+    return paths
+
+
+def main() -> None:
+    args = parse_args()
+    paths = plot_transition(args.run_dir, args.outdir)
+    for path in paths:
+        print(f"plot: {path}")
+
+
+if __name__ == "__main__":
+    main()