Reset rholrate per iteration to match Fortran

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

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")

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

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