Repository navigation
Interferometer + Delaunay sparse-operator curvature math is wrong (silent ~34% error) #314
Description
Summary
InversionInterferometerSparse computes a curvature matrix that disagrees with the mapping path by ~34% Frobenius norm when paired with a Delaunay-mesh mapper. The regularized matrix then loses positive-definiteness, raising numpy.linalg.LinAlgError: Matrix is not positive definite deep inside Inversion.log_det_curvature_reg_matrix_term's Cholesky call.
PR https://github.com/PyAutoLabs/PyAutoArray/pull/_TBD_ adds a defensive NotImplementedError at the entry point. This issue tracks the actual math fix.
Reproducer
Realistic ALMA setup (Hannah Stacey's profiler): Interferometer.from_fits(transformer_class=TransformerDFT) at 16984 visibilities × (40, 40) real-space mask × Delaunay(pixels=578, zeroed_pixels=30) + ConstantSplit regularization. Then dataset.apply_sparse_operator(use_jax=True) and try fit.figure_of_merit — fails.
Minimal in-Python compare against the mapping path (one-off script, since deleted; reproducible from the public profiler):
=== curvature_matrix comparison ===
mapping shape: (578, 578)
sparse shape: (578, 578)
||sparse - mapping||_F = 20.7711
||mapping||_F = 60.3154
on-diagonal disagreement |diag(sparse - mapping)|.max = 4.33995
off-diagonal disagreement max = 1.80408
=== curvature_reg_matrix_reduced eigenvalues ===
mapping eigvals: min = 0.664671 max = 11.447
sparse eigvals: min = -4.82475 max = 12.6188
=== Cholesky attempts ===
mapping Cholesky: PASS
sparse Cholesky: FAIL (Matrix is not positive definite)
What is/isn't the trigger
- Not
zeroed_pixels: same exact failure (Frobenius diff 20.77, smallest eigvalue -4.82) withzeroed_pixels=0. - Not
use_jax:apply_sparse_operator(use_jax=False)matchesuse_jax=Trueto 5.5e-14 — both produce the same wrong matrix. - Not Rectangular mapping:
rectangular_sparse.pyin autolens_workspace_test works; the existing unit tests attest_abstract.py::test__curvature_matrix__via_sparse_operator__identical_to_mappingpass. - Triggered by Delaunay (Pmax = 3 source pixels per image pixel via barycentric interpolation, weights summing to 1). Top-10 most-disagreeing rows are all on active vertices — not in the zeroed-edge block — and the disagreement is mostly off-diagonal, so source-pixel cross-coupling is structurally wrong.
Root cause hypothesis (unverified)
The implementation in autoarray/inversion/inversion/interferometer/inversion_interferometer_util.py:677–839 (InterferometerSparseOperator.curvature_matrix_via_sparse_operator_from) was written assuming a single source pixel per image pixel (Pmax=1, weight=1, the Rectangular case). For Delaunay each image pixel has three source pixels and weights summing to 1; the COO accumulation + FFT-operator-application + segment_sum cols aggregation likely doesn't compose correctly when multiple source pixels contribute to the same image pixel.
The Imaging sparse path (inversion_imaging_util.curvature_matrix_diag_from) handles Delaunay correctly — see test_abstract.py::test__curvature_matrix_via_sparse_operator__includes_source_interpolation__identical_to_mapping. It uses mapper.sparse_triplets_curvature (flat COO via rows, cols, vals) instead of the (M_masked, Pmax) shape the interferometer path uses. The interferometer path may need rewriting to mirror the imaging path's structure.
Test coverage gap
Before this fix, no test combined Interferometer + Delaunay + sparse:
test_autoarray/inversion/inversion/test_abstract.py:101and:160useaa.Imaging(imaging sparse path) — different code.autolens_workspace_test/scripts/jax_likelihood_functions/interferometer/rectangular_sparse.pyis the only end-to-end interferometer sparse test, Rectangular-only.autolens_workspace_test/scripts/jax_likelihood_functions/imaging/delaunay_mge.py:123has# dataset = dataset.apply_sparse_operator(batch_size=128)commented out — never activated.
Recommended fix scope
Rewrite InterferometerSparseOperator.curvature_matrix_via_sparse_operator_from to mirror the structure of the imaging sparse path. Add a sibling regression test test__curvature_matrix__interferometer_sparse_operator__delaunay__identical_to_mapping, modeled on the existing imaging interpolation test, that verifies sparse-vs-mapping agreement at rtol=1e-4 for a Delaunay interferometer dataset. Once that's in, the guard added by the linked PR can be removed.
Workaround
Use the plain DFT or NUFFT path (i.e. don't call apply_sparse_operator) for Delaunay interferometer fits. The just-shipped JAX-traceable DFT path at autolens_workspace_developer/jax_profiling/jit/{interferometer,datacube}/delaunay.py is the reference setup.
🤖 Filed via Claude Code
Activity
Jammy2211 commentedon May 15, 2026
Defensive guard shipping in #315. The real math fix (the contents of this issue) remains open.
Merge pull request #315 from PyAutoLabs/feature/datacube-sparse-operator
Jammy2211 commentedon May 16, 2026
Library PR Created — 2026-05-16
PR
Status
Workspace changes needed (one call-site update in
autolens_workspace_test/scripts/jax_assertions/sparse_operators.pyplus the stalerectangular_sparse.pyliteral). Next:/ship_workspace.What the fix does
We replaced the
InterferometerSparseOperatorcurvature path with one that mirrorsImagingSparseOperator.curvature_matrix_diag_from(rows, cols, vals, *, S)and uses extent-flat row indices (newMask2D.extent_index_for_masked_pixelproperty). The old path used native-flat indices, so for any real_space_mask withshape_native_masked_pixels < shape_native(i.e. every non-trivial circular mask) most of the JAX.at[rows, bc].add(v)writes silently fell out-of-bounds and were dropped by JAX.That same indexing bug also explains the previously-documented ~0.4% Pmax = 1 "numerical reformulation" gap on
rectangular_sparse.py. With the fix, all three paths (apply_sparse_operator(use_jax=True),TransformerDFT-no-sparse,TransformerNUFFT-no-sparse) agree to ~1e-13.Test coverage
test__apply_sparse_operator__delaunay_mapper__raises_not_implementedis converted totest__curvature_matrix__interferometer_sparse_operator__delaunay__identical_to_mapping, asserting sparse-vs-mapping parity atrtol=1e-4on a Delaunay mapper. This is the missing regression that PR fix(interferometer-sparse): guard against Delaunay mappers (issue #314) #315's guard was a placeholder for.PR #316 supersedes the defensive guard from #315.
Jammy2211 commentedon May 16, 2026
Workspace PR Created — 2026-05-16
PR
Status
Both PRs are now up. Workspace PR is gated on #316 — must wait for the library PR to land on
mainbefore merging the workspace PR.Jammy2211 commentedon May 16, 2026
Smoke Test Results — 2026-05-16
All 11 active smoke tests passed across autolens_workspace_test (99s wall-clock). 1 skipped:
database/scrape/general(NEEDS_FIX 2026-04-27 — pre-existing skip, unrelated to this PR).Jammy2211 commentedon May 16, 2026
Shipped — 2026-05-16
PRs
Summary
We replaced the InterferometerSparseOperator curvature path with a structure mirroring ImagingSparseOperator.curvature_matrix_diag_from(rows, cols, vals, *, S) and routed it through a new Mask2D.extent_index_for_masked_pixel property so the rows land in the operator's actual extent-flat scatter buffer. The new test__curvature_matrix__interferometer_sparse_operator__delaunay__identical_to_mapping confirms sparse-vs-mapping parity at rtol=1e-4 on a Delaunay mapper, and the rectangular_sparse.py end-to-end test shows all three paths (sparse, DFT-no-sparse, NUFFT-no-sparse) now agree to ~1e-13. The defensive NotImplementedError guard from #315 is removed.
Session Notes
Closing.