Add AAAI depth experiments and diagnostic figures

20 files changed, 2490 insertions, 9 deletions

diff --git a/.gitignore b/.gitignore
index 12c5ea2..148c1cd 100644
--- a/.gitignore
+++ b/.gitignore
@@ -40,6 +40,9 @@ Thumbs.db
 # Logs
 *.log
 
+# Local handoff bundles
+facap_aaai_*.zip
+
 # Paper build artifacts
 paper/*.pdf
 paper/*.aux
diff --git a/README.md b/README.md
index 171bd24..2bdd779 100644
--- a/README.md
+++ b/README.md
@@ -1,9 +1,9 @@
 # facap — What Does Random Feedback Cost?
 
-Theory + experiments quantifying the training cost of **feedback alignment**
-(FA: fixed random backward matrices instead of transposed weights) relative to
-backpropagation in MLPs. Target venue: **AAAI-27** (abstract 2026-07-20,
-paper 2026-07-27).
+Theory + experiments quantifying the training cost of **fixed random
+feedback** relative to backpropagation.  Exact initialization results cover
+FA and DFA; finite-time analysis and experiments focus on FA. Target venue:
+**AAAI-27** (abstract 2026-07-20, paper 2026-07-27).
 
 ## Results in one table
 
@@ -11,17 +11,19 @@ paper 2026-07-27).
 |---|---|
 | static alignment is an exact Beta(1/2,(D-1)/2) law; log-volume cost adds to Theta(Ln^2) | 700k samples, KS p 0.17-0.99; Gamma(L,1) max KS 0.0042 |
 | no prior-free feedback initialization beats isotropic (minimax 1/D) | exact theorem + lam_min table |
-| initial FA/BP speed gap = hidden-layer BP speed share, exactly, any depth | synthetic 0.4378 vs 0.4403; MNIST diff <= 0.0089 over share 0.39-0.82 |
+| exact expected initial deficit = nonoutput BP decrease share, for FA and DFA | depth 1/2/3/4/6: max row error 0.0092; CE CNN max error 0.0083 |
 | one-hidden-layer e0 is exactly Gaussian | KS p 0.13-0.999, real-backward sampler |
-| the finite-time gap is a closed-form **soft ramp** — no phase transition | dense T=30000: log-linear R^2 0.993; closed form corr 0.933 matched / 0.982 ensemble; lam_min corr 0.987 |
-| a two-snapshot early-velocity estimator predicts the gap with no fit | 256 traj: MAE 0.0019, corr 0.999; stress grid 0.994-0.99989; MNIST MAE 0.0020 (3.9% rel) |
+| random-direction deletion is not an FA cost model | up to 250x finite-time overprediction |
+| matrix mismatch does not set FA's finite-time cost | depth 1->6: joint alignment proxy drops 18.1 orders while width-64 gap grows 0.0421->0.3461; frozen MAE 0.0669 vs two-snapshot 0.0104/0.0139 |
+| initialization-rate compensation is material but incomplete | gap reduction 29-51% across depths; independently tuned BP/FA retain a gap |
 | optimization cost is invariant; test-side effect is task-dependent | teacher task: test gap sign flip (-0.120+/-0.013 at w=32) |
 
 ## Layout
 
-- `notes/` — project memory. Entry points: `notes/41_paper_plan.md` (plan),
+- `notes/` — project memory. Entry points: `notes/43_aaai_accept_bar_experiments.md`
+  (new experiments), `notes/41_paper_plan.md` (older writing plan),
   `notes/36_evidence_ledger.md` (numbers), `notes/40_reproduction_manifest.md`
-  (runs). `notes/README.md` indexes all 41 notes with status.
+  (runs). `notes/README.md` indexes all 43 notes with status.
 - `scripts/` — all experiments (python, float64, CPU; deterministic seeds).
   `scripts/README.md` has the status table and recipes.
 - `outputs/` — gitignored, regenerable. Paper-grade runs are exactly those in
diff --git a/notes/40_reproduction_manifest.md b/notes/40_reproduction_manifest.md
index 009f058..d7390cb 100644
--- a/notes/40_reproduction_manifest.md
+++ b/notes/40_reproduction_manifest.md
@@ -54,10 +54,23 @@ deterministic given the seeds shown. Long runs use `nohup`.
 | kernel code self-test vs autograd (<= 5.3e-15) | `python scripts/real_data_validation.py --self-test` | stdout | 38 |
 | teacher task: train gap 0.249->0.0069 soft; test gap sign flip (-0.120+/-0.013 at w=32) | `python scripts/downstream_capacity_sweep.py --task mlp --teacher-width 64 --teacher-hidden-layers 2 --normalize-targets --widths 8 12 16 24 32 48 64 96 --train-samples 256 --test-samples 2048 --probe-samples 32 --steps 8000 --lr 0.01 --optimizer sgd --init-seeds 3 --feedback-seeds 8 --data-seed 7 --skip-jacobian --plot --outdir outputs/teacher_test_gap_sgd_T8000` then `python scripts/plot_teacher_test_gap.py` | `teacher_test_gap_sgd_T8000/` | 39 |
 
+## Focused-paper controls added 2026-07-21
+
+| claim | command | outputs | note |
+|---|---|---|---|
+| FA and DFA exact initialization calibration through depth 6: 60 rows, max abs error 0.009164 | `python scripts/aaai_depth_experiments.py --part init --depths 1 2 3 4 6 --width 64 --train-samples 64 --init-seeds 6 --feedback-draws 256 --torch-threads 16 --outdir outputs/aaai_depth_experiments` | `aaai_depth_experiments/initialization_rows.csv` | 43 |
+| FA finite-time depth x width grid: matrix-alignment proxy drops 18.06 orders while width-64 gap grows 0.0421 -> 0.3461; frozen MAE 0.06693; two-snapshot MAE 0.01043 / 0.01394 | `python scripts/aaai_depth_experiments.py --part finite --depths 1 2 3 4 6 --finite-widths 32 64 96 --train-samples 64 --finite-init-seeds 4 --finite-feedback-seeds 8 --lr 0.001 --horizon 50 --early-step 20 --torch-threads 16 --outdir outputs/aaai_depth_experiments`, then `--part summarize` | `aaai_depth_experiments/finite_time_*`, `mismatch_vs_cost_by_depth.*` | 43 |
+| parameter-matched depth control around P=5376: widths 269/64/47/39/31; gap is nonmonotone and peaks at depth 4 | `python scripts/parameter_matched_depth_control.py --depths 1 2 3 4 6 --reference-depth 2 --reference-width 64 --train-samples 64 --init-seeds 4 --feedback-seeds 8 --lr 0.001 --horizon 50 --early-step 20 --torch-threads 16 --outdir outputs/parameter_matched_depth_control` | `parameter_matched_depth_control/` | 43 |
+| initialization-rate compensation reduces but does not remove the gap (29-51% reduction); separate stability-confirmed tuning preserves a gap | `python scripts/learning_rate_compensation.py --depths 1 2 3 4 6 --width 64 --train-samples 128 --steps 200 --base-lr 0.001 --tune-init-seeds 4 --tune-feedback-seeds 4 --eval-init-seeds 6 --eval-feedback-seeds 8 --torch-threads 16 --outdir outputs/learning_rate_compensation` | `learning_rate_compensation/` | 43 |
+| cross-entropy CNN initialization identity for FA and DFA: max abs error 0.008309, max standardized error 1.308 | `python scripts/cnn_initialization_validation.py --train-samples 128 --image-size 8 --channels 8 12 --init-seeds 4 --feedback-draws 2048 --torch-threads 8 --outdir outputs/cnn_initialization_validation` | `cnn_initialization_validation/` | 43 |
+| two focused three-panel figures | `python scripts/plot_aaai_main_figures.py --outdir outputs/aaai_main_figures` | `aaai_main_figures/` | 43 |
+
 ## Figure -> file map (current best versions)
 
 | paper figure (plan in note 41) | file |
 |---|---|
+| focused Figure 1: mismatch vs exact initial cost | `aaai_main_figures/figure1_mismatch_vs_initial_cost.{png,pdf}` |
+| focused Figure 2: later cost and operator evolution | `aaai_main_figures/figure2_training_dynamics.{png,pdf}` |
 | soft ramp + closed form + ramp law | `closed_form_soft_ramp/closed_form_vs_measured_ramp.png` |
 | e0 moment calibration (synthetic) | `actual_fa_initial_operator_moments/predicted_vs_empirical_initial_erosion_mean.png` |
 | e0 Gaussian panels | `actual_fa_initial_erosion_distribution/e0_distribution_theory_vs_empirical.png` |
diff --git a/notes/41_paper_plan.md b/notes/41_paper_plan.md
index d3bd12f..41f6a05 100644
--- a/notes/41_paper_plan.md
+++ b/notes/41_paper_plan.md
@@ -1,5 +1,11 @@
 # AAAI-27 Paper Plan
 
+> **Superseded framing (2026-07-21).** The original four-figure soft-ramp plan
+> below predates the focused matrix-mismatch-versus-optimization-cost review.
+> New experiment results and the two replacement composite figures are in
+> note 43.  Keep this file for venue constraints and historical scheduling,
+> not as the current narrative or figure plan.
+
 The single planning document for the writing phase. Inputs: evidence ledger
 (note 36), reproduction manifest (note 40), corrections (note 37). Supersedes
 note 03 and the figure list in note 35.
diff --git a/notes/42_master_technical_reference.md b/notes/42_master_technical_reference.md
index 82c436e..65a2ae5 100644
--- a/notes/42_master_technical_reference.md
+++ b/notes/42_master_technical_reference.md
@@ -5,6 +5,12 @@ derivations, exact numbers, configs, and negative results. Built for writing
 the paper by hand: you should not need any other file open. Cross-checked
 against the raw CSVs on 2026-06-13.
 
+> **2026-07-21 experiment addendum:** note 43 adds the FA/DFA depth-1-to-6
+> calibration, FA finite-time depth/width grid, learning-rate controls, and
+> cross-entropy CNN validation.  Its results supersede this file's earlier
+> statement that DFA and non-MLP validation were absent; the theorem/proof
+> material below remains the technical source.
+
 Conventions in this file:
 - "no fit" = no post-hoc scalar tuned to the reported data; predictions come
   from initialization or two early kernel snapshots.
diff --git a/notes/43_aaai_accept_bar_experiments.md b/notes/43_aaai_accept_bar_experiments.md
new file mode 100644
index 0000000..11a981d
--- /dev/null
+++ b/notes/43_aaai_accept_bar_experiments.md
@@ -0,0 +1,272 @@
+# AAAI Accept-Bar Experiments
+
+Status: **complete, paper-grade local outputs** (2026-07-21).  This note
+records the experiments added after the focused-paper review.  It supplements
+the pre-existing master reference (note 42); it does not rewrite the paper.
+
+## 1. What was added
+
+Four missing controls are now implemented and run:
+
+1. FA and DFA initialization-cost calibration at hidden depths `{1,2,3,4,6}`.
+2. FA finite-time gap and operator-predictor errors over depth
+   `{1,2,3,4,6}` and width `{32,64,96}`.
+3. Same-rate, exact-initialization-rate-compensated, and independently tuned
+   BP/FA learning-rate controls.
+4. A non-MLP initialization validation: a two-convolution-layer CNN with
+   cross-entropy loss on an MNIST subset, for both FA and DFA.
+
+Two three-panel composite figures are generated from these results plus the
+existing random-direction-deletion negative control.
+
+## 2. Shared implementation checks
+
+New shared implementation: `scripts/feedback_rules.py`.
+
+- Arbitrary-depth manual BP agrees with autograd to `5.55e-17` max error.
+- FA and DFA output-layer gradients agree with BP exactly.
+- A 4096-feedback-draw numerical moment test for each rule lands within the
+  Monte Carlo interval predicted by the exact initialization result.
+- The CNN manual BP agrees with autograd to `5.55e-17`; FA and DFA preserve
+  the exact output gradient.
+
+Commands:
+
+```bash
+python scripts/test_feedback_rules.py
+python scripts/cnn_initialization_validation.py --self-test
+```
+
+## 3. FA/DFA initialization calibration across depth
+
+Configuration:
+
+- synthetic random regression, `N=64`, input dimension 16, output dimension 4;
+- fixed width 64;
+- hidden depth `{1,2,3,4,6}`;
+- 6 fixed forward initializations per depth;
+- 256 independent feedback draws per forward initialization and rule;
+- real FA/DFA backward pass for every draw (15,360 feedback evaluations).
+
+The predicted quantity is the relative deficit in BP's first-order loss
+decrease:
+
+```text
+1 - ||g_output^BP||^2 / sum_l ||g_l^BP||^2.
+```
+
+Depth-aggregated results:
+
+| depth | predicted | FA measured | DFA measured |
+|---:|---:|---:|---:|
+| 1 | 0.050380 | 0.050200 | 0.050578 |
+| 2 | 0.276966 | 0.275983 | 0.278078 |
+| 3 | 0.402635 | 0.402218 | 0.402282 |
+| 4 | 0.458542 | 0.457015 | 0.456779 |
+| 6 | 0.623193 | 0.624234 | 0.623791 |
+
+Across all 60 rule/depth/init rows, max absolute calibration error is
+`0.009164`; max standardized error is `2.412` Monte Carlo standard errors.
+The same exact prediction covers FA and DFA even though their feedback paths
+differ.
+
+Output: `outputs/aaai_depth_experiments/initialization_rows.csv` and
+`initial_deficit_calibration.png`.
+
+## 4. Finite-time FA cost across depth and width
+
+Configuration:
+
+- same synthetic task family, `N=64`, output dimension 4;
+- hidden depth `{1,2,3,4,6}`, width `{32,64,96}`;
+- 4 forward initializations x 8 FA feedback initializations per cell;
+- full-batch SGD, fixed learning rate `1e-3`, horizon `T=50`;
+- early operator snapshot `s=20`;
+- 480 FA trajectories plus cached matched BP trajectories.
+
+Curve uncertainty is computed across forward-initialization means (feedback
+draws are nested within each initialization), rather than treating all 32
+feedback trajectories in a cell as independent.
+
+At width 64, the measured mean gap rises with depth:
+
+| depth | measured | frozen `K(0)` | linear velocity | early retangent |
+|---:|---:|---:|---:|---:|
+| 1 | 0.042081 | 0.046148 | 0.041747 | 0.043214 |
+| 2 | 0.152697 | 0.176001 | 0.150869 | 0.159079 |
+| 3 | 0.252548 | 0.301208 | 0.246641 | 0.263831 |
+| 4 | 0.281210 | 0.360586 | 0.274093 | 0.302059 |
+| 6 | 0.346053 | 0.518519 | 0.312131 | 0.371649 |
+
+For the same fixed-width architecture, the independent per-block alignment
+proxy `prod_l 1/D_l` falls from `3.91e-3` at depth 1 to `3.39e-21` at depth 6,
+an 18.06-order-of-magnitude change.  Over that sweep, the measured initial
+deficit changes from `0.0504` to `0.6240`, and the measured finite-time gap
+from `0.0421` to `0.3461`.  Thus depth makes the matrix-level mismatch
+astronomically worse without producing a commensurate optimization-cost
+explosion.  The direct comparison is saved in
+`mismatch_vs_cost_by_depth.{csv,png}`.
+
+Across the full 480-row grid:
+
+| predictor | MAE | relative MAE vs mean measured gap | correlation |
+|---|---:|---:|---:|
+| frozen `K(0)` | 0.06693 | 28.6% | 0.97166 |
+| linear operator velocity | 0.01043 | 4.45% | 0.99492 |
+| early retangent | 0.01394 | 5.95% | 0.99830 |
+
+The frozen operator increasingly underpredicts progress at larger depth.  Two
+snapshots recover most of the magnitude, supporting the scoped statement that
+the later cost depends on operator evolution.  The finite-time result remains
+FA-only.
+
+Output: `outputs/aaai_depth_experiments/finite_time_rows.csv`,
+`finite_time_metrics.csv`, and the corresponding plots.
+
+### Parameter-matched supplement
+
+The fixed-width result is accompanied by an equal-parameter-count control.
+Using depth-2 width-64 as the target (`P=5376`, no biases), the closest
+uniform widths are `{269,64,47,39,31}` for depths `{1,2,3,4,6}`, with all
+parameter counts within 49 of the target.  Each cell contains 4 forward
+initializations x 8 feedback initializations; error bars again use the four
+forward-initialization means.
+
+| depth | width | parameters | measured gap |
+|---:|---:|---:|---:|
+| 1 | 269 | 5380 | 0.02695 |
+| 2 | 64 | 5376 | 0.15270 |
+| 3 | 47 | 5358 | 0.18857 |
+| 4 | 39 | 5343 | 0.37599 |
+| 6 | 31 | 5425 | 0.21024 |
+
+The gap is not monotone in depth once parameter count is fixed: it peaks at
+depth 4 and falls at depth 6.  This is another reason not to infer optimization
+cost directly from the number of independently mismatched matrices.  Output:
+`outputs/parameter_matched_depth_control/`.
+
+## 5. Learning-rate alternative explanation
+
+Configuration:
+
+- fixed width 64, depth `{1,2,3,4,6}`, `N=128`, 200 steps;
+- baseline BP/FA learning rate `1e-3`;
+- compensated FA rate `1e-3 / rho`, with `rho` read from one BP backward pass;
+- independent BP and FA rates selected using 4 tuning initializations x 4 FA
+  feedback seeds from
+  `{0.00025,0.0005,0.001,0.002,0.004,0.008,0.016,0.032,0.064,0.128,0.256}`;
+- stability-confirmed evaluation on 6 new forward initializations x 8 feedback
+  initializations; all 720 reported evaluation rows stable.
+
+Learning-rate curve uncertainty is likewise computed across the six forward
+initializations, with eight feedback draws nested within each one.
+
+| depth | same-rate gap | compensated gap | reduction | independently tuned gap |
+|---:|---:|---:|---:|---:|
+| 1 | 0.05906 | 0.04053 | 31.4% | 0.56940 |
+| 2 | 0.19916 | 0.09729 | 51.2% | 0.64210 |
+| 3 | 0.27346 | 0.16263 | 40.5% | 1.13002 |
+| 4 | 0.29340 | 0.18442 | 37.1% | 1.30812 |
+| 6 | 0.40664 | 0.28850 | 29.1% | 1.75739 |
+
+Initialization compensation removes a material part of the same-rate gap but
+does not remove it at any tested depth.  Independent tuning does not rescue
+FA: the selected BP rates are `0.064-0.128`, whereas the selected FA rates
+fall from `0.064` at depth 1 to stability-confirmed `0.004` at depth 6.  BP therefore makes much
+more progress in the same 200-step budget.  This result should be described as
+a fixed-horizon optimization comparison, not a claim about asymptotic
+convergence.
+
+Output: `outputs/learning_rate_compensation/`.
+
+## 6. Cross-entropy CNN initialization validation
+
+Configuration:
+
+- MNIST subset, 128 examples resized to `8x8`;
+- `Conv(1,8,3x3) -> ReLU -> Conv(8,12,3x3) -> ReLU -> global average pool -> linear`;
+- cross-entropy loss;
+- 4 forward initializations;
+- 2048 real-backward feedback draws per initialization and rule (FA and DFA).
+
+The predicted deficit spans `0.592-0.816`.  Across the eight rows, max
+absolute calibration error is `0.008309`, and max standardized error is only
+`1.308` Monte Carlo standard errors.  This validates the initialization
+identity outside squared-loss MLPs without making a finite-time CNN claim.
+
+Output: `outputs/cnn_initialization_validation/`.
+
+## 7. Figures
+
+- `outputs/aaai_main_figures/figure1_mismatch_vs_initial_cost.{png,pdf}`:
+  BP/FA/DFA schematic; minimax mismatch plus the 250x null-model
+  overprediction; joint initialization calibration across rule, depth, MNIST
+  MLP, and cross-entropy CNN.
+- `outputs/aaai_main_figures/figure2_training_dynamics.{png,pdf}`:
+  measured finite-time cost over depth/width; frozen-operator calibration;
+  two-snapshot calibration.
+- Supplemental learning-rate figures:
+  `outputs/learning_rate_compensation/learning_rate_control_by_depth.png` and
+  `learning_rate_tuning_curves.png`.
+- Direct depth contrast:
+  `outputs/aaai_depth_experiments/mismatch_vs_cost_by_depth.{csv,png}`.
+- Parameter-matched supplement:
+  `outputs/parameter_matched_depth_control/parameter_matched_gap_by_depth.png`.
+
+## 8. Claim boundary after these experiments
+
+Supported:
+
+- The exact expected initialization result applies to both FA and DFA under
+  independent zero-mean feedback, and is validated through depth 6.
+- It also holds for a cross-entropy CNN in the tested convolutional feedback
+  implementation.
+- Simple initial-rate compensation reduces but does not eliminate the
+  finite-time FA/BP gap in the tested MLP setting.
+- FA's later gap and frozen-model error grow with depth; two early operator
+  snapshots remain predictive.
+
+Not supported / not claimed:
+
+- DFA finite-time dynamics.
+- CNN finite-time training equivalence or modern large-scale scaling.
+- A universal absence of width/depth thresholds.
+- A per-eigenmode erosion theorem.
+
+## 9. Reproduction commands
+
+```bash
+python scripts/aaai_depth_experiments.py --part init \
+  --depths 1 2 3 4 6 --width 64 --train-samples 64 \
+  --init-seeds 6 --feedback-draws 256 --torch-threads 16 \
+  --outdir outputs/aaai_depth_experiments
+
+python scripts/aaai_depth_experiments.py --part finite \
+  --depths 1 2 3 4 6 --finite-widths 32 64 96 --train-samples 64 \
+  --finite-init-seeds 4 --finite-feedback-seeds 8 \
+  --lr 0.001 --horizon 50 --early-step 20 --torch-threads 16 \
+  --outdir outputs/aaai_depth_experiments
+
+python scripts/aaai_depth_experiments.py --part summarize \
+  --outdir outputs/aaai_depth_experiments
+
+python scripts/parameter_matched_depth_control.py \
+  --depths 1 2 3 4 6 --reference-depth 2 --reference-width 64 \
+  --train-samples 64 --init-seeds 4 --feedback-seeds 8 \
+  --lr 0.001 --horizon 50 --early-step 20 --torch-threads 16 \
+  --outdir outputs/parameter_matched_depth_control
+
+python scripts/learning_rate_compensation.py \
+  --depths 1 2 3 4 6 --width 64 --train-samples 128 --steps 200 \
+  --base-lr 0.001 --tune-init-seeds 4 --tune-feedback-seeds 4 \
+  --eval-init-seeds 6 --eval-feedback-seeds 8 --torch-threads 16 \
+  --outdir outputs/learning_rate_compensation
+
+python scripts/cnn_initialization_validation.py \
+  --train-samples 128 --image-size 8 --channels 8 12 \
+  --init-seeds 4 --feedback-draws 2048 --torch-threads 8 \
+  --outdir outputs/cnn_initialization_validation
+
+python scripts/plot_aaai_main_figures.py \
+  --outdir outputs/aaai_main_figures
+```
diff --git a/notes/README.md b/notes/README.md
index 853c122..021b966 100644
--- a/notes/README.md
+++ b/notes/README.md
@@ -7,6 +7,7 @@ capacity & operator-gap analysis)** — target AAAI-27.
 
 | doc | role |
 |---|---|
+| [43_aaai_accept_bar_experiments.md](43_aaai_accept_bar_experiments.md) | **new focused-paper experiments**: FA/DFA depth calibration, finite-time depth grid, learning-rate controls, CE CNN validation, and two composite figures |
 | [42_master_technical_reference.md](42_master_technical_reference.md) | **write the paper from this**: every theorem with full statement + proof, every result with exact numbers/configs, derivations, negative results, honesty boundaries -- self-contained |
 | [41_paper_plan.md](41_paper_plan.md) | **the plan**: venue rules, narrative, section/page budget, theorem & figure numbering, related-work list, schedule, risks |
 | [36_evidence_ledger.md](36_evidence_ledger.md) | **the numbers, condensed**: every claim with its corrected headline figure, claim discipline, booby traps |
@@ -73,6 +74,11 @@ Phase 6 — wrap-up for AAAI-27 (2026-06-09..13):
   test gap**; 40 reproduction manifest; 41 paper plan; 42 **master technical
   reference** (full proofs + numbers, the write-from doc).
 
+Phase 7 — focused accept-bar controls (2026-07-21):
+- 43 **FA/DFA depth calibration, FA finite-time depth/width grid, learning-rate
+  compensation and independent tuning, cross-entropy CNN initialization
+  validation, and focused Figures 1-2**.
+
 ## Note conventions
 
 - Prefer formulas and precise claims over prose-only ideas.
diff --git a/paper/figures/aaai/cnn_initialization_calibration.png b/paper/figures/aaai/cnn_initialization_calibration.png
new file mode 100644
index 0000000..7312fe6
--- /dev/null
+++ b/paper/figures/aaai/cnn_initialization_calibration.png
diff --git a/paper/figures/aaai/figure1_mismatch_vs_initial_cost.png b/paper/figures/aaai/figure1_mismatch_vs_initial_cost.png
new file mode 100644
index 0000000..1195734
--- /dev/null
+++ b/paper/figures/aaai/figure1_mismatch_vs_initial_cost.png
diff --git a/paper/figures/aaai/figure2_training_dynamics.png b/paper/figures/aaai/figure2_training_dynamics.png
new file mode 100644
index 0000000..ddaf7d4
--- /dev/null
+++ b/paper/figures/aaai/figure2_training_dynamics.png
diff --git a/paper/figures/aaai/learning_rate_control_by_depth.png b/paper/figures/aaai/learning_rate_control_by_depth.png
new file mode 100644
index 0000000..d363922
--- /dev/null
+++ b/paper/figures/aaai/learning_rate_control_by_depth.png
diff --git a/paper/figures/aaai/parameter_matched_gap_by_depth.png b/paper/figures/aaai/parameter_matched_gap_by_depth.png
new file mode 100644
index 0000000..954423f
--- /dev/null
+++ b/paper/figures/aaai/parameter_matched_gap_by_depth.png
diff --git a/scripts/README.md b/scripts/README.md
index 51dcbe3..c7c5ead 100644
--- a/scripts/README.md
+++ b/scripts/README.md
@@ -20,6 +20,12 @@ the recipes further down this file cover the static-geometry scripts).
 | `finite_time_kernel_diagnostic.py` | oracle time-varying K(t) vs fixed K(0) | 13 |
 | `real_data_validation.py` | MNIST: T3 + estimator (has `--self-test`) | 38 |
 | `plot_teacher_test_gap.py` | teacher train-vs-test gap figure | 39 |
+| `feedback_rules.py`, `test_feedback_rules.py` | shared arbitrary-depth BP/FA/DFA implementation + numerical checks | 43 |
+| `aaai_depth_experiments.py` | FA/DFA initialization depth calibration + FA finite-time depth/width grid | 43 |
+| `parameter_matched_depth_control.py` | equal-parameter-count finite-time depth control for the supplement | 43 |
+| `learning_rate_compensation.py` | same-rate, initialization-compensated, and independently tuned controls | 43 |
+| `cnn_initialization_validation.py` | FA/DFA initialization theorem on a cross-entropy CNN | 43 |
+| `plot_aaai_main_figures.py` | two focused three-panel main figures | 43 |
 
 **Diagnostics / negative results (kept for the record):**
 `fa_tangent_kernel_capacity.py` (kernel capacity probe, note 07),
diff --git a/scripts/aaai_depth_experiments.py b/scripts/aaai_depth_experiments.py
new file mode 100644
index 0000000..bb52f85
--- /dev/null
+++ b/scripts/aaai_depth_experiments.py
@@ -0,0 +1,633 @@
+#!/usr/bin/env python3
+"""AAAI depth experiments for fixed random feedback.
+
+Part ``init`` validates the exact expected initial first-order deficit for both
+FA and DFA at hidden depths 1, 2, 3, 4, and 6.  Part ``finite`` measures the
+FA/BP finite-time loss gap and evaluates the frozen-initialization, linear
+operator-velocity, and early-retangent predictors over depth and width.
+"""
+
+from __future__ import annotations
+
+import argparse
+import csv
+import json
+import math
+import sys
+import time
+from dataclasses import asdict, dataclass
+from pathlib import Path
+
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+import torch
+
+SCRIPT_DIR = Path(__file__).resolve().parent
+if str(SCRIPT_DIR) not in sys.path:
+    sys.path.insert(0, str(SCRIPT_DIR))
+
+import feedback_rules as fr  # noqa: E402
+import real_data_validation as rdv  # noqa: E402
+
+
+@dataclass(frozen=True)
+class InitRow:
+    rule: str
+    depth: int
+    width: int
+    init_seed: int
+    feedback_draws: int
+    bp_speed: float
+    output_speed: float
+    predicted_deficit: float
+    empirical_deficit: float
+    empirical_stderr: float
+    empirical_std: float
+    calibration_error: float
+
+
+@dataclass(frozen=True)
+class FiniteRow:
+    depth: int
+    width: int
+    train_samples: int
+    init_seed: int
+    feedback_seed: int
+    lr: float
+    horizon: int
+    early_step: int
+    initial_deficit_prediction: float
+    empirical_bp_loss: float
+    empirical_fa_loss: float
+    empirical_gap: float
+    fixed_gap: float
+    velocity_gap: float
+    retangent_gap: float
+    fixed_error: float
+    velocity_error: float
+    retangent_error: float
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument(
+        "--part", choices=["init", "finite", "both", "summarize"], default="both"
+    )
+    parser.add_argument("--depths", type=int, nargs="+", default=[1, 2, 3, 4, 6])
+    parser.add_argument("--width", type=int, default=64, help="fixed width for initialization")
+    parser.add_argument("--finite-widths", type=int, nargs="+", default=[32, 64, 96])
+    parser.add_argument("--input-dim", type=int, default=16)
+    parser.add_argument("--output-dim", type=int, default=4)
+    parser.add_argument("--train-samples", type=int, default=64)
+    parser.add_argument("--data-seed", type=int, default=2027)
+    parser.add_argument("--init-seeds", type=int, default=6)
+    parser.add_argument("--feedback-draws", type=int, default=256)
+    parser.add_argument("--finite-init-seeds", type=int, default=3)
+    parser.add_argument("--finite-feedback-seeds", type=int, default=4)
+    parser.add_argument("--lr", type=float, default=1e-3)
+    parser.add_argument("--horizon", type=int, default=50)
+    parser.add_argument("--early-step", type=int, default=20)
+    parser.add_argument("--feedback-scale", choices=["relu", "fan-in", "unit"], default="relu")
+    parser.add_argument("--torch-threads", type=int, default=16)
+    parser.add_argument("--outdir", type=Path, default=Path("outputs/aaai_depth_experiments"))
+    return parser.parse_args()
+
+
+def make_data(args: argparse.Namespace) -> tuple[torch.Tensor, torch.Tensor]:
+    generator = torch.Generator().manual_seed(args.data_seed)
+    x = torch.randn(args.train_samples, args.input_dim, generator=generator, dtype=torch.float64)
+    y = torch.randn(args.train_samples, args.output_dim, generator=generator, dtype=torch.float64)
+    return x, y
+
+
+def write_dataclass_csv(path: Path, rows: list[object]) -> None:
+    if not rows:
+        return
+    path.parent.mkdir(parents=True, exist_ok=True)
+    first = asdict(rows[0])  # type: ignore[arg-type]
+    with path.open("w", newline="") as handle:
+        writer = csv.DictWriter(handle, fieldnames=list(first))
+        writer.writeheader()
+        for row in rows:
+            writer.writerow(asdict(row))  # type: ignore[arg-type]
+
+
+def run_initialization(args: argparse.Namespace, x: torch.Tensor, y: torch.Tensor) -> list[InitRow]:
+    rows: list[InitRow] = []
+    for depth in args.depths:
+        dims = [args.input_dim, *([args.width] * depth), args.output_dim]
+        for init_index in range(args.init_seeds):
+            init_seed = 10_000 + init_index
+            weights = fr.initialize_mlp(dims, seed=init_seed)
+            bp_grads = fr.gradients(weights, x, y)
+            bp_speed = fr.squared_norm(bp_grads)
+            output_speed = fr.squared_norm([bp_grads[-1]])
+            predicted = 1.0 - output_speed / bp_speed
+            for rule_index, rule in enumerate(("fa", "dfa")):
+                deficits: list[float] = []
+                for draw in range(args.feedback_draws):
+                    feedback_seed = 100_000_000 + 10_000 * depth + 1_000 * init_index + 2 * draw + rule_index
+                    feedback = fr.init_feedback(
+                        dims,
+                        seed=feedback_seed,
+                        rule=rule,
+                        mode=args.feedback_scale,
+                    )
+                    rule_grads = fr.gradients(weights, x, y, rule=rule, feedback=feedback)
+                    deficits.append(1.0 - fr.inner_product(bp_grads, rule_grads) / bp_speed)
+                values = np.asarray(deficits)
+                empirical = float(values.mean())
+                row = InitRow(
+                    rule=rule.upper(),
+                    depth=depth,
+                    width=args.width,
+                    init_seed=init_seed,
+                    feedback_draws=args.feedback_draws,
+                    bp_speed=bp_speed,
+                    output_speed=output_speed,
+                    predicted_deficit=predicted,
+                    empirical_deficit=empirical,
+                    empirical_stderr=float(values.std(ddof=1) / math.sqrt(len(values))),
+                    empirical_std=float(values.std(ddof=1)),
+                    calibration_error=empirical - predicted,
+                )
+                rows.append(row)
+                print(
+                    f"[init] {rule.upper()} depth={depth} init={init_index}: "
+                    f"prediction={predicted:.5f}, measured={empirical:.5f} "
+                    f"+/- {row.empirical_stderr:.5f}",
+                    flush=True,
+                )
+    return rows
+
+
+def train_with_snapshots(
+    weights0: list[torch.Tensor],
+    x: torch.Tensor,
+    y: torch.Tensor,
+    lr: float,
+    horizon: int,
+    feedback: list[torch.Tensor] | None,
+    early_step: int,
+) -> tuple[list[torch.Tensor], list[torch.Tensor]]:
+    final, snapshots = rdv.train_with_snapshots(
+        weights0,
+        x,
+        y,
+        lr,
+        horizon,
+        feedback,
+        {early_step},
+    )
+    return final, snapshots[early_step]
+
+
+def retangent_loss(
+    kernel_s: torch.Tensor,
+    residual_s: torch.Tensor,
+    lr: float,
+    samples: int,
+    horizon: int,
+    early_step: int,
+) -> float:
+    return rdv.rollout_loss(
+        kernel_s,
+        None,
+        residual_s,
+        lr,
+        samples,
+        horizon - early_step,
+        early_step,
+    )
+
+
+def run_finite_time(args: argparse.Namespace, x: torch.Tensor, y: torch.Tensor) -> list[FiniteRow]:
+    rows: list[FiniteRow] = []
+    n = x.shape[0]
+    for depth in args.depths:
+        for width in args.finite_widths:
+            dims = [args.input_dim, *([width] * depth), args.output_dim]
+            for init_index in range(args.finite_init_seeds):
+                block_start = time.time()
+                init_seed = 20_000 + init_index
+                weights0 = fr.initialize_mlp(dims, seed=init_seed)
+                r0 = (fr.predict(weights0, x) - y).reshape(-1)
+                bp_grads = fr.gradients(weights0, x, y)
+                initial_prediction = 1.0 - fr.squared_norm([bp_grads[-1]]) / fr.squared_norm(bp_grads)
+
+                k_bp0 = rdv.tangent_kernel(weights0, x, None)
+                bp_final, bp_early = train_with_snapshots(
+                    weights0, x, y, args.lr, args.horizon, None, args.early_step
+                )
+                k_bps = rdv.tangent_kernel(bp_early, x, None)
+                bp_emp = fr.mse(bp_final, x, y)
+                bp_fixed = rdv.rollout_loss(
+                    k_bp0, None, r0, args.lr, n, args.horizon, args.early_step
+                )
+                bp_velocity = rdv.rollout_loss(
+                    k_bp0,
+                    k_bps - k_bp0,
+                    r0,
+                    args.lr,
+                    n,
+                    args.horizon,
+                    args.early_step,
+                )
+                bp_residual_s = (fr.predict(bp_early, x) - y).reshape(-1)
+                bp_retangent = retangent_loss(
+                    k_bps, bp_residual_s, args.lr, n, args.horizon, args.early_step
+                )
+
+                for feedback_index in range(args.finite_feedback_seeds):
+                    feedback_seed = 200_000 + 1_000 * init_index + feedback_index
+                    feedback = fr.init_feedback(
+                        dims,
+                        seed=feedback_seed,
+                        rule="fa",
+                        mode=args.feedback_scale,
+                    )
+                    k_fa0 = rdv.tangent_kernel(weights0, x, feedback)
+                    fa_final, fa_early = train_with_snapshots(
+                        weights0,
+                        x,
+                        y,
+                        args.lr,
+                        args.horizon,
+                        feedback,
+                        args.early_step,
+                    )
+                    k_fas = rdv.tangent_kernel(fa_early, x, feedback)
+                    fa_emp = fr.mse(fa_final, x, y)
+                    fa_fixed = rdv.rollout_loss(
+                        k_fa0, None, r0, args.lr, n, args.horizon, args.early_step
+                    )
+                    fa_velocity = rdv.rollout_loss(
+                        k_fa0,
+                        k_fas - k_fa0,
+                        r0,
+                        args.lr,
+                        n,
+                        args.horizon,
+                        args.early_step,
+                    )
+                    fa_residual_s = (fr.predict(fa_early, x) - y).reshape(-1)
+                    fa_retangent = retangent_loss(
+                        k_fas, fa_residual_s, args.lr, n, args.horizon, args.early_step
+                    )
+
+                    empirical_gap = fa_emp - bp_emp
+                    fixed_gap = fa_fixed - bp_fixed
+                    velocity_gap = fa_velocity - bp_velocity
+                    retangent_gap = fa_retangent - bp_retangent
+                    rows.append(
+                        FiniteRow(
+                            depth=depth,
+                            width=width,
+                            train_samples=n,
+                            init_seed=init_seed,
+                            feedback_seed=feedback_seed,
+                            lr=args.lr,
+                            horizon=args.horizon,
+                            early_step=args.early_step,
+                            initial_deficit_prediction=initial_prediction,
+                            empirical_bp_loss=bp_emp,
+                            empirical_fa_loss=fa_emp,
+                            empirical_gap=empirical_gap,
+                            fixed_gap=fixed_gap,
+                            velocity_gap=velocity_gap,
+                            retangent_gap=retangent_gap,
+                            fixed_error=fixed_gap - empirical_gap,
+                            velocity_error=velocity_gap - empirical_gap,
+                            retangent_error=retangent_gap - empirical_gap,
+                        )
+                    )
+                print(
+                    f"[finite] depth={depth} width={width} init={init_index}: "
+                    f"{args.finite_feedback_seeds} feedbacks in {time.time() - block_start:.1f}s",
+                    flush=True,
+                )
+    return rows
+
+
+def safe_corr(left: np.ndarray, right: np.ndarray) -> float:
+    if len(left) < 2 or float(np.std(left) * np.std(right)) == 0.0:
+        return math.nan
+    return float(np.corrcoef(left, right)[0, 1])
+
+
+def summarize_finite(rows: list[FiniteRow]) -> list[dict[str, float | int | str]]:
+    metrics: list[dict[str, float | int | str]] = []
+    groups = sorted({(row.depth, row.width) for row in rows})
+    for depth, width in groups:
+        subset = [row for row in rows if row.depth == depth and row.width == width]
+        empirical = np.asarray([row.empirical_gap for row in subset])
+        for predictor, field in (
+            ("frozen K(0)", "fixed_gap"),
+            ("linear velocity", "velocity_gap"),
+            ("early retangent", "retangent_gap"),
+        ):
+            predicted = np.asarray([getattr(row, field) for row in subset])
+            metrics.append(
+                {
+                    "depth": depth,
+                    "width": width,
+                    "predictor": predictor,
+                    "rows": len(subset),
+                    "empirical_gap_mean": float(empirical.mean()),
+                    "empirical_gap_std": float(empirical.std(ddof=1)),
+                    "mae": float(np.mean(np.abs(predicted - empirical))),
+                    "bias": float(np.mean(predicted - empirical)),
+                    "corr": safe_corr(predicted, empirical),
+                }
+            )
+    return metrics
+
+
+def write_dict_csv(path: Path, rows: list[dict[str, object]]) -> None:
+    if not rows:
+        return
+    with path.open("w", newline="") as handle:
+        writer = csv.DictWriter(handle, fieldnames=list(rows[0]))
+        writer.writeheader()
+        writer.writerows(rows)
+
+
+def plot_initial(rows: list[InitRow], outdir: Path) -> None:
+    fig, ax = plt.subplots(figsize=(5.7, 5.2), dpi=180)
+    colors = plt.cm.viridis(np.linspace(0.1, 0.9, len(sorted({row.depth for row in rows}))))
+    for color, depth in zip(colors, sorted({row.depth for row in rows})):
+        for rule, marker in (("FA", "o"), ("DFA", "s")):
+            subset = [row for row in rows if row.depth == depth and row.rule == rule]
+            ax.errorbar(
+                [row.predicted_deficit for row in subset],
+                [row.empirical_deficit for row in subset],
+                yerr=[2 * row.empirical_stderr for row in subset],
+                fmt=marker,
+                ms=4.5,
+                color=color,
+                alpha=0.78,
+                capsize=1.5,
+                lw=0.8,
+                label=f"depth {depth}, {rule}",
+            )
+    values = [value for row in rows for value in (row.predicted_deficit, row.empirical_deficit)]
+    lo, hi = min(values), max(values)
+    pad = 0.04 * (hi - lo + 1e-12)
+    ax.plot([lo - pad, hi + pad], [lo - pad, hi + pad], color="black", lw=1)
+    ax.set_xlabel("exact expected initial deficit")
+    ax.set_ylabel("measured mean initial deficit")
+    ax.set_title("FA and DFA share the exact initialization cost")
+    ax.grid(alpha=0.16)
+    ax.legend(fontsize=6.5, ncols=2)
+    fig.tight_layout()
+    fig.savefig(outdir / "initial_deficit_calibration.png", bbox_inches="tight")
+    plt.close(fig)
+
+
+def plot_finite(rows: list[FiniteRow], metrics: list[dict[str, object]], outdir: Path) -> None:
+    fig, ax = plt.subplots(figsize=(6.4, 4.7), dpi=180)
+    for width, marker in zip(sorted({row.width for row in rows}), ("o", "s", "^", "D")):
+        depths = sorted({row.depth for row in rows if row.width == width})
+        means = [np.mean([row.empirical_gap for row in rows if row.width == width and row.depth == d]) for d in depths]
+        sems = []
+        for depth in depths:
+            cell = [row for row in rows if row.width == width and row.depth == depth]
+            init_means = [
+                np.mean([row.empirical_gap for row in cell if row.init_seed == init_seed])
+                for init_seed in sorted({row.init_seed for row in cell})
+            ]
+            sems.append(float(np.std(init_means, ddof=1) / math.sqrt(len(init_means))))
+        ax.errorbar(depths, means, yerr=sems, marker=marker, capsize=3, label=f"width {width}")
+    ax.set_xlabel("hidden-layer depth")
+    ax.set_ylabel("FA loss - BP loss at fixed horizon")
+    ax.set_title("Finite-time optimization cost across depth and width")
+    ax.grid(alpha=0.16)
+    ax.legend()
+    fig.tight_layout()
+    fig.savefig(outdir / "finite_gap_by_depth_width.png", bbox_inches="tight")
+    plt.close(fig)
+
+    fig, ax = plt.subplots(figsize=(7.0, 4.8), dpi=180)
+    predictors = ["frozen K(0)", "linear velocity", "early retangent"]
+    for predictor, marker in zip(predictors, ("^", "o", "s")):
+        subset = [row for row in metrics if row["predictor"] == predictor and row["width"] == 64]
+        ax.plot(
+            [int(row["depth"]) for row in subset],
+            [float(row["mae"]) for row in subset],
+            marker=marker,
+            label=predictor,
+        )
+    ax.set_xlabel("hidden-layer depth")
+    ax.set_ylabel("gap prediction MAE (width 64)")
+    ax.set_title("Operator drift, not initialization mismatch, controls later error")
+    ax.set_yscale("log")
+    ax.grid(alpha=0.16, which="both")
+    ax.legend()
+    fig.tight_layout()
+    fig.savefig(outdir / "predictor_mae_by_depth.png", bbox_inches="tight")
+    plt.close(fig)
+
+
+def summarize_mismatch_vs_cost(
+    init_rows: list[dict[str, str]],
+    finite_rows: list[dict[str, str]],
+    width: int,
+    output_dim: int,
+    outdir: Path,
+) -> list[dict[str, object]]:
+    """Contrast exponentially compounding matrix mismatch with measured costs."""
+    depths = sorted({int(row["depth"]) for row in init_rows})
+    records: list[dict[str, object]] = []
+    for depth in depths:
+        init_subset = [row for row in init_rows if int(row["depth"]) == depth]
+        finite_subset = [
+            row
+            for row in finite_rows
+            if int(row["depth"]) == depth and int(row["width"]) == width
+        ]
+        # One width x output block plus depth-1 width x width FA blocks.
+        joint_alignment_proxy = (1.0 / (width * output_dim)) * (
+            1.0 / (width * width)
+        ) ** (depth - 1)
+        records.append(
+            {
+                "depth": depth,
+                "width": width,
+                "joint_alignment_proxy": joint_alignment_proxy,
+                "negative_log10_joint_alignment": -math.log10(joint_alignment_proxy),
+                "initial_deficit_mean": float(
+                    np.mean([float(row["empirical_deficit"]) for row in init_subset])
+                ),
+                "finite_gap_mean": float(
+                    np.mean([float(row["empirical_gap"]) for row in finite_subset])
+                ),
+            }
+        )
+
+    write_dict_csv(outdir / "mismatch_vs_cost_by_depth.csv", records)
+    fig, ax = plt.subplots(figsize=(7.0, 4.8), dpi=180)
+    ax.plot(
+        [int(row["depth"]) for row in records],
+        [float(row["joint_alignment_proxy"]) for row in records],
+        "o--",
+        color="black",
+        label=r"joint matrix-alignment proxy $\prod_l 1/D_l$",
+    )
+    ax.set_yscale("log")
+    ax.set_xlabel("hidden-layer depth")
+    ax.set_ylabel("joint squared-alignment proxy")
+    ax.grid(alpha=0.16, which="both")
+    ax2 = ax.twinx()
+    ax2.plot(
+        [int(row["depth"]) for row in records],
+        [float(row["initial_deficit_mean"]) for row in records],
+        "s-",
+        color="#2f6f9f",
+        label="measured initial deficit",
+    )
+    ax2.plot(
+        [int(row["depth"]) for row in records],
+        [float(row["finite_gap_mean"]) for row in records],
+        "^-",
+        color="#c65f16",
+        label="measured finite-time gap",
+    )
+    ax2.set_ylabel("optimization cost")
+    handles1, labels1 = ax.get_legend_handles_labels()
+    handles2, labels2 = ax2.get_legend_handles_labels()
+    ax2.legend(handles1 + handles2, labels1 + labels2, fontsize=8, loc="center right")
+    ax.set_title("Matrix mismatch compounds exponentially; optimization cost does not")
+    fig.tight_layout()
+    fig.savefig(outdir / "mismatch_vs_cost_by_depth.png", bbox_inches="tight")
+    plt.close(fig)
+    return records
+
+
+def main() -> None:
+    args = parse_args()
+    if args.early_step <= 0 or args.early_step >= args.horizon:
+        raise ValueError("early-step must lie strictly between 0 and horizon")
+    torch.set_num_threads(args.torch_threads)
+    args.outdir.mkdir(parents=True, exist_ok=True)
+    x, y = make_data(args)
+    summary: dict[str, object] = {
+        "config": {key: str(value) if isinstance(value, Path) else value for key, value in vars(args).items()}
+    }
+
+    if args.part == "summarize":
+        with (args.outdir / "initialization_rows.csv").open() as handle:
+            existing_init = list(csv.DictReader(handle))
+        with (args.outdir / "finite_time_rows.csv").open() as handle:
+            existing_finite = list(csv.DictReader(handle))
+        init_errors = np.asarray([float(row["calibration_error"]) for row in existing_init])
+        init_stderr = np.asarray([float(row["empirical_stderr"]) for row in existing_init])
+        summary["initialization"] = {
+            "rows": len(existing_init),
+            "max_abs_calibration_error": float(np.max(np.abs(init_errors))),
+            "max_standardized_error": float(np.max(np.abs(init_errors) / init_stderr)),
+            "rules": sorted({row["rule"] for row in existing_init}),
+            "depths": sorted({int(row["depth"]) for row in existing_init}),
+        }
+        empirical = np.asarray([float(row["empirical_gap"]) for row in existing_finite])
+        summary["finite_time"] = {
+            "rows": len(existing_finite),
+            "depths": sorted({int(row["depth"]) for row in existing_finite}),
+            "widths": sorted({int(row["width"]) for row in existing_finite}),
+            "overall": {
+                predictor: {
+                    "mae": float(np.mean(np.abs(predicted - empirical))),
+                    "corr": safe_corr(predicted, empirical),
+                }
+                for predictor, predicted in (
+                    (
+                        "frozen",
+                        np.asarray([float(row["fixed_gap"]) for row in existing_finite]),
+                    ),
+                    (
+                        "linear_velocity",
+                        np.asarray([float(row["velocity_gap"]) for row in existing_finite]),
+                    ),
+                    (
+                        "early_retangent",
+                        np.asarray([float(row["retangent_gap"]) for row in existing_finite]),
+                    ),
+                )
+            },
+        }
+        mismatch_records = summarize_mismatch_vs_cost(
+            existing_init,
+            existing_finite,
+            width=args.width,
+            output_dim=args.output_dim,
+            outdir=args.outdir,
+        )
+        summary["mismatch_vs_cost"] = {
+            "joint_alignment_proxy_depth_1": mismatch_records[0]["joint_alignment_proxy"],
+            "joint_alignment_proxy_depth_6": mismatch_records[-1]["joint_alignment_proxy"],
+            "orders_of_magnitude_drop": float(
+                math.log10(
+                    float(mismatch_records[0]["joint_alignment_proxy"])
+                    / float(mismatch_records[-1]["joint_alignment_proxy"])
+                )
+            ),
+        }
+        (args.outdir / "summary.json").write_text(json.dumps(summary, indent=2) + "\n")
+        print(f"summary: {args.outdir / 'summary.json'}")
+        return
+
+    if args.part in {"init", "both"}:
+        init_rows = run_initialization(args, x, y)
+        write_dataclass_csv(args.outdir / "initialization_rows.csv", init_rows)
+        plot_initial(init_rows, args.outdir)
+        summary["initialization"] = {
+            "rows": len(init_rows),
+            "max_abs_calibration_error": max(abs(row.calibration_error) for row in init_rows),
+            "max_standardized_error": max(
+                abs(row.calibration_error) / row.empirical_stderr for row in init_rows
+            ),
+            "rules": sorted({row.rule for row in init_rows}),
+            "depths": sorted({row.depth for row in init_rows}),
+        }
+
+    if args.part in {"finite", "both"}:
+        finite_rows = run_finite_time(args, x, y)
+        write_dataclass_csv(args.outdir / "finite_time_rows.csv", finite_rows)
+        metrics = summarize_finite(finite_rows)
+        write_dict_csv(args.outdir / "finite_time_metrics.csv", metrics)  # type: ignore[arg-type]
+        plot_finite(finite_rows, metrics, args.outdir)
+        summary["finite_time"] = {
+            "rows": len(finite_rows),
+            "depths": sorted({row.depth for row in finite_rows}),
+            "widths": sorted({row.width for row in finite_rows}),
+            "overall": {
+                predictor: {
+                    "mae": float(
+                        np.mean(
+                            [
+                                abs(getattr(row, field) - row.empirical_gap)
+                                for row in finite_rows
+                            ]
+                        )
+                    ),
+                    "corr": safe_corr(
+                        np.asarray([getattr(row, field) for row in finite_rows]),
+                        np.asarray([row.empirical_gap for row in finite_rows]),
+                    ),
+                }
+                for predictor, field in (
+                    ("frozen", "fixed_gap"),
+                    ("linear_velocity", "velocity_gap"),
+                    ("early_retangent", "retangent_gap"),
+                )
+            },
+        }
+
+    (args.outdir / "summary.json").write_text(json.dumps(summary, indent=2) + "\n")
+    print(f"summary: {args.outdir / 'summary.json'}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/scripts/cnn_initialization_validation.py b/scripts/cnn_initialization_validation.py
new file mode 100644
index 0000000..dcf4f92
--- /dev/null
+++ b/scripts/cnn_initialization_validation.py
@@ -0,0 +1,320 @@
+#!/usr/bin/env python3
+"""Cross-entropy CNN validation of the exact FA/DFA initialization cost.
+
+The model has two 3x3 convolutional ReLU layers, global average pooling, and a
+linear classifier.  The experiment runs the real convolutional FA and DFA
+backward rules for every feedback draw and compares their mean first-order
+loss-decrease deficit with the hidden-parameter share from one BP backward
+pass.  This directly tests that the initialization result is not a
+squared-loss MLP identity.
+"""
+
+from __future__ import annotations
+
+import argparse
+import csv
+import json
+import math
+from dataclasses import asdict, dataclass
+from pathlib import Path
+
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+import torch
+import torch.nn.functional as functional
+
+
+Tensor = torch.Tensor
+
+
+@dataclass(frozen=True)
+class CNNRow:
+    architecture: str
+    loss: str
+    rule: str
+    init_seed: int
+    feedback_draws: int
+    bp_speed: float
+    output_speed: float
+    predicted_deficit: float
+    empirical_deficit: float
+    empirical_stderr: float
+    empirical_std: float
+    calibration_error: float
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--data-root", type=Path, default=Path("data"))
+    parser.add_argument("--train-samples", type=int, default=128)
+    parser.add_argument("--image-size", type=int, default=8)
+    parser.add_argument("--channels", type=int, nargs=2, default=[8, 12])
+    parser.add_argument("--init-seeds", type=int, default=4)
+    parser.add_argument("--feedback-draws", type=int, default=512)
+    parser.add_argument("--torch-threads", type=int, default=16)
+    parser.add_argument("--self-test", action="store_true")
+    parser.add_argument(
+        "--outdir", type=Path, default=Path("outputs/cnn_initialization_validation")
+    )
+    return parser.parse_args()
+
+
+def load_mnist(root: Path, samples: int, image_size: int) -> tuple[Tensor, Tensor]:
+    from torchvision import datasets
+
+    dataset = datasets.MNIST(root=str(root), train=True, download=True)
+    generator = torch.Generator().manual_seed(2027)
+    indices = torch.randperm(len(dataset), generator=generator)[:samples]
+    images = dataset.data[indices].to(torch.float64).unsqueeze(1) / 255.0
+    images = functional.interpolate(
+        images, size=(image_size, image_size), mode="bilinear", align_corners=False
+    )
+    images = (images - images.mean()) / (images.std() + 1e-12)
+    labels = dataset.targets[indices]
+    return images, labels
+
+
+def initialize_weights(channels: tuple[int, int], classes: int, seed: int) -> list[Tensor]:
+    c1, c2 = channels
+    generator = torch.Generator().manual_seed(seed)
+    w1 = torch.randn(c1, 1, 3, 3, generator=generator, dtype=torch.float64) * math.sqrt(2 / 9)
+    w2 = torch.randn(c2, c1, 3, 3, generator=generator, dtype=torch.float64) * math.sqrt(
+        2 / (9 * c1)
+    )
+    w3 = torch.randn(classes, c2, generator=generator, dtype=torch.float64) / math.sqrt(c2)
+    return [w1, w2, w3]
+
+
+def forward(weights: list[Tensor], x: Tensor) -> tuple[Tensor, tuple[Tensor, ...]]:
+    w1, w2, w3 = weights
+    z1 = functional.conv2d(x, w1, padding=1)
+    h1 = torch.relu(z1)
+    z2 = functional.conv2d(h1, w2, padding=1)
+    h2 = torch.relu(z2)
+    pooled = h2.mean(dim=(2, 3))
+    logits = pooled @ w3.T
+    return logits, (z1, h1, z2, h2, pooled)
+
+
+def autograd_gradients(weights: list[Tensor], x: Tensor, labels: Tensor) -> list[Tensor]:
+    leaves = [weight.clone().detach().requires_grad_(True) for weight in weights]
+    logits, _ = forward(leaves, x)
+    loss = functional.cross_entropy(logits, labels)
+    return list(torch.autograd.grad(loss, leaves))
+
+
+def init_feedback(
+    weights: list[Tensor], classes: int, seed: int, rule: str
+) -> list[Tensor]:
+    w1, w2, _w3 = weights
+    c1, c2 = w1.shape[0], w2.shape[0]
+    generator = torch.Generator().manual_seed(seed)
+    if rule == "fa":
+        conv_map = torch.randn(w2.shape, generator=generator, dtype=torch.float64) * math.sqrt(
+            2 / (9 * c1)
+        )
+        output_map = torch.randn(
+            classes, c2, generator=generator, dtype=torch.float64
+        ) * math.sqrt(2 / c2)
+        return [conv_map, output_map]
+    if rule == "dfa":
+        direct_first = torch.randn(
+            classes, c1, generator=generator, dtype=torch.float64
+        ) * math.sqrt(2 / c1)
+        direct_second = torch.randn(
+            classes, c2, generator=generator, dtype=torch.float64
+        ) * math.sqrt(2 / c2)
+        return [direct_first, direct_second]
+    raise ValueError("rule must be fa or dfa")
+
+
+def manual_gradients(
+    weights: list[Tensor],
+    x: Tensor,
+    labels: Tensor,
+    rule: str,
+    feedback: list[Tensor] | None = None,
+) -> list[Tensor]:
+    if rule not in {"bp", "fa", "dfa"}:
+        raise ValueError("rule must be bp, fa, or dfa")
+    if rule != "bp" and feedback is None:
+        raise ValueError("FA/DFA require feedback maps")
+
+    w1, w2, w3 = weights
+    logits, (z1, h1, z2, h2, pooled) = forward(weights, x)
+    batch, height, width = h2.shape[0], h2.shape[2], h2.shape[3]
+    probabilities = torch.softmax(logits, dim=1)
+    probabilities[torch.arange(batch), labels] -= 1.0
+    delta_output = probabilities / batch
+    grad_w3 = delta_output.T @ pooled
+
+    if rule == "bp":
+        delta_pooled = delta_output @ w3
+        delta_z2 = delta_pooled[:, :, None, None].expand_as(h2) / (height * width)
+        delta_z2 = delta_z2 * (z2 > 0)
+        delta_h1 = functional.conv_transpose2d(delta_z2, w2, padding=1)
+        delta_z1 = delta_h1 * (z1 > 0)
+    elif rule == "fa":
+        assert feedback is not None
+        conv_map, output_map = feedback
+        delta_pooled = delta_output @ output_map
+        delta_z2 = delta_pooled[:, :, None, None].expand_as(h2) / (height * width)
+        delta_z2 = delta_z2 * (z2 > 0)
+        delta_h1 = functional.conv_transpose2d(delta_z2, conv_map, padding=1)
+        delta_z1 = delta_h1 * (z1 > 0)
+    else:
+        assert feedback is not None
+        direct_first, direct_second = feedback
+        direct_z2 = delta_output @ direct_second
+        delta_z2 = direct_z2[:, :, None, None].expand_as(h2) / (height * width)
+        delta_z2 = delta_z2 * (z2 > 0)
+        direct_z1 = delta_output @ direct_first
+        delta_z1 = direct_z1[:, :, None, None].expand_as(z1) / (height * width)
+        delta_z1 = delta_z1 * (z1 > 0)
+
+    grad_w2 = torch.nn.grad.conv2d_weight(h1, w2.shape, delta_z2, padding=1)
+    grad_w1 = torch.nn.grad.conv2d_weight(x, w1.shape, delta_z1, padding=1)
+    return [grad_w1, grad_w2, grad_w3]
+
+
+def squared_norm(grads: list[Tensor]) -> float:
+    return float(sum(torch.sum(grad * grad) for grad in grads))
+
+
+def inner_product(left: list[Tensor], right: list[Tensor]) -> float:
+    return float(sum(torch.sum(a * b) for a, b in zip(left, right)))
+
+
+def self_test() -> None:
+    generator = torch.Generator().manual_seed(31)
+    x = torch.randn(5, 1, 6, 6, generator=generator, dtype=torch.float64)
+    labels = torch.randint(0, 3, (5,), generator=generator)
+    weights = initialize_weights((4, 5), classes=3, seed=41)
+    automatic = autograd_gradients(weights, x, labels)
+    manual = manual_gradients(weights, x, labels, rule="bp")
+    error = max(float((a - b).abs().max()) for a, b in zip(automatic, manual))
+    print(f"CNN BP manual/autograd max error: {error:.3e}")
+    assert error < 1e-11
+    for rule in ("fa", "dfa"):
+        feedback = init_feedback(weights, classes=3, seed=51, rule=rule)
+        rule_grads = manual_gradients(weights, x, labels, rule=rule, feedback=feedback)
+        output_error = float((rule_grads[-1] - manual[-1]).abs().max())
+        print(f"CNN {rule.upper()} output-gradient error: {output_error:.3e}")
+        assert output_error < 1e-12
+    print("CNN initialization self-test PASSED")
+
+
+def run(args: argparse.Namespace, x: Tensor, labels: Tensor) -> list[CNNRow]:
+    rows: list[CNNRow] = []
+    classes = 10
+    channels = (args.channels[0], args.channels[1])
+    for init_index in range(args.init_seeds):
+        init_seed = 10_000 + init_index
+        weights = initialize_weights(channels, classes, init_seed)
+        bp = manual_gradients(weights, x, labels, rule="bp")
+        bp_speed = squared_norm(bp)
+        output_speed = squared_norm([bp[-1]])
+        prediction = 1.0 - output_speed / bp_speed
+        for rule_index, rule in enumerate(("fa", "dfa")):
+            deficits: list[float] = []
+            for draw in range(args.feedback_draws):
+                feedback = init_feedback(
+                    weights,
+                    classes,
+                    seed=100_000 + 10_000 * init_index + 2 * draw + rule_index,
+                    rule=rule,
+                )
+                rule_grads = manual_gradients(weights, x, labels, rule, feedback)
+                deficits.append(1.0 - inner_product(bp, rule_grads) / bp_speed)
+            values = np.asarray(deficits)
+            empirical = float(values.mean())
+            row = CNNRow(
+                architecture=f"Conv({channels[0]},{channels[1]})-GAP-linear",
+                loss="cross_entropy",
+                rule=rule.upper(),
+                init_seed=init_seed,
+                feedback_draws=args.feedback_draws,
+                bp_speed=bp_speed,
+                output_speed=output_speed,
+                predicted_deficit=prediction,
+                empirical_deficit=empirical,
+                empirical_stderr=float(values.std(ddof=1) / math.sqrt(len(values))),
+                empirical_std=float(values.std(ddof=1)),
+                calibration_error=empirical - prediction,
+            )
+            rows.append(row)
+            print(
+                f"[CNN] {rule.upper()} init={init_index}: prediction={prediction:.5f}, "
+                f"measured={empirical:.5f} +/- {row.empirical_stderr:.5f}",
+                flush=True,
+            )
+    return rows
+
+
+def write_rows(path: Path, rows: list[CNNRow]) -> None:
+    with path.open("w", newline="") as handle:
+        writer = csv.DictWriter(handle, fieldnames=list(CNNRow.__annotations__))
+        writer.writeheader()
+        for row in rows:
+            writer.writerow(asdict(row))
+
+
+def plot(rows: list[CNNRow], outdir: Path) -> None:
+    fig, ax = plt.subplots(figsize=(5.5, 5.0), dpi=180)
+    for rule, marker, color in (("FA", "o", "#2f6f9f"), ("DFA", "s", "#c65f16")):
+        subset = [row for row in rows if row.rule == rule]
+        ax.errorbar(
+            [row.predicted_deficit for row in subset],
+            [row.empirical_deficit for row in subset],
+            yerr=[2 * row.empirical_stderr for row in subset],
+            fmt=marker,
+            color=color,
+            capsize=2,
+            label=rule,
+        )
+    values = [value for row in rows for value in (row.predicted_deficit, row.empirical_deficit)]
+    lo, hi = min(values), max(values)
+    pad = 0.05 * (hi - lo + 1e-12)
+    ax.plot([lo - pad, hi + pad], [lo - pad, hi + pad], color="black", lw=1)
+    ax.set_xlabel("exact expected initial deficit")
+    ax.set_ylabel("measured mean initial deficit")
+    ax.set_title("Cross-entropy CNN: FA and DFA initialization cost")
+    ax.grid(alpha=0.16)
+    ax.legend()
+    fig.tight_layout()
+    fig.savefig(outdir / "cnn_initialization_calibration.png", bbox_inches="tight")
+    plt.close(fig)
+
+
+def main() -> None:
+    args = parse_args()
+    torch.set_num_threads(args.torch_threads)
+    if args.self_test:
+        self_test()
+        return
+    args.outdir.mkdir(parents=True, exist_ok=True)
+    x, labels = load_mnist(args.data_root, args.train_samples, args.image_size)
+    self_test()
+    rows = run(args, x, labels)
+    write_rows(args.outdir / "cnn_initialization_rows.csv", rows)
+    plot(rows, args.outdir)
+    payload = {
+        "config": {key: str(value) if isinstance(value, Path) else value for key, value in vars(args).items()},
+        "rows": len(rows),
+        "max_abs_calibration_error": max(abs(row.calibration_error) for row in rows),
+        "max_standardized_error": max(
+            abs(row.calibration_error) / row.empirical_stderr for row in rows
+        ),
+        "rules": sorted({row.rule for row in rows}),
+        "architecture": rows[0].architecture,
+        "loss": rows[0].loss,
+    }
+    (args.outdir / "summary.json").write_text(json.dumps(payload, indent=2) + "\n")
+    print(f"summary: {args.outdir / 'summary.json'}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/scripts/feedback_rules.py b/scripts/feedback_rules.py
new file mode 100644
index 0000000..b658d1d
--- /dev/null
+++ b/scripts/feedback_rules.py
@@ -0,0 +1,187 @@
+#!/usr/bin/env python3
+"""Shared BP, FA, and DFA utilities for arbitrary-depth MLP experiments.
+
+The older experiment scripts predate the DFA validation and expose only the
+layerwise-FA rule.  This module keeps the new paper-critical experiments on a
+single, explicitly tested implementation.  All losses use
+``0.5 * mean(sum((prediction - target)**2, dim=1))``.
+"""
+
+from __future__ import annotations
+
+import math
+from collections.abc import Sequence
+
+import torch
+
+
+Tensor = torch.Tensor
+Rule = str
+
+
+def initialize_mlp(dims: Sequence[int], seed: int, device: str = "cpu") -> list[Tensor]:
+    """He-initialize hidden layers and fan-in initialize the linear output."""
+    if len(dims) < 3:
+        raise ValueError("dims must contain input, at least one hidden layer, and output")
+    generator = torch.Generator(device=device)
+    generator.manual_seed(seed)
+    weights: list[Tensor] = []
+    for layer, (fan_in, fan_out) in enumerate(zip(dims[:-1], dims[1:])):
+        scale = math.sqrt(2.0 / fan_in) if layer < len(dims) - 2 else 1.0 / math.sqrt(fan_in)
+        weights.append(
+            torch.randn(
+                fan_out,
+                fan_in,
+                generator=generator,
+                dtype=torch.float64,
+                device=device,
+            )
+            * scale
+        )
+    return weights
+
+
+def clone_weights(weights: Sequence[Tensor]) -> list[Tensor]:
+    return [weight.clone() for weight in weights]
+
+
+def forward(weights: Sequence[Tensor], x: Tensor) -> tuple[list[Tensor], list[Tensor]]:
+    activations = [x]
+    preacts: list[Tensor] = []
+    current = x
+    for layer, weight in enumerate(weights):
+        current = current @ weight.T
+        preacts.append(current)
+        if layer < len(weights) - 1:
+            current = torch.relu(current)
+        activations.append(current)
+    return activations, preacts
+
+
+def predict(weights: Sequence[Tensor], x: Tensor) -> Tensor:
+    return forward(weights, x)[0][-1]
+
+
+def mse(weights: Sequence[Tensor], x: Tensor, y: Tensor) -> float:
+    error = predict(weights, x) - y
+    return float(0.5 * torch.mean(torch.sum(error * error, dim=1)))
+
+
+def _feedback_scale(rows: int, mode: str) -> float:
+    if mode == "relu":
+        return math.sqrt(2.0 / rows)
+    if mode == "fan-in":
+        return math.sqrt(1.0 / rows)
+    if mode == "unit":
+        return 1.0
+    raise ValueError(f"unknown feedback scale: {mode}")
+
+
+def init_feedback(
+    dims: Sequence[int],
+    seed: int,
+    rule: Rule,
+    mode: str = "relu",
+    device: str = "cpu",
+) -> list[Tensor]:
+    """Initialize independent zero-mean feedback maps for FA or DFA.
+
+    FA map ``l`` has shape ``(dims[l+1], dims[l+2])`` and replaces the
+    transpose action of the next forward matrix.  DFA map ``l`` has shape
+    ``(dims[l+1], dims[-1])`` and maps output error directly to hidden layer
+    ``l``.  The output layer never receives a feedback matrix.
+    """
+    if rule not in {"fa", "dfa"}:
+        raise ValueError("feedback rule must be 'fa' or 'dfa'")
+    generator = torch.Generator(device=device)
+    generator.manual_seed(seed)
+    shapes = (
+        [(dims[layer + 1], dims[layer + 2]) for layer in range(len(dims) - 2)]
+        if rule == "fa"
+        else [(dims[layer + 1], dims[-1]) for layer in range(len(dims) - 2)]
+    )
+    return [
+        torch.randn(rows, cols, generator=generator, dtype=torch.float64, device=device)
+        * _feedback_scale(rows, mode)
+        for rows, cols in shapes
+    ]
+
+
+def gradients(
+    weights: Sequence[Tensor],
+    x: Tensor,
+    y: Tensor,
+    rule: Rule = "bp",
+    feedback: Sequence[Tensor] | None = None,
+) -> list[Tensor]:
+    """Return BP gradients or FA/DFA pseudo-gradients for an MLP."""
+    if rule not in {"bp", "fa", "dfa"}:
+        raise ValueError("rule must be one of: bp, fa, dfa")
+    if rule == "bp" and feedback is not None:
+        raise ValueError("BP does not take feedback maps")
+    if rule != "bp" and (feedback is None or len(feedback) != len(weights) - 1):
+        raise ValueError(f"{rule.upper()} needs one feedback map per hidden layer")
+
+    activations, preacts = forward(weights, x)
+    batch = x.shape[0]
+    deltas: list[Tensor] = [torch.empty(0, dtype=x.dtype, device=x.device) for _ in weights]
+    deltas[-1] = (activations[-1] - y) / batch
+
+    if rule == "dfa":
+        assert feedback is not None
+        output_delta = deltas[-1]
+        for layer in range(len(weights) - 2, -1, -1):
+            deltas[layer] = (output_delta @ feedback[layer].T) * (preacts[layer] > 0)
+    else:
+        for layer in range(len(weights) - 2, -1, -1):
+            if rule == "bp":
+                back = deltas[layer + 1] @ weights[layer + 1]
+            else:
+                assert feedback is not None
+                back = deltas[layer + 1] @ feedback[layer].T
+            deltas[layer] = back * (preacts[layer] > 0)
+
+    return [delta.T @ activations[layer] for layer, delta in enumerate(deltas)]
+
+
+def train(
+    weights0: Sequence[Tensor],
+    x: Tensor,
+    y: Tensor,
+    lr: float,
+    steps: int,
+    rule: Rule = "bp",
+    feedback: Sequence[Tensor] | None = None,
+) -> list[Tensor]:
+    weights = clone_weights(weights0)
+    for _ in range(steps):
+        grads = gradients(weights, x, y, rule=rule, feedback=feedback)
+        weights = [weight - lr * grad for weight, grad in zip(weights, grads)]
+        if not all(bool(torch.isfinite(weight).all()) for weight in weights):
+            break
+    return weights
+
+
+def squared_norm(grads: Sequence[Tensor]) -> float:
+    return float(sum(torch.sum(grad * grad) for grad in grads))
+
+
+def inner_product(left: Sequence[Tensor], right: Sequence[Tensor]) -> float:
+    return float(sum(torch.sum(a * b) for a, b in zip(left, right)))
+
+
+def first_order_statistics(
+    bp_grads: Sequence[Tensor], rule_grads: Sequence[Tensor]
+) -> tuple[float, float, float, float]:
+    """Return BP speed, output speed, rule/BP decrease ratio, and deficit."""
+    bp_speed = squared_norm(bp_grads)
+    output_speed = squared_norm([bp_grads[-1]])
+    ratio = inner_product(bp_grads, rule_grads) / bp_speed
+    return bp_speed, output_speed, ratio, 1.0 - ratio
+
+
+def autograd_bp_gradients(weights: Sequence[Tensor], x: Tensor, y: Tensor) -> list[Tensor]:
+    leaves = [weight.clone().detach().requires_grad_(True) for weight in weights]
+    error = predict(leaves, x) - y
+    loss = 0.5 * torch.mean(torch.sum(error * error, dim=1))
+    return list(torch.autograd.grad(loss, leaves))
diff --git a/scripts/learning_rate_compensation.py b/scripts/learning_rate_compensation.py
new file mode 100644
index 0000000..10e1a64
--- /dev/null
+++ b/scripts/learning_rate_compensation.py
@@ -0,0 +1,488 @@
+#!/usr/bin/env python3
+"""Learning-rate controls for the finite-time FA/BP optimization gap.
+
+For every depth, the script first tunes BP and FA learning rates on held-out
+initialization/feedback seeds.  It then evaluates three FA conditions against
+BP at its tuned rate:
+
+1. the same learning rate as BP;
+2. initialization compensation ``eta_FA = eta_BP / rho``, where ``rho`` is
+   the exact expected fraction of BP's first-order decrease retained by FA;
+3. FA's independently tuned stable learning rate.
+"""
+
+from __future__ import annotations
+
+import argparse
+import csv
+import json
+import math
+import sys
+from dataclasses import asdict, dataclass
+from pathlib import Path
+
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+import torch
+
+SCRIPT_DIR = Path(__file__).resolve().parent
+if str(SCRIPT_DIR) not in sys.path:
+    sys.path.insert(0, str(SCRIPT_DIR))
+
+import feedback_rules as fr  # noqa: E402
+
+
+@dataclass(frozen=True)
+class TuneRow:
+    depth: int
+    rule: str
+    lr: float
+    trials: int
+    stable_trials: int
+    mean_final_loss: float
+    median_final_loss: float
+
+
+@dataclass(frozen=True)
+class EvaluationRow:
+    depth: int
+    width: int
+    init_seed: int
+    feedback_seed: int
+    regime: str
+    steps: int
+    bp_lr: float
+    fa_lr: float
+    rho: float
+    initial_loss: float
+    bp_final_loss: float
+    fa_final_loss: float
+    gap_to_bp: float
+    stable: bool
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--depths", type=int, nargs="+", default=[1, 2, 3, 4, 6])
+    parser.add_argument("--width", type=int, default=64)
+    parser.add_argument("--input-dim", type=int, default=16)
+    parser.add_argument("--output-dim", type=int, default=4)
+    parser.add_argument("--train-samples", type=int, default=128)
+    parser.add_argument("--steps", type=int, default=200)
+    parser.add_argument("--base-lr", type=float, default=1e-3)
+    parser.add_argument(
+        "--lr-grid",
+        type=float,
+        nargs="+",
+        default=[
+            0.00025,
+            0.0005,
+            0.001,
+            0.002,
+            0.004,
+            0.008,
+            0.016,
+            0.032,
+            0.064,
+            0.128,
+            0.256,
+        ],
+    )
+    parser.add_argument("--tune-init-seeds", type=int, default=2)
+    parser.add_argument("--tune-feedback-seeds", type=int, default=2)
+    parser.add_argument("--eval-init-seeds", type=int, default=4)
+    parser.add_argument("--eval-feedback-seeds", type=int, default=4)
+    parser.add_argument("--data-seed", type=int, default=4242)
+    parser.add_argument("--feedback-scale", choices=["relu", "fan-in", "unit"], default="relu")
+    parser.add_argument("--stability-factor", type=float, default=10.0)
+    parser.add_argument("--torch-threads", type=int, default=16)
+    parser.add_argument("--postprocess-only", action="store_true")
+    parser.add_argument("--outdir", type=Path, default=Path("outputs/learning_rate_compensation"))
+    return parser.parse_args()
+
+
+def make_data(args: argparse.Namespace) -> tuple[torch.Tensor, torch.Tensor]:
+    generator = torch.Generator().manual_seed(args.data_seed)
+    x = torch.randn(args.train_samples, args.input_dim, generator=generator, dtype=torch.float64)
+    y = torch.randn(args.train_samples, args.output_dim, generator=generator, dtype=torch.float64)
+    return x, y
+
+
+def stable_loss(
+    weights: list[torch.Tensor],
+    x: torch.Tensor,
+    y: torch.Tensor,
+    initial_loss: float,
+    stability_factor: float,
+) -> tuple[float, bool]:
+    loss = fr.mse(weights, x, y)
+    stable = math.isfinite(loss) and loss <= stability_factor * initial_loss
+    return loss, stable
+
+
+def tune_depth(
+    args: argparse.Namespace,
+    depth: int,
+    x: torch.Tensor,
+    y: torch.Tensor,
+) -> tuple[float, float, list[TuneRow]]:
+    dims = [args.input_dim, *([args.width] * depth), args.output_dim]
+    records: list[TuneRow] = []
+    for rule in ("bp", "fa"):
+        for lr in args.lr_grid:
+            losses: list[float] = []
+            stable_flags: list[bool] = []
+            for init_index in range(args.tune_init_seeds):
+                init_seed = 700_000 + init_index
+                weights0 = fr.initialize_mlp(dims, init_seed)
+                initial_loss = fr.mse(weights0, x, y)
+                feedback_count = 1 if rule == "bp" else args.tune_feedback_seeds
+                for feedback_index in range(feedback_count):
+                    feedback = None
+                    if rule == "fa":
+                        feedback = fr.init_feedback(
+                            dims,
+                            710_000 + 1_000 * init_index + feedback_index,
+                            rule="fa",
+                            mode=args.feedback_scale,
+                        )
+                    final = fr.train(
+                        weights0,
+                        x,
+                        y,
+                        lr,
+                        args.steps,
+                        rule=rule,
+                        feedback=feedback,
+                    )
+                    loss, stable = stable_loss(
+                        final, x, y, initial_loss, args.stability_factor
+                    )
+                    losses.append(loss if math.isfinite(loss) else math.inf)
+                    stable_flags.append(stable)
+            stable_losses = [loss for loss, stable in zip(losses, stable_flags) if stable]
+            records.append(
+                TuneRow(
+                    depth=depth,
+                    rule=rule.upper(),
+                    lr=lr,
+                    trials=len(losses),
+                    stable_trials=sum(stable_flags),
+                    mean_final_loss=float(np.mean(stable_losses)) if stable_losses else math.inf,
+                    median_final_loss=float(np.median(stable_losses)) if stable_losses else math.inf,
+                )
+            )
+            print(
+                f"[tune] depth={depth} {rule.upper()} lr={lr:g}: "
+                f"stable={sum(stable_flags)}/{len(losses)}, "
+                f"mean={records[-1].mean_final_loss:.5g}",
+                flush=True,
+            )
+
+    best: dict[str, float] = {}
+    for rule in ("BP", "FA"):
+        candidates = [
+            row
+            for row in records
+            if row.rule == rule and row.stable_trials == row.trials and math.isfinite(row.mean_final_loss)
+        ]
+        if not candidates:
+            raise RuntimeError(f"no fully stable {rule} learning rate at depth {depth}")
+        best[rule] = min(candidates, key=lambda row: row.mean_final_loss).lr
+    print(
+        f"[tune] depth={depth}: selected BP={best['BP']:g}, FA={best['FA']:g}",
+        flush=True,
+    )
+    return best["BP"], best["FA"], records
+
+
+def evaluate_depth(
+    args: argparse.Namespace,
+    depth: int,
+    bp_lr: float,
+    fa_lr: float,
+    x: torch.Tensor,
+    y: torch.Tensor,
+) -> list[EvaluationRow]:
+    dims = [args.input_dim, *([args.width] * depth), args.output_dim]
+    rows: list[EvaluationRow] = []
+    for init_index in range(args.eval_init_seeds):
+        init_seed = 10_000 + init_index
+        weights0 = fr.initialize_mlp(dims, init_seed)
+        initial_loss = fr.mse(weights0, x, y)
+        bp_grads = fr.gradients(weights0, x, y)
+        rho = fr.squared_norm([bp_grads[-1]]) / fr.squared_norm(bp_grads)
+        bp_losses: dict[str, tuple[float, bool]] = {}
+        for name, rate in (("base", args.base_lr), ("tuned", bp_lr)):
+            bp_final = fr.train(weights0, x, y, rate, args.steps, rule="bp")
+            bp_losses[name] = stable_loss(
+                bp_final, x, y, initial_loss, args.stability_factor
+            )
+            if not bp_losses[name][1]:
+                raise RuntimeError(
+                    f"{name} BP rate became unstable at depth={depth}, init={init_seed}"
+                )
+
+        regimes = (
+            ("same_lr", args.base_lr, args.base_lr, "base"),
+            ("initial_compensation", args.base_lr, args.base_lr / rho, "base"),
+            ("independently_tuned", bp_lr, fa_lr, "tuned"),
+        )
+        for feedback_index in range(args.eval_feedback_seeds):
+            feedback_seed = 100_000 + 1_000 * init_index + feedback_index
+            feedback = fr.init_feedback(
+                dims,
+                feedback_seed,
+                rule="fa",
+                mode=args.feedback_scale,
+            )
+            for regime, bp_rate, fa_rate, bp_key in regimes:
+                fa_final = fr.train(
+                    weights0,
+                    x,
+                    y,
+                    fa_rate,
+                    args.steps,
+                    rule="fa",
+                    feedback=feedback,
+                )
+                fa_final_loss, stable = stable_loss(
+                    fa_final, x, y, initial_loss, args.stability_factor
+                )
+                rows.append(
+                    EvaluationRow(
+                        depth=depth,
+                        width=args.width,
+                        init_seed=init_seed,
+                        feedback_seed=feedback_seed,
+                        regime=regime,
+                        steps=args.steps,
+                        bp_lr=bp_rate,
+                        fa_lr=fa_rate,
+                        rho=rho,
+                        initial_loss=initial_loss,
+                        bp_final_loss=bp_losses[bp_key][0],
+                        fa_final_loss=fa_final_loss,
+                        gap_to_bp=fa_final_loss - bp_losses[bp_key][0],
+                        stable=stable,
+                    )
+                )
+        print(
+            f"[eval] depth={depth} init={init_index}: rho={rho:.4f}, "
+            f"compensation={args.base_lr / rho:.5g}",
+            flush=True,
+        )
+    return rows
+
+
+def write_dataclass_csv(path: Path, rows: list[object]) -> None:
+    if not rows:
+        return
+    first = asdict(rows[0])  # type: ignore[arg-type]
+    with path.open("w", newline="") as handle:
+        writer = csv.DictWriter(handle, fieldnames=list(first))
+        writer.writeheader()
+        for row in rows:
+            writer.writerow(asdict(row))  # type: ignore[arg-type]
+
+
+def summarize(rows: list[EvaluationRow]) -> list[dict[str, object]]:
+    result: list[dict[str, object]] = []
+    for depth in sorted({row.depth for row in rows}):
+        for regime in ("same_lr", "initial_compensation", "independently_tuned"):
+            subset = [row for row in rows if row.depth == depth and row.regime == regime]
+            gaps = np.asarray([row.gap_to_bp for row in subset])
+            fa_losses = np.asarray([row.fa_final_loss for row in subset])
+            init_gap_means = np.asarray(
+                [
+                    np.mean(
+                        [row.gap_to_bp for row in subset if row.init_seed == init_seed]
+                    )
+                    for init_seed in sorted({row.init_seed for row in subset})
+                ]
+            )
+            result.append(
+                {
+                    "depth": depth,
+                    "regime": regime,
+                    "rows": len(subset),
+                    "stable_rows": sum(row.stable for row in subset),
+                    "bp_lr": subset[0].bp_lr,
+                    "fa_lr_mean": float(np.mean([row.fa_lr for row in subset])),
+                    "rho_mean": float(np.mean([row.rho for row in subset])),
+                    "bp_loss_mean": float(np.mean([row.bp_final_loss for row in subset])),
+                    "fa_loss_mean": float(fa_losses.mean()),
+                    "gap_mean": float(gaps.mean()),
+                    "gap_std": float(init_gap_means.std(ddof=1)),
+                    "gap_sem": float(
+                        init_gap_means.std(ddof=1) / math.sqrt(len(init_gap_means))
+                    ),
+                }
+            )
+    return result
+
+
+def write_dict_csv(path: Path, rows: list[dict[str, object]]) -> None:
+    if not rows:
+        return
+    with path.open("w", newline="") as handle:
+        writer = csv.DictWriter(handle, fieldnames=list(rows[0]))
+        writer.writeheader()
+        writer.writerows(rows)
+
+
+def plot_results(summary: list[dict[str, object]], tune_rows: list[TuneRow], outdir: Path) -> None:
+    fig, ax = plt.subplots(figsize=(7.0, 4.8), dpi=180)
+    styles = {
+        "same_lr": ("o", "same learning rate"),
+        "initial_compensation": ("s", "initial-rate compensation"),
+        "independently_tuned": ("^", "independently tuned"),
+    }
+    for regime, (marker, label) in styles.items():
+        subset = [row for row in summary if row["regime"] == regime]
+        ax.errorbar(
+            [int(row["depth"]) for row in subset],
+            [float(row["gap_mean"]) for row in subset],
+            yerr=[float(row["gap_sem"]) for row in subset],
+            marker=marker,
+            capsize=3,
+            label=label,
+        )
+    ax.axhline(0, color="black", lw=1)
+    ax.set_xlabel("hidden-layer depth")
+    ax.set_ylabel("FA loss - BP loss at fixed steps")
+    ax.set_title("Finite-time gap under three learning-rate controls")
+    ax.grid(alpha=0.16)
+    ax.legend()
+    fig.tight_layout()
+    fig.savefig(outdir / "learning_rate_control_by_depth.png", bbox_inches="tight")
+    plt.close(fig)
+
+    fig, axes = plt.subplots(1, 2, figsize=(10.5, 4.2), dpi=180, sharey=True)
+    for ax, rule in zip(axes, ("BP", "FA")):
+        for depth in sorted({row.depth for row in tune_rows}):
+            subset = [row for row in tune_rows if row.rule == rule and row.depth == depth]
+            ax.plot(
+                [row.lr for row in subset],
+                [row.mean_final_loss for row in subset],
+                marker="o",
+                ms=3,
+                label=f"depth {depth}",
+            )
+        ax.set_xscale("log")
+        ax.set_yscale("log")
+        ax.set_xlabel("learning rate")
+        ax.set_title(rule)
+        ax.grid(alpha=0.16, which="both")
+    axes[0].set_ylabel("held-out tuning loss")
+    axes[1].legend(fontsize=7, ncols=2)
+    fig.tight_layout()
+    fig.savefig(outdir / "learning_rate_tuning_curves.png", bbox_inches="tight")
+    plt.close(fig)
+
+
+def main() -> None:
+    args = parse_args()
+    torch.set_num_threads(args.torch_threads)
+    args.outdir.mkdir(parents=True, exist_ok=True)
+    if args.postprocess_only:
+        with (args.outdir / "tuning_rows.csv").open() as handle:
+            tune_dicts = list(csv.DictReader(handle))
+        with (args.outdir / "evaluation_rows.csv").open() as handle:
+            eval_dicts = list(csv.DictReader(handle))
+        tune_rows = [
+            TuneRow(
+                depth=int(row["depth"]),
+                rule=row["rule"],
+                lr=float(row["lr"]),
+                trials=int(row["trials"]),
+                stable_trials=int(row["stable_trials"]),
+                mean_final_loss=float(row["mean_final_loss"]),
+                median_final_loss=float(row["median_final_loss"]),
+            )
+            for row in tune_dicts
+        ]
+        eval_rows = [
+            EvaluationRow(
+                depth=int(row["depth"]),
+                width=int(row["width"]),
+                init_seed=int(row["init_seed"]),
+                feedback_seed=int(row["feedback_seed"]),
+                regime=row["regime"],
+                steps=int(row["steps"]),
+                bp_lr=float(row["bp_lr"]),
+                fa_lr=float(row["fa_lr"]),
+                rho=float(row["rho"]),
+                initial_loss=float(row["initial_loss"]),
+                bp_final_loss=float(row["bp_final_loss"]),
+                fa_final_loss=float(row["fa_final_loss"]),
+                gap_to_bp=float(row["gap_to_bp"]),
+                stable=row["stable"] == "True",
+            )
+            for row in eval_dicts
+        ]
+        summary_rows = summarize(eval_rows)
+        write_dict_csv(args.outdir / "regime_summary.csv", summary_rows)
+        plot_results(summary_rows, tune_rows, args.outdir)
+        summary_path = args.outdir / "summary.json"
+        payload = json.loads(summary_path.read_text())
+        payload["regime_summary"] = summary_rows
+        payload["uncertainty"] = "SEM and SD are computed across forward-initialization means."
+        summary_path.write_text(json.dumps(payload, indent=2) + "\n")
+        print(f"summary: {summary_path}")
+        return
+    x, y = make_data(args)
+    all_tune: list[TuneRow] = []
+    all_eval: list[EvaluationRow] = []
+    selected: dict[str, dict[str, float]] = {}
+    for depth in args.depths:
+        bp_lr, fa_lr, tune_rows = tune_depth(args, depth, x, y)
+        all_tune.extend(tune_rows)
+        while True:
+            evaluation = evaluate_depth(args, depth, bp_lr, fa_lr, x, y)
+            tuned_rows = [row for row in evaluation if row.regime == "independently_tuned"]
+            if all(row.stable for row in tuned_rows):
+                break
+            lower_candidates = [
+                row
+                for row in tune_rows
+                if row.rule == "FA"
+                and row.lr < fa_lr
+                and row.stable_trials == row.trials
+                and math.isfinite(row.mean_final_loss)
+            ]
+            if not lower_candidates:
+                raise RuntimeError(
+                    f"no lower FA rate passed stability confirmation at depth {depth}"
+                )
+            previous = fa_lr
+            fa_lr = min(lower_candidates, key=lambda row: row.mean_final_loss).lr
+            print(
+                f"[stability] depth={depth}: FA lr {previous:g} failed on held-out "
+                f"evaluation seeds; retrying {fa_lr:g}",
+                flush=True,
+            )
+        selected[str(depth)] = {"bp_lr": bp_lr, "fa_lr": fa_lr}
+        all_eval.extend(evaluation)
+
+    write_dataclass_csv(args.outdir / "tuning_rows.csv", all_tune)
+    write_dataclass_csv(args.outdir / "evaluation_rows.csv", all_eval)
+    summary_rows = summarize(all_eval)
+    write_dict_csv(args.outdir / "regime_summary.csv", summary_rows)
+    plot_results(summary_rows, all_tune, args.outdir)
+    payload = {
+        "config": {key: str(value) if isinstance(value, Path) else value for key, value in vars(args).items()},
+        "selected_learning_rates": selected,
+        "rows": len(all_eval),
+        "regime_summary": summary_rows,
+        "uncertainty": "SEM and SD are computed across forward-initialization means.",
+    }
+    (args.outdir / "summary.json").write_text(json.dumps(payload, indent=2) + "\n")
+    print(f"summary: {args.outdir / 'summary.json'}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/scripts/parameter_matched_depth_control.py b/scripts/parameter_matched_depth_control.py
new file mode 100644
index 0000000..24f1f2b
--- /dev/null
+++ b/scripts/parameter_matched_depth_control.py
@@ -0,0 +1,174 @@
+#!/usr/bin/env python3
+"""Parameter-count-matched supplement for the FA finite-time depth sweep."""
+
+from __future__ import annotations
+
+import argparse
+import csv
+import json
+import math
+import sys
+from dataclasses import asdict
+from pathlib import Path
+from types import SimpleNamespace
+
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+import torch
+
+SCRIPT_DIR = Path(__file__).resolve().parent
+if str(SCRIPT_DIR) not in sys.path:
+    sys.path.insert(0, str(SCRIPT_DIR))
+
+import aaai_depth_experiments as depth_exp  # noqa: E402
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--depths", type=int, nargs="+", default=[1, 2, 3, 4, 6])
+    parser.add_argument("--reference-depth", type=int, default=2)
+    parser.add_argument("--reference-width", type=int, default=64)
+    parser.add_argument("--input-dim", type=int, default=16)
+    parser.add_argument("--output-dim", type=int, default=4)
+    parser.add_argument("--train-samples", type=int, default=64)
+    parser.add_argument("--data-seed", type=int, default=2027)
+    parser.add_argument("--init-seeds", type=int, default=4)
+    parser.add_argument("--feedback-seeds", type=int, default=8)
+    parser.add_argument("--lr", type=float, default=1e-3)
+    parser.add_argument("--horizon", type=int, default=50)
+    parser.add_argument("--early-step", type=int, default=20)
+    parser.add_argument("--feedback-scale", choices=["relu", "fan-in", "unit"], default="relu")
+    parser.add_argument("--torch-threads", type=int, default=16)
+    parser.add_argument(
+        "--outdir", type=Path, default=Path("outputs/parameter_matched_depth_control")
+    )
+    return parser.parse_args()
+
+
+def parameter_count(input_dim: int, output_dim: int, depth: int, width: int) -> int:
+    return input_dim * width + (depth - 1) * width * width + width * output_dim
+
+
+def closest_width(input_dim: int, output_dim: int, depth: int, target: int) -> int:
+    upper = max(2, math.ceil(target / (input_dim + output_dim)) + 2)
+    return min(
+        range(1, upper + 1),
+        key=lambda width: abs(parameter_count(input_dim, output_dim, depth, width) - target),
+    )
+
+
+def write_rows(path: Path, rows: list[depth_exp.FiniteRow]) -> None:
+    with path.open("w", newline="") as handle:
+        writer = csv.DictWriter(handle, fieldnames=list(depth_exp.FiniteRow.__annotations__))
+        writer.writeheader()
+        for row in rows:
+            writer.writerow(asdict(row))
+
+
+def main() -> None:
+    args = parse_args()
+    torch.set_num_threads(args.torch_threads)
+    args.outdir.mkdir(parents=True, exist_ok=True)
+    target = parameter_count(
+        args.input_dim, args.output_dim, args.reference_depth, args.reference_width
+    )
+    widths = {
+        depth: closest_width(args.input_dim, args.output_dim, depth, target)
+        for depth in args.depths
+    }
+    generator = torch.Generator().manual_seed(args.data_seed)
+    x = torch.randn(args.train_samples, args.input_dim, generator=generator, dtype=torch.float64)
+    y = torch.randn(args.train_samples, args.output_dim, generator=generator, dtype=torch.float64)
+
+    rows: list[depth_exp.FiniteRow] = []
+    for depth in args.depths:
+        run_args = SimpleNamespace(
+            depths=[depth],
+            finite_widths=[widths[depth]],
+            input_dim=args.input_dim,
+            output_dim=args.output_dim,
+            finite_init_seeds=args.init_seeds,
+            finite_feedback_seeds=args.feedback_seeds,
+            lr=args.lr,
+            horizon=args.horizon,
+            early_step=args.early_step,
+            feedback_scale=args.feedback_scale,
+        )
+        rows.extend(depth_exp.run_finite_time(run_args, x, y))
+
+    write_rows(args.outdir / "parameter_matched_rows.csv", rows)
+    metrics = depth_exp.summarize_finite(rows)
+    depth_exp.write_dict_csv(args.outdir / "parameter_matched_metrics.csv", metrics)  # type: ignore[arg-type]
+
+    summaries: list[dict[str, object]] = []
+    for depth in args.depths:
+        subset = [row for row in rows if row.depth == depth]
+        values = np.asarray([row.empirical_gap for row in subset])
+        init_means = np.asarray(
+            [
+                np.mean([row.empirical_gap for row in subset if row.init_seed == init_seed])
+                for init_seed in sorted({row.init_seed for row in subset})
+            ]
+        )
+        summaries.append(
+            {
+                "depth": depth,
+                "width": widths[depth],
+                "parameter_count": parameter_count(
+                    args.input_dim, args.output_dim, depth, widths[depth]
+                ),
+                "parameter_count_error": parameter_count(
+                    args.input_dim, args.output_dim, depth, widths[depth]
+                )
+                - target,
+                "rows": len(subset),
+                "empirical_gap_mean": float(values.mean()),
+                "empirical_gap_sem": float(
+                    init_means.std(ddof=1) / math.sqrt(len(init_means))
+                ),
+            }
+        )
+    depth_exp.write_dict_csv(args.outdir / "parameter_matched_summary.csv", summaries)
+
+    fig, ax = plt.subplots(figsize=(6.6, 4.7), dpi=180)
+    ax.errorbar(
+        [int(row["depth"]) for row in summaries],
+        [float(row["empirical_gap_mean"]) for row in summaries],
+        yerr=[float(row["empirical_gap_sem"]) for row in summaries],
+        marker="o",
+        capsize=3,
+        color="#6f4c9b",
+    )
+    for row in summaries:
+        ax.annotate(
+            f"w={row['width']}",
+            (int(row["depth"]), float(row["empirical_gap_mean"])),
+            xytext=(0, 9),
+            textcoords="offset points",
+            ha="center",
+            fontsize=8,
+        )
+    ax.set_xlabel("hidden-layer depth")
+    ax.set_ylabel("FA loss - BP loss at fixed horizon")
+    ax.set_title(f"Parameter-matched depth control (target P={target})")
+    ax.grid(alpha=0.16)
+    fig.tight_layout()
+    fig.savefig(args.outdir / "parameter_matched_gap_by_depth.png", bbox_inches="tight")
+    plt.close(fig)
+
+    payload = {
+        "config": {key: str(value) if isinstance(value, Path) else value for key, value in vars(args).items()},
+        "target_parameter_count": target,
+        "widths": {str(depth): width for depth, width in widths.items()},
+        "rows": len(rows),
+        "summary": summaries,
+    }
+    (args.outdir / "summary.json").write_text(json.dumps(payload, indent=2) + "\n")
+    print(f"summary: {args.outdir / 'summary.json'}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/scripts/plot_aaai_main_figures.py b/scripts/plot_aaai_main_figures.py
new file mode 100644
index 0000000..cd12011
--- /dev/null
+++ b/scripts/plot_aaai_main_figures.py
@@ -0,0 +1,311 @@
+#!/usr/bin/env python3
+"""Render the two three-panel composite figures for the focused AAAI story."""
+
+from __future__ import annotations
+
+import argparse
+import csv
+import math
+from pathlib import Path
+
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--depth-dir", type=Path, default=Path("outputs/aaai_depth_experiments"))
+    parser.add_argument(
+        "--cnn-dir", type=Path, default=Path("outputs/cnn_initialization_validation")
+    )
+    parser.add_argument(
+        "--mnist-dir", type=Path, default=Path("outputs/real_data_validation_mnist")
+    )
+    parser.add_argument(
+        "--null-dir", type=Path, default=Path("outputs/soft_erosion_theory_vs_empirical")
+    )
+    parser.add_argument("--outdir", type=Path, default=Path("outputs/aaai_main_figures"))
+    return parser.parse_args()
+
+
+def read_csv(path: Path) -> list[dict[str, str]]:
+    if not path.exists():
+        raise FileNotFoundError(path)
+    with path.open() as handle:
+        return list(csv.DictReader(handle))
+
+
+def panel_label(ax: plt.Axes, label: str) -> None:
+    ax.text(-0.13, 1.06, label, transform=ax.transAxes, fontsize=12, fontweight="bold")
+
+
+def draw_schematic(ax: plt.Axes) -> None:
+    ax.set_xlim(0, 1)
+    ax.set_ylim(0, 1)
+    ax.axis("off")
+    rows = [(0.82, "BP", "#2f6f9f"), (0.50, "FA", "#c65f16"), (0.18, "DFA", "#2b8a3e")]
+    x_positions = [0.18, 0.42, 0.66, 0.90]
+    labels = ["$h_1$", "$h_2$", "$h_3$", "error"]
+    for y, name, color in rows:
+        ax.text(0.01, y, name, va="center", ha="left", fontweight="bold", color=color)
+        for x, node in zip(x_positions, labels):
+            ax.add_patch(plt.Circle((x, y), 0.035, facecolor="white", edgecolor=color, lw=1.5))
+            ax.text(x, y - 0.075, node, ha="center", va="top", fontsize=8)
+        if name in {"BP", "FA"}:
+            arrow_label = r"$W^\top$" if name == "BP" else "$B$"
+            for right, left in zip(x_positions[:0:-1], x_positions[-2::-1]):
+                ax.annotate(
+                    "",
+                    xy=(left + 0.04, y),
+                    xytext=(right - 0.04, y),
+                    arrowprops=dict(arrowstyle="->", color=color, lw=1.5),
+                )
+                ax.text((left + right) / 2, y + 0.035, arrow_label, ha="center", fontsize=8)
+        else:
+            for left in x_positions[:-1]:
+                ax.annotate(
+                    "",
+                    xy=(left + 0.04, y),
+                    xytext=(x_positions[-1] - 0.04, y),
+                    arrowprops=dict(arrowstyle="->", color=color, lw=1.2, alpha=0.9),
+                )
+            ax.text(0.56, y + 0.055, "direct random maps", ha="center", fontsize=8)
+    ax.set_title("Fixed random feedback rules", fontsize=10)
+
+
+def figure_one(args: argparse.Namespace) -> Path:
+    init_rows = read_csv(args.depth_dir / "initialization_rows.csv")
+    cnn_rows = read_csv(args.cnn_dir / "cnn_initialization_rows.csv")
+    mnist_rows = read_csv(args.mnist_dir / "e0_mnist_rows.csv")
+    theory_rows = read_csv(args.null_dir / "theory_gap_summary.csv")
+    empirical_rows = read_csv(args.null_dir / "empirical_gap_summary.csv")
+
+    fig, axes = plt.subplots(1, 3, figsize=(14.2, 4.1), dpi=220)
+    draw_schematic(axes[0])
+    panel_label(axes[0], "a")
+
+    ax = axes[1]
+    theory_by_width = {int(row["width"]): float(row["mean"]) for row in theory_rows}
+    empirical_by_width = {int(row["width"]): float(row["mean"]) for row in empirical_rows}
+    widths = sorted(set(theory_by_width) & set(empirical_by_width))
+    ax.plot(
+        widths,
+        [theory_by_width[width] for width in widths],
+        "o-",
+        color="#8c6d31",
+        label="random-direction deletion",
+    )
+    ax.plot(
+        widths,
+        [empirical_by_width[width] for width in widths],
+        "s-",
+        color="#c65f16",
+        label="measured FA cost",
+    )
+    ax.set_yscale("log")
+    ax.set_xlabel("width")
+    ax.set_ylabel("finite-time loss gap")
+    ax.grid(alpha=0.16, which="both")
+    ax2 = ax.twinx()
+    ax2.plot(widths, [1.0 / (width * width) for width in widths], "k--", lw=1.4, label="$1/D$")
+    ax2.set_yscale("log")
+    ax2.set_ylabel("best guaranteed squared alignment")
+    handles1, labels1 = ax.get_legend_handles_labels()
+    handles2, labels2 = ax2.get_legend_handles_labels()
+    ax.legend(handles1 + handles2, labels1 + labels2, fontsize=7, loc="upper right")
+    worst_ratio = max(theory_by_width[width] / empirical_by_width[width] for width in widths)
+    ax.text(0.04, 0.05, f"null overestimate: up to {worst_ratio:.0f}×", transform=ax.transAxes, fontsize=8)
+    ax.set_title("Matrix mismatch is not optimization cost", fontsize=10)
+    panel_label(ax, "b")
+
+    ax = axes[2]
+    depths = sorted({int(row["depth"]) for row in init_rows})
+    cmap = plt.cm.viridis
+    depth_color = {depth: cmap(index / max(len(depths) - 1, 1)) for index, depth in enumerate(depths)}
+    for row in init_rows:
+        depth = int(row["depth"])
+        marker = "o" if row["rule"] == "FA" else "s"
+        ax.errorbar(
+            float(row["predicted_deficit"]),
+            float(row["empirical_deficit"]),
+            yerr=2 * float(row["empirical_stderr"]),
+            fmt=marker,
+            ms=3.8,
+            color=depth_color[depth],
+            alpha=0.68,
+            lw=0.7,
+            capsize=1,
+        )
+    for row in mnist_rows:
+        ax.plot(float(row["hidden_share"]), float(row["empirical_mean"]), "^", color="#7b6fa6", ms=4, alpha=0.65)
+    for row in cnn_rows:
+        marker = "P" if row["rule"] == "FA" else "X"
+        ax.errorbar(
+            float(row["predicted_deficit"]),
+            float(row["empirical_deficit"]),
+            yerr=2 * float(row["empirical_stderr"]),
+            fmt=marker,
+            color="#b23a48",
+            ms=5,
+            capsize=1,
+        )
+    values = [
+        float(row[key])
+        for row in init_rows
+        for key in ("predicted_deficit", "empirical_deficit")
+    ]
+    values += [
+        float(row[key])
+        for row in cnn_rows
+        for key in ("predicted_deficit", "empirical_deficit")
+    ]
+    values += [float(row[key]) for row in mnist_rows for key in ("hidden_share", "empirical_mean")]
+    lo, hi = min(values), max(values)
+    pad = 0.04 * (hi - lo + 1e-12)
+    ax.plot([lo - pad, hi + pad], [lo - pad, hi + pad], "k-", lw=1)
+    legend_handles = [
+        plt.Line2D([], [], marker="o", ls="", color="0.35", label="synthetic FA"),
+        plt.Line2D([], [], marker="s", ls="", color="0.35", label="synthetic DFA"),
+        plt.Line2D([], [], marker="^", ls="", color="#7b6fa6", label="MNIST MLP"),
+        plt.Line2D([], [], marker="P", ls="", color="#b23a48", label="CE CNN FA/DFA"),
+    ]
+    shape_legend = ax.legend(handles=legend_handles, fontsize=7, loc="lower right")
+    ax.add_artist(shape_legend)
+    depth_handles = [
+        plt.Line2D([], [], marker="o", ls="", color=depth_color[depth], label=f"depth {depth}")
+        for depth in depths
+    ]
+    ax.legend(handles=depth_handles, fontsize=6, ncols=2, loc="upper left")
+    ax.set_xlabel("exact expected initial deficit")
+    ax.set_ylabel("measured mean deficit")
+    ax.set_title("One calibration across rule, depth, and data", fontsize=10)
+    ax.grid(alpha=0.16)
+    panel_label(ax, "c")
+
+    fig.tight_layout(w_pad=2.0)
+    path = args.outdir / "figure1_mismatch_vs_initial_cost.png"
+    fig.savefig(path, bbox_inches="tight")
+    fig.savefig(path.with_suffix(".pdf"), bbox_inches="tight")
+    plt.close(fig)
+    return path
+
+
+def scatter_predictor(
+    ax: plt.Axes,
+    rows: list[dict[str, str]],
+    fields: list[tuple[str, str, str]],
+    title: str,
+) -> None:
+    depths = sorted({int(row["depth"]) for row in rows})
+    colors = {depth: plt.cm.viridis(i / max(len(depths) - 1, 1)) for i, depth in enumerate(depths)}
+    values: list[float] = []
+    for field, marker, label in fields:
+        for row in rows:
+            x = float(row[field])
+            y = float(row["empirical_gap"])
+            values.extend([x, y])
+            ax.scatter(x, y, s=13, marker=marker, color=colors[int(row["depth"])], alpha=0.60)
+        ax.scatter([], [], marker=marker, color="0.35", label=label)
+    lo, hi = min(values), max(values)
+    pad = 0.04 * (hi - lo + 1e-12)
+    ax.plot([lo - pad, hi + pad], [lo - pad, hi + pad], "k-", lw=1)
+    ax.set_xlabel("predicted gap")
+    ax.set_ylabel("measured gap")
+    ax.set_title(title, fontsize=10)
+    ax.grid(alpha=0.16)
+    ax.legend(fontsize=7)
+
+
+def figure_two(args: argparse.Namespace) -> Path:
+    rows = read_csv(args.depth_dir / "finite_time_rows.csv")
+    fig, axes = plt.subplots(1, 3, figsize=(14.0, 4.1), dpi=220)
+
+    ax = axes[0]
+    widths = sorted({int(row["width"]) for row in rows})
+    markers = ["o", "s", "^"]
+    for width, marker in zip(widths, markers):
+        depths = sorted({int(row["depth"]) for row in rows if int(row["width"]) == width})
+        means = []
+        sems = []
+        for depth in depths:
+            cell = [
+                row
+                for row in rows
+                if int(row["width"]) == width and int(row["depth"]) == depth
+            ]
+            values = np.asarray([float(row["empirical_gap"]) for row in cell])
+            init_means = np.asarray(
+                [
+                    np.mean(
+                        [
+                            float(row["empirical_gap"])
+                            for row in cell
+                            if int(row["init_seed"]) == init_seed
+                        ]
+                    )
+                    for init_seed in sorted({int(row["init_seed"]) for row in cell})
+                ]
+            )
+            means.append(float(values.mean()))
+            sems.append(float(init_means.std(ddof=1) / math.sqrt(len(init_means))))
+        ax.errorbar(depths, means, yerr=sems, marker=marker, capsize=2, label=f"width {width}")
+    ax.set_xlabel("hidden-layer depth")
+    ax.set_ylabel("FA loss - BP loss")
+    ax.set_title("Later cost changes with training dynamics", fontsize=10)
+    ax.grid(alpha=0.16)
+    ax.legend(fontsize=7)
+    panel_label(ax, "a")
+
+    scatter_predictor(axes[1], rows, [("fixed_gap", "^", "frozen $K(0)$")], "Frozen initialization misses drift")
+    panel_label(axes[1], "b")
+
+    scatter_predictor(
+        axes[2],
+        rows,
+        [
+            ("velocity_gap", "o", "linear velocity"),
+            ("retangent_gap", "s", "early retangent"),
+        ],
+        "Two snapshots recover the later cost",
+    )
+    panel_label(axes[2], "c")
+
+    predictor_legend = axes[2].get_legend()
+    if predictor_legend is not None:
+        axes[2].add_artist(predictor_legend)
+
+    depth_handles = [
+        plt.Line2D(
+            [],
+            [],
+            marker="o",
+            ls="",
+            color=plt.cm.viridis(i / max(len(sorted({int(row['depth']) for row in rows})) - 1, 1)),
+            label=f"depth {depth}",
+        )
+        for i, depth in enumerate(sorted({int(row["depth"]) for row in rows}))
+    ]
+    axes[2].legend(handles=depth_handles, fontsize=6, ncols=2, loc="lower right")
+
+    fig.tight_layout(w_pad=2.0)
+    path = args.outdir / "figure2_training_dynamics.png"
+    fig.savefig(path, bbox_inches="tight")
+    fig.savefig(path.with_suffix(".pdf"), bbox_inches="tight")
+    plt.close(fig)
+    return path
+
+
+def main() -> None:
+    args = parse_args()
+    args.outdir.mkdir(parents=True, exist_ok=True)
+    path1 = figure_one(args)
+    path2 = figure_two(args)
+    print(f"figure: {path1}")
+    print(f"figure: {path2}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/scripts/test_feedback_rules.py b/scripts/test_feedback_rules.py
new file mode 100644
index 0000000..9957315
--- /dev/null
+++ b/scripts/test_feedback_rules.py
@@ -0,0 +1,54 @@
+#!/usr/bin/env python3
+"""Numerical tests for the shared BP/FA/DFA implementation."""
+
+from __future__ import annotations
+
+import torch
+
+import feedback_rules as fr
+
+
+def main() -> None:
+    torch.set_num_threads(2)
+    dims = [5, 7, 6, 3]
+    weights = fr.initialize_mlp(dims, seed=11)
+    generator = torch.Generator().manual_seed(7)
+    x = torch.randn(9, dims[0], generator=generator, dtype=torch.float64)
+    y = torch.randn(9, dims[-1], generator=generator, dtype=torch.float64)
+
+    manual = fr.gradients(weights, x, y)
+    auto = fr.autograd_bp_gradients(weights, x, y)
+    bp_error = max(float((left - right).abs().max()) for left, right in zip(manual, auto))
+
+    fa = fr.init_feedback(dims, seed=21, rule="fa")
+    dfa = fr.init_feedback(dims, seed=22, rule="dfa")
+    fa_grads = fr.gradients(weights, x, y, rule="fa", feedback=fa)
+    dfa_grads = fr.gradients(weights, x, y, rule="dfa", feedback=dfa)
+    output_fa_error = float((manual[-1] - fa_grads[-1]).abs().max())
+    output_dfa_error = float((manual[-1] - dfa_grads[-1]).abs().max())
+
+    bp_speed = fr.squared_norm(manual)
+    output_share = fr.squared_norm([manual[-1]]) / bp_speed
+    for rule in ("fa", "dfa"):
+        ratios = []
+        for draw in range(4096):
+            feedback = fr.init_feedback(dims, seed=100_000 + draw, rule=rule)
+            rule_grads = fr.gradients(weights, x, y, rule=rule, feedback=feedback)
+            ratios.append(fr.inner_product(manual, rule_grads) / bp_speed)
+        mean_ratio = sum(ratios) / len(ratios)
+        stderr = torch.tensor(ratios, dtype=torch.float64).std(unbiased=True).item() / len(ratios) ** 0.5
+        error = abs(mean_ratio - output_share)
+        print(
+            f"{rule.upper()}: mean ratio={mean_ratio:.6f}, prediction={output_share:.6f}, "
+            f"error={error:.3g}, stderr={stderr:.3g}"
+        )
+        assert error <= 4.5 * stderr + 1e-4
+
+    print(f"BP manual/autograd max error: {bp_error:.3e}")
+    print(f"FA/DFA output-gradient errors: {output_fa_error:.3e}, {output_dfa_error:.3e}")
+    assert max(bp_error, output_fa_error, output_dfa_error) < 1e-11
+    print("feedback_rules self-test PASSED")
+
+
+if __name__ == "__main__":
+    main()