Fix convergence and precision bugs in the LM and VarPro minimizers (#388)

## Fix convergence and precision bugs in the LM and VarPro minimizers

Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

RELEASE.md

@@ -36,7 +36,10 @@
   analytic model derivatives and needs far fewer function evaluations.
   ``Minimizer.LevenbergMarquardt`` is a built-in, dependency-free
   Levenberg-Marquardt solver with the analytic Jacobian; it reflects steps at
-  parameter limits and fits data cubes without per-pixel allocations. With
+  parameter limits and fits data cubes without per-pixel allocations. A fit
+  that stalls away from a minimum, or whose width collapses below the spacing
+  of the scan points, is reported as failed (all zeros) rather than returned
+  as a result, so noisy or signal-free pixels may fail. With
   ``compute_errors``, Fumili and LevenbergMarquardt report parameter errors
   from their linearised covariance instead of running Hesse. Pixels for which
   Fumili does not reach a valid minimum, which Minuit2's implementation cannot
@@ -56,7 +59,10 @@
   derivatives by the tests). Each trial point solves those exactly and the
   Levenberg-Marquardt iteration runs over the nonlinear parameters only, so
   the fit needs fewer evaluations and does not depend on start values of the
-  linear parameters. Limits on linear parameters, pixels that do not converge
+  linear parameters. When the basis functions are nearly collinear, as for a
+  polynomial on a narrow range far from x = 0, one more pass over the data
+  refines the linear solution to the precision of the other minimizers.
+  Limits on linear parameters, pixels that do not converge
   and models without the separable structure fall back to
   ``LevenbergMarquardt``. Parameter errors are the same Gauss-Newton
   estimates. The Gaussian, plateau, charge-sharing and S-curve models also
@@ -141,6 +147,14 @@
 ### Bugfixes:
 - Fitting a 3D data cube in Python without ``y_err`` now returns ``par_err``
   when ``compute_errors`` is set, as fitting a single pixel already did.
+- Samples with a zero error are ignored by every minimizer, as documented,
+  also when they hold NaN. They are now removed before the start values are
+  estimated; previously a masked NaN at a sample used by the estimate, such
+  as the first or last point of ``Pol1``, made even Migrad fail. A pixel with
+  fewer remaining samples than free parameters is reported as failed.
+- Fitting a single pixel checks that ``x``, ``y`` and ``y_err`` have the same
+  size and raises ``RuntimeError`` otherwise. Previously a shorter ``y`` or
+  ``y_err`` was read past its end.
 - The Python ``Cluster`` constructors validate that the data array holds
   exactly one value per pixel and raise ``ValueError`` otherwise. Previously a
   longer array wrote past the end of the cluster data.

docs/src/pyFit.rst

@@ -130,7 +130,11 @@ control the minimizer:
 - ``tolerance`` is the EDM tolerance in Minuit's convention. Migrad,
   LevenbergMarquardt and VarPro stop when the estimated distance
   to the minimum is below ``0.002 * tolerance``, Fumili below
-  ``1e-4 * tolerance``.
+  ``1e-4 * tolerance``. The EDM is measured in units of ``chi2``, not
+  relative to the data: without ``y_err`` every sample counts with an error
+  of 1, so a fit of small values (amplitudes well below 1) meets the
+  tolerance while still far from the minimum. Pass ``y_err`` or lower
+  ``tolerance`` for such data.
 - ``strategy`` is the Minuit2 strategy (0 is fast, 1 is Minuit2's default).
   LevenbergMarquardt and VarPro ignore it.
 - ``compute_errors`` adds ``par_err`` to the result. Fixed parameters and

include/aare/FitModel.hpp

@@ -43,7 +43,9 @@ namespace aare {
  * report an error of 0 and a failed fit returns zeros. `tolerance` is the
  * EDM tolerance in Minuit's convention (Migrad, LevenbergMarquardt and
  * VarPro stop below 0.002 * tolerance, Fumili below
- * 1e-4 * tolerance). `max_calls` bounds the work per pixel in the
+ * 1e-4 * tolerance). The EDM is in units of chi2, so an unweighted fit of
+ * small values meets it far from the minimum; pass errors or a lower
+ * tolerance for such data. `max_calls` bounds the work per pixel in the
  * minimizer's own units (Minuit2 function calls, or model evaluations for
  * LevenbergMarquardt and VarPro, which spend one evaluation per
  * iteration and two when a step crosses a limit). `strategy` only affects

python/tests/test_Fit.py

@@ -175,6 +175,15 @@ def test_invalid_parameter_access_raises():
         model.SetParLimits("mu", 1.0, 1.0)
 
 
+def test_mismatched_array_sizes_raise(gaussian_data, minimizer):
+    x, y, y_err = gaussian_data
+    model = aare.Gaussian(minimizer=minimizer)
+    with pytest.raises(RuntimeError):
+        model.fit(x, y[:-1])
+    with pytest.raises(RuntimeError):
+        model.fit(x, y, y_err[:-1])
+
+
 @pytest.mark.parametrize("cls", [aare.RisingScurve, aare.FallingScurve])
 def test_scurve_fit(cls, minimizer):
     truth = np.array([10.0, 0.1, 50.0, 3.0, 1000.0, 2.0])
@@ -239,3 +248,154 @@ def test_all_minimizers_agree_on_3d_data_with_errors():
         np.testing.assert_allclose(result["par"], np.broadcast_to(expected, (2, 3, 3)), atol=2e-3)
         np.testing.assert_allclose(result["chi2"], 0.0, atol=1e-3)
         np.testing.assert_allclose(result["par_err"], reference["par_err"], rtol=2e-2)
