You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
* Docs: comprehensive validation evidence (#108)
Restructure guides/validation.md into a full evidence page: at-a-glance
summary table, source-density bit-exactness table, cross-platform
device/precision invariance (phase3 equivalence matrix + IC topomaps),
EEGLAB drop-in round-trip, and an other-validated-behaviors table.
Collapse the 40-run per-run table for scannability. Amari-distance
figures cross-checked against amari_summary.json. Builds under
mkdocs build --strict.
* Docs: address PR review findings on validation page
Scope source-density bit-exactness to the non-GG families (the default
GG has no literal Fortran oracle; validated by single-model parity),
fix two test-path citations to include the pyAMICA/ prefix, reconcile
the two k=30/70ch rows (different backend sets), define IC/LL before
first use, expand ZCA, and drop the unused EM abbreviation.
**Correctness in pyAMICA is defined as parity with the reference Fortran binary,
4
4
not merely as convergence.** A run is correct when it reproduces the Fortran output within numerical tolerance.
5
+
This page collects the full verification evidence: bit-exact score functions, single-model parity,
6
+
multi-model distributional equivalence, cross-platform device and precision invariance,
7
+
the EEGLAB drop-in round-trip, and the remaining validated behaviors.
8
+
Every result uses the bundled real sample EEG and the reference Fortran binary; none uses synthetic data.
9
+
Throughout, IC abbreviates independent component and LL log-likelihood.
10
+
11
+
## Validation at a glance
12
+
13
+
| What was checked | How | Result |
14
+
|---|---|---|
15
+
| Source-density score and log-density (non-GG families) | vs the literal `amica15.f90` expressions | bit-exact ($<10^{-12}$) |
16
+
| Per-block sufficient statistics and one M-step | vs Fortran | bit-exact ($\sim\!10^{-15}$) |
17
+
| Single-model solution | log-likelihood, component correlation, Amari distance vs Fortran | LL within ~0.005 of $-3.4018$; correlation 0.997; Amari 0.006 |
18
+
| Multi-model solution | distributional equivalence over 20-run ensembles | indistinguishable from Fortran's own run-to-run spread ($p = 0.96$) |
19
+
| Device and precision invariance | same independent components across CPU/CUDA/MPS/MLX, float32/float64, Linux/macOS | identical (1.000) across all eight torch/MLX combinations |
20
+
| Cross-backend log-likelihood | converged LL across every backend | agree to ~3 significant digits (max pairwise ~0.003) |
21
+
| EEGLAB output |`write_amica_output` round-trip through `loadmodout15`| single-model bytes are an exact serialization; loads with correct layout |
22
+
| Degenerate fits | NaN or singular log-likelihood | refused, never returned as NaN sources |
5
23
6
24
## The validation harness
7
25
@@ -23,6 +41,31 @@ On real sample EEG the natural-gradient backend reaches Fortran's solution:
23
41
The fixed source-density families are bit-exact against the literal Fortran score/derivative expressions (~1e-15),
24
42
and the backend converges to the binary's solution within ~0.005 log-likelihood.
25
43
44
+
### Source-density families are bit-exact
45
+
46
+
AMICA models each source with one of the reference's five `pdftype` density families.
47
+
For every family other than the default generalized Gaussian, the vectorized log-density and score reproduce the literal `amica15.f90` expressions
48
+
to float64 precision (test bound $<10^{-12}$, observed $\sim\!10^{-15}$):
49
+
the source model is not an approximation of the Fortran one, it is the same function.
50
+
The generalized Gaussian has no closed-form literal to compare against, since its score depends on the adaptive shape $\rho$,
51
+
so the default family is validated by the single-model parity above instead.
52
+
The oracle column below records which check applies to each family:
53
+
54
+
|`pdftype`| Density family | Score $f_p(y)$ | Bit-exact oracle |
`pdftype=0` is the default and is byte-for-byte the pre-family implementation.
63
+
The `pdftype=1` extended-Infomax switcher flips each source between the super-Gaussian (code 1) and
64
+
sub-Gaussian (code 4) densities on a kurtosis schedule; its dynamic switch has no bit-exact oracle
65
+
(the reference's `do_choose_pdfs` accumulator is dead code in the binary), so it is validated by real-data
66
+
log-likelihood instead. Each fixed family converges within ~0.005 LL of the binary at a matched Newton budget.
67
+
See `pyAMICA/tests/torch_tests/test_ng_pdf_families.py` and ADR 0002.
68
+
26
69
## Multi-model equivalence
27
70
28
71
Multi-model AMICA is not partition-identifiable, so exact partition parity with Fortran is the wrong acceptance bar.
@@ -71,54 +114,81 @@ This pairing correction is not free: on this ensemble it lowers the reported dis
71
114
72
115
The same run-level permutation test (20000 permutations, intact 40-run units) finds no evidence that between-implementation agreement is worse than Fortran's own run-to-run agreement ($p > 0.999$), agreeing with the correlation-based conclusion above.
73
116
74
-
#### Per-run detail
117
+
??? note "Per-run detail (all 40 runs, both metrics)"
118
+
119
+
Table 1 in the paper and the group summaries above report distribution means;
120
+
the table below gives each of the 40 runs' own mean agreement to its own group's other 19 runs (`within`) and to all 20 opposite-implementation runs (`between`), for both metrics.
121
+
Regenerate with `uv run python .context/issue-27/amari_distance.py`, which writes `.context/issue-27/per_run_detail.csv`.
122
+
123
+
| implementation | run | corr within | corr between | Amari within | Amari between |
The strongest reassurance that pyAMICA is a single, well-defined implementation is that it recovers the *same*
169
+
independent components no matter where or how it runs. Fitting the same real EEG (ds002718 sub-002, 147,000 frames, 70 channels, 2000 iterations)
170
+
on every backend and Hungarian-matching the unmixing components across them:
171
+
172
+
{ width=680 }
173
+
/// caption
174
+
Mean Hungarian-matched \|correlation\| of the recovered components between every pair of backends.
175
+
The eight torch/MLX combinations, CPU, CUDA, MPS, and MLX, at both float32 and float64, on macOS-arm64 and Linux-x86_64, all agree at **1.000**.
176
+
///
75
177
76
-
Table 1 in the paper and the group summaries above report distribution means;
77
-
the table below gives each of the 40 runs' own mean agreement to its own group's other 19 runs (`within`) and to all 20 opposite-implementation runs (`between`), for both metrics.
78
-
Regenerate with `uv run python .context/issue-27/amari_distance.py`, which writes `.context/issue-27/per_run_detail.csv`.
178
+
**Every torch/MLX backend is identical to every other at 1.000**: the same decomposition on any device, at any
179
+
precision, on either operating system. This is the definitive "float32 == float64" and "GPU == CPU" result.
180
+
The two native-Fortran builds (macOS-arm64 and Linux-x86_64) agree with each other at 0.972 and with the
181
+
torch/MLX cluster at ~0.90. That residual is not a backend defect: Fortran does not even reach 1.000 against
182
+
*itself* across platforms, because each native run is seeded from the clock and settles into a different (equally valid)
183
+
local optimum on the weakly-determined components only. As the decomposition becomes better-determined that gap closes
184
+
(next section), and the component maps are visibly the same down every row:
79
185
80
-
| implementation | run | corr within | corr between | Amari within | Amari between |
As `k` grows, cross-backend component equivalence rises toward 1.0;
130
200
at the rule-of-thumb minimum (`k` around 20-30) only the strongest components are backend-reproducible,
131
201
while the rest are under-determined and settle into different but equally valid local optima (AMICA is non-convex).
202
+
This is why the native-Fortran rows above sit at ~0.90 (70 channels, `k` = 30) rather than 1.0. Two independent sweeps confirm it.
203
+
204
+
### Sweeping channels at fixed frames
205
+
206
+
Holding frames at 147,000 and increasing the channel count lowers `k`.
207
+
Comparing MLX-float32 against the independent native-Fortran-float64 build (the hardest cross-implementation pair):
208
+
209
+
| channels | k = frames/ch² | mean matched \|corr\|| components > 0.95 |
210
+
|---:|---:|---:|---:|
211
+
| 16 | 574 |**0.997**| 16/16 |
212
+
| 32 | 144 | 0.974 | 27/32 |
213
+
| 48 | 64 | 0.954 | 34/48 |
214
+
| 70 | 30 | 0.898 | 20/70 |
215
+
216
+
At high `k`**every backend, including the independent Fortran build, recovers identical components** (0.997, all 16/16 at `k` = 574).
217
+
At `k` = 30, the rule-of-thumb minimum, only the strongest ~20/70 components are reproducible; the rest are under-determined.
132
218
133
-
### Data-size sweep: equivalence versus k
219
+
### Sweeping frames at fixed channels
134
220
135
-
Holding channels fixed at 70 and increasing the number of frames (so `k` rises),
136
-
on real EEG (ds002718 sub-002),
137
-
cross-backend IC equivalence climbs sharply and then saturates once the decomposition is well-determined:
221
+
Holding channels at 70 and increasing frames raises `k`, on the same real EEG:
138
222
139
223
{ width=640 }
140
224
/// caption
141
-
Mean Hungarian-matched cross-backend |correlation| versus $k = \text{frames} / \text{channels}^2$ (70 channels, 2000 iterations,
225
+
Mean Hungarian-matched cross-backend \|correlation\| versus $k = \text{frames} / \text{channels}^2$ (70 channels, 2000 iterations,
142
226
native-Fortran and PyTorch-CUDA float64/float32 backends).
143
227
Equivalence saturates at ~0.98 once $k \geq 60$.
144
228
///
@@ -151,6 +235,9 @@ Equivalence saturates at ~0.98 once $k \geq 60$.
151
235
| 490,000 | 100 | 0.983 | 94.8% |
152
236
| 747,750 | 152 | 0.982 | 92.4% |
153
237
238
+
The `k` = 30 row here (0.929) sits above the 0.898 at `k` = 30 in the channel-sweep table because the two average different backend sets:
239
+
the channel sweep reports only the hardest MLX-versus-native-Fortran pair, while this frame sweep averages over the native-Fortran and PyTorch-CUDA float64/float32 cluster.
240
+
154
241
!!! note "The threshold is data-specific"
155
242
For this recording the equivalence knee falls **between k=30 and k=60**; below
156
243
it the backends settle into different (equally valid) local optima, above it
@@ -172,6 +259,27 @@ The two double-precision implementations, an independent native Fortran binary a
172
259
173
260
This is cross-*implementation* agreement, not just cross-device. The residual gap is dominated by the float32 path (rounding accumulated over 2000 iterations, plus an early stop when the natural-gradient learning rate hit its floor), which is a convergence/precision effect rather than a backend defect.
174
261
262
+
## EEGLAB drop-in round-trip
263
+
264
+
pyAMICA writes the same on-disk format EEGLAB's AMICA plugin reads, so a fit is a drop-in replacement:
265
+
no re-sorting, sign-flipping, or reformatting. After a fit, `write_amica_output(dir)` writes the raw binary
Beyond the core parity results, the following AMICA features are implemented and validated. Where the reference
361
+
binary contains a bit-exact oracle it is used; where the reference code path is unrunnable (declared but never
362
+
allocated, so it cannot be exercised even in Fortran) the feature is behavior-validated on real EEG instead, and
363
+
guarded to a no-op so the parity results above stay byte-for-byte unchanged.
364
+
365
+
| Behavior | Status | Validation |
366
+
|---|---|---|
367
+
| Best-iterate safeguard (`keep_best`, #51) | on by default | returns the highest-LL iterate; cuts multi-model LL sd from 12.7x to 2.0x Fortran's. Single-model parity stays bit-exact (monotone, no restore). ADR 0003 |
368
+
| Per-model bias `c` update (#27) | on for `n_models>1`| Fortran `update_c`; per-block stats bit-exact; no-op for `n_models=1`|
369
+
| Component sharing (`share_comps`, #60) | off by default | Fortran `identify_shared_comps` ported; no bit-exact oracle (`Spinv2` is unrunnable), behavior-validated; byte-identical when unshared |
370
+
| Outlier rejection (`do_reject`, #123) | off by default |`good_idx` mechanism in both backends; NumPy port validated vs the PyTorch backend |
371
+
| Degenerate-fit contract (#50) | always | a NaN or singular fit is refused (`converged_` / `stop_reason_`); `transform`/`get_*`/`save`/`write_amica_output` raise instead of returning NaN sources |
372
+
373
+
Tests live under `pyAMICA/tests/`: `torch_tests/test_ng_backend.py`, `torch_tests/test_ng_sharing.py`, `torch_tests/test_amica_ng_wrapper.py`, and `test_numpy_reject.py`.
374
+
375
+
## Reproducing these results
376
+
377
+
Everything on this page runs from the bundled sample data and the reference binary, with no external download:
378
+
379
+
```bash
380
+
uv run python validate_implementations.py # single- and multi-model parity report
381
+
uv run pytest # the full parity/behavior test suite
382
+
```
383
+
384
+
The multi-model ensemble and Amari detail regenerate from saved fits (no re-fitting) with
385
+
`uv run python .context/issue-27/amari_distance.py`. The cross-platform benchmark and equivalence
386
+
figures are produced by `benchmarks/benchmark_decompose.py` (and the sweep scripts alongside it);
387
+
the underlying findings are in `.context/issue-84/` and `.context/issue-90/`.
0 commit comments