@@ -1,8 +1,10 @@
 """
-Parity tests for the matrix-free interferometer utilities (phase 1 of PyAutoArray #638).
+Parity tests for the matrix-free interferometer utilities (phases 1-3 of PyAutoArray #638).
 
 The matrix-free sparse mapping matrix, curvature matrix action, data vector and PCG solve are held
-against the dense quantities of `InversionInterferometerSparse` on a small synthetic dataset.
+against the dense quantities of `InversionInterferometerSparse` on a small synthetic dataset, and
+(phase 3) the exact preconditioner log-determinant, the stochastic Lanczos quadrature estimate of
+`ln det(F + H)` and the assembled log evidence against their dense counterparts.
 
 These tests deliberately carry no `importorskip("jax")`: the module is NumPy / scipy only, so the
 `unittest-nojax` CI leg runs them.
@@ -26,10 +28,16 @@
 )
 
 
-def _dataset_sparse_from(seed=1, n_visibilities=40):
+def _dataset_sparse_from(seed=1, n_visibilities=40, uv_scale=1.0e4):
     """
     A circular-masked dataset, so the slim grid is a strict subset of its unmasked extent and the
     scatter / gather between them is exercised, with random visibilities of unit noise.
+
+    `uv_scale` sets the baseline lengths in wavelengths. At the default `1e4` the fringes are far
+    broader than the mask, so `W~` is close to a constant and nowhere near the Powell preconditioner's
+    `κ₀ I`; with many visibilities on long baselines (`500`, `1e6`) the dirty beam is compact, as for
+    real data, and `P⁻¹(F + H)` has eigenvalues near one, the regime the stochastic log-determinant
+    is built for.
     """
     mask = aa.Mask2D.circular(shape_native=(12, 12), pixel_scales=0.5, radius=2.5)
 
@@ -38,7 +46,7 @@ def _dataset_sparse_from(seed=1, n_visibilities=40):
     dataset = aa.Interferometer(
         data=aa.Visibilities(visibilities=rng.normal(size=(n_visibilities, 2))),
         noise_map=aa.VisibilitiesNoiseMap(visibilities=np.ones((n_visibilities, 2))),
-        uv_wavelengths=rng.normal(size=(n_visibilities, 2)) * 1.0e4,
+        uv_wavelengths=rng.normal(size=(n_visibilities, 2)) * uv_scale,
         real_space_mask=mask,
         transformer_class=aa.TransformerDFT,
     )
@@ -307,11 +315,240 @@ def test__powell_preconditioner__equals_kappa0_mtm_plus_coefficient_regularizati
     )
 
 
+def _matern_sparse_system_from(coefficient=1.0):
+    """
+    The phase-3 case: a `MaternSparse`-regularized rectangular mesh on the data-dominated dataset
+    (500 visibilities on long baselines), returning the dense-path inversion and the matrix-free
+    operator, `H` and its factor `R` with `H = RᵀR` (the coefficient folded into `R`).
+    """
+    mask, dataset_sparse = _dataset_sparse_from(n_visibilities=500, uv_scale=1.0e6)
+
+    regularization = aa.reg.MaternSparse(coefficient=coefficient, scale=1.0)
+
+    mapper = _rectangular_mapper_from(mask=mask, regularization=regularization)
+    inversion = _inversion_from(dataset_sparse=dataset_sparse, mapper=mapper)
+
+    operator = mf.MatrixFreeDataOperator.from_inversion(
+        inversion=inversion, mapper=mapper
+    )
+
+    regularization_matrix = regularization.regularization_matrix_sparse_from(mapper)
+    regularization_factor = np.sqrt(
+        regularization.coefficient
+    ) * regularization.regularization_factor_sparse_from(mapper)
+
+    return (
+        dataset_sparse,
+        regularization,
+        mapper,
+        inversion,
+        operator,
+        regularization_matrix,
+        regularization_factor,
+    )
+
+
+@pytest.mark.parametrize("coefficient", [1.0, 0.01])
+def test__log_det_preconditioner_from__equals_dense_slogdet(coefficient):
+    _, _, _, _, operator, regularization_matrix, _ = _matern_sparse_system_from(
+        coefficient=coefficient
+    )
+
+    preconditioner, factor = mf.powell_preconditioner_from(
+        operator=operator,
+        regularization_matrix_sparse=regularization_matrix,
+        coefficient=1.0,
+    )
+
+    sign, log_det = np.linalg.slogdet(preconditioner.toarray())
+
+    assert sign == 1.0
+    assert mf.log_det_preconditioner_from(factor) == pytest.approx(log_det, abs=1.0e-10)
+
+
+@pytest.mark.parametrize("coefficient", [1.0, 0.01])
+def test__log_det_curvature_reg_via_slq__within_half_percent_of_dense_and_deterministic(
+    coefficient,
+):
+    (
+        _,
+        _,
+        _,
+        inversion,
+        operator,
+        regularization_matrix,
+        regularization_factor,
+    ) = _matern_sparse_system_from(coefficient=coefficient)
+
+    sign, log_det_dense = np.linalg.slogdet(np.asarray(inversion.curvature_reg_matrix))
+
+    assert sign == 1.0
+
+    estimate, log_det_probes = mf.log_det_curvature_reg_via_slq_from(
+        operator=operator,
+        regularization_matrix_sparse=regularization_matrix,
+        regularization_factor_sparse=regularization_factor,
+        coefficient=1.0,
+        n_probes=10,
+        n_steps=40,
+    )
+
+    assert log_det_probes.shape == (10,)
+    assert estimate == pytest.approx(np.mean(log_det_probes), 1.0e-14)
+    assert abs(estimate - log_det_dense) <= 0.005 * abs(log_det_dense)
+
+    estimate_repeat, log_det_probes_repeat = mf.log_det_curvature_reg_via_slq_from(
+        operator=operator,
+        regularization_matrix_sparse=regularization_matrix,
+        regularization_factor_sparse=regularization_factor,
+        coefficient=1.0,
+        n_probes=10,
+        n_steps=40,
+    )
+
+    assert estimate_repeat == estimate
+    assert np.array_equal(log_det_probes_repeat, log_det_probes)
+
+
+def test__log_det_curvature_reg_via_slq__error_shrinks_with_probes():
+    """
+    Loose convergence check: over a handful of seeds the mean error at 10 probes is below that at 2,
+    and the Monte-Carlo standard error of the mean shrinks.
+    """
+    (
+        _,
+        _,
+        _,
+        inversion,
+        operator,
+        regularization_matrix,
+        regularization_factor,
+    ) = _matern_sparse_system_from()
+
+    _, log_det_dense = np.linalg.slogdet(np.asarray(inversion.curvature_reg_matrix))
+
+    def errors_and_standard_errors(n_probes):
+        errors, standard_errors = [], []
+
+        for seed in range(5):
+            estimate, log_det_probes = mf.log_det_curvature_reg_via_slq_from(
+                operator=operator,
+                regularization_matrix_sparse=regularization_matrix,
+                regularization_factor_sparse=regularization_factor,
+                coefficient=1.0,
+                n_probes=n_probes,
+                seed=seed,
+            )
+            errors.append(abs(estimate - log_det_dense))
+            standard_errors.append(np.std(log_det_probes) / np.sqrt(n_probes))
+
+        return np.mean(errors), np.mean(standard_errors)
+
+    error_2, standard_error_2 = errors_and_standard_errors(n_probes=2)
+    error_10, standard_error_10 = errors_and_standard_errors(n_probes=10)
+
+    assert error_10 < error_2
+    assert standard_error_10 < standard_error_2
+
+
+def test__log_det_curvature_reg_via_slq__factor_not_matching_regularization_matrix__raises():
+    _, regularization, mapper, _, operator, regularization_matrix, _ = (
+        _matern_sparse_system_from(coefficient=0.01)
+    )
+
+    with pytest.raises(ValueError):
+        mf.log_det_curvature_reg_via_slq_from(
+            operator=operator,
+            regularization_matrix_sparse=regularization_matrix,
+            regularization_factor_sparse=regularization.regularization_factor_sparse_from(
+                mapper
+            ),
+            coefficient=1.0,
+        )
+
+
+@pytest.mark.parametrize("coefficient", [1.0, 0.01])
+def test__log_evidence_matrix_free_from__matches_dense_sparse_path_evidence(
+    coefficient,
+):
+    """
+    The dense reference is `FitInterferometer.log_evidence`'s assembly, `fit_util.log_evidence_from`,
+    of the unconstrained sparse-path inversion's terms. With the dense `ln det(F + H)` substituted the
+    other four terms agree to round-off; with the SLQ estimate the total agrees within the SLQ
+    tolerance (0.5% of `ln det(F + H)`).
+    """
+    from autoarray.fit import fit_util
+
+    (
+        dataset_sparse,
+        regularization,
+        mapper,
+        inversion,
+        operator,
+        regularization_matrix,
+        regularization_factor,
+    ) = _matern_sparse_system_from(coefficient=coefficient)
+
+    noise_normalization = dataset_sparse.sparse_operator.noise_normalization
+
+    log_evidence_dense = fit_util.log_evidence_from(
+        chi_squared=inversion.fast_chi_squared,
+        regularization_term=inversion.regularization_term,
+        log_curvature_regularization_term=inversion.log_det_curvature_reg_matrix_term,
+        log_regularization_term=inversion.log_det_regularization_matrix_term,
+        noise_normalization=noise_normalization,
+    )
+
+    reconstruction, info = mf.reconstruction_via_pcg_from(
+        operator=operator,
+        regularization_matrix_sparse=regularization_matrix,
+        coefficient=1.0,
+        rtol=1.0e-12,
+    )
+
+    assert info == 0
+
+    log_det_regularization = regularization.log_det_regularization_matrix_term_from(
+        mapper
+    )
+
+    assert log_det_regularization == pytest.approx(
+        inversion.log_det_regularization_matrix_term, abs=1.0e-10
+    )
+
+    kwargs = dict(
+        operator=operator,
+        regularization_matrix_sparse=regularization_matrix,
+        regularization_factor_sparse=regularization_factor,
+        coefficient=1.0,
+        reconstruction=reconstruction,
+        log_det_regularization=log_det_regularization,
+        data_term=dataset_sparse.sparse_operator.data_term,
+        noise_normalization=noise_normalization,
+    )
+
+    log_evidence_exact_log_det = mf.log_evidence_matrix_free_from(
+        **kwargs, log_det_curvature_reg=inversion.log_det_curvature_reg_matrix_term
+    )
+
+    assert log_evidence_exact_log_det == pytest.approx(
+        log_evidence_dense, abs=1.0e-10 * abs(log_evidence_dense)
+    )
+
+    log_evidence_slq = mf.log_evidence_matrix_free_from(**kwargs)
+
+    _, log_det_dense = np.linalg.slogdet(np.asarray(inversion.curvature_reg_matrix))
+
+    assert abs(log_evidence_slq - log_evidence_dense) <= 0.5 * 0.005 * abs(
+        log_det_dense
+    )
+
+
 def test__matrix_free_util__imports_and_runs_with_jax_blocked():
     """
     The module must work in the `unittest-nojax` leg. JAX is in `sys.modules` for the rest of a
     normal test session, so a subprocess hides `jax`, `jaxlib` and `nufftax` from the import system,
