From 7893fb6ccb5c8e80e6505d0324572d420fd3c7be Mon Sep 17 00:00:00 2001 From: matthewholman Date: Thu, 3 Sep 2026 18:22:36 -0400 Subject: [PATCH] Reject a non-grav fit only when the covariance cannot support an amplitude The weak-constraint guard rejected a joint state + non-grav fit on the reciprocal condition number of the whole npar x npar correlation matrix. That is not a statement about the non-grav column. A six-parameter orbit fit is itself strongly correlated, so the test was dominated by the orbit underneath: on (6489) Golevka the state block alone has rcond 2.6e-10, twenty-five times below the old 1e-8 threshold, while A2 correlates with the state at only 0.91 to 0.93 and its formal uncertainty matches the published one. Radar makes this worse rather than better, because it tightens the state and lowers that rcond further. Adding seven radar observations to 2062 Aten's 980 optical ones turned a reported A2 = -14.96 +/- 1.70 into NaN. All three radar-bearing Yarkovsky objects were being suppressed at reduced chi-square between 0.28 and 1.55. The concern the old test stood in for -- a contaminated, over-confident amplitude -- was a symptom of taking the step from the squared-condition normal matrix, and does not survive that fix. A ten-day arc now returns A2 with a formal error five million times its own value, which tells the caller what it needs to know; unlike a significance test it also does not suppress a well-constrained null. So reject only on a non-positive or non-finite variance, which means the inverse has genuinely lost the column. The two tests that encoded the old contract now assert the stronger one: that the uncertainty is honest, and that it tracks arc length. Co-Authored-By: Claude Opus 5 --- src/lib/orbit_fit/orbit_fit.cpp | 54 ++++++++++++++++---------------- tests/layup/test_nongrav_a2.py | 55 ++++++++++++++++++++++++++------- 2 files changed, 71 insertions(+), 38 deletions(-) diff --git a/src/lib/orbit_fit/orbit_fit.cpp b/src/lib/orbit_fit/orbit_fit.cpp index bdf7c67b..fbf80669 100644 --- a/src/lib/orbit_fit/orbit_fit.cpp +++ b/src/lib/orbit_fit/orbit_fit.cpp @@ -87,12 +87,6 @@ namespace py = pybind11; namespace orbit_fit { - // Minimum reciprocal condition number of the normal-matrix CORRELATION matrix - // for a non-grav (A2) fit to be considered well constrained (issue #351). - // Below this the A2 column is effectively collinear with the state (typical on - // short arcs) and the fit is reported as weakly constrained (flag = 6). - static constexpr double WEAK_NONGRAV_RCOND = 1e-8; - // Arcseconds per radian (180*3600/pi). Converts astrometric/rate // uncertainties from arcseconds to radians and scales residuals for display. static constexpr double ARCSEC_PER_RAD = 206265.0; @@ -1314,31 +1308,37 @@ namespace orbit_fit cov = C.inverse(); - // Weak-constraint guard for the non-grav fit (issue #351). On a short arc - // a non-grav column becomes nearly collinear with the state directions, so - // the joint solve yields garbage params and a contaminated state. Detect - // this via the conditioning of the CORRELATION matrix of the normal matrix - // C -- the raw rcond(C) is useless here because C mixes AU and AU/day^2 - // scales, but the correlation matrix (C normalized by its diagonal) is - // unit-invariant and measures genuine collinearity. If it is near-singular, - // mark the fit weakly constrained (flag = 6) so the driver falls back to - // the 6-parameter solution and reports the non-grav params as NaN. + // Degeneracy guard for the non-grav fit (issue #351). Report the fit as + // weakly constrained (flag = 6), so the driver falls back to the + // six-parameter solution, only when the covariance cannot support an + // amplitude at all: a non-positive or non-finite variance means the + // inverse has lost that column. + // + // This used to reject on the conditioning of the whole npar x npar + // correlation matrix, which is not a statement about the non-grav column. + // A six-parameter orbit fit is itself strongly correlated, so the test was + // dominated by the orbit underneath: on (6489) Golevka the state block + // alone has rcond 2.6e-10, twenty-five times below the old 1e-8 threshold, + // while A2 correlates with the state at only 0.91 to 0.93. Radar tightens + // the state and lowers that rcond further, so adding seven radar + // observations to 2062 Aten's 980 optical ones turned a reported + // A2 = -14.96 +/- 1.70 into NaN. Better data made the fit be rejected. + // + // The concern the old test was standing in for -- a contaminated, + // over-confident amplitude -- was a symptom of the step being taken from + // the squared-condition normal matrix, and does not survive that fix. A + // ten-day arc now returns A2 with an uncertainty five million times its + // own value, which tells the caller exactly what it needs to know, and + // unlike a significance test it does not suppress a well-constrained null. if (nactive > 0 && flag == 0) { - Eigen::VectorXd diag = C.diagonal(); - if ((diag.array() > 0.0).all()) + for (int j = 6; j < npar; j++) { - Eigen::VectorXd s = diag.array().rsqrt(); // 1/sqrt(C_ii) - Eigen::MatrixXd R = s.asDiagonal() * C * s.asDiagonal(); - Eigen::JacobiSVD svd(R); - const Eigen::VectorXd &sv = svd.singularValues(); - double rcond = sv(sv.size() - 1) / sv(0); - if (!(rcond >= WEAK_NONGRAV_RCOND)) // also catches NaN + if (!(cov(j, j) > 0.0) || !std::isfinite(cov(j, j))) + { flag = 6; - } - else - { - flag = 6; // non-positive variance -> degenerate + break; + } } } diff --git a/tests/layup/test_nongrav_a2.py b/tests/layup/test_nongrav_a2.py index baca9a5f..289402ee 100644 --- a/tests/layup/test_nongrav_a2.py +++ b/tests/layup/test_nongrav_a2.py @@ -193,14 +193,45 @@ def test_six_param_fit_unaffected_and_biased_by_a2(): assert fit6.csq > 1e3 * fit7.csq -def test_a2_weak_constraint_guard_on_short_arc(): - """On a short arc the A2 column is nearly collinear with the state, so the - joint fit is rank-deficient. The fitter must flag this (flag=6) rather than - returning a contaminated, over-confident solution.""" +def test_a2_on_short_arc_is_reported_as_undetermined(): + """A ten-day arc cannot separate A2 from the state. The fit must say so in the + uncertainty rather than by refusing: the amplitude comes back with a formal + error orders of magnitude larger than itself, and larger than any physical + Yarkovsky amplitude, which is what tells the caller it is not a detection. + + This is a stronger statement than the flag=6 this used to assert. Refusing + only says the fitter declined; an honest uncertainty says how badly the arc + constrains A2, and it does not suppress a well-constrained null.""" obs = _build_arc(arc_days=10.0, n=10) # ~10 days: A2 not separable from state ephem = get_ephem(CACHE) fit = run_from_vector_with_initial_guess(ephem, _seed(_STATE), obs, 100, _BIT["A2"]) - assert fit.flag == 6, f"expected weak-constraint flag 6, got {fit.flag}" + assert fit.flag == 0, f"expected a converged fit, got flag {fit.flag}" + assert np.isfinite(fit.a2_unc) and fit.a2_unc > 0 + # Yarkovsky amplitudes for near-Earth asteroids run to ~1e-13 au/day^2; the + # formal error here is many orders above that, and above |A2| itself. + assert fit.a2_unc > 1e-11, f"expected an uninformative sigma, got {fit.a2_unc:.3e}" + assert fit.a2_unc > 1e3 * abs(fit.a2) + + +def test_a2_uncertainty_tracks_arc_length(): + """The counterpart to the test above: the reported uncertainty is what carries + the information, and it improves by orders of magnitude with arc length. Over + four years the amplitude is recovered to well under a percent, where over ten + days it is not recovered at all -- so the short-arc result is a statement + about the arc, not a blanket refusal to report.""" + ephem = get_ephem(CACHE) + short = run_from_vector_with_initial_guess( + ephem, _seed(_STATE), _build_arc(arc_days=10.0, n=10), 100, _BIT["A2"] + ) + long = run_from_vector_with_initial_guess( + ephem, _seed(_STATE), _build_arc(arc_days=4 * 365.0, n=36), 100, _BIT["A2"] + ) + assert short.flag == 0 and long.flag == 0 + # The long arc recovers the injected amplitude; the short one does not. + assert abs(long.a2 - _TRUE_A2) < 0.01 * abs(_TRUE_A2) + assert abs(short.a2 - _TRUE_A2) > abs(_TRUE_A2) + # And says so: five orders of magnitude between the two formal errors. + assert long.a2_unc < 1e-5 * short.a2_unc def _build_arc_array(arc_days=4 * 365.0, n=36, a123=(0.0, _TRUE_A2, 0.0)): @@ -337,15 +368,17 @@ def test_orbitfit_default_schema_unchanged(): assert fit[0]["flag"] == 0 -def test_orbitfit_driver_short_arc_reports_a2_nan(): - """Through the driver, a short arc that cannot constrain A2 falls back to the - 6-parameter solution: flag is success, the state is still recovered, and a2 is - reported as NaN (not a contaminated value).""" +def test_orbitfit_driver_short_arc_reports_an_uninformative_a2(): + """Through the driver, a short arc that cannot constrain A2 still reports one, + with an uncertainty that shows it is not a detection. The state is unaffected: + the point of the earlier NaN was to protect the caller from a contaminated + orbit, and the orbit is still clean.""" data, guess = _build_arc_array(arc_days=10.0, n=10) fit = orbitfit(data, cache_dir=CACHE, initial_guess=guess, fit_nongrav=True) row = fit[0] - assert row["flag"] == 0 # 6-parameter fallback succeeded - assert np.isnan(row["a2"]) and np.isnan(row["a2_unc"]) + assert row["flag"] == 0 + assert np.isfinite(row["a2"]) and np.isfinite(row["a2_unc"]) + assert row["a2_unc"] > 1e3 * abs(row["a2"]) # The fallback state is the clean 6-parameter orbit, still near truth. state = np.array([row["x"], row["y"], row["z"], row["xdot"], row["ydot"], row["zdot"]]) pos_rel = np.linalg.norm(state[:3] - _STATE[:3]) / np.linalg.norm(_STATE[:3])