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

`reconstruction_noise_map_with_covariance` applied np.sqrt elementwise to the
whole inverse of curvature_reg_matrix. The off-diagonals of a covariance
matrix are covariances and are routinely negative, so every one was NaN by
construction -- for any matrix, however well-conditioned -- and each call
emitted "RuntimeWarning: invalid value encountered in sqrt". The property's
docstring promised a matrix accounting for "the covariance of the noise
between pixels"; the entries carrying that covariance were the broken ones.

The 1D reconstruction_noise_map was NOT affected: np.sqrt is elementwise, so
it commutes with taking the diagonal.

Replaced with reconstruction_covariance_matrix, which returns C directly and
forms it from a Cholesky factorization. An A/B against np.linalg.inv refuted
two hypotheses and confirmed one:

  - inv gives negative diagonals on well-formed SPD matrices: 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
    (regularization keeps the matrix PD, cond ~6.8e7 at duplicated columns)
  - inv returns asymmetric output: CONFIRMED
    (5.2e-7 at cond 1e12, vs 2.6e-16 for cho_solve)
  - inv silently succeeds on indefinite matrices: CONFIRMED

The last is the substantive one. cho_factor raises LinAlgError on a matrix
with a negative eigenvalue; inv succeeds and returns a plausible-looking
covariance. At eigenvalue -1e-8 all 300 diagonals came back negative; 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. Both existing call sites already catch LinAlgError
(inversion_plots.py:169 and :395), verified rather than assumed.

reconstruction_noise_map now computes sqrt(diag(C)) directly. It was correct
only incidentally before, via the elementwise sqrt; the invariant is now
stated, and it cannot silently become a variance if the matrix changes.

The old name stays as a deprecated alias. Its values do change -- diagonal
from std-dev to variance, off-diagonals from NaN to covariances -- so the
warning says so explicitly.

Tests: the only prior assertion checked [0, 0], a diagonal element, which is
why this shipped. Added off-diagonal finiteness under -W error::RuntimeWarning,
exact symmetry, the sqrt(diag(C)) invariant, the indefinite-matrix raise, and
the deprecation. Repointed two plotter monkeypatch sites to the new property.

Full suite: 1150 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

@@ -1,4 +1,5 @@
 import copy
+import warnings
 
 import numpy as np
 from typing import Dict, List, Optional, Type, Union
@@ -836,27 +837,78 @@ def log_det_regularization_matrix_term(self) -> float:
         return self._log_det_symmetric_from(self.regularization_matrix_reduced)
 
     @property
-    def reconstruction_noise_map_with_covariance(self) -> np.ndarray:
+    def reconstruction_covariance_matrix(self) -> np.ndarray:
         """
-        Returns the noise-map of the reconstruction as a two dimension matrix which accounts for the covariance
-        of the noise between pixels.
+        Returns the covariance matrix of the reconstruction, ``C = [F + reg_coeff*H]^-1``.
+
+        This is the inverse of the curvature matrix with regularization -- the same matrix used to solve for the
+        reconstruction via the linear inversion. Its diagonal holds the variance of each reconstructed pixel and
+        its off-diagonal entries the covariances between pixels; the off-diagonals are routinely negative, which is
+        a property of the covariance and not an error.
+
+        For the RMS standard deviation of each pixel (the quantity used for scientific analysis) use
+        `reconstruction_noise_map`, which takes the square root of this matrix's diagonal.
+
+        The inverse is formed from a Cholesky factorization rather than `np.linalg.inv`, for two reasons:
+
+        - `cho_factor` raises `LinAlgError` when the matrix is not positive-definite. `np.linalg.inv` raises only
+          on an exactly singular matrix, so an indefinite `curvature_reg_matrix` -- which the inversion does
+          encounter, hence `Settings.no_regularization_add_to_curvature_diag_value` -- previously returned a
+          plausible-looking covariance with no error and no warning. A covariance is only defined for a
+          positive-definite matrix, so failing here is correct and is handled by the callers in
+          `inversion/plot/inversion_plots.py`.
+        - `np.linalg.inv` is LU-based and exploits neither the symmetry nor the positive-definiteness this matrix
+          has. Its output drifts out of symmetry as conditioning worsens (measured at ~5e-7 absolute at
+          `cond ~ 1e12`, against ~3e-16 for the Cholesky solve), while a covariance matrix is symmetric by
+          definition.
+
+        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.
+
+        Returns
+        -------
+        The covariance matrix of the reconstruction, of shape [total_params, total_params].
+        """
+        from scipy.linalg import cho_factor, cho_solve
+
+        matrix = np.asarray(self.curvature_reg_matrix)
+
+        covariance = cho_solve(
+            cho_factor(matrix), np.eye(matrix.shape[0], dtype=matrix.dtype)
+        )
+
+        # cho_solve is accurate but not bitwise symmetric; a covariance matrix is symmetric by definition.
+        return 0.5 * (covariance + covariance.T)
 
