Covariance is formed by inverting the normal matrix, so it comes back indefinite on ill-conditioned short arcs

A converged fit can return a covariance whose position block has a negative trace, with flag = 0 and no indication that anything is wrong. In a re-fit of the MPC catalog this happens on about 5 per cent of the converged fits in one short-arc population (40 of 768), at N = 4-60 and a median reduced chi-square of 1.40 — they look like ordinary good fits.

The cause is that the step and the covariance take different routes through the same matrix.

The step avoids squaring the condition number, which is what #536 was for. orbit_fit.cpp:908-916 records why, on (6489) Golevka: cond(sqrt(W) B) = 2.6e10 becomes cond(B^T W B) = 9.2e19, past what double precision resolves, and the solve silently rank-truncates. So the step is solved as a damped least-squares problem by Householder QR:

orbit_fit.cpp:938    dX = A.householderQr().solve(rhs);      // A = [sqrt(W) B ; sqrt(lambda) I]

The covariance is not:

orbit_fit.cpp:887    C = Bt * W * B;
orbit_fit.cpp:1318   cov = C.inverse();

C.inverse() is the squared-condition-number path the QR was introduced to avoid. Where cond(sqrt(W) B) reaches ~1e10 the inverse is dominated by rounding and is no longer positive definite, so cov is emitted with negative position variances.

The marginal cases show it is conditioning rather than a code path that writes bad values. One object fits identically on two runs — same N, n_rej, ndof, niter, arc — with chi-square differing by 3e-6 relative, and the sign of the position trace flips between them. Repeated fits on one node are bit-identical; the variation is between nodes at the 1e-6 level.

Two fixes, and they are independent:

1. Form the covariance from the same decomposition. With sqrt(W) B = QR (undamped), B^T W B = R^T R, so cov = R^-1 R^-T by two triangular solves, at unsquared conditioning. The R already computed is of the damped system, whose sqrt(lambda) I rows flatter it precisely when lambda is large, so this needs one QR of the undamped sqrt(W) B at convergence.

2. Regardless of (1), test the result before emitting it. An LDLT that fails, or any non-positive diagonal entry, should set a flag rather than return silently. A caller cannot detect this: the matrix is full and finite, and only the eigenvalues give it away.

Reported from outside the package, where a downstream sqrt(max(trace, 0)) had been masking it as a zero covariance.

Drafted with Claude Opus 5.

Start in orbit_fit.cpp at lines 887, 908-916, 938, and 1318, tracing how the converged fit builds C and cov versus the QR step. Reproduce the short-arc cases described for (6489) Golevka and inspect covariance validity under repeated fits. Done means ill-conditioned results cannot silently emit an indefinite covariance and converged-fit behavior is covered by regression checks.