Kernel regularization: opt-in Cholesky-solve evidence term + relative jitter

[Jammy2211]

Closes leg 2 of PyAutoMind draft/feature/autoarray/regularization_jax_gradient_gaps.md.

Both changes are opt-in and default-off, so default evidence values are unchanged and archived results stay comparable — the constraint the #391 adversarial probe verdict imposed.

Leg 2a — s^T H s via one Cholesky solve

The kernel schemes define H = coefficient * C^-1 and were forming that inverse explicitly. Forming an inverse is slower and less accurate than solving: the error is amplified by cond(C), which reaches ~1e9 on the clustered traced vertices of the kNN mesh families, putting a noise floor on the likelihood the sampler sees.

regularization_term now computes coefficient * s^T C^-1 s from a single Cholesky solve against s. Shape mirrors the log-det shortcut shipped in #391 exactly: a regularization_term_from hook on AbstractRegularization returning None by default, overridden by the four kernel schemes, consulted behind a new Settings.regularization_term_method ("matmul" default, "cho_solve" opt-in).

Measured on a clustered fixture at cond(C) = 3.2e9, graded against an exactly-known reference (choosing s = Cv makes C^-1 s = v, so the true form is s^T v and needs no inverse anywhere):

	relative error
explicit inverse	5.99e-08
Cholesky solve	2.93e-16

Why a separate setting from log_det_method

Deliberate. The two evidence terms can be moved onto their exact factorizations independently, which is what makes an evidence shift attributable to one term rather than both. Gating s^T H s on a setting named for log-determinants would also have been a lie in the name.

Scope limit worth stating

This does not remove the explicit inverse from the inversion as a whole. curvature_reg_matrix is a dense F + H feeding the dense solve for the reconstruction, so H is still formed there. Only the evidence terms avoid it. Doing better needs an iterative-solver design — a different piece of work.

Leg 2b — relative jitter

The covariance jitter was a fixed absolute 1e-8 * I, which only makes sense when diag(C) ~ 1. Measured: true for the three unweighted kernels (K(0) == 1), false for MaternAdaptKernel, whose C_ii = w_i^2 spans the adaptive-weight dynamic range.

Distortion of the faintest pixel, 40-pixel fixture:

inner/outer	faintest C_ii	distortion
1.0 / 1.0	1.0e+00	1.0e-08
0.5 / 4.0	3.9e-03	2.6e-06
1.0 / 20.0	6.3e-06	1.6e-03
0.1 / 100.0	1.0e-08	1.0e+00

inner_coefficient/outer_coefficient are free model parameters, so a sampler can reach that bottom row mid-fit — silently, with no exception and no NaN, just wrong smoothing.

jitter_relative=True applies C_ii *= (1 + jitter). Since C = D^(1/2) R D^(1/2), that is exactly D^(1/2) (R + jitter I) D^(1/2) — the jitter lands on the correlation matrix, so every pixel gets the same relative protection whatever its scale.

A rejected alternative, recorded so it isn't retried: jitter = N * eps * max(diag) ("as small as possible above the round-off floor") fixes the distortion but lets cond(C) reach 3.2e15, at the edge of float64 — reintroducing exactly the noise this branch exists to remove. The correlation-relative rule leaves conditioning unchanged (3.16e9 both conventions, measured).

Two traps found while implementing, both pinned by tests

• MaternAdaptKernel passes coefficient=0.0 to MaternKernel.__init__ (its weights live inside C_w). Inheriting the parent's term would have silently zeroed it — wrong but finite, the worst failure mode. It needs its own override.
• GaussianKernel's formed matrix carries a symmetrisation plus trace-scaled jitter that the analytic shortcut excludes — consistent with how its log-det shortcut already behaves, since both exist only to guard a factorization these paths avoid.

Schemes with no factorization (Constant, Adapt, the split families) return None, and one None falls the whole computation back to the formed matrix, so mixed inversions stay correct.

Testing

29 new tests across three files:

• test_kernel_regularization_term.py (9) — shortcut vs formed matrix for all four schemes, the gate in both directions, the fallback, the MaternAdapt zero-coefficient trap, the ill-conditioning gain.
• test_kernel_jitter_relative.py (7) — default byte-identical to previous behaviour, relative jitter as a pure diagonal rescaling, the adaptive-weight bug, conditioning preserved on clustered vertices, flag threaded to all three covariance call sites per scheme.
• test_kernel_jax_gradients.py (13) — the JAX leg of the gate.

JAX certification:

check	result
FD certification of d(s^T C^-1 s)/d(scale)	4.6e-09
eager vs jit	3.6e-15
implicit vs explicit inverse (JAX)	1.1e-13
numpy vs jax parity	8.4e-12

Full suite: 936 passed, 52 skipped. Three test_transformer.py failures are pre-existing on main and unrelated (missing optional pynufft).

Matérn's JAX path is not covered by the new JAX tests — it needs the modified Bessel from tfp-nightly, an optional-of-an-optional. Gaussian and Exponential exercise the same shared apply_jitter / quadratic_form_via_cholesky code without it. The workspace-level autolens_workspace_test/scripts/imaging/jax_grad/regularization.py still covers the Matérn end-to-end and is unaffected by this branch's defaults.

---

Generated by Claude Code

[Jammy2211]

CI is red on unittest (3.12) and unittest (3.13), but the failure is pre-existing on main and not caused by this branch.

Both jobs fail on a single test:

FAILED test_autoarray/util/test_cholesky_degenerate.py::test__fnnls_cholesky__never_returns_a_non_finite_solution[0.0] - assert 0 > 0

The same test fails identically on main at efaf3041 — this PR's base commit — in run 31219374045 (2026-08-07):

	result
main @ efaf3041	1 failed, 959 passed, 1 skipped
this PR @ dd2cd72	1 failed, 989 passed, 1 skipped

Same test, same assertion, +30 passing tests from this branch and no new failures.

Why it is not this branch. The assertion is the test's own sanity check at test_cholesky_degenerate.py:130:

# Sanity: the degenerate band must actually be exercising the guard,
# otherwise the assertion above is passing vacuously.
assert raised > 0

It requires fnnls_cholesky to raise LinAlgError for at least one of 40 seeds at jitter=0.0. Whether a near-degenerate matrix raises depends on the LAPACK/BLAS build, so this is environment-sensitive by construction. It passes locally (32 passed) and fails on the CI runners.

fnnls_cholesky lives in autoarray/util/fnnls.py, which this branch does not touch. The diff is confined to autoarray/inversion/regularization/*, inversion/abstract.py, settings.py, config/general.yaml, and three new test files — there is no code path from any of those to fnnls.

Not fixing it here. Deciding what that sanity assertion should be (drop it, seed-tune it, or make the guard deterministic) is a judgement about a different subsystem and belongs in its own change rather than being folded into this one. Happy to file it separately.

---

Generated by Claude Code