+
+
+# Regressions of the built-in solvers: fits that used to be reported as valid
+# although they had stopped away from the minimum or returned garbage.
+
+
+def test_gaussian_fit_to_noise_fails_or_stays_bounded(minimizer):
+    # A width collapsing between two scan points let the amplitude diverge
+    # (up to 1e159) in a fit that was still reported as valid.
+    x = np.linspace(-6.0, 6.0, 61)
+    y = np.random.default_rng(0).normal(0.0, 1.0, (20, 20, x.size))
+    result = aare.Gaussian(minimizer=minimizer).fit(x, y, n_threads=2)
+
+    par = result["par"]
+    failed = np.all(par == 0.0, axis=-1)
+    bounded = np.all(np.isfinite(par) & (np.abs(par) < 1e6), axis=-1)
+    assert np.all(failed | bounded)
+    assert np.all(np.isfinite(result["chi2"]))
+
+
+def test_step_reflected_onto_start_is_not_convergence(minimizer):
+    # The step of +2 reflects at the limit 1 back onto the start 0; the
+    # truncated step must still be tried.
+    model = aare.Pol1(minimizer=minimizer)
+    model.SetParameter("p0", 0.0)
+    model.FixParameter("p1", 0.0)
+    model.SetParLimits("p0", -1.0, 1.0)
+    result = model.fit(np.arange(4.0), np.full(4, 2.2))
+
+    assert result["par"][0] == pytest.approx(1.0, abs=1e-6)
+    assert result["chi2"][0] == pytest.approx(5.76, rel=1e-6)
+
+
+@pytest.mark.parametrize(
+    "amplitude, offset", [(1e8, 0.0), (1e10, 0.0), (120.0, 1e9)]
+)
+def test_large_parameter_does_not_stop_the_fit(minimizer, amplitude, offset):
+    # One relative step norm over all parameters let a large amplitude or
+    # position stop the fit at the start estimate of the width.
+    x = np.linspace(-6.0, 6.0, 61)
+    y = gaussian(x, amplitude, 0.8, 1.3)
+    y_err = np.sqrt(y) + 1.0
+    result = aare.Gaussian(minimizer=minimizer, max_calls=500).fit(
+        x + offset, y, y_err
+    )
+
+    par = result["par"]
+    assert par[0] == pytest.approx(amplitude, rel=1e-3)
+    assert par[1] - offset == pytest.approx(0.8, abs=1e-3)
+    assert par[2] == pytest.approx(1.3, rel=1e-3)
+
+
+def test_rejected_trial_on_limit_does_not_freeze_interior_parameter(minimizer):
+    # A rejected trial truncated onto the lower limit froze mu at its start.
+    x = np.linspace(-0.1, 0.1, 11)
+    mu = np.sqrt(2.0 * np.log(2.0))
+    y = gaussian(x, 10.0, mu, 1.0)
+    model = aare.Gaussian(minimizer=minimizer)
+    model.FixParameter("A", 10.0)
+    model.FixParameter("sigma", 1.0)
+    model.SetParameter("mu", 3.0)
+    model.SetParLimits("mu", 0.0, 5.0)
+    result = model.fit(x, y)
+
+    assert result["par"][1] == pytest.approx(mu, abs=1e-3)
+    assert result["chi2"][0] < 1e-4
+
+
+# Samples with a zero error are ignored, also when they hold NaN (a masked
+# dead frame): a zero weight does not remove a NaN from a sum.
+
+
+@pytest.mark.parametrize("masked", [4, 0, 9])
+def test_masked_nan_sample_is_ignored(minimizer, masked):
+    # The start estimate of Pol1 uses the first and the last sample.
+    x = np.arange(10.0)
+    y = 2.0 + 0.5 * x
+    y_err = np.ones_like(y)
+    y[masked] = np.nan
+    y_err[masked] = 0.0
+    result = aare.Pol1(minimizer=minimizer, compute_errors=True).fit(x, y, y_err)
+
+    np.testing.assert_allclose(result["par"], [2.0, 0.5], atol=1e-5)
+    assert np.all(result["par_err"] > 0.0)
+    assert result["chi2"][0] == pytest.approx(0.0, abs=1e-6)
+
+
+def test_too_few_unmasked_samples_fail(minimizer):
+    x = np.arange(10.0)
+    y = np.full_like(x, np.nan)
+    y_err = np.zeros_like(x)
+    y[0], y_err[0] = 2.0, 1.0  # one sample left for two free parameters
+    result = aare.Pol1(minimizer=minimizer).fit(x, y, y_err)
+    assert np.all(result["par"] == 0.0)
+    assert np.all(result["chi2"] == 0.0)
+
+
+def test_masked_samples_in_cube_match_1d_fits(minimizer):
+    x = np.linspace(-6.0, 6.0, 61)
+    y = np.stack([gaussian(x, 80.0 + 10.0 * c, 0.8, 1.2) for c in range(4)])
+    y = y.reshape(1, 4, x.size)
+    y_err = np.ones_like(y)
+    for col, masked in [(1, 10), (2, 0), (3, 34)]:
+        y[0, col, masked] = np.nan
+        y_err[0, col, masked] = 0.0
+    model = aare.Gaussian(minimizer=minimizer, compute_errors=True)
+    cube = model.fit(x, y, y_err, n_threads=1)
+
+    for col in range(4):
+        single = model.fit(x, y[0, col], y_err[0, col])
+        np.testing.assert_allclose(cube["par"][0, col], single["par"], atol=1e-12)
+        np.testing.assert_allclose(cube["par_err"][0, col], single["par_err"], atol=1e-12)
+        np.testing.assert_allclose(
+            cube["par"][0, col], [80.0 + 10.0 * col, 0.8, 1.2], rtol=1e-3
+        )
+
+
+def test_parameter_whose_minimum_lies_on_its_limit_reports_no_error(minimizer):
+    # At a minimum that lies exactly on the limit the gradient is zero.
+    x = np.array([-1.0, -1.0, 1.0, 1.0])
+    model = aare.Pol1(minimizer=minimizer, compute_errors=True)
+    model.SetParLimits("p0", 0.0, 1.0)
+    result = model.fit(x, 1.0 + 2.0 * x)
+
+    assert result["par"][0] == pytest.approx(1.0, abs=1e-6)
+    assert result["par_err"][0] == 0.0
+    assert result["par_err"][1] == pytest.approx(0.5, rel=1e-6)
+
+
+@pytest.mark.parametrize("amplitude", [0.1, 1e-4])
+def test_fit_does_not_stop_at_start_without_a_step(minimizer, amplitude):
+    # The EDM tolerance is absolute in chi2 units: an unweighted fit of small
+    # data passes it at the start estimate (sigma 1.3617) already.
+    x = np.linspace(-6.0, 6.0, 61)
+    y = gaussian(x, amplitude, 0.8, 1.3)
+    result = aare.Gaussian(minimizer=minimizer).fit(x, y)
+    assert abs(result["par"][2] - 1.3) < 0.03  # Migrad ends at 1.3205
+
+    one_call = aare.Gaussian(minimizer=minimizer, max_calls=1).fit(x, y)
+    assert np.all(one_call["par"] == 0.0)
+
+
+@pytest.mark.parametrize("lo, span", [(8000, 1000), (20000, 500), (60000, 200)])
+def test_polynomial_far_from_zero_keeps_its_precision(minimizer, lo, span):
+    # Nearly collinear columns 1, x, x^2: VarPro's normal equations lost six
+    # digits at x = 20000.
+    x = np.linspace(lo, lo + span, 40)
+    c = np.array([1.5 - 2e-3 * lo + 3e-7 * lo**2, 2e-3 - 6e-7 * lo, 3e-7])
+    y = c[0] + c[1] * x + c[2] * x**2
+    result = aare.Pol2(minimizer=minimizer).fit(x, y)
+    np.testing.assert_allclose(result["par"], c, rtol=1e-9)

src/fit/DampedGaussNewton.hpp

@@ -154,15 +154,19 @@ bool gauss_newton_errors(const Mat<N> &H, const std::array<bool, N> &skip,
  *   points outward is frozen for that iteration.
  * - Convergence in Minuit's convention: the iteration stops when the
  *   estimated distance to the minimum of chi2 (EDM, from the Gauss-Newton
- *   Hessian) drops below 0.002 * tolerance, or when the step becomes
- *   negligible.
+ *   Hessian) drops below 0.002 * tolerance. A start point that passes
+ *   already still takes one undamped step, see converged(). When the damped
+ *   step becomes negligible for every parameter before that, the iteration
+ *   has stalled and the fit is valid only if the EDM test passes at that
+ *   point.
  * - max_calls bounds the number of evaluations: the initial point and every
  *   trial point count, so an iteration costs one evaluation, or two when a
  *   limit is crossed.
  * - With reject_degenerate, a trial point at which the model has lost its
  *   sensitivity to a free parameter (its Hessian diagonal collapsed by more
- *   than eight orders of magnitude) is rejected like an inadmissible point,
- *   so the damping grows and a shorter step is tried instead.
+ *   than eight orders of magnitude relative to the current point or to the
+ *   start point) is rejected like an inadmissible point, so the damping
+ *   grows and a shorter step is tried instead.
  * - With polish, a converged fit takes one more undamped Gauss-Newton step
  *   and keeps it when it lowers chi2, see polish_step().
  *
