fix: raise LinAlgError on non-finite input, symmetrize before Cholesky

Addresses an adversarial review of the previous commit.

1. MAJOR regression, now fixed. scipy's cho_factor/cho_solve default to
   check_finite=True and raise ValueError -- not LinAlgError -- on a NaN/inf
   matrix. Both callers (inversion_plots.py:169 and :397) catch only
   LinAlgError, and the CSV writer's docstring explicitly promises that a
   failure there may not abort the enclosing model-fit. So a NaN-contaminated
   curvature matrix would have killed a fit that previously wrote a nan column
   and continued. Verified before and after: the old code returned [nan, nan];
   the new code raised builtins.ValueError past the guards.

   Fixed by checking finiteness explicitly and raising LinAlgError with a
   message naming the cause, rather than by broadening the callers' except
   clauses -- the property's contract is now "raises LinAlgError for any input
   that has no covariance", which holds downstream too. check_finite=False is
   then passed to scipy, so the explicit check costs nothing net.

2. cho_factor reads only the upper triangle, so an asymmetric input was
   silently inverted as though its lower triangle matched. Verified: for
   [[2.0, 0.5], [0.1, 2.0]] it gave diag 0.5333 against the true inverse's
   0.5063. The output symmetrization does not address this -- it symmetrizes
   the result, not the input -- so the input is now symmetrized too.
   curvature_reg_matrix is F + H and symmetric by construction, so this is
   defensive, but the silent wrong answer is not an acceptable failure mode.

3. Softened the "no value change" claim for reconstruction_noise_map. The two
   forms are algebraically identical but only numerically equivalent, since the
   covariance is now formed by Cholesky rather than LU: ~7e-15 relative at
   cond 1e3 rising to ~4e-5 at cond 1e13. Neither is the more correct result.

4. Tests. The symmetry test was tautological -- 0.5 * (C + C.T) is bitwise
   symmetric for any C, so it could only fail if the line were deleted. It now
   also asserts accuracy against an exactly-constructed ground truth
   (A = Q diag(w) Q.T), with the tolerance set from the measured error (~3e-9
   at cond 1e9) and a comment recording that np.linalg.inv is marginally more
   accurate on that matrix, not less. Added tests for the non-finite raise and
   the asymmetric-input symmetrization. Fixed a stale test name that said
   "asymmetric_curvature_reg_matrix" over a symmetric matrix.

Downstream check for the value-changing deprecated alias: no consumer of
reconstruction_noise_map_with_covariance exists in PyAutoGalaxy (3ca31bf),
PyAutoLens (87e5827) or autolens_workspace -- zero hits for "with_covariance"
in any of them. reconstruction_noise_map IS used, in autolens_workspace
source_science.py scripts computing signal_to_noise_map = reconstruction /
reconstruction_noise_map, and its values are unchanged.

Full suite: 1152 passed, 1 skipped. The 3 test_transformer.py pynufft
failures are pre-existing on a6b07cd and unrelated.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0133X4XhMV91SFjzV2mK4Ejh

Author: claude

autoarray/inversion/inversion/abstract.py

@@ -862,6 +862,14 @@ def reconstruction_covariance_matrix(self) -> np.ndarray:
           `cond ~ 1e12`, against ~3e-16 for the Cholesky solve), while a covariance matrix is symmetric by
           definition.
 
+        Every failure mode raises `LinAlgError`, including a non-finite `curvature_reg_matrix`. That case is
+        checked explicitly because scipy would otherwise raise `ValueError`, which the callers -- here and
+        downstream -- do not catch; the previous implementation returned a silently all-`NaN` matrix instead.
+
+        The matrix is symmetrized on input. `cho_factor` reads only the upper triangle, so an asymmetric input
+        would otherwise be inverted as though its lower triangle matched, silently and with no diagnostic.
+        `curvature_reg_matrix` is `F + H` and symmetric by construction, so this is defensive only.
+
         This property is NumPy-only: the input is coerced with `np.asarray`, so a JAX `curvature_reg_matrix`
         forces a device-to-host transfer. It is a post-fit diagnostic, not part of the likelihood, so it is not on
         the JIT path.
@@ -874,8 +882,21 @@ def reconstruction_covariance_matrix(self) -> np.ndarray:
 
         matrix = np.asarray(self.curvature_reg_matrix)
 
+        if not np.isfinite(matrix).all():
+            raise np.linalg.LinAlgError(
+                "The curvature_reg_matrix contains non-finite entries (NaN or inf), so the reconstruction "
+                "covariance is undefined. Raised as LinAlgError so the plotting and CSV callers, which guard "
+                "on LinAlgError, degrade gracefully rather than aborting the model-fit."
+            )
+
+        # cho_factor reads only the upper triangle; symmetrize so an asymmetric input cannot be silently
+        # inverted as though its lower triangle matched its upper.
+        matrix = 0.5 * (matrix + matrix.T)
+
         covariance = cho_solve(
-            cho_factor(matrix), np.eye(matrix.shape[0], dtype=matrix.dtype)
+            cho_factor(matrix, check_finite=False),
+            np.eye(matrix.shape[0], dtype=matrix.dtype),
+            check_finite=False,
         )
 
         # cho_solve is accurate but not bitwise symmetric; a covariance matrix is symmetric by definition.
