From f1817af0ae23acf72b160d173c0e5b00c57eba03 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Fri, 10 Apr 2026 17:29:13 +0100 Subject: [PATCH] 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 --- .../inversion/regularization/gaussian_kernel.py | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/autoarray/inversion/regularization/gaussian_kernel.py b/autoarray/inversion/regularization/gaussian_kernel.py index b44900c5b..bb979971d 100644 --- a/autoarray/inversion/regularization/gaussian_kernel.py +++ b/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