-    imports the module and runs a PCG solve.
+    imports the module and runs a PCG solve and the phase-3 log-determinant / evidence functions.
     """
     script = textwrap.dedent(
         """
@@ -391,6 +628,34 @@ def _find_spec_blocking_jax(cls, name, path=None, target=None):
         assert info == 0
         assert np.allclose(reconstruction, np.asarray(inversion.reconstruction))
 
+        regularization_matrix = scipy.sparse.csr_matrix(
+            np.asarray(inversion.regularization_matrix)
+        )
+        regularization_factor = scipy.sparse.csr_matrix(
+            np.linalg.cholesky(regularization_matrix.toarray()).T
+        )
+
+        _, factor = mf.powell_preconditioner_from(operator, regularization_matrix, 1.0)
+        assert np.isfinite(mf.log_det_preconditioner_from(factor))
+
+        log_det, _ = mf.log_det_curvature_reg_via_slq_from(
+            operator, regularization_matrix, regularization_factor, 1.0, n_probes=2
+        )
+        assert np.isfinite(log_det)
+
+        log_evidence = mf.log_evidence_matrix_free_from(
+            operator,
+            regularization_matrix,
+            regularization_factor,
+            1.0,
+            reconstruction,
+            np.linalg.slogdet(regularization_matrix.toarray())[1],
+            1.0,
+            1.0,
+            n_probes=2,
+        )
+        assert np.isfinite(log_evidence)
+
         assert not any(
             name.split(".")[0] in ("jax", "jaxlib", "nufftax") for name in sys.modules
         )