@@ -925,6 +946,12 @@ def reconstruction_noise_map(self):
         It is computed as the square root of the diagonal of `reconstruction_covariance_matrix`, which is the
         inverse of the same matrix used to solve for the reconstruction via the linear inversion.
 
+        This previously took the diagonal of an elementwise-square-rooted matrix. The two are algebraically
+        identical -- `np.sqrt` is elementwise, so it commutes with taking the diagonal -- but only numerically
+        equivalent, since the covariance is now formed by Cholesky rather than LU. The difference is
+        conditioning-limited roundoff, measured at ~7e-15 relative at `cond ~ 1e3` rising to ~4e-5 at
+        `cond ~ 1e13`; neither result is the more correct one.
+
         Returns
         -------
         The noise-map of the reconstruction as a one dimensional ndarray, which does not account for the covariance

test_autoarray/inversion/inversion/test_abstract.py

@@ -678,7 +678,7 @@ def test__log_det_method__slogdet_is_finite_where_cholesky_fails_on_non_positive
     assert result == pytest.approx(np.linalg.slogdet(matrix)[1], 1.0e-8)
 
 
-def test__reconstruction_noise_map__asymmetric_curvature_reg_matrix__correct_diagonal_noise_values():
+def test__reconstruction_noise_map__correct_diagonal_noise_values():
     curvature_reg_matrix = np.array([[1.0, 1.0, 1.0], [1.0, 2.0, 1.0], [1.0, 1.0, 3.0]])
 
     inversion = aa.m.MockInversion(curvature_reg_matrix=curvature_reg_matrix)
@@ -713,20 +713,61 @@ def test__reconstruction_covariance_matrix__off_diagonals_are_finite_and_negativ
     assert covariance == pytest.approx(np.linalg.inv(curvature_reg_matrix), 1.0e-8)
 
 
-def test__reconstruction_covariance_matrix__is_exactly_symmetric():
-    """A covariance matrix is symmetric by definition; `np.linalg.inv` drifts out of symmetry."""
+def test__reconstruction_covariance_matrix__is_accurate_and_symmetric_when_ill_conditioned():
+    """
+    Ground truth is exact by construction: for `A = Q diag(w) Q.T` the inverse is `Q diag(1/w) Q.T`.
+
+    The symmetry half only guards the symmetrization line -- `0.5 * (C + C.T)` is bitwise symmetric for any C --
+    so the accuracy assertion against the constructed truth is what tests the factorization itself.
+    """
     rng = np.random.default_rng(1234)
     q, _ = np.linalg.qr(rng.standard_normal((25, 25)))
-    curvature_reg_matrix = (q * np.logspace(0, 9, 25)) @ q.T
+    eigenvalues = np.logspace(0, 9, 25)
+
+    curvature_reg_matrix = (q * eigenvalues) @ q.T
     curvature_reg_matrix = 0.5 * (curvature_reg_matrix + curvature_reg_matrix.T)
 
+    covariance_true = (q * (1.0 / eigenvalues)) @ q.T
+
     inversion = aa.m.MockInversion(curvature_reg_matrix=curvature_reg_matrix)
 
     covariance = inversion.reconstruction_covariance_matrix
 
+    # cond ~ 1e9, so the achievable accuracy is eps * cond ~ 2e-7; the measured error is ~3e-9. This is not a
+    # claim that Cholesky beats LU here -- it does not, `np.linalg.inv` measures ~7e-10 on this matrix.
+    assert covariance == pytest.approx(covariance_true, abs=1.0e-7)
     assert covariance == pytest.approx(covariance.T, abs=1.0e-15)
 
 
+def test__reconstruction_covariance_matrix__asymmetric_input_is_symmetrized_not_silently_upper_triangle():
+    """
+    `cho_factor` reads only the upper triangle, so an asymmetric input would be inverted as though its lower
+    triangle matched its upper -- silently, and differing from the true inverse.
+    """
+    curvature_reg_matrix = np.array([[2.0, 0.5], [0.1, 2.0]])
+    symmetrized = 0.5 * (curvature_reg_matrix + curvature_reg_matrix.T)
+
+    inversion = aa.m.MockInversion(curvature_reg_matrix=curvature_reg_matrix)
+
+    assert inversion.reconstruction_covariance_matrix == pytest.approx(
+        np.linalg.inv(symmetrized), 1.0e-8
+    )
+
+
+def test__reconstruction_covariance_matrix__non_finite_matrix_raises_lin_alg_error():
+    """
+    scipy raises `ValueError` on a non-finite matrix, which the plotting and CSV callers do not catch -- they
+    guard on `LinAlgError`. The CSV writer explicitly promises not to abort the enclosing model-fit, so the
+    non-finite case is converted rather than allowed to escape.
+    """
+    curvature_reg_matrix = np.array([[1.0, np.nan], [np.nan, 2.0]])
+
+    inversion = aa.m.MockInversion(curvature_reg_matrix=curvature_reg_matrix)
+
+    with pytest.raises(np.linalg.LinAlgError, match="non-finite"):
+        inversion.reconstruction_covariance_matrix
+
+
 def test__reconstruction_noise_map__is_sqrt_of_covariance_diagonal():
     """
     The invariant, asserted directly rather than via hand-computed values.