Cholesky uses an absolute pivot threshold and `Factor()` masks failure

[gabrielfrasantos]

Severity: medium
Domain: math
Status: VERIFIED — read against a210d34 on 2026-08-10
Suggested labels: bug, numerical-robustness, math

Summary

Two related weaknesses in CholeskyDecomposition:

1. The positive-definiteness test uses a fixed absolute threshold, so a well-conditioned but
   small-magnitude SPD matrix is falsely rejected.
2. Factor() converts a detected failure into a zero matrix via value_or, destroying the status
   that TryFactor() correctly produced.

Location

numerical/math/CholeskyDecomposition.hpp

Evidence

const float sum = ToFloat(a.at(i, j)) - detail::CholeskyInnerProduct(l, i, j);

if (i != j)
    l.at(i, j) = T(sum / ToFloat(l.at(j, j)));
else if (sum < 1e-10f)              // absolute, scale-dependent
    return std::nullopt;
else
    l.at(i, j) = T(math::Sqrt(sum));

template<typename T, std::size_t N>
constexpr SquareMatrix<T, N> CholeskyDecomposition<T, N>::Factor(const SquareMatrix<T, N>& a)
{
    return TryFactor(a).value_or(SquareMatrix<T, N>{});   // silent zero matrix
}

Consequences

False rejection. A = 1e-6 · I is SPD with condition number exactly 1, but
sum = 1e-6 ≥ 1e-10 passes while A = 1e-11 · I (also perfectly conditioned) is rejected. Whether
a matrix factorises depends on its physical units, not its conditioning.

Failure masking. A caller of Factor() cannot distinguish "Cholesky of the zero matrix" from
"factorisation failed". Callers that then solve with this factor divide by zero.

Suggested fix

1. Scale the threshold relative to the pivot magnitude, e.g. reject when
   sum <= eps * a.at(i,i) (or track max|diag| and compare against eps * maxDiag).
2. Deprecate or remove Factor(); if it must stay, have it assert rather than silently return a
   zero matrix.

Notes

TryFactor/Solve returning std::optional is the correct pattern and is used properly by
ConsistencyMetrics. This issue is about the threshold and the value_or escape hatch only.

Grep for other value_or( uses on factorisation results — the same masking pattern may exist
elsewhere.