Guarantee GaussianKernel regularization matrix is PD

Symmetrise and add trace-scaled jitter to the inverted covariance
matrix so that cholesky(regularization_matrix_reduced) in the log
evidence calculation cannot fail on floating-point noise when the
Gaussian scale is large relative to pixel spacing.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

Author: 2 people

autoarray/inversion/regularization/gaussian_kernel.py

@@ -115,4 +115,19 @@ def regularization_matrix_from(self, linear_obj: LinearObj, xp=np) -> np.ndarray
             xp=xp,
         )
 
-        return self.coefficient * xp.linalg.inv(covariance_matrix)
+        regularization_matrix = self.coefficient * xp.linalg.inv(covariance_matrix)
+
+        # inv() loses exact symmetry and can introduce tiny negative eigenvalues
+        # when the covariance matrix is near-singular (e.g. scale >> pixel
+        # spacing). Symmetrise and add a trace-scaled diagonal jitter so the
+        # downstream cholesky in log_det_regularization_matrix_term cannot
+        # fail on floating-point noise.
+        regularization_matrix = 0.5 * (regularization_matrix + regularization_matrix.T)
+        N = regularization_matrix.shape[0]
+        diag_mean = xp.mean(xp.diag(regularization_matrix))
+        jitter = 1e-8 * xp.abs(diag_mean)
+        regularization_matrix = regularization_matrix + xp.eye(
+            N, dtype=regularization_matrix.dtype
+        ) * jitter
+
+        return regularization_matrix