perf: apply the interferometer sparse operator with rfft2/irfft2 - #540
Merged
Merged
Conversation
`InterferometerSparseOperator` zero-padded a real batch to (2y, 2x) and applied the operator with the complex `fft2`/`ifft2` pair, taking `Re(...)` of an already-real result; `Khat` was likewise the complex transform of a real preload. Because the preload and the batch are both real, the real-transform pair is exact here, not an approximation: the discarded half of the spectrum is the conjugate mirror of the half that is kept, and `irfft2` reconstructs it. - `from_nufft_precision_operator` now caches `Khat = rfft2(preload)`, shape (2y, x + 1) rather than (2y, 2x) complex. - `apply_operator` does `rfft2` -> multiply -> `irfft2(s=(2y, 2x))` -> crop, dropping the `jnp.real` call. - Pinned against the complex `fft2`/`ifft2` algorithm at `rtol=1e-10`, `atol=1e-10 * max|G|` on the 7x7 / K=5 fixture and a seeded 12x12 / K=64 case (max abs error 2.8e-14). A loosened tolerance would hide a bug, so the pin stays at round-off. - Measured 1.27-1.61x faster on every geometry and mesh (autolens_profiling#226, `results/notes/numba_interferometer_verdict.md`). No API change: the four consumers only call `apply_operator`, `Khat` is not serialised and the workspace `.npy` caches store the raw preload. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018hLF3ZAcz5MmaSJBEcLkvF
This was referenced Sep 8, 2026
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
InterferometerSparseOperatorzero-padded a real batch to(2y, 2x)and applied theoperator with the complex
fft2/ifft2pair, takingRe(...)of an already-realresult.
Khatwas likewise the complex transform of a real preload. Because both thepreload and the batch are real, the real-transform pair
rfft2/irfft2computes thesame result — exactly, not approximately: the discarded half of the spectrum is the
conjugate mirror of the half that is kept, and
irfft2reconstructs it.from_nufft_precision_operatorcachesKhat = jnp.fft.rfft2(preload)— shape(2y, x + 1)instead of(2y, 2x).apply_operatordoesrfft2→ multiply →irfft2(s=(2y, 2x))→ crop, dropping thenow-redundant
jnp.real.Closes #538. Follow-up from
autolens_profiling#226 (phase 2 ofnumba-interferometer-revisit); numbers inautolens_profiling/results/notes/numba_interferometer_verdict.mdandresults/breakdown/interferometer/bakeoff_v2026.8.17.1.json.Measured saving
Single thread, i9-10885H,
autoarray2026.8.17.1, same synthetic inputs, batch size 128,scipy stack — so the comparison isolates the transform:
fft2rfft2rfft2on the NumPy stack also beats the library's own JAX-CPU route by 1.10-1.30× atevery one of those cells — i.e. the transform change is worth more than the backend.
Exactness pin
test__interferometer_sparse_operator__apply_operator__rfft2_matches_complex_fft2_referencebuilds the operator from the shared 7×7 / K=5 fixture and a seeded 12×12 / K=64 case,
and compares
apply_operator(np.eye(M))against an inline NumPy complexfft2/ifft2reference at
rtol=1e-10,atol=1e-10·max|G|— max absolute error 2.8e-14. The changeis exact, so a loosened tolerance would be hiding a bug; the pin stays at round-off. It also
asserts the new
Khat.shape == (2y, x + 1).API Changes
None — internal changes only. The
Khatfield's shape and meaning change(
(2y, 2x)complexfft2→(2y, x + 1)rfft2), but it is an internal cache field:the four consumers (
curvature_matrix_diag_from,_off_diag_from,_func_list_from,operated_matrix_slim_from) only callapply_operator, and no public signature moves.Downstream / workspace impact
None. No API change and no numerical result changes (round-off only), so no workspace
migration or new demo is needed:
grepacrossautofit_workspace/,autogalaxy_workspace/,autolens_workspace/,autolens_workspace_test/,euclid_strong_lens_modeling_pipeline/andHowToLens/finds exactly one consumer,
autolens_workspace_test/scripts/misc/jax_assertions/sparse_operators.py, which callsInterferometerSparseOperator.from_nufft_precision_operator(...)with an unchangedsignature and pins the result at
rtol=1e-4— well outside the 2.8e-14 shift.Khatis not serialised: the class has no__getstate__and rides as pytree auxdata (
simulator.py:56 no_flatten). The workspace.npypreload caches store the rawpreload, not
Khat, so no cache is invalidated by the shape change.PyAutoHands/autohands/config/no_run.yamlis not present in this checkout, so nodisabled-script hidden risk could be checked there.
/prm precondition — A100 confirmation
The prompt asks: "GPU: the saving should hold or improve (half the FFT work and half the
intermediate memory); confirm on the A100 before merging rather than assuming it."
No GPU was available in the session that opened this PR, so this is shipped on CPU
parity plus the phase-2 CPU measurements above, and the GPU check is left as a
/prmprecondition for the human to require or waive.
Recipe, if required:
/mnt/ral/jnightin/PyAutoinstall — it lagsmainand carries live jobs) and put itfirst on
PYTHONPATH; check outmainfor the "old" leg and this branch for the "new"leg.
jaxmicrobenchmark: build a random real(560, 560)preload(a 280×280 extent →
(2y, 2x)),InterferometerSparseOperator.from_nufft_precision_operator,then time
apply_operatoron a(M, 128)batch —block_until_ready(), one warm-upcall for the jit, then ~10 timed repeats,
x64on.hpc/sync push-submit gpu <script>(SLURMgpupartition,JAX picks up the GPU automatically), then
hpc/sync jobs/tail gpu/pull.irfft2'sXLA lowering is worse than the complex pair there, and is the one outcome the CPU
numbers cannot rule out.
Test Plan
pytest test_autoarray/inversion test_autoarray/dataset/interferometer test_autoarray/operators -q→ 604 passed (61 s), CPU JAX x64.rtol=1e-10(max abs err 2.8e-14).test_interferometer.pysparse-vs-mapping end-to-end comparisons (1e-4) unchanged.Full API Changes (for automation & release notes)
Changed Behaviour
autoarray.InterferometerSparseOperator.from_nufft_precision_operator(...)— the cachedKhatfield is nowrfft2(preload)of shape(2y, x + 1)rather thanfft2(preload)of shape
(2y, 2x). Signature and return type unchanged.autoarray.InterferometerSparseOperator.apply_operator(Fbatch_flat)— computed withrfft2/irfft2; results identical to the previous complex route to round-off(max abs error 2.8e-14 on the test fixtures). Signature and return shape unchanged.
Removed / Added / Renamed / Changed Signature
Migration
Generated by the PyAutoLabs agent workflow.
🤖 Generated with Claude Code
https://claude.ai/code/session_018hLF3ZAcz5MmaSJBEcLkvF