-        The diagonal of this matrix is the noise-map of the reconstruction, which can be used for analysing the
-        reconstruction with noise properties that are representative of the fit and therefore should be used
-        for any scientific analysis (e.g. source reconstructions of strong lenses).
+    @property
+    def reconstruction_noise_map_with_covariance(self) -> np.ndarray:
+        """
+        Deprecated alias of `reconstruction_covariance_matrix`.
 
-        This noise-map is defined as the RMS standard deviation of the noise in every pixel of the reconstruction.
-        This definition is identical to the `noise_map` attributes of dataset objects.
+        This property previously returned ``np.sqrt(np.linalg.inv(curvature_reg_matrix))`` -- an elementwise square
+        root of the whole covariance matrix. Because the off-diagonal entries of a covariance matrix are routinely
+        negative, every such entry was `NaN` by construction, for any input matrix however well-conditioned, and
+        each call emitted `RuntimeWarning: invalid value encountered in sqrt`.
 
-        It is computed as the square root of the inverse of the curvature matrix with regularization, which is the
-        same matrix used to solve for the reconstruction via the linear inversion.
+        It now returns the covariance matrix itself, so the values differ: the diagonal holds variances rather
+        than standard deviations, and the off-diagonals hold covariances rather than `NaN`.
 
         Returns
         -------
-        The noise-map of the reconstruction as a two dimension matrix which accounts for the covariance of the noise
-        between pixels.
+        The covariance matrix of the reconstruction (see `reconstruction_covariance_matrix`).
         """
-        return np.sqrt(np.linalg.inv(self.curvature_reg_matrix))
+        warnings.warn(
+            "`reconstruction_noise_map_with_covariance` is deprecated; use "
+            "`reconstruction_covariance_matrix` instead. Note the values have changed: it now returns the "
+            "covariance matrix, so its diagonal holds variances rather than standard deviations (the previous "
+            "elementwise square root made every off-diagonal NaN). For the RMS noise of each pixel use "
+            "`reconstruction_noise_map`.",
+            DeprecationWarning,
+            stacklevel=2,
+        )
+
+        return self.reconstruction_covariance_matrix
 
     @property
     def reconstruction_noise_map(self):
@@ -870,15 +922,15 @@ def reconstruction_noise_map(self):
         The noise-map of the reconstruction is the RMS standard deviation of the noise in every pixel of the
         reconstruction. This definition is identical to the `noise_map` attributes of dataset objects.
 
