feat: kernel-regularization linear algebra — cholesky inverses, jitter kwarg, slogdet-gated exact log-det

Two tiers, split by whether likelihood values move:

No likelihood change at defaults (always on):
- ExponentialKernel and GaussianKernel invert their SPD kernel covariance
  via inv_via_cholesky (as MaternKernel) instead of xp.linalg.inv — the
  identical quantity with better accuracy and symmetry; GaussianKernel
  keeps its trace-scaled stabilisation jitter on the formed matrix.
- All four kernel schemes (Matern/MaternAdapt/Gaussian/Exponential) expose
  the covariance diagonal jitter as a constructor kwarg. The default None
  resolves to the historical 1e-8 (byte-identical matrices, asserted by a
  new test) and keeps af.Model prior counts unchanged.

Gated on Settings.log_det_method == "slogdet" (opt-in, default off):
- New AbstractRegularization.log_det_regularization_matrix_term_from hook
  (base returns None). The kernel schemes override it with the analytically
  exact log det H = pixels*log(coeff) - log det C from a single Cholesky of
  their covariance C (MaternAdaptKernel omits the coeff term — its adaptive
  weights live inside C), avoiding the explicit inverse whose round-off
  (amplified by cond(C) ~ 1e9 on clustered traced mesh vertices) puts a
  ~1e-6 absolute noise floor on the evidence. C's conditioning does not
  depend on the regularization coefficient, so the term is also finite at
  any coefficient.
- AbstractInversion.log_det_regularization_matrix_term consults the hook
  only under "slogdet" with every linear object regularized (per-object
  sum = block-diagonal log det), falling back wholesale to the existing
  slogdet path when any scheme lacks a shortcut. The default "cholesky"
  evidence path never calls it — default values are unchanged (asserted).

Extends the PyAutoArray#391 opt-in rather than adding a new toggle; Settings
docstring and packaged general.yaml comment updated. Six new unit tests in
test_kernel_log_det.py; regularization + inversion-abstract suites green.

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

Author: claude

autoarray/config/general.yaml

@@ -9,7 +9,7 @@ inversion:
   nnls_jacobi_preconditioning: true   # If True (default), the curvature matrix passed to jaxnnls.solve_nnls_primal is Jacobi-preconditioned (D Q D y = D q, x = D y). Fixes NaN backward-pass gradients on ill-conditioned Q and roughly halves forward solve time. Set False to restore the raw unpreconditioned solve.
   nnls_target_kappa: 1.0e-11          # Central-path relaxation parameter passed to jaxnnls.solve_nnls_primal. Larger values smooth the relaxed-KKT backward pass and prevent NaN gradients on ill-conditioned Q; smaller values tighten the primal solve. Verified finite gradients across all MGE/rectangular/delaunay pipelines (imaging + interferometer) with scale invariance over 5 orders of magnitude in noise. jaxnnls's own default (1e-3) is too aggressive for the backward pass.
   reconstruction_vmax_factor: 0.5     # Plots of an Inversion's reconstruction use the reconstructed data's bright value multiplied by this factor.
-  log_det_method: cholesky            # How the Bayesian-evidence log-determinant terms are computed. "cholesky" (default) is the historical 2*sum(log(diag(cholesky(M)))); "slogdet" uses logabsdet of slogdet(M), which is identical where M is positive-definite but finite (not NaN) where the Cholesky fails, for gradient-based searches (opt-in, non-default; does not change the default evidence). See PyAutoArray#391.
+  log_det_method: cholesky            # How the Bayesian-evidence log-determinant terms are computed. "cholesky" (default) is the historical 2*sum(log(diag(cholesky(M)))); "slogdet" uses logabsdet of slogdet(M), which is identical where M is positive-definite but finite (not NaN) where the Cholesky fails, for gradient-based searches (opt-in, non-default; does not change the default evidence). Under "slogdet" the kernel regularization schemes (Matern/Gaussian/Exponential) also compute the regularization log-det analytically from a Cholesky of their covariance instead of factorizing the formed inverse. See PyAutoArray#391.
 numba:
   use_numba: true
   cache: true

