diff --git a/pamica/numpy_impl/core.py b/pamica/numpy_impl/core.py index 23dcb8b..6dbe3c6 100644 --- a/pamica/numpy_impl/core.py +++ b/pamica/numpy_impl/core.py @@ -628,6 +628,9 @@ def _reinitialize_for_restart(self): self.good_idx = saved_good_idx self.num_good_samples = saved_num_good self.lrate = self.lrate0 + # rholrate is the (maxdecs-ratcheted) rho-rate ceiling; restore it to the + # pristine rholrate0 so a re-fit starts fresh (issue #193). + self.rholrate = self.rholrate0 self.ll = [] self.nd = [] @@ -1468,10 +1471,13 @@ def _check_convergence( Mirrors Fortran's per-iteration convergence handling (amica17.f90:1062-1103). On a likelihood decrease Fortran does NOT stop - at ``maxdecs``; it lowers the learning-rate *ceiling* (``lrate0``, and - ``newtrate``/``rholrate0`` under Newton) and continues, which is what - keeps a long run from oscillating and drifting past its converged - solution (issue #41). The updated ``numdecs``/``numincs`` counters are + at ``maxdecs``; it lowers the learning-rate *ceilings* (``lrate0``; the + rho rate once ``iter > newt_start``; ``newtrate`` under Newton) and + continues, which is what keeps a long run from oscillating and drifting + past its converged solution (issue #41). The rho ceiling is + ``self.rholrate`` here (reset to ``rholrate0`` each fit), ratcheted only + at ``maxdecs`` -- never per-decrease (issue #193). The updated + ``numdecs``/``numincs`` counters are returned so they accumulate across iterations (they previously did not). Parameters @@ -1517,8 +1523,8 @@ def _check_convergence( grad_norm = self.nd[-1] if len(self.nd) > 0 else None # Likelihood decrease (Fortran amica17.f90:1062-1083): reduce the current - # lrate, and once maxdecs decreases have accrued, ratchet the ceiling - # (lrate0, plus newtrate/rholrate0 under Newton) down and continue -- + # lrate, and once maxdecs decreases have accrued, ratchet the ceilings + # (lrate0; the rho rate; newtrate under Newton) down and continue -- # NOT stop. Only a lrate/gradient floor terminates on a decrease. if self.ll[-1] < self.ll[-2]: if self.lrate <= self.minlrate or ( @@ -1531,12 +1537,17 @@ def _check_convergence( numincs, ) self.lrate *= self.lratefact - self.rholrate *= self.rholratefact numdecs += 1 if numdecs >= self.max_decs: self.lrate0 *= self.lratefact if self.iter > self.newt_start: - self.rholrate0 *= self.rholratefact + # rho rate is a ceiling reset to rholrate0 each iteration + # (Fortran amica15.f90:1788); it ratchets ONLY here at maxdecs + # (amica15.f90:1050), never per LL-decrease. The old per-decrease + # self.rholrate *= rholratefact was a monotone decay with no + # reset that collapsed the rho rate and froze rho at a stale + # shape (issue #193). + self.rholrate *= self.rholratefact if self.do_newton and self.iter > self.newt_start: self.newtrate *= self.lratefact numdecs = 0 diff --git a/pamica/tests/test_sample_data.py b/pamica/tests/test_sample_data.py index 9b120b8..b4f971a 100644 --- a/pamica/tests/test_sample_data.py +++ b/pamica/tests/test_sample_data.py @@ -665,9 +665,14 @@ def _get_updates_and_likelihood(self): def test_check_convergence_ratchets_lrate_on_decrease(): """After max_decs consecutive LL decreases, _check_convergence lowers the - learning-rate ceiling (lrate0, and newtrate under Newton) and resets numdecs - -- it does NOT stop -- matching Fortran's ratchet (#41). Drives the - convergence handler directly with a decreasing LL history. + learning-rate ceilings (lrate0; the rho rate once iter>newt_start; newtrate + under Newton) and resets numdecs -- it does NOT stop -- matching Fortran's + ratchet (#41). Drives the convergence handler directly with a decreasing LL + history. + + Also pins the #193 fix: the rho rate is a maxdecs-ratcheted CEILING, NOT a + per-decrease monotone decay. rholrate must stay untouched on a decrease below + max_decs and only ratchet when numdecs hits max_decs. """ model = AMICA( num_models=1, @@ -679,16 +684,19 @@ def test_check_convergence_ratchets_lrate_on_decrease(): use_grad_norm=False, use_min_dll=False, ) - model.iter = 100 # past newt_start, so the Newton-rate ratchet applies + model.iter = 100 # past newt_start, so the Newton/rho-rate ratchets apply model.nd = [1.0] # gradient norm well above min_grad_norm (no floor stop) lrate0_before = model.lrate0 newtrate_before = model.newtrate + rholrate_before = model.rholrate - # First decrease: numdecs -> 1 (below max_decs), so just reduce lrate. + # First decrease: numdecs -> 1 (below max_decs), so just reduce lrate. The rho + # ceiling must NOT move here -- pre-#193 it decayed on every decrease. model.ll = [-3.0, -3.1] conv, _, numdecs, numincs = model._check_convergence(0, 0) assert conv is False assert numdecs == 1 + assert model.rholrate == rholrate_before # Second consecutive decrease: numdecs hits max_decs=2 -> ratchet ceilings, # reset numdecs, and CONTINUE (not converged). @@ -698,6 +706,34 @@ def test_check_convergence_ratchets_lrate_on_decrease(): assert numdecs == 0 assert model.lrate0 == lrate0_before * 0.5 assert model.newtrate == newtrate_before * 0.5 + assert model.rholrate == pytest.approx(rholrate_before * model.rholratefact) + + +def test_reinitialize_for_restart_resets_rho_ceiling(): + """Issue #193: a mid-fit restart (non-finite LL, Fortran restartiter path) + resets the rho-rate ceiling to rholrate0, exactly like the lrate reset -- a + previously-ratcheted rholrate must not carry across the restart. + """ + model = AMICA(num_models=1, num_mix=3, seed=0) + model.data_dim = 4 + model.num_samples = 16 + model.num_comps = 4 + model.good_idx = None + model.num_good_samples = 16 + model._initialize_parameters() + + # Simulate a fit that ratcheted both ceilings, then hit a non-finite LL. + model.rholrate = model.rholrate0 * 0.25 + model.lrate = model.lrate0 * 0.25 + model.ll = [-3.0, -3.1] + model.nd = [0.5] + + model._reinitialize_for_restart() + + # Both ceilings restored to pristine; history cleared for a fresh judgment. + assert model.rholrate == model.rholrate0 + assert model.lrate == model.lrate0 + assert model.ll == [] @pytest.mark.skipif(not op.exists(eeglab_data_file), reason="sample data missing") diff --git a/pamica/tests/torch_tests/test_ng_backend.py b/pamica/tests/torch_tests/test_ng_backend.py index 590a35c..7a842f0 100644 --- a/pamica/tests/torch_tests/test_ng_backend.py +++ b/pamica/tests/torch_tests/test_ng_backend.py @@ -1380,3 +1380,43 @@ def test_keep_best_inactive_under_reject(): m.fit(data, max_iter=12, verbose=False) assert m.numrej >= 1 # rejection actually fired, so the good set changed assert m.final_ll_ == m.ll_history[-1] # no best-iterate restore under reject + + +@pytest.mark.slow +@pytest.mark.skipif(not DATA_FILE.exists(), reason="sample data missing") +def test_rholrate_ratchets_at_maxdecs_not_per_decrease(): + """Issue #193: the rho learning rate is a maxdecs-ratcheted *ceiling*, not a + per-LL-decrease monotone decay. + + Fortran resets ``rholrate = rholrate0`` every iteration before the rho update + (amica15.f90:1788) and only tightens the ceiling at ``maxdecs`` + (amica15.f90:1050, gated on ``iter > newt_start``). torch previously decayed + ``rholrate`` on EVERY LL decrease with no reset, collapsing it to ~1e-5 within + a few hundred iterations and freezing rho at a stale shape. + + On the real sample data a long-enough Newton run overshoots and triggers + several LL decreases. The surviving ``rholrate`` must have ratcheted exactly + as often as ``newtrate`` (both gated on ``iter > newt_start`` at ``maxdecs``, + matched 0.5 factor here), NOT once per decrease. + """ + data = _load_real_data() + m = _fresh_ng( + block_size=512, do_newton=True, newt_start=50, newtrate=1.0, lrate=0.05, + lratefact=0.5, rholrate=0.05, rholratefact=0.5, maxdecs=3, + ) # fmt: skip + m.fit(data, max_iter=300, verbose=False) + + ll = m.ll_history + n_dec = sum(1 for i in range(1, len(ll)) if ll[i] < ll[i - 1]) + # The run must actually exercise the schedule (Newton overshoots past + # newt_start), or the assertions below prove nothing. + assert m.newtrate < m.newtrate0, "no newtrate ratchet fired; budget too short" + assert n_dec > m.maxdecs, "too few LL decreases to distinguish the schedules" + + # rholrate and newtrate share the maxdecs ratchet schedule, so the surviving + # rho ceiling ratcheted the same number of times as newtrate. + assert m.rholrate == pytest.approx(m.rholrate0 * (m.newtrate / m.newtrate0)) + # The old per-decrease decay (rholrate0 * rholratefact**n_dec) sits orders of + # magnitude below the fixed ceiling; guard against a regression to it. + buggy = m.rholrate0 * (m.rholratefact**n_dec) + assert m.rholrate > buggy * 10 diff --git a/pamica/torch_impl/core.py b/pamica/torch_impl/core.py index 8f4f89d..e51d87e 100644 --- a/pamica/torch_impl/core.py +++ b/pamica/torch_impl/core.py @@ -1824,12 +1824,24 @@ def fit( # Learning-rate control, ported from Fortran (amica17.f90:1062-1108). # Natural-gradient/Newton ascent is not monotonic at a fixed rate: - # when the log-likelihood decreases, anneal the working rates - # (lrate, rholrate). If decreases persist for maxdecs iterations, - # ratchet the *ceilings* down (lrate_cap, and newtrate once Newton - # is running) so the per-iteration ramp can no longer re-inflate - # lrate back to the overshooting value -- without this the ramp and - # a one-shot halving just oscillate and the LL drifts down. + # when the log-likelihood decreases, anneal the working lrate. If + # decreases persist for maxdecs iterations, ratchet the *ceilings* + # down (lrate_cap; newtrate once Newton is running; and the rho rate) + # so the per-iteration ramp can no longer re-inflate lrate back to the + # overshooting value -- without this the ramp and a one-shot halving + # just oscillate and the LL drifts down. + # + # rholrate is a CEILING here, not a per-decrease-annealed working + # rate. Fortran resets rholrate=rholrate0 each iteration before the + # rho update (amica15.f90:1788) and only tightens the rholrate0 + # ceiling at maxdecs (amica15.f90:1050, gated on iter>newt_start), so + # its per-decrease rholrate*=rholratefact (:1045) is always overwritten + # by the reset and never reaches the rho update. rho has no ramp, so + # self.rholrate carries that ceiling directly (reset to rholrate0 each + # fit, nothing re-inflates it) and must ratchet ONLY at maxdecs. The + # previous per-decrease self.rholrate*=rholratefact was a monotone + # decay with no reset that collapsed the rho rate to ~1e-5 within a few + # hundred iterations and froze rho at a stale shape (issue #193). if len(self.ll_history) > 1 and ll < self.ll_history[-2]: if self.lrate <= self.minlrate: logger.warning( @@ -1840,10 +1852,11 @@ def fit( self.stop_reason = "lrate_floor" break self.lrate *= self.lratefact - self.rholrate *= self.rholratefact numdecs += 1 if numdecs >= self.maxdecs: self.lrate_cap *= self.lratefact + if it > self.newt_start: + self.rholrate *= self.rholratefact if self.do_newton and it > self.newt_start: self.newtrate *= self.lratefact numdecs = 0