@@ -194,6 +198,7 @@ template <int NP> class DampedGaussNewton {
                         const Vec<NP> &upper, const std::array<bool, NP> &fixed,
                         const Options &opt) {
         nfree_ = 0;
+        accepted_ = 0;
         for (int k = 0; k < NP; ++k) {
             lower_[k] = lower[k];
             upper_[k] = upper[k];
@@ -214,6 +219,8 @@ template <int NP> class DampedGaussNewton {
         if (!evaluate(eval, q_, cur_))
             return res;
         res.calls = 1;
+        for (int a = 0; a < nfree_; ++a)
+            start_diag_[a] = cur_.H[free_[a]][free_[a]];
         const double edm_target = 0.002 * opt.tolerance;
         const int max_calls = opt.max_calls;
         const bool reject_degenerate = opt.reject_degenerate;
@@ -224,7 +231,9 @@ template <int NP> class DampedGaussNewton {
             extract_free();
             mark_active();
             if (gmax_ < 1e-12) {
-                res.valid = true;
+                edm(); // the undamped step for converged()
+                res.valid =
+                    converged(eval, res, max_calls, reject_degenerate, false);
                 break;
             }
             // Solve for the damped step first. The solve yields
@@ -237,10 +246,9 @@ template <int NP> class DampedGaussNewton {
                 if (e < edm_target) {
                     // A point already deep inside the tolerance needs no
                     // polishing.
-                    if (opt.polish && e >= 0.01 * edm_target &&
-                        res.calls < max_calls)
-                        polish_step(eval, res, reject_degenerate);
-                    res.valid = true;
+                    res.valid =
+                        converged(eval, res, max_calls, reject_degenerate,
+                                  opt.polish && e >= 0.01 * edm_target);
                     break;
                 }
             }
@@ -253,8 +261,16 @@ template <int NP> class DampedGaussNewton {
                     break;
                 continue;
             }
-            double step2 = 0.0;
-            double q2 = 0.0;
+            // A stalled iteration is not a minimum: repeated rejections
+            // shrink the damped step until it is negligible anywhere. The
+            // test looks at the step before it is reflected or truncated at
+            // a limit, so a step that a limit folds back onto the current
+            // point still gets its truncated alternative evaluated.
+            if (negligible_step()) {
+                res.valid = edm() < edm_target;
+                break;
+            }
+            bool moved = false;
             bool crossed = false;
             for (int a = 0; a < nfree_; ++a) {
                 const int k = free_[a];
@@ -271,17 +287,16 @@ template <int NP> class DampedGaussNewton {
                 }
                 q_new_[a] = reflected;
                 q_alt_[a] = truncated;
-                const double d = q_new_[a] - q_[a];
-                step2 += d * d;
-                q2 += q_[a] * q_[a];
+                moved = moved || reflected != q_[a];
             }
-            if (std::sqrt(step2) <= 1e-8 * (std::sqrt(q2) + 1e-8)) {
-                res.valid = true;
-                break;
+            // A reflection that lands on the current point needs no
+            // evaluation; the truncated step below is the only candidate.
+            bool ok = false;
+            if (moved) {
+                ++res.calls;
+                ok = evaluate(eval, q_new_, trial_) &&
+                     !(reject_degenerate && degenerate(trial_));
             }
-            ++res.calls;
-            bool ok = evaluate(eval, q_new_, trial_) &&
-                      !(reject_degenerate && degenerate(trial_));
             if (crossed && res.calls < max_calls) {
                 // Also try the step truncated at the limit, keep the better.
                 ++res.calls;
@@ -300,6 +315,7 @@ template <int NP> class DampedGaussNewton {
                 std::swap(cur_, trial_);
                 for (int a = 0; a < nfree_; ++a)
                     p_[free_[a]] = q_[a];
+                ++accepted_;
                 const double t = 2.0 * rho - 1.0;
                 u *= std::max(1.0 / 3.0, 1.0 - t * t * t);
                 nu = 2.0;
@@ -369,31 +385,56 @@ template <int NP> class DampedGaussNewton {
     }
 
     // A parameter sitting on a limit with the descent direction pointing
-    // outward is frozen for this iteration. The test deliberately looks at
-    // the last evaluated point p_, which after a rejected step is the trial
-    // rather than the current point q_: a rejected trial that came off the
-    // limit unfreezes the parameter, and the reflected step can then leave
-    // the limit. Testing q_ instead was measured to leave more fits stuck on
-    // a limit with a higher chi2 (FallingScurve and GaussianChargeSharing
-    // cubes), although it helped Pol2.
+    // outward is frozen for this iteration. Position and gradient are both
+    // taken at the current point q_, so a parameter stays frozen at the end
+    // only where the constrained minimum lies on its limit. An earlier
+    // version tested the last evaluated point p_, a rejected trial, so that a
+    // trial that came off a limit unfroze the parameter: testing q_ had left
+    // more FallingScurve and GaussianChargeSharing fits stuck with the width
+    // on its lower limit and a higher chi2. But pairing a trial position with
+    // the current gradient could freeze an interior parameter and stop the
+    // fit there. The width no longer reaches its limit that way, because
+    // LevenbergMarquardt now rejects degenerate trial points as
+    // VariableProjection does, see degenerate().
     void mark_active() {
         for (int a = 0; a < nfree_; ++a) {
             const int k = free_[a];
-            active_[a] = (p_[k] <= lower_[k] && g_[a] > 0.0) ||
-                         (p_[k] >= upper_[k] && g_[a] < 0.0);
+            active_[a] = (q_[a] <= lower_[k] && g_[a] > 0.0) ||
+                         (q_[a] >= upper_[k] && g_[a] < 0.0);
         }
     }
 
+    // True when the damped step is negligible for every parameter that is
+    // not frozen: below 1e-8 of its value plus its own scale 1/sqrt(H_aa),
+    // the parameter change that raises chi2 by about one. Judging each
+    // parameter on its own scale keeps a large parameter (an amplitude of
+    // 1e10, a peak position at 1e9) from hiding the steps of the others.
+    bool negligible_step() const {
+        for (int a = 0; a < nfree_; ++a) {
+            if (active_[a])
+                continue;
+            const double scale =
+                std::abs(q_[a]) + 1.0 / std::sqrt(std::max(jtj_[a][a], 1e-300));
+            if (!(std::abs(delta_[a]) <= 1e-8 * scale))
+                return false;
+        }
+        return true;
+    }
+
     // A trial point where the Hessian diagonal of a free parameter collapsed
-    // relative to the current point: the model no longer responds to that
-    // parameter there, so the iteration could stop at a spurious stationary
-    // point (a step function with the width driven to zero, for example).
+    // relative to the current point or to the start point: the model no
+    // longer responds to that parameter there, so the iteration could stop
+    // at a spurious stationary point (a step function with the width driven
+    // to zero, for example). The start point catches a collapse spread over
+    // many accepted steps, each of which passes the test against its
+    // predecessor.
     bool degenerate(const Normal &trial) const {
         for (int a = 0; a < nfree_; ++a) {
             if (active_[a])
                 continue;
             const int k = free_[a];
-            if (!(trial.H[k][k] >= 1e-8 * cur_.H[k][k]))
+            const double ref = std::max(cur_.H[k][k], start_diag_[a]);
+            if (!(trial.H[k][k] >= 1e-8 * ref))
                 return true;
         }
         return false;
@@ -425,6 +466,25 @@ template <int NP> class DampedGaussNewton {
         return e;
     }
 
+    // The convergence test passed. A point that passes it before any step
+    // was accepted may merely lie within the tolerance of chi2's scale, as
+    // the start estimate of an unweighted fit of data of small magnitude
+    // does, so it still takes the undamped step of polish_step(), as Migrad
+    // always iterates at least once. Without the budget for that step the
+    // fit has not converged. With polish, a converged fit takes the same
+    // step when the budget allows. Returns whether the fit is valid.
+    template <typename Eval>
+    bool converged(Eval &eval, MinimizerResult &res, int max_calls,
+                   bool reject_degenerate, bool polish) {
+        const bool first = accepted_ == 0 && nfree_ > 0;
+        if (!first && !polish)
+            return true;
+        if (res.calls >= max_calls)
+            return !first;
+        polish_step(eval, res, reject_degenerate);
+        return true;
+    }
+
     // One undamped Gauss-Newton step from a point that passed the EDM test,
     // kept when it lowers F. The test stops the iteration within the
     // tolerance, but a solver that converges in few iterations stops while
@@ -507,6 +567,7 @@ template <int NP> class DampedGaussNewton {
     }
 
     int nfree_ = 0;
+    int accepted_ = 0; // trial points accepted in this fit
     std::array<int, NP> free_{};
     std::array<int, NP> free_pos_{};
     std::array<bool, NP> fixed_{};
@@ -520,6 +581,7 @@ template <int NP> class DampedGaussNewton {
     Vec<NP> g_{};
     Vec<NP> delta_{};
     Vec<NP> gn_step_{};
+    Vec<NP> start_diag_{}; // Hessian diagonal at the start point
     Mat<NP> jtj_{};
     // The two solves that bound the EDM differ by rounding only when the
     // damping is negligible; the margin keeps the stop decision exact.

src/fit/Fit.cpp

@@ -7,7 +7,9 @@
 #include "aare/Models.hpp"
 #include "aare/utils/par.hpp"
 #include "aare/utils/task.hpp"
+#include <algorithm>
 #include <array>
+#include <cmath>
 #include <cstddef>
 #include <stdexcept>
 #include <string>
@@ -141,6 +143,34 @@ template <typename Model> class PixelFitter {
     bool fit(const FitModel<Model> &model, NDView<double, 1> x,
              NDView<double, 1> y, NDView<double, 1> y_err, double *par_out,
              double *err_out, double &chi2) {
+        // Samples with a zero error are ignored. They are dropped here,
+        // before the start estimate and the minimizer: a masked sample may
+        // hold NaN, and a zero weight does not remove a NaN from a sum.
+        if (y_err.size() > 0 &&
+            std::find(y_err.begin(), y_err.end(), 0.0) != y_err.end()) {
+            keep_x_.clear();
+            keep_y_.clear();
+            keep_err_.clear();
+            for (ssize_t i = 0; i < y.size(); ++i) {
+                if (y_err[i] == 0.0)
+                    continue;
+                keep_x_.push_back(x[i]);
+                keep_y_.push_back(y[i]);
+                keep_err_.push_back(y_err[i]);
+            }
+            const auto n = static_cast<ssize_t>(keep_x_.size());
+            x = NDView<double, 1>(keep_x_.data(), {n});
+            y = NDView<double, 1>(keep_y_.data(), {n});
+            y_err = NDView<double, 1>(keep_err_.data(), {n});
+        }
+        // Fewer samples than free parameters leave the fit undetermined,
+        // and the start estimates read at least two samples.
+        ssize_t n_free = 0;
+        for (std::size_t k = 0; k < Model::npar; ++k)
+            n_free += model.is_user_fixed(static_cast<unsigned int>(k)) ? 0 : 1;
+        if (y.size() < std::max<ssize_t>(n_free, 2))
+            return fail(par_out, err_out, chi2);
+
         const auto start = start_values(model, x, y);
 
         if (model.minimizer() == Minimizer::VarPro) {
@@ -150,7 +180,7 @@ template <typename Model> class PixelFitter {
             if constexpr (model::is_separable<Model>::value) {
                 const auto res =
                     vp_.fit(model, x, y, y_err, start, par_out, err_out);
-                if (res.valid) {
+                if (res.valid && finite(par_out, err_out, res.chi2)) {
                     chi2 = res.chi2;
                     return true;
                 }
@@ -164,15 +194,40 @@ template <typename Model> class PixelFitter {
     }
 
   private:
+    // The built-in solvers do not guard every operation against NaN or
+    // infinity, so a result is checked before it is reported as a valid fit.
+    static bool finite(const double *par, const double *err, double chi2) {
+        bool ok = std::isfinite(chi2);
+        for (std::size_t k = 0; k < Model::npar; ++k)
+            ok = ok && std::isfinite(par[k]) && (!err || std::isfinite(err[k]));
+        return ok;
+    }
+
     bool fit_lm(const FitModel<Model> &model, NDView<double, 1> x,
                 NDView<double, 1> y, NDView<double, 1> y_err,
                 const std::array<double, Model::npar> &start, double *par_out,
                 double *err_out, double &chi2) {
         const auto res = lm_.fit(model, x, y, y_err, start, par_out, err_out);
-        chi2 = res.valid ? res.chi2 : 0.0;
-        return res.valid;
+        if (res.valid && finite(par_out, err_out, res.chi2)) {
+            chi2 = res.chi2;
+            return true;
+        }
+        return fail(par_out, err_out, chi2);
     }
 
+    // A failed fit reports zeros.
+    static bool fail(double *par_out, double *err_out, double &chi2) {
+        std::fill(par_out, par_out + Model::npar, 0.0);
+        if (err_out)
+            std::fill(err_out, err_out + Model::npar, 0.0);
+        chi2 = 0.0;
+        return false;
+    }
+
+    // The samples of the current pixel that have a nonzero error, see fit().
+    std::vector<double> keep_x_;
+    std::vector<double> keep_y_;
+    std::vector<double> keep_err_;
     LevenbergMarquardt<Model> lm_;
     std::conditional_t<model::is_separable<Model>::value,
                        VariableProjection<Model>, NotSeparable>
@@ -185,6 +240,7 @@ template <typename Model>
 NDArray<double, 1> fit_pixel(const FitModel<Model> &model, NDView<double, 1> x,
                              NDView<double, 1> y, NDView<double, 1> y_err) {
     constexpr std::size_t npar = Model::npar;
+    detail::check_fit_pixel_shapes(x, y, y_err);
     const bool want_errors = model.compute_errors();
     const auto result_size =
         static_cast<ssize_t>(want_errors ? (2 * npar + 1) : (npar + 1));

src/fit/Fit.test.cpp

@@ -9,6 +9,7 @@
 #include <catch2/generators/catch_generators.hpp>
 #include <cmath>
 #include <limits>
+#include <random>
 #include <stdexcept>
 #include <string>
 #include <vector>
@@ -439,6 +440,21 @@ TEST_CASE("FitModel validates parameter names, indices and limits", "[fit]") {
     CHECK(model.minimizer() == Minimizer::LevenbergMarquardt);
 }
 
+TEST_CASE("fit_pixel rejects arrays of different sizes", "[fit]") {
+    const auto minimizer =
+        GENERATE(Minimizer::Migrad, Minimizer::Fumili,
+                 Minimizer::LevenbergMarquardt, Minimizer::VarPro);
+    INFO("minimizer " << name(minimizer));
+    auto d = gaussian_data();
+    const auto model = make_model<aare::model::Gaussian>(minimizer);
+    NDView<double, 1> short_y(d.y.data(), {n_points - 1});
+    NDView<double, 1> short_err(d.y_err.data(), {n_points - 1});
+    CHECK_THROWS_AS(aare::fit_pixel(model, d.x.view(), short_y),
+                    std::runtime_error);
+    CHECK_THROWS_AS(aare::fit_pixel(model, d.x.view(), d.y.view(), short_err),
+                    std::runtime_error);
+}
+
 TEST_CASE("Copying a FitModel keeps the minimizer", "[fit]") {
     const aare::FitModel<aare::model::Pol1> original(
         0, 100, 0.5, false, Minimizer::LevenbergMarquardt);
@@ -925,3 +941,368 @@ TEST_CASE("A parameter ending exactly on a limit reports an error of 0",
     CHECK(res(2) == 0.0);
     CHECK(res(3) == 0.0);
 }
+
+// ---------------------------------------------------------------------------
+// Regressions of the built-in solvers: fits that used to be reported as valid
+// although they had stopped away from the minimum or returned garbage.
+// ---------------------------------------------------------------------------
+
+namespace {
+
+// Standard normal noise from std::mt19937, whose output the standard fixes,
+// and a Box-Muller transform written out here: std::normal_distribution
+// differs between standard libraries.
+std::vector<double> normal_noise(std::size_t n, unsigned seed) {
+    std::mt19937 rng(seed);
+    const auto uniform = [&rng] {
+        return (static_cast<double>(rng()) + 0.5) / 4294967296.0;
+    };
+    std::vector<double> out(n);
+    for (std::size_t i = 0; i < n; i += 2) {
+        const double r = std::sqrt(-2.0 * std::log(uniform()));
+        const double phi = 2.0 * std::acos(-1.0) * uniform();
+        out[i] = r * std::cos(phi);
+        if (i + 1 < n)
+            out[i + 1] = r * std::sin(phi);
+    }
+    return out;
+}
+
+} // namespace
+
+TEST_CASE("A Gaussian fit to pure noise fails or stays finite and bounded",
+          "[fit]") {
+    // Fitting noise, a narrow width can collapse between two scan points;
+    // the amplitude that fits a single point then diverges. Such a fit must
+    // be reported as failed, not as a valid result with |A| up to 1e159.
+    const auto minimizer =
+        GENERATE(Minimizer::Migrad, Minimizer::Fumili,
+                 Minimizer::LevenbergMarquardt, Minimizer::VarPro);
+    INFO("minimizer " << name(minimizer));
+    constexpr ssize_t rows = 20;
+    constexpr ssize_t cols = 20;
+    NDArray<double, 1> x({n_points});
+    for (ssize_t i = 0; i < n_points; ++i)
+        x(i) = -6.0 + 0.2 * static_cast<double>(i);
+    NDArray<double, 3> y({rows, cols, n_points});
+    const auto noise = normal_noise(static_cast<std::size_t>(y.size()), 1);
+    std::copy(noise.begin(), noise.end(), y.begin());
+
+    NDArray<double, 3> par({rows, cols, 3});
+    NDArray<double, 2> chi2({rows, cols});
+    const auto model = make_model<aare::model::Gaussian>(minimizer);
+    aare::fit_3d(model, x.view(), y.view(), NDView<double, 3>{}, par.view(),
+                 NDView<double, 3>{}, chi2.view(), 1);
+
+    int absurd = 0;
+    for (ssize_t row = 0; row < rows; ++row) {
+        for (ssize_t col = 0; col < cols; ++col) {
+            bool ok = std::isfinite(chi2(row, col));
+            for (ssize_t k = 0; k < 3; ++k)
+                ok = ok && std::isfinite(par(row, col, k)) &&
+                     std::abs(par(row, col, k)) < 1e6;
+            absurd += ok ? 0 : 1;
+        }
+    }
+    CHECK(absurd == 0);
+}
+
+TEST_CASE("A step reflected back onto the start point is not convergence",
+          "[fit]") {
+    // p0 starts at 0 with limits [-1, 1]; the damped step of +2 reflects at
+    // the upper limit back onto 0. The step truncated at the limit must still
+    // be tried: the constrained minimum is p0 = 1.
+    const auto minimizer =
+        GENERATE(Minimizer::Migrad, Minimizer::Fumili,
+                 Minimizer::LevenbergMarquardt, Minimizer::VarPro);
+    INFO("minimizer " << name(minimizer));
+    NDArray<double, 1> x({4});
+    NDArray<double, 1> y({4}, 2.2);
+    for (ssize_t i = 0; i < 4; ++i)
+        x(i) = static_cast<double>(i);
+
+    auto model = make_model<aare::model::Pol1>(minimizer);
+    model.SetParameter("p0", 0.0);
+    model.FixParameter("p1", 0.0);
+    model.SetParLimits("p0", -1.0, 1.0);
+    const auto res = aare::fit_pixel(model, x.view(), y.view());
+    CHECK(res(0) == Approx(1.0).margin(limit_margin(minimizer)));
+    CHECK(res(2) == Approx(5.76).epsilon(1e-6));
+}
+
+TEST_CASE("A large parameter does not hide the steps of the others", "[fit]") {
+    // The negligible-step test used one relative norm over all parameters,
+    // so an amplitude of 1e10 or a position of 1e9 stopped the fit at the
+    // start estimate of the width.
+    const auto minimizer =
+        GENERATE(Minimizer::Migrad, Minimizer::Fumili,
+                 Minimizer::LevenbergMarquardt, Minimizer::VarPro);
+    INFO("minimizer " << name(minimizer));
+    const auto [amplitude, offset] = GENERATE(
+        std::pair{1e8, 0.0}, std::pair{1e10, 0.0}, std::pair{120.0, 1e9});
+    INFO("amplitude " << amplitude << ", offset " << offset);
+
+    Data d{NDArray<double, 1>({n_points}), NDArray<double, 1>({n_points}),
+           NDArray<double, 1>({n_points})};
+    fill_gaussian(d.x, d.y, amplitude, 0.8, 1.3);
+    for (ssize_t i = 0; i < n_points; ++i) {
+        d.x(i) += offset;
+        d.y_err(i) = std::sqrt(d.y(i)) + 1.0;
+    }
+    const auto model = make_model<aare::model::Gaussian>(minimizer, false, 500);
+    const auto res =
+        aare::fit_pixel(model, d.x.view(), d.y.view(), d.y_err.view());
+    // The EDM tolerance leaves a few 1e-4 relative, far below the errors;
+    // the stuck fits were 0.7 % (A = 1e8) and 5 % (start estimate) off.
+    CHECK(res(0) == Approx(amplitude).epsilon(1e-3));
+    CHECK(res(1) - offset == Approx(0.8).margin(1e-3));
+    CHECK(res(2) == Approx(1.3).epsilon(1e-3));
+}
+
+TEST_CASE("A rejected trial on a limit does not freeze an interior parameter",
+          "[fit]") {
+    // mu starts at 3 inside [0, 5]. A trial truncated onto the lower limit is
+    // rejected; testing that trial's position against the gradient of the
+    // accepted point froze mu, and the fit stopped at the start value.
+    const auto minimizer =
+        GENERATE(Minimizer::Migrad, Minimizer::Fumili,
+                 Minimizer::LevenbergMarquardt, Minimizer::VarPro);
+    INFO("minimizer " << name(minimizer));
+    constexpr ssize_t n = 11;
+    const double mu = std::sqrt(2.0 * std::log(2.0));
+    NDArray<double, 1> x({n});
+    NDArray<double, 1> y({n});
+    for (ssize_t i = 0; i < n; ++i) {
+        x(i) = -0.1 + 0.02 * static_cast<double>(i);
+        y(i) = 10.0 * std::exp(-0.5 * (x(i) - mu) * (x(i) - mu));
+    }
+    auto model = make_model<aare::model::Gaussian>(minimizer);
+    model.FixParameter("A", 10.0);
+    model.FixParameter("sigma", 1.0);
+    model.SetParameter("mu", 3.0);
+    model.SetParLimits("mu", 0.0, 5.0);
+    const auto res = aare::fit_pixel(model, x.view(), y.view());
+    CHECK(res(1) == Approx(mu).margin(1e-3));
+    CHECK(res(3) < 1e-4);
+}
+
+// ---------------------------------------------------------------------------
+// Samples with a zero error are ignored. They may hold NaN (a masked dead
+// frame), which a zero weight does not remove from a sum, so they must not
+// reach the start estimate or the minimizer at all.
+// ---------------------------------------------------------------------------
+
+TEST_CASE("Samples with a zero error are ignored even when NaN", "[fit]") {
+    const auto minimizer =
+        GENERATE(Minimizer::Migrad, Minimizer::Fumili,
+                 Minimizer::LevenbergMarquardt, Minimizer::VarPro);
+    INFO("minimizer " << name(minimizer));
+    const double nan = std::numeric_limits<double>::quiet_NaN();
+
+    constexpr ssize_t n = 10;
+    NDArray<double, 1> x({n});
+    NDArray<double, 1> y({n});
+    NDArray<double, 1> y_err({n}, 1.0);
+    for (ssize_t i = 0; i < n; ++i) {
+        x(i) = static_cast<double>(i);
+        y(i) = 2.0 + 0.5 * x(i);
+    }
+    const auto line = make_model<aare::model::Pol1>(minimizer, true);
+
+    // The start estimate of Pol1 uses the first and the last sample.
+    const auto masked = GENERATE(as<ssize_t>{}, 4, 0, n - 1);
+    SECTION("masked NaN on a line") {
+        INFO("masked sample " << masked);
+        y(masked) = nan;
+        y_err(masked) = 0.0;
+        const auto res =
+            aare::fit_pixel(line, x.view(), y.view(), y_err.view());
+        CHECK(res(0) == Approx(2.0).epsilon(1e-5).margin(1e-6));
+        CHECK(res(1) == Approx(0.5).epsilon(1e-5).margin(1e-6));
+        CHECK(res(2) > 0.0);
+        CHECK(res(4) == Approx(0.0).margin(1e-6));
+    }
+
+    SECTION("masked NaN at the peak of a Gaussian") {
+        auto d = gaussian_data();
+        constexpr ssize_t peak = 34; // x = 0.8
+        d.y(peak) = nan;
+        d.y_err(peak) = 0.0;
+        const auto model = make_model<aare::model::Gaussian>(minimizer);
+        const auto res =
+            aare::fit_pixel(model, d.x.view(), d.y.view(), d.y_err.view());
+        const double tol = truth_tolerance(minimizer);
+        CHECK(res(0) == Approx(120.0).epsilon(tol));
+        CHECK(res(1) == Approx(0.8).epsilon(tol));
+        CHECK(res(2) == Approx(1.3).epsilon(tol));
+    }
+
+    SECTION("too few samples left fail cleanly") {
+        // One sample left for two free parameters, then none at all.
+        for (ssize_t i = 1; i < n; ++i) {
+            y(i) = nan;
+            y_err(i) = 0.0;
+        }
+        auto res = aare::fit_pixel(line, x.view(), y.view(), y_err.view());
+        for (ssize_t k = 0; k < res.size(); ++k)
+            CHECK(res(k) == 0.0);
+        y_err(0) = 0.0;
+        res = aare::fit_pixel(line, x.view(), y.view(), y_err.view());
+        for (ssize_t k = 0; k < res.size(); ++k)
+            CHECK(res(k) == 0.0);
+    }
+}
+
+TEST_CASE("Masked samples in a data cube match per-pixel fits", "[fit]") {
+    // The kept samples of each pixel go through reused buffers, which must
+    // not carry anything over from one pixel to the next.
+    const auto minimizer =
+        GENERATE(Minimizer::Migrad, Minimizer::Fumili,
+                 Minimizer::LevenbergMarquardt, Minimizer::VarPro);
+    INFO("minimizer " << name(minimizer));
+    constexpr ssize_t rows = 2;
+    constexpr ssize_t cols = 3;
+    // masked sample per pixel, -1 for none
+    constexpr ssize_t masked[rows][cols] = {{-1, 10, 0},
+                                            {34, -1, n_points - 1}};
+    NDArray<double, 1> x({n_points});
+    NDArray<double, 1> tmp({n_points});
+    NDArray<double, 3> y({rows, cols, n_points});
+    NDArray<double, 3> y_err({rows, cols, n_points}, 1.0);
+    for (ssize_t row = 0; row < rows; ++row) {
+        for (ssize_t col = 0; col < cols; ++col) {
+            fill_gaussian(x, tmp, 80.0 + 10.0 * static_cast<double>(col), 0.8,
+                          1.0 + 0.2 * static_cast<double>(row));
+            for (ssize_t i = 0; i < n_points; ++i)
+                y(row, col, i) = tmp(i);
+            if (masked[row][col] >= 0) {
+                y(row, col, masked[row][col]) =
+                    std::numeric_limits<double>::quiet_NaN();
+                y_err(row, col, masked[row][col]) = 0.0;
+            }
+        }
+    }
+
+    NDArray<double, 3> par({rows, cols, 3});
+    NDArray<double, 3> err({rows, cols, 3});
+    NDArray<double, 2> chi2({rows, cols});
+    const auto model = make_model<aare::model::Gaussian>(minimizer, true);
+    aare::fit_3d(model, x.view(), y.view(), y_err.view(), par.view(),
+                 err.view(), chi2.view(), 1);
+
+    for (ssize_t row = 0; row < rows; ++row) {
+        for (ssize_t col = 0; col < cols; ++col) {
+            INFO("pixel " << row << ", " << col);
+            NDView<double, 1> yv(&y(row, col, 0), {n_points});
+            NDView<double, 1> ev(&y_err(row, col, 0), {n_points});
+            const auto single = aare::fit_pixel(model, x.view(), yv, ev);
+            for (ssize_t k = 0; k < 3; ++k) {
+                CHECK(par(row, col, k) == Approx(single(k)).margin(1e-12));
+                CHECK(err(row, col, k) == Approx(single(3 + k)).margin(1e-12));
+            }
+            CHECK(chi2(row, col) == Approx(single(6)).margin(1e-12));
+            // Close to the truth within the EDM tolerance; the exact check
+            // is the agreement with the per-pixel fit above.
+            const double tol = 1e-3;
+            CHECK(par(row, col, 0) ==
+                  Approx(80.0 + 10.0 * static_cast<double>(col)).epsilon(tol));
+            CHECK(par(row, col, 1) == Approx(0.8).epsilon(tol));
+            CHECK(par(row, col, 2) ==
+                  Approx(1.0 + 0.2 * static_cast<double>(row)).epsilon(tol));
+        }
+    }
+}
+
+TEST_CASE("A parameter whose minimum lies on its limit reports no error",
+          "[fit]") {
+    // At a minimum that lies exactly on a limit the gradient is zero, so a
+    // test for an outward gradient does not see the parameter as bound.
+    const auto minimizer =
+        GENERATE(Minimizer::Migrad, Minimizer::Fumili,
+                 Minimizer::LevenbergMarquardt, Minimizer::VarPro);
+    INFO("minimizer " << name(minimizer));
+
+    SECTION("linear parameter") {
+        NDArray<double, 1> x({4});
+        NDArray<double, 1> y({4});
+        const double xs[4] = {-1.0, -1.0, 1.0, 1.0};
+        for (ssize_t i = 0; i < 4; ++i) {
+            x(i) = xs[i];
+            y(i) = 1.0 + 2.0 * xs[i];
+        }
+        auto model = make_model<aare::model::Pol1>(minimizer, true);
+        model.SetParLimits("p0", 0.0, 1.0);
+        const auto res = aare::fit_pixel(model, x.view(), y.view());
+        CHECK(res(0) == Approx(1.0).margin(limit_margin(minimizer)));
+        CHECK(res(2) == 0.0);
+        CHECK(res(3) == Approx(0.5).epsilon(1e-6)); // 1 / sqrt(sum x^2)
+    }
+
+    SECTION("nonlinear parameter") {
+        auto d = gaussian_data();
+        auto model = make_model<aare::model::Gaussian>(minimizer, true);
+        model.SetParLimits("mu", -5.0, 0.8);
+        const auto res =
+            aare::fit_pixel(model, d.x.view(), d.y.view(), d.y_err.view());
+        CHECK(res(1) == Approx(0.8).margin(1e-6));
+        CHECK(res(4) == 0.0);
+        CHECK(res(3) > 0.0);
+        CHECK(res(5) > 0.0);
+    }
+}
+
+TEST_CASE("A fit does not stop at its start point without taking a step",
+          "[fit]") {
+    // The EDM tolerance is absolute in chi2 units. An unweighted fit of data
+    // of small magnitude passes it at the start estimate already; the fit
+    // must still take a step, as Migrad always iterates at least once.
+    const auto minimizer =
+        GENERATE(Minimizer::Migrad, Minimizer::Fumili,
+                 Minimizer::LevenbergMarquardt, Minimizer::VarPro);
+    INFO("minimizer " << name(minimizer));
+    const double amplitude = GENERATE(0.1, 1e-4);
+    INFO("amplitude " << amplitude);
+    Data d{NDArray<double, 1>({n_points}), NDArray<double, 1>({n_points}),
+           NDArray<double, 1>({n_points}, 1.0)};
+    fill_gaussian(d.x, d.y, amplitude, 0.8, 1.3);
+
+    // The start estimate of sigma is 1.3617; Migrad ends at 1.3205.
+    const auto model = make_model<aare::model::Gaussian>(minimizer);
+    const auto res = aare::fit_pixel(model, d.x.view(), d.y.view());
+    CHECK(std::abs(res(2) - 1.3) < 0.03);
+
+    // Without the budget for a single step the fit has not converged.
+    const auto one_call =
+        make_model<aare::model::Gaussian>(minimizer, false, 1);
+    const auto none = aare::fit_pixel(one_call, d.x.view(), d.y.view());
+    for (ssize_t k = 0; k < none.size(); ++k)
+        CHECK(none(k) == 0.0);
+}
+
+TEST_CASE("A polynomial far from x = 0 keeps its precision", "[fit]") {
+    // On a narrow range far from 0 the columns 1, x, x^2 are nearly
+    // collinear. VarPro solves them from the normal equations and lost six
+    // digits at x = 20000; from about x = 60000 it falls back to LM.
+    const auto minimizer =
+        GENERATE(Minimizer::Migrad, Minimizer::Fumili,
+                 Minimizer::LevenbergMarquardt, Minimizer::VarPro);
+    INFO("minimizer " << name(minimizer));
+    const auto [lo, span] =
+        GENERATE(std::pair{8000.0, 1000.0}, std::pair{20000.0, 500.0},
+                 std::pair{60000.0, 200.0});
+    INFO("x in [" << lo << ", " << lo + span << "]");
+    // 1.5 + 2e-3 (x - lo) + 3e-7 (x - lo)^2 expanded in powers of x.
+    const double c[3] = {1.5 - 2e-3 * lo + 3e-7 * lo * lo, 2e-3 - 6e-7 * lo,
+                         3e-7};
+    constexpr ssize_t n = 40;
+    NDArray<double, 1> x({n});
+    NDArray<double, 1> y({n});
+    for (ssize_t i = 0; i < n; ++i) {
+        x(i) = lo + span * static_cast<double>(i) / static_cast<double>(n - 1);
+        y(i) = c[0] + c[1] * x(i) + c[2] * x(i) * x(i);
+    }
+    const auto model = make_model<aare::model::Pol2>(minimizer);
+    const auto res = aare::fit_pixel(model, x.view(), y.view());
+    for (ssize_t k = 0; k < 3; ++k)
+        CHECK(res(k) == Approx(c[k]).epsilon(1e-9));
+}

src/fit/FitHelpers.hpp

@@ -35,6 +35,16 @@ std::array<double, Model::npar> start_values(const FitModel<Model> &model,
     return start;
 }
 
+/** @brief Validate the array sizes passed to fit_pixel. */
+inline void check_fit_pixel_shapes(NDView<double, 1> x, NDView<double, 1> y,
+                                   NDView<double, 1> y_err) {
+    if (x.size() != y.size())
+        throw std::runtime_error("fit_pixel: x.size() must match y.size().");
+    if (y_err.size() > 0 && y_err.size() != y.size())
+        throw std::runtime_error(
+            "fit_pixel: y and y_err must have identical size.");
+}
+
 /** @brief Validate the array shapes passed to fit_3d. */
 template <typename Model>
 void check_fit_3d_shapes(NDView<double, 1> x, NDView<double, 3> y,

src/fit/LevenbergMarquardt.hpp

@@ -26,8 +26,9 @@ namespace aare::detail {
  * equations g = J^T r and J^T J, which DampedGaussNewton iterates on: it
  * handles fixed parameters, limits, damping and the convergence test, see
  * there. Rejected trials are rare with its damping rule, and a rejected one
- * only wastes a single pass. As in Minuit2, max_calls = 0 selects the
- * default budget 200 + 100 * npar + 5 * npar^2.
+ * only wastes a single pass. Trial points at which the model loses its
+ * sensitivity to a free parameter are rejected. As in Minuit2, max_calls = 0
+ * selects the default budget 200 + 100 * npar + 5 * npar^2.
  *
  * Errors are Gauss-Newton estimates, sqrt(diag((J^T J)^-1)), over the free
  * parameters that do not sit on a limit. Fixed and limit-bound parameters
@@ -95,6 +96,11 @@ template <typename Model> class LevenbergMarquardt {
         // Start cautiously: the estimated start values of the linear
         // parameters can be far off, and the first steps are then wild.
         opt.damping0 = 0.1;
+        // A width driven onto its lower limit turns a model into a step
+        // function whose chi2 no longer responds to the width and jumps
+        // with the position; the iteration then stalls there. Such trial
+        // points are rejected, as in VariableProjection.
+        opt.reject_degenerate = true;
         const Result res = driver_.fit(eval, start, lower, upper, fixed, opt);
         if (!res.valid) {
             std::fill(par_out, par_out + npar, 0.0);

src/fit/VariableProjection.hpp

@@ -101,6 +101,7 @@ template <typename Model> class VariableProjection {
         }
         pvec_.assign(start.begin(), start.end());
         have_last_ = false;
+        have_ref_ = false;
         nfree_lin_ = 0;
         for (int j = 0; j < nlin; ++j) {
             const auto k = static_cast<unsigned int>(Traits::linear_index(j));
@@ -195,8 +196,9 @@ template <typename Model> class VariableProjection {
   private:
     // F = chi2 / 2, gradient and Kaufman Hessian of the projected residual
     // at the nonlinear point alpha, and the linear solution in out.aux.
-    // Returns false when the model rejects the point or the basis functions
-    // are linearly dependent there.
+    // Returns false when the model rejects the point, or when the basis
+    // functions are numerically linearly dependent there or one of them has
+    // vanished on the scan points, see reduce().
     bool evaluate(const Vec<ndrv> &alpha, Normal &out) {
         for (int k = 0; k < nnl; ++k)
             pvec_[Traits::nonlinear_par[k]] = alpha[k];
@@ -293,13 +295,64 @@ template <typename Model> class VariableProjection {
             for (int b = 0; b < nfree_lin_; ++b)
                 A[a][b] = M[ja][free_lin_[b]];
         }
+        // A free basis function that has all but vanished on the scan points
+        // (a Gaussian narrower than the point spacing, centred between two
+        // points) leaves its coefficient undetermined, and the solve returns
+        // an absurd one that fits a single point. Its squared norm is compared
+        // with that at the pixel's first evaluation, the start point.
+        if (!have_ref_) {
+            for (int j = 0; j < nlin; ++j)
+                ref_norm2_[j] = M[j][j];
+            have_ref_ = true;
+        }
+        for (int a = 0; a < nfree_lin_; ++a) {
+            const int j = free_lin_[a];
+            if (!(M[j][j] >= vanished_basis * ref_norm2_[j]))
+                return false;
+        }
         Mat<nlin> L{};
         if (!cholesky_factor<nlin>(nfree_lin_, A, L))
             return false;
+        // Nearly collinear basis functions: a pivot of the factorisation far
+        // below its diagonal element means that the linear solution loses
+        // about log10(diagonal / pivot) digits, beyond which it is noise.
+        bool refine = false;
+        for (int a = 0; a < nfree_lin_; ++a) {
+            if (!(L[a][a] * L[a][a] >= collinear_basis * A[a][a]))
+                return false;
+            refine = refine || L[a][a] * L[a][a] < refine_basis * A[a][a];
+        }
         Vec<nlin> beta_free{};
         cholesky_substitute<nlin>(nfree_lin_, L, v, beta_free);
         for (int a = 0; a < nfree_lin_; ++a)
             pvec_[Traits::linear_index(free_lin_[a])] = beta_free[a];
+        // Short of that, the digits lost to the normal equations are
+        // recovered by one step of iterative refinement on the data (the
+        // corrected semi-normal equations): the correction solves the same
+        // equations for Phi^T W r of the residual of the first solution. A
+        // polynomial fitted far from x = 0 needs it; x in [20000, 20500]
+        // otherwise costs a Pol2 six digits.
+        if (refine) {
+            Vec<nlin> beta;
+            for (int j = 0; j < nlin; ++j)
+                beta[j] = pvec_[Traits::linear_index(j)];
+            Vec<nlin> s{};
+            for (ssize_t i = 0; i < n; ++i) {
+                double r = yd[i];
+                for (int j = 0; j < nlin; ++j)
+                    r -= beta[j] * pc[j][i];
+                const double rw = Weighted ? r * w2[i] : r;
+                for (int j = 0; j < nlin; ++j)
+                    s[j] += pc[j][i] * rw;
+            }
+            Vec<nlin> sf{};
+            for (int a = 0; a < nlin; ++a)
+                sf[a] = a < nfree_lin_ ? s[free_lin_[a]] : 0.0;
+            Vec<nlin> delta{};
+            cholesky_substitute<nlin>(nfree_lin_, L, sf, delta);
+            for (int a = 0; a < nfree_lin_; ++a)
+                pvec_[Traits::linear_index(free_lin_[a])] += delta[a];
+        }
         Vec<nlin> beta;
         for (int j = 0; j < nlin; ++j) {
             beta[j] = pvec_[Traits::linear_index(j)];
@@ -438,6 +491,14 @@ template <typename Model> class VariableProjection {
     Mat<nlin> last_M_{};
     std::array<Vec<nlin>, ndrv> last_Q_{};
     Mat<ndrv> last_N_{};
+    // Squared norms of the basis functions at the start point, see reduce().
+    bool have_ref_ = false;
+    Vec<nlin> ref_norm2_{};
+    static constexpr double vanished_basis = 1e-8;
+    static constexpr double collinear_basis = 1e-10;
+    // Pivot ratio below which the linear solution is refined, about four
+    // lost digits.
+    static constexpr double refine_basis = 1e-4;
     int nfree_lin_ = 0;
     std::array<int, nlin> free_lin_{};
     std::array<bool, nlin> lin_fixed_{};