autoarray/inversion/inversion/abstract.py

@@ -758,6 +758,17 @@ def log_det_regularization_matrix_term(self) -> float:
         used for the source reconstruction, this uses scipy sparse linear algebra to solve the determinant efficiently
         (or ``slogdet`` when ``Settings.log_det_method == "slogdet"`` — see :meth:`_log_det_symmetric_from`).
 
+        Under ``log_det_method == "slogdet"`` (opt-in, default off — PyAutoArray#391), regularization schemes
+        which know a factorization of their own matrix may additionally supply the term directly via
+        :meth:`AbstractRegularization.log_det_regularization_matrix_term_from` — the kernel schemes
+        (``MaternKernel`` etc.) return the analytically exact ``pixels * log(coeff) - log det C`` from a single
+        Cholesky of their covariance ``C``, avoiding the round-off of factorizing the explicitly formed inverse
+        (which reaches ~1e-6 absolute in the evidence at cond(C) ~ 1e9 on clustered traced mesh vertices).
+        Because ``regularization_matrix_reduced`` is the block diagonal of the per-object matrices when every
+        linear object is regularized, the term is the sum of the per-object terms; if any scheme has no
+        shortcut (returns ``None``) the whole computation falls back to ``slogdet`` of the formed matrix.
+        The default ``"cholesky"`` path never consults the shortcut, so default evidence values are unchanged.
+
         Returns
         -------
         float
@@ -766,6 +777,22 @@ def log_det_regularization_matrix_term(self) -> float:
         if not self.has(cls=AbstractRegularization):
             return 0.0
 
+        if (
+            self.settings.log_det_method == "slogdet"
+            and self.all_linear_obj_have_regularization
+        ):
+            term_list = [
+                regularization.log_det_regularization_matrix_term_from(
+                    linear_obj=linear_obj, xp=self._xp
+                )
+                for linear_obj, regularization in zip(
+                    self.linear_obj_list, self.regularization_list
+                )
+            ]
+
+            if all(term is not None for term in term_list):
+                return sum(term_list)
+
         return self._log_det_symmetric_from(self.regularization_matrix_reduced)
 
     @property

autoarray/inversion/regularization/abstract.py

@@ -1,6 +1,6 @@
 from __future__ import annotations
 import numpy as np
-from typing import TYPE_CHECKING
+from typing import Optional, TYPE_CHECKING
 
 if TYPE_CHECKING:
     from autoarray.inversion.linear_obj.linear_obj import LinearObj
@@ -167,3 +167,38 @@ def regularization_matrix_from(self, linear_obj: LinearObj, xp=np) -> np.ndarray
         The regularization matrix.
         """
         raise NotImplementedError
+
+    def log_det_regularization_matrix_term_from(
+        self, linear_obj: LinearObj, xp=np
+    ) -> Optional[float]:
+        """
+        Returns ``log det H`` of this scheme's regularization matrix computed from a
+        factorization the scheme itself knows about, or ``None`` when no such shortcut
+        exists (the default).
+
+        The kernel regularization schemes build ``H = coefficient * C^-1`` from a dense
+        covariance ``C`` and can therefore return the analytically exact
+        ``pixels * log(coefficient) - log det C`` from a single Cholesky of ``C`` —
+        avoiding the explicit inverse, whose round-off (amplified by ``cond(C)``, which
+        reaches ~1e9 on clustered traced mesh vertices) otherwise leaks into the
+        evidence. ``C``'s conditioning does not depend on the regularization
+        coefficient, so this path is also finite where factorizing the formed ``H``
+        fails at extreme coefficients.
+
+        The inversion consumes this ONLY when ``Settings.log_det_method == "slogdet"``
+        (the opt-in gradient-safe log-det, PyAutoArray#391) — the default
+        ``"cholesky"`` evidence path never calls it, so default likelihood values are
+        unchanged. See ``AbstractInversion.log_det_regularization_matrix_term``.
+
+        Parameters
+        ----------
+        linear_obj
+            The linear object (e.g. a ``Mapper``) whose regularization matrix the
+            log-determinant is of.
+
+        Returns
+        -------
+        The log determinant of the regularization matrix, or ``None`` when this scheme
+        has no factorization-aware shortcut.
+        """
+        return None

autoarray/inversion/regularization/exponential_kernel.py

@@ -1,6 +1,6 @@
 from __future__ import annotations
 import numpy as np
-from typing import TYPE_CHECKING
+from typing import Optional, TYPE_CHECKING
 
 if TYPE_CHECKING:
     from autoarray.inversion.linear_obj.linear_obj import LinearObj
@@ -11,6 +11,7 @@
 def exp_cov_matrix_from(
     scale: float,
     pixel_points: np.ndarray,  # shape (N, 2)
+    jitter: float = 1e-8,
     xp=np,
 ) -> np.ndarray:  # shape (N, N)
     """
@@ -53,13 +54,18 @@ def exp_cov_matrix_from(
 
     # add a small jitter on the diagonal
     N = pts.shape[0]
-    cov = cov + xp.eye(N, dtype=cov.dtype) * 1e-8
+    cov = cov + xp.eye(N, dtype=cov.dtype) * jitter
 
     return cov
 
 
 class ExponentialKernel(AbstractRegularization):
-    def __init__(self, coefficient: float = 1.0, scale: float = 1.0):
+    def __init__(
+        self,
+        coefficient: float = 1.0,
+        scale: float = 1.0,
+        jitter: Optional[float] = None,
+    ):
         """
         Regularization which uses an Exponential smoothing kernel to regularize the solution.
 
@@ -88,12 +94,22 @@ def __init__(self, coefficient: float = 1.0, scale: float = 1.0):
             The regularization coefficient which controls the degree of smooth of the inversion reconstruction.
         scale
             The typical scale of the exponential regularization pattern.
+        jitter
+            The small value added to the covariance diagonal for numerical stability.
+            ``None`` (default) uses the historical value 1e-8 — behaviour is identical
+            to not having this parameter (it is a fixed setting, not a free model
+            parameter, hence the ``None`` default).
         """
         self.coefficient = coefficient
         self.scale = scale
+        self.jitter = jitter
 
         super().__init__()
 
+    @property
+    def jitter_value(self) -> float:
+        return 1e-8 if self.jitter is None else self.jitter
+
     def regularization_weights_from(self, linear_obj: LinearObj, xp=np) -> np.ndarray:
         """
         Returns the regularization weights of this regularization scheme.
@@ -128,10 +144,40 @@ def regularization_matrix_from(self, linear_obj: LinearObj, xp=np) -> np.ndarray
         -------
         The regularization matrix.
         """
+        from autoarray.inversion.regularization.matern_kernel import inv_via_cholesky
+
+        covariance_matrix = exp_cov_matrix_from(
+            scale=self.scale,
+            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            jitter=self.jitter_value,
+            xp=xp,
+        )
+
+        # The SPD inverse via Cholesky (as MaternKernel) — identical quantity to
+        # inv(), with better accuracy and symmetry on the SPD covariance.
+        return self.coefficient * inv_via_cholesky(covariance_matrix, xp=xp)
+
+    def log_det_regularization_matrix_term_from(
+        self, linear_obj: LinearObj, xp=np
+    ) -> float:
+        """
+        The analytically exact ``log det H`` from a single Cholesky of the kernel
+        covariance: ``H = coefficient * C^-1``, so
+        ``log det H = pixels * log(coefficient) - 2 * sum(log(diag(cholesky(C))))``.
+
+        Consumed by the inversion only when ``Settings.log_det_method == "slogdet"``
+        (see :meth:`AbstractRegularization.log_det_regularization_matrix_term_from`);
+        the default evidence path factorizes the formed ``H`` and is unchanged.
+        """
         covariance_matrix = exp_cov_matrix_from(
             scale=self.scale,
             pixel_points=linear_obj.source_plane_mesh_grid.array,
+            jitter=self.jitter_value,
             xp=xp,
         )
 
-        return self.coefficient * xp.linalg.inv(covariance_matrix)
+        log_det_covariance = 2.0 * xp.sum(
+            xp.log(xp.diag(xp.linalg.cholesky(covariance_matrix)))
+        )
+
+        return linear_obj.params * np.log(self.coefficient) - log_det_covariance

autoarray/inversion/regularization/gaussian_kernel.py

@@ -1,6 +1,6 @@
 from __future__ import annotations
 import numpy as np
-from typing import TYPE_CHECKING
+from typing import Optional, TYPE_CHECKING
 
 if TYPE_CHECKING:
     from autoarray.inversion.linear_obj.linear_obj import LinearObj
@@ -9,7 +9,10 @@
 
 
 def gauss_cov_matrix_from(
-    scale: float, pixel_points: np.ndarray, xp=np  # shape (N, 2)
+    scale: float,
+    pixel_points: np.ndarray,  # shape (N, 2)
+    jitter: float = 1e-8,
+    xp=np,
 ) -> np.ndarray:
     """
     Construct the source‐pixel Gaussian covariance matrix for regularization.
@@ -43,13 +46,18 @@ def gauss_cov_matrix_from(
 
     # Add tiny jitter on the diagonal
     N = pts.shape[0]
-    cov = cov + xp.eye(N, dtype=cov.dtype) * 1e-8
+    cov = cov + xp.eye(N, dtype=cov.dtype) * jitter
 
     return cov
 
 
 class GaussianKernel(AbstractRegularization):
-    def __init__(self, coefficient: float = 1.0, scale: float = 1.0):
+    def __init__(
+        self,
+        coefficient: float = 1.0,
+        scale: float = 1.0,
+        jitter: Optional[float] = None,
+    ):
         """
         Regularization which uses a Gaussian smoothing kernel to regularize the solution.
 
@@ -76,11 +84,21 @@ def __init__(self, coefficient: float = 1.0, scale: float = 1.0):
             The regularization coefficient which controls the degree of smooth of the inversion reconstruction.
         scale
             The typical scale of the exponential regularization pattern.
+        jitter
+            The small value added to the covariance diagonal for numerical stability.
+            ``None`` (default) uses the historical value 1e-8 — behaviour is identical
+            to not having this parameter (it is a fixed setting, not a free model
+            parameter, hence the ``None`` default).
         """
         self.coefficient = coefficient
         self.scale = scale
+        self.jitter = jitter
         super().__init__()
 
+    @property
+    def jitter_value(self) -> float:
+        return 1e-8 if self.jitter is None else self.jitter
+
     def regularization_weights_from(self, linear_obj: LinearObj, xp=np) -> np.ndarray:
         """
         Returns the regularization weights of this regularization scheme.
@@ -115,25 +133,63 @@ def regularization_matrix_from(self, linear_obj: LinearObj, xp=np) -> np.ndarray
         -------
         The regularization matrix.
         """
+        from autoarray.inversion.regularization.matern_kernel import inv_via_cholesky
+
         covariance_matrix = gauss_cov_matrix_from(
             scale=self.scale,
             pixel_points=linear_obj.source_plane_mesh_grid.array,
+            jitter=self.jitter_value,
             xp=xp,
         )
 
-        regularization_matrix = self.coefficient * xp.linalg.inv(covariance_matrix)
+        # The SPD inverse via Cholesky (as MaternKernel) — identical quantity to
+        # inv(), with better accuracy and symmetry on the SPD covariance.
+        regularization_matrix = self.coefficient * inv_via_cholesky(
+            covariance_matrix, xp=xp
+        )
 
-        # 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.
+        # The inverse can still lose exact symmetry and 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)
+        h_jitter = 1e-8 * xp.abs(diag_mean)
         regularization_matrix = regularization_matrix + xp.eye(
             N, dtype=regularization_matrix.dtype
-        ) * jitter
+        ) * h_jitter
 
         return regularization_matrix
+
+    def log_det_regularization_matrix_term_from(
+        self, linear_obj: LinearObj, xp=np
+    ) -> float:
+        """
+        The analytically exact ``log det H`` from a single Cholesky of the kernel
+        covariance: ``H = coefficient * C^-1``, so
+        ``log det H = pixels * log(coefficient) - 2 * sum(log(diag(cholesky(C))))``.
+
+        This is the log-determinant of the analytic ``coefficient * C^-1`` — it
+        deliberately excludes the trace-scaled stabilisation jitter that
+        :meth:`regularization_matrix_from` adds to the formed matrix (that jitter
+        exists only to guard the factorization of the explicit inverse, which this
+        shortcut avoids entirely).
+
+        Consumed by the inversion only when ``Settings.log_det_method == "slogdet"``
+        (see :meth:`AbstractRegularization.log_det_regularization_matrix_term_from`);
+        the default evidence path factorizes the formed ``H`` and is unchanged.
+        """
+        covariance_matrix = gauss_cov_matrix_from(
+            scale=self.scale,
+            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            jitter=self.jitter_value,
+            xp=xp,
+        )
+
+        log_det_covariance = 2.0 * xp.sum(
+            xp.log(xp.diag(xp.linalg.cholesky(covariance_matrix)))
+        )
+
+        return linear_obj.params * np.log(self.coefficient) - log_det_covariance

autoarray/inversion/regularization/matern_adapt_kernel.py

@@ -1,6 +1,6 @@
 from __future__ import annotations
 import numpy as np
-from typing import TYPE_CHECKING
+from typing import Optional, TYPE_CHECKING
 
 from autoarray.inversion.regularization.matern_kernel import MaternKernel
 
@@ -21,6 +21,7 @@ def __init__(
         inner_coefficient: float = 1.0,
         outer_coefficient: float = 1.0,
         signal_scale: float = 1.0,
+        jitter: Optional[float] = None,
     ):
         """
         Regularization which uses a Matern smoothing kernel to regularize the solution with regularization weights
@@ -63,7 +64,7 @@ def __init__(
             receive significantly higher weights (and faint pixels lower weights), while smaller values produce a
             more uniform weighting. Typical values are of order unity (e.g. 0.5–2.0).
         """
-        super().__init__(coefficient=0.0, scale=scale, nu=nu)
+        super().__init__(coefficient=0.0, scale=scale, nu=nu, jitter=jitter)
         self.inner_coefficient = inner_coefficient
         self.outer_coefficient = outer_coefficient
         self.signal_scale = signal_scale
@@ -109,7 +110,35 @@ def regularization_matrix_from(self, linear_obj: LinearObj, xp=np) -> np.ndarray
             pixel_points=pixel_points,
             nu=self.nu,
             weights=kernel_weights,
+            jitter=self.jitter_value,
             xp=xp,
         )
 
         return inv_via_cholesky(covariance_matrix, xp=xp)
+
+    def log_det_regularization_matrix_term_from(
+        self, linear_obj: LinearObj, xp=np
+    ) -> float:
+        """
+        The analytically exact ``log det H`` from a single Cholesky of the weighted
+        kernel covariance: this scheme's ``H = C_w^-1`` (no coefficient scaling — the
+        adaptive weights are inside ``C_w``), so ``log det H = -log det C_w``.
+
+        Consumed by the inversion only when ``Settings.log_det_method == "slogdet"``
+        (see :meth:`AbstractRegularization.log_det_regularization_matrix_term_from`);
+        the default evidence path factorizes the formed ``H`` and is unchanged.
+        """
+        kernel_weights = 1.0 / self.regularization_weights_from(
+            linear_obj=linear_obj, xp=xp
+        )
+
+        covariance_matrix = matern_cov_matrix_from(
+            scale=self.scale,
+            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            nu=self.nu,
+            weights=kernel_weights,
+            jitter=self.jitter_value,
+            xp=xp,
+        )
+
+        return -2.0 * xp.sum(xp.log(xp.diag(xp.linalg.cholesky(covariance_matrix))))