-        It is computed as the square root of the diagonal of the `reconstruction_noise_map_with_covariance` matrix,
-        which is the same matrix used to solve for the reconstruction via the linear inversion.
+        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.
 
         Returns
         -------
         The noise-map of the reconstruction as a one dimensional ndarray, which does not account for the covariance
         of the noise between pixels.
         """
-        return np.diagonal(self.reconstruction_noise_map_with_covariance)
+        return np.sqrt(np.diag(self.reconstruction_covariance_matrix))
 
     @property
     def reconstruction_noise_map_dict(self) -> Dict[LinearObj, np.ndarray]:

test_autoarray/inversion/inversion/test_abstract.py

@@ -1,3 +1,5 @@
+import warnings
+
 import numpy as np
 import pytest
 
@@ -681,14 +683,97 @@ def test__reconstruction_noise_map__asymmetric_curvature_reg_matrix__correct_dia
 
     inversion = aa.m.MockInversion(curvature_reg_matrix=curvature_reg_matrix)
 
-    assert inversion.reconstruction_noise_map_with_covariance[0, 0] == pytest.approx(
-        np.sqrt(2.5), 1.0e-2
-    )
+    assert inversion.reconstruction_covariance_matrix[0, 0] == pytest.approx(2.5, 1.0e-2)
     assert inversion.reconstruction_noise_map == pytest.approx(
         np.sqrt(np.array([2.5, 1.0, 0.5])), 1.0e-3
     )
 
 
+def test__reconstruction_covariance_matrix__off_diagonals_are_finite_and_negative():
+    """
+    The off-diagonal entries of a covariance matrix are covariances and are routinely negative.
+
+    `reconstruction_covariance_matrix` previously applied `np.sqrt` elementwise to the whole inverse, so every
+    negative off-diagonal became NaN by construction -- for any matrix, however well-conditioned -- while
+    emitting `RuntimeWarning: invalid value encountered in sqrt`. Only the [0, 0] diagonal element was asserted,
+    so nothing caught it.
+    """
+    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)
+
+    with warnings.catch_warnings():
+        warnings.simplefilter("error", RuntimeWarning)
+        covariance = inversion.reconstruction_covariance_matrix
+
+    assert np.all(np.isfinite(covariance))
+
+    # this matrix has anti-correlated pixels, so the off-diagonals are genuinely negative
+    assert covariance[0, 1] < 0.0
+    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."""
+    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
+    curvature_reg_matrix = 0.5 * (curvature_reg_matrix + curvature_reg_matrix.T)
+
+    inversion = aa.m.MockInversion(curvature_reg_matrix=curvature_reg_matrix)
+
+    covariance = inversion.reconstruction_covariance_matrix
+
+    assert covariance == pytest.approx(covariance.T, abs=1.0e-15)
+
+
+def test__reconstruction_noise_map__is_sqrt_of_covariance_diagonal():
+    """
+    The invariant, asserted directly rather than via hand-computed values.
+
+    `reconstruction_noise_map` used to be `np.diagonal(...)` of an already-square-rooted matrix, which was
+    correct only incidentally -- because `np.sqrt` is elementwise. It now takes the square root of the
+    covariance diagonal itself, so the relationship is stated rather than emergent.
+    """
+    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)
+
+    assert inversion.reconstruction_noise_map == pytest.approx(
+        np.sqrt(np.diag(inversion.reconstruction_covariance_matrix)), 1.0e-12
+    )
+
+
+def test__reconstruction_covariance_matrix__raises_on_a_non_positive_definite_matrix():
+    """
+    A covariance is only defined for a positive-definite matrix.
+
+    `np.linalg.inv` raises only on an exactly singular matrix, so an indefinite `curvature_reg_matrix` returned
+    a plausible-looking covariance with no error and no warning. The Cholesky factorization rejects it, and the
+    plotting and CSV callers already catch `LinAlgError`.
+    """
+    # symmetric, non-singular, but indefinite (eigenvalues +1 and -1)
+    curvature_reg_matrix = np.array([[0.0, 1.0], [1.0, 0.0]])
+
+    inversion = aa.m.MockInversion(curvature_reg_matrix=curvature_reg_matrix)
+
+    assert np.isfinite(np.linalg.inv(curvature_reg_matrix)).all()  # inv is silent here
+
+    with pytest.raises(np.linalg.LinAlgError):
+        inversion.reconstruction_covariance_matrix
+
+
+def test__reconstruction_noise_map_with_covariance__is_deprecated_alias():
+    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)
+
+    with pytest.warns(DeprecationWarning, match="reconstruction_covariance_matrix"):
+        covariance = inversion.reconstruction_noise_map_with_covariance
+
+    assert covariance == pytest.approx(inversion.reconstruction_covariance_matrix, 1.0e-12)
+
+
 def test__max_pixel_list_from_and_centre__returns_top_pixels_and_brightest_centre():
 
     source_plane_mesh_grid = aa.Grid2DIrregular(

test_autoarray/inversion/plot/test_inversion_plotters.py

@@ -79,8 +79,8 @@ def test__inversion_subplot_of_mapper__singular_curvature_reg_matrix(
 
     monkeypatch.setattr(
         type(inversion),
-        "reconstruction_noise_map_with_covariance",
-        property(lambda self: np.sqrt(np.linalg.inv(np.zeros((params, params))))),
+        "reconstruction_covariance_matrix",
+        property(lambda self: np.linalg.inv(np.zeros((params, params)))),
     )
 
     with pytest.raises(np.linalg.LinAlgError):
@@ -107,8 +107,8 @@ def test__save_reconstruction_csv__singular_curvature_reg_matrix(
 
     monkeypatch.setattr(
         type(inversion),
-        "reconstruction_noise_map_with_covariance",
-        property(lambda self: np.sqrt(np.linalg.inv(np.zeros((params, params))))),
+        "reconstruction_covariance_matrix",
+        property(lambda self: np.linalg.inv(np.zeros((params, params)))),
     )
 
     with pytest.raises(np.linalg.LinAlgError):