fix: form the reconstruction covariance via Cholesky, not an elementwise sqrt

[Jammy2211]

Closes #468.

The defect

AbstractInversion.reconstruction_noise_map_with_covariance was:

return np.sqrt(np.linalg.inv(self.curvature_reg_matrix))

np.sqrt is applied elementwise to the whole inverse. That inverse is the reconstruction covariance matrix, whose off-diagonals are covariances and are routinely negative — so every one was NaN by construction, for any input matrix however well-conditioned, with a RuntimeWarning on every call. The docstring promised a matrix that "accounts for the covariance of the noise between pixels"; the entries carrying that covariance were exactly the broken ones.

Found during the reproduction gate for #467.

This did not affect the 1D reconstruction_noise_map. np.sqrt is elementwise, so it commutes with taking the diagonal (diagonal(sqrt(C))[i] == sqrt(C[i,i])). Only the covariance-aware consumer and the warning spam were hit.

Changes

• new reconstruction_covariance_matrix — returns C itself, formed by Cholesky (cho_factor/cho_solve), input and output symmetrized, with an explicit finiteness guard.
• reconstruction_noise_map — decoupled, now sqrt(diag(C)) directly. It was previously correct only incidentally, via the elementwise sqrt; the invariant is now stated, so it cannot silently become a variance if the matrix changes.
• reconstruction_noise_map_with_covariance — deprecated alias whose warning states the value change explicitly.

Why Cholesky — on measured evidence, not assumption

An A/B was run before writing the fix, and it refuted two of the arguments originally made for this change:

Claim	Verdict
inv gives negative diagonals on well-formed SPD	Refuted — 0 across cond 1e3–1e15, n=400, 20 trials each
inv is materially less accurate on the diagonal	Refuted — matches cho_solve; at cond 1e15 inv was marginally better
Near-coincident mesh vertices degrade the inverse	Refuted — with regularization the matrix stays PD (cond ~6.8e7 even at exactly duplicated columns)
inv returns asymmetric output	Confirmed — 5.2e-7 at cond 1e12, vs 2.6e-16
inv silently succeeds on indefinite matrices	Confirmed

The case for Cholesky is detection, not accuracy. cho_factor raises LinAlgError on a negative eigenvalue; np.linalg.inv raises only on an exactly singular matrix and otherwise returns a plausible-looking covariance. At eigenvalue -1e-8 all 300 diagonals came back negative (whole noise map NaN); at -1.0, zero did — no NaN, no warning, no error, and wrong numbers. A covariance is only defined for a positive-definite matrix, so this adds a definiteness check the code lacked entirely.

Not established: that a real curvature_reg_matrix is indefinite in a converged fit. Settings.no_regularization_add_to_curvature_diag_value and the curvature_matrix_with_added_to_diag_from docstring say non-PD matrices happen, but no fit was instrumented to confirm it.

Review fixes (second commit)

An adversarial review of the first commit found a regression it had introduced:

• Non-finite input raised ValueError, not LinAlgError. scipy defaults to check_finite=True. Both callers (inversion_plots.py:169, :397) catch only LinAlgError, and the CSV writer's docstring explicitly promises a failure there may not abort the enclosing model-fit. Verified: old code returned [nan, nan]; first fix raised ValueError past the guards. Now raises LinAlgError with a message naming the cause, so the contract holds for downstream callers too. check_finite=False is then passed to scipy, so the explicit check costs nothing net.
• cho_factor reads only the upper triangle, so an asymmetric input was silently inverted as though its lower triangle matched — [[2.0, 0.5], [0.1, 2.0]] gave diag 0.5333 vs the true inverse's 0.5063. Output symmetrization does not address this; the input is now symmetrized too.

API and behaviour changes

reconstruction_noise_map_with_covariance values change — diagonal from std-dev to variance, off-diagonals from NaN to covariances. Grepped rather than assumed safe:

Repo	with_covariance	reconstruction_noise_map
PyAutoGalaxy 3ca31bf	none	none
PyAutoLens 87e5827	none	none
autolens_workspace	none	4 scripts + notebooks

Control greps confirm the checkouts are real (394 / 237 / 465 .py files), so the nulls are genuine. No consumer of the alias exists in the stack. reconstruction_noise_map is used — source_science.py computes signal_to_noise_map = reconstruction / reconstruction_noise_map — and its values are unchanged. Not checked: autogalaxy_workspace, HowTo repos, external user code.

reconstruction_noise_map is numerically, not bitwise, equivalent. Algebraically identical, but the covariance is now formed by Cholesky rather than LU: ~7e-15 relative at cond 1e3 rising to ~4e-5 at cond 1e13. Neither result is the more correct one.

reconstruction_covariance_matrix now raises on indefinite/non-finite input where inv silently returned numbers. Both existing call sites already catch LinAlgError — verified, not assumed.

Tests

The only prior assertion on this property checked [0, 0], a diagonal element — which is why the off-diagonal NaNs shipped. Added:

• off-diagonals finite and genuinely negative, under -W error::RuntimeWarning (fails on the old implementation)
• accuracy against an exactly-constructed ground truth A = Q diag(w) Q.T, plus symmetry
• reconstruction_noise_map == sqrt(diag(C)) as an explicit invariant
• raises on an indefinite matrix, asserting inv is silent on the same input
• raises LinAlgError (not ValueError) on a non-finite matrix
• asymmetric input is symmetrized rather than silently upper-triangled
• the deprecated alias warns and matches the new property

Two plotter monkeypatch sites repointed to the new property. Fixed a stale test name asserting asymmetric_curvature_reg_matrix over a symmetric matrix.

Local run: 1152 passed, 1 skipped. Three test_transformer.py pynufft failures are pre-existing on a6b07cd (verified by stashing) and unrelated — a sandbox dependency issue, expected to be green in CI.

Out of scope

The estimator-level defects are tracked separately and deliberately excluded: the covariance describes the unconstrained Warren & Dye solve while use_positive_only_solver: true (NNLS) is the shipped default; the noise map ignores zeroed_ids_to_keep under use_edge_zeroed_pixels: true; and that setting is silently ignored when the positive-only solver is off. Those need a statistical decision, not a code repair.

---

🤖 Generated with Claude Code

https://claude.ai/code/session_0133X4XhMV91SFjzV2mK4Ejh

---

Generated by Claude Code