Skip to content

Return best iterate from AMICATorchNG.fit (#51) - #53

Merged
neuromechanist merged 2 commits into
mainfrom
51-multi-model-ng-log-likelihood-is-002-lower-and-more-variable-than-fortran
Jul 7, 2026
Merged

Return best iterate from AMICATorchNG.fit (#51)#53
neuromechanist merged 2 commits into
mainfrom
51-multi-model-ng-log-likelihood-is-002-lower-and-more-variable-than-fortran

Conversation

@neuromechanist

Copy link
Copy Markdown
Member

Closes #51.

Problem

Multi-model NG log-likelihood was ~0.02 lower and ~13x more variable than Fortran (n_models=2, 100 iters, real sample EEG), even though the M-step is bit-exact vs Fortran (#27).

Root cause: return-last, not a bad basin

fit() returned the last EM iterate under the deliberately non-monotone lrate schedule (both NG and Fortran anneal the rate only after an LL decrease, amica15.f90:1038-1058). An NG-only 20-seed sweep localized the variance to a late overshoot: the sole variance-driving seed climbed to LL -3.3573 (in the pack) by iter 97, then crashed to -3.5452 in its final ~2 iterations after Newton went non-positive-definite and the lrate ramp re-inflated the step. 9/20 "good" seeds also ended below their own peak.

Fix

AMICATorchNG.fit returns the highest-LL iterate:

  • keep_best: bool = True restores the best iterate when the run ends more than _KEEP_BEST_TOL (1e-9) below its peak. A monotone single-model run has best == last, so no restore fires and issue Verify NG init-basin vs Newton bit-parity via Fortran load_* init-matching #24 parity stays byte-for-byte identical (verified: max param diff 0.0 with keep_best on vs off). Inactive under do_reject (the good-sample set, hence the LL normalization, changes across iterations, making per-iteration LLs incomparable).
  • final_ll_ reports the returned iterate's LL; ll_history stays the true trajectory. The AMICA wrapper, validate_implementations.py, and both ensemble scripts read final_ll_.
  • state_dict format bumped 2 → 3 (adds keep_best, final_ll).

Results (ensemble_ll.py, N=20, real EEG + amica15mac, matched schedule)

max_iter=100 mean LL sd mean gap sd ratio
Fortran -3.3541 0.0031 -- --
NG return-last -3.3738 0.0399 -0.0197 12.7x
NG keep_best -3.3634 0.0064 -0.0093 2.0x

keep_best cuts the LL variance from 12.7x → 2.0x Fortran's sd (the headline "~13x more variable" defect).

The residual mean gap is convergence speed, not a worse optimum

NG keep_best mean LL vs iteration budget (8 seeds): -3.3639 @100, -3.3541 (== Fortran) @200, -3.3523 @300. NG needs ~2x iterations for the same LL but converges to the identical optimum (the M-step is bit-exact vs Fortran). Not a correctness issue.

Tests

  • test_keep_best_single_model_is_bit_exact (issue Verify NG init-basin vs Newton bit-parity via Fortran load_* init-matching #24 guard: identical params on vs off)
  • test_keep_best_snapshot_restore_roundtrip (snapshot clones, restore reverts)
  • test_keep_best_returns_within_tol_of_peak (contract: returned LL within tol of peak, params consistent)
  • Full torch suite: 69 passed, 5 skipped. Ruff clean; ty no new diagnostics.

Real sample data + Fortran binary only (NO MOCK). ADR .context/decisions/0003-best-iterate-safeguard.md; analysis + reproduction in .context/issue-51/.

Multi-model NG log-likelihood was ~0.02 lower and ~13x more variable than
Fortran because fit() returned the LAST EM iterate under the deliberately
non-monotone lrate schedule (both NG and Fortran anneal only after an LL
decrease). The variance was driven by late Newton-fallback overshoots: on
the sample EEG one seed peaked at LL -3.357 then crashed to -3.545 in its
final iterations. It was return-last, not a bad basin.

fit() now tracks and returns the highest-LL iterate:
- keep_best (default True) restores the best iterate when the run ends more
  than _KEEP_BEST_TOL (1e-9) below its peak; a monotone single-model run has
  best == last, so no restore fires and issue #24 parity stays byte-for-byte
  identical (verified: max param diff 0.0). Inactive under do_reject (the
  good-sample set, hence the LL normalization, changes across iterations).
- final_ll_ reports the returned iterate's LL; ll_history stays the true
  trajectory. The AMICA wrapper, validate_implementations, and the ensemble
  scripts read final_ll_.
- state_dict format bumped 2 -> 3 (adds keep_best, final_ll).

At matched 100-iter budget this cuts the LL sd from 12.7x to 2.0x Fortran's.
The residual ~0.009 mean gap is convergence speed, not a worse optimum: NG
reaches Fortran's exact mean (-3.3541) at 200 iters and exceeds it at 300.

Tested: 69 torch tests pass (+3 keep_best tests); ruff clean; ty no new
diagnostics. ADR 0003, .context/issue-51/.
PR review (4 Sonnet reviewers) findings:
- final_ll_ was left at the last finite ll_history value after a degenerate
  stop while the model held diverged params; now set to NaN there
  (silent-failure review).
- Roll back n_kurt_done with pdtype in the best-iterate snapshot so a restored
  model's adaptive-PDF switch count stays consistent with its pdtype.
- Warn once when keep_best is requested together with do_reject (the safeguard
  is silently inactive there).
- Fix the restore-skip comment (a singular_ll stop leaves A/W finite, not
  non-finite; the skip is a scope decision deferred to #50) and qualify the
  _KEEP_BEST_TOL min_dll scale reference (Fortran's, not the legacy NumPy one).

Tests added/expanded:
- keep_best inactive under do_reject (final_ll_ == ll_history[-1]).
- keep_best False-vs-True differential on an overshooting multi-model run, with
  an explicit skip when the run is monotone so restore-branch coverage is
  visible, not vacuous.
- wrapper final_ll_ survives save/load.

69 -> 72 torch tests pass; ruff clean.
@neuromechanist

Copy link
Copy Markdown
Member Author

Review response (4 Sonnet reviewers: code, silent-failure, tests, comments)

code-reviewer: no ≥80-confidence issues (ran the new tests + the #24 parity guard, confirmed invariants, 0 new ty diagnostics). Its two sub-threshold notes are addressed below.

Fixed (commit 4e063cf)

  • final_ll_ stale on a degenerate stop (silent-failure, critical — backend part). It was left at the last finite ll_history value while the model held diverged parameters. Now set to NaN when stop_reason is degenerate.
  • n_kurt_done not rolled back with pdtype (silent-failure medium + code-reviewer nit). The best-iterate snapshot now captures n_kurt_done, so a restored/saved adaptive-PDF model's switch count stays consistent with its pdtype.
  • do_reject silently disables the default-on safeguard (silent-failure medium). fit now emits a one-time logger.warning when keep_best and do_reject.
  • Restore-skip comment was wrong (comment-analyzer, critical). A singular_ll stop leaves A/W finite-but-singular, not non-finite; the comment now states the real reason (salvaging would pre-empt AMICA wrapper marks is_fitted_ even on a degenerate (nan_ll/singular_ll) fit #50's degenerate-fit contract; state_dict already refuses on stop_reason).
  • _KEEP_BEST_TOL scale reference ambiguous (comment-analyzer). Qualified to Fortran's min_dll (normalized), explicitly not the legacy NumPy min_dll (un-normalized).
  • Tests (test-analyzer gaps): added test_keep_best_inactive_under_reject (final_ll_ == ll_history[-1] under rejection); expanded the peak test with the keep_best=False differential and an explicit pytest.skip when the run is monotone so restore-branch coverage is visible rather than vacuous; the wrapper save/load test now checks final_ll_. 72 torch tests pass.

Deferred to #50 (with rationale)

Declined (with rationale)

  • INFO → WARNING for the restore log (silent-failure, high). A restore is the routine, correct action of a default-on safeguard and fires on most multi-model fits (never on the monotone single-model default). Promoting it to WARNING would be alarm-fatigue noise on normal runs, and the outcome is fully visible via final_ll_/ll_history. The sibling nan_ll/lrate_floor warnings are genuine abnormal stops — a different category. Kept at INFO.
  • Dedicated degenerate-stop test not added. The backend is robust on the clean sample EEG (even lrate=50 does not diverge to nan_ll/singular_ll), so a natural degenerate stop is not deterministically inducible without rank-deficient/synthetic input, which conflicts with the NO-MOCK policy. The degenerate branches (final_ll_=NaN, restore-skip) are simple guards; the degenerate stop_reason path is exercised by test_state_dict_refuses_degenerate_model.

@neuromechanist
neuromechanist merged commit 2b49ac6 into main Jul 7, 2026
5 checks passed
@neuromechanist
neuromechanist deleted the 51-multi-model-ng-log-likelihood-is-002-lower-and-more-variable-than-fortran branch July 7, 2026 05:17
neuromechanist added a commit that referenced this pull request Aug 10, 2026
Added tests that the five new config keys (use_min_dll/min_dll/
maxincs/use_grad_norm/min_nd) round-trip through state_dict()/
from_state_dict(), and that a simulated pre-#207 payload (format_
version 3, missing those keys) still loads and falls back to the
Fortran defaults. Also documented, at the format_version check itself,
why it deliberately was not bumped for this change (prior precedent
#52/#53 bumped it; the additive-only new keys don't need to).

PR #213 review finding 6.
neuromechanist added a commit that referenced this pull request Aug 10, 2026
Added tests that the five new config keys (use_min_dll/min_dll/
maxincs/use_grad_norm/min_nd) round-trip through state_dict()/
from_state_dict(), and that a simulated pre-#207 payload (format_
version 3, missing those keys) still loads and falls back to the
Fortran defaults. Also documented, at the format_version check itself,
why it deliberately was not bumped for this change (prior precedent
#52/#53 bumped it; the additive-only new keys don't need to).

PR #213 review finding 6.
neuromechanist added a commit that referenced this pull request Aug 10, 2026
* feat: add NG convergence stops (issue #207)

AMICATorchNG was missing three Fortran convergence criteria
(amica15.f90): use_min_dll/maxincs (consecutive small-gain stop),
use_grad_norm/min_nd (weight-gradient-norm stop), and the
decrease-branch's ".or. ndtmpsum <= min_nd" half -- the last is
why lrate could sit at newtrate and oscillate under do_newton=True
without ever stopping. All three are Fortran-faithful defaults
(True/1e-9/5/True/1e-7, matching amica15_header.f90).

ndtmpsum is now computed every iteration (Fortran-faithful,
including the comp_used mask): the direction/dAk computation in
_update_parameters was decoupled from the share_comps A-freeze
gate, since Fortran computes dAk/ndtmpsum unconditionally in
accum_updates_and_likelihood, strictly before the separately-gated
update_A step. Default (non-sharing) path is unaffected.

New stop_reason values (min_dll, grad_norm, grad_norm_floor) are
converged, not degenerate; plumbed through AMICA (**kwargs) and
state_dict()/from_state_dict(). Corrected a misleading comment
claiming only amica17 normalizes LL before the min_dll comparison
-- amica15 (the actual reference binary's source) normalizes
identically; the real divergence is numpy_impl's un-normalized
raw-sum LL comparison, a separate pre-existing gap this does not
touch.

Tested: full torch suite green (188 passed, 5 pre-existing skips);
validate_implementations.py unchanged at max-iter 100 and 2000
(matches sample_params.json's budget) -- PyTorch LL/iteration count
identical before/after, neither stop fires on the bundled 32-channel
sample within that budget.

* test: add NG convergence-stop suite (issue #207)

Real bundled sample EEG only. Covers: each stop firing with the
right stop_reason (min_dll, grad_norm, grad_norm_floor, and the
pre-existing lrate_floor unshadowed by the new check); the maxincs
consecutive-count rule including reset-on-larger-gain, verified
against an independent reimplementation applied to a stops-disabled
reference trajectory rather than hardcoded iteration numbers; the
have_prev guard (never fires before two LL values exist); the
share_comps freeze window still computing a fresh (non-stale)
ndtmpsum every iteration; a converged stop leaving transform/
state_dict/AMICA.save usable; keep_best and do_reject interactions;
and that both stops disabled reproduces pre-#207 behavior (never
emits the three new stop_reason values).

14/14 pass; ruff and ty clean.

* test: drop slow marker from non-Fortran convergence tests

Three tests in test_ng_convergence.py used only bundled EEG and pure
PyTorch (no Fortran binary) and ran in seconds, but were marked
@pytest.mark.slow. CI runs pytest -m "not slow", documented as
excluding tests that invoke the macOS-only Fortran reference binary,
so these three never ran in CI. One of them,
test_a_frozen_window_still_computes_fresh_grad_norm, is the only test
of the dAk/A-freeze decoupling (issue #207), so that change had zero
CI coverage.

PR #213 review finding 1.

* fix: document and test stop_reason shadowing (issue #207)

None of the three fit()-loop stop blocks (decrease branch; min_dll;
grad_norm) short-circuits on an earlier one having already fired the
same iteration, matching Fortran's independent leave=.true. structure
(not a fidelity bug). But the standalone grad_norm check runs
unconditionally after the decrease branch, so under the shipped
use_grad_norm=True default it always wins: "grad_norm_floor" is
unreachable as a final stop_reason, and the use_grad_norm docstring
wrongly implied it was the fix for the reported CUDA case. Corrected
the docstring/comments in torch_impl/core.py and amica.py, and added a
test proving the shadowing under shipped True/True defaults.

PR #213 review finding 2.

* test: rename mislabeled do_reject test, add missing coverage

test_do_reject_interaction_min_dll_stop_leaves_good_idx_usable set
use_min_dll=False, so it actually exercised grad_norm_floor, not
min_dll; renamed to match. Added the two genuinely missing do_reject
combinations: the standalone min_dll and grad_norm stops.

PR #213 review finding 3.

* test: exercise a genuine keep_best overshoot restore

The old test's trajectory was monotonically increasing, so
final_ll_ == max(ll_history) == ll_history[-1] held whether the
restore logic worked or was a no-op. Reworked it around the known
non-monotone recipe from test_write_amica_output_ll_matches_kept_iterate
(#92), combined with a loosened min_dll so the run stops a few
iterations past its peak via the new min_dll stop_reason, and assert
final_ll_ != ll_history[-1] to prove the restore branch actually ran.

PR #213 review finding 4.

* test: exercise AMICA.save/load in convergence-stop tests

The wrapper usability test's own docstring claimed "transform()/
save() usable" but never called AMICA.save() anywhere in the file.
Added a real save()/load() round trip and confirmed the reloaded
model reports the same stop_reason_ and reproduces transform()
exactly.

PR #213 review finding 5.

* test: cover issue #207 config persistence round-trip

Added tests that the five new config keys (use_min_dll/min_dll/
maxincs/use_grad_norm/min_nd) round-trip through state_dict()/
from_state_dict(), and that a simulated pre-#207 payload (format_
version 3, missing those keys) still loads and falls back to the
Fortran defaults. Also documented, at the format_version check itself,
why it deliberately was not bumped for this change (prior precedent
#52/#53 bumped it; the additive-only new keys don't need to).

PR #213 review finding 6.

* test: add stop reachability at literal shipped default thresholds

Every prior min_dll/grad_norm test loosens the threshold by 5-6 orders
of magnitude to force a fast stop, which would not catch a scale bug
in the comparison itself (issue #212 found exactly that in the
numpy_impl backend). Added a fast (a few seconds), non-slow test at
the literal shipped defaults (min_dll=1e-9, maxincs=5,
use_grad_norm=True, min_nd=1e-7, none overridden) that reaches
min_dll via an early Newton start on a small real-data subset.

PR #213 review finding 7.

* test: cover mir_history_ vs keep_best restore and save/load

Issue #161 flagged two documented-but-untested mir_history_ claims:
that a keep_best restore (#51) does not rewrite it (so its last entry
can be from discarded, pre-restore parameters, distinct from
model.mir(X) on the returned ones), and that it comes back empty after
a save/load round trip (not persisted in state_dict()). Both verified
true on real data; also updated the module docstring to summarize the
full set of PR #213 review additions in this file.

Folds in issue #161.

* docs: add changelog entry for issue #207 convergence stops

* fix: correct timing claim in threshold-reachability docstring

The test is actually the slowest in the file (~8s from the 326-iteration
fit itself), not "well under 3s" as originally written; corrected to
state the real number and why it's still not a slow-marker candidate.

* Give the shipped-default reachability test real headroom

The iteration at which min_dll fires is BLAS-dependent: 326 on
macOS-arm64, 412 on Linux-x86_64 with a CUDA torch build, and past 500
on the GitHub Linux runner, where CI failed with stop_reason=max_iter.
At max_iter=500 the test consumed 82 percent of its budget on the
fastest platform, so any numerical variation tipped it over.

The claim under test is that the default threshold is reachable at all,
not that it is reached by a given iteration, so raise the budget well
above the observed spread and record the spread in the docstring.

Refs #207

* Assert convergence behaviour, not the iteration it happens on

CI failed twice on iteration-count assumptions. First the stop reason,
because max_iter=500 left no headroom; then a leftover
len(ll_history) < 500 bound after the budget was raised. The stop fires
at 326 on macOS-arm64, 412 on Linux-x86_64 with a CUDA torch build, and
1076 on the GitHub runner, so any constant fitted to one machine is a
trap.

Both bounds now track the budget. The unrelated len(ll_history) == 23
golden value in the lrate_floor test is loosened for the same reason,
before it fails the same way.

Refs #207

* Rebuild paper.pdf [skip ci]

---------

Co-authored-by: github-actions[bot] <github-actions[bot]@users.noreply.github.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Multi-model NG log-likelihood is ~0.02 lower and more variable than Fortran

1 participant