perf(numba): HST curvature matrix phase 2 — mapper×mapper block + MGE operated matrix

[Jammy2211]

Overview

Successor to #505. #505 took the HST rectangular numba likelihood evaluation from 1.56 s to 0.60 s by
routing the mapper x linear-func block of the curvature matrix F through the batched FFT Convolver.
The post-#505 split is now mapper x mapper 0.277 s (46%), MGE operated mapping matrix ~0.22 s
(37%), everything else <= 0.05 s. Phase 2 attacks those two blocks, targeting <= ~0.35 s/eval
with the log-likelihood unchanged to pinned tolerance.

Blocked until #505 merges — phase 2 edits the same files (sparse.py,
inversion_imaging_numba_util.py, the two autolens_profiling breakdown harness scripts, the results
note). Branch from main once #505 is in. Registered in PyAutoMind/planned.md as blocked.

Plan

• Step 0 — instrument first (checkpoint gates everything after). Split the MGE operated mapping matrix harness row into three timed sub-rows and record the geometry constants into the JSON, so the
  lever is chosen from a measurement rather than an estimate.
• Step 1 — mapper x mapper, bit-identical restructuring. Write the oracle test the kernel does not
  have today, then hoist the re-gathered (pix, weight) pairs into locals. Estimated 10-20%.
• Step 2 — mapper x mapper, two-stage reformulation. Per-data-pixel dense source-space accumulator
  followed by contiguous AXPYs, replacing ~2e8 irregular read-modify-writes into a 4.9 MB matrix with
  ~5e7 L1 scatters plus vectorisable dense adds. Estimated 1.7-2.5x on the block (0.277 -> 0.11-0.16 s),
  at rtol=1e-6.
• Step 3 — MGE operated matrix (gated by the step-0 split). A small bit-identical OverSampler
  divisor cache in PyAutoArray, plus a shared-geometry fast path in PyAutoGalaxy for
  LightProfileLinearObjFuncList — the 60 MGE Gaussians share centre/ell_comps and differ only in
  sigma, so the transform and eccentric-radius grid are recomputed 120x per evaluation. Estimated
  MGE 0.22 -> ~0.10-0.13 s.
• Step 4 — measure and ship. Paired same-session before/after (alternate B/A/B/A; host variance is
  20-30%), three pins at explicit rtol=1e-6 plus a curvature_matrix array comparison, a Nautilus
  pool run, note appended in autolens_profiling, then ship_library PyAutoArray ->
  ship_library PyAutoGalaxy -> ship_workspace autolens_profiling.

Two findings from the planning read-through shaped this: the MGE half of the cost lives in
PyAutoGalaxy, not PyAutoArray (the Convolver batching from #505 is already applied there); and two
of the four mapper x mapper candidate levers are dead on arrival — upper-triangle symmetry is already
exploited (the sparse preload stores ip1 >= ip0 and the kernel folds A + A^T) and unique-mappings
compression is already what the kernel iterates over. The live lever is the two-stage reformulation.

Detailed implementation plan

Affected Repositories

• PyAutoArray (primary)
• PyAutoGalaxy
• autolens_profiling

Branch Survey

Repository	Current Branch	Dirty?
./PyAutoArray	main	clean
./PyAutoGalaxy	main	clean
./autolens_profiling	main	clean

worktree_check_conflict reports PyAutoArray and autolens_profiling are both claimed by the in-flight
task numba-hst-curvature-matrix-speedup (#505), and autolens_profiling additionally by
nuts-warm-start-driver-and-a100-probe. Hence: queued in planned.md, not started.

Suggested branch: feature/numba-hst-curvature-matrix-phase2
Worktree: ~/Code/PyAutoLabs-wt/numba-hst-curvature-matrix-phase2/

Implementation Steps

Step 0 — instrument (checkpoint gates everything after)

Files: autolens_profiling/scripts/imaging/likelihood_breakdown/{pixelization_numba,delaunay_numba}.py
(keep STEP_ACCESSORS in sync as #505 did).

1. Split the MGE operated mapping matrix row into 3 timed sub-rows — 60x image_2d_from(grid),
   60x image_2d_from(blurring_grid), 1x psf.convolved_mapping_matrix_via_real_space_np_from(...) —
   with the existing row becoming the residual (same scheme as the F rows; the override is a
   cached_property on the func object so sub-timings do not prime it; the four rows must sum).
2. Record geometry constants once into the JSON: sparse_operator.lengths.sum()/mean()/max(),
   unique_mappings.pix_lengths.mean()/max(), data_to_pix_unique.shape, mapper.params,
   n_image_pixels, PSF shape.
3. Do not split mapper x mapper further (single mapper -> one kernel call; fold loops are
   n_source^2, noise).
4. Run hst/euclid rectangular + hst Delaunay-1250 (OMP_NUM_THREADS=1 AUTOARRAY_NUMBA_OPERATED_MEMO=0, n_repeats 10); post the split to this issue.
   Checkpoint: confirms (a) profile-eval vs convolution share of MGE, (b) the mapper x mapper op count.

Step 1 — mapper x mapper, bit-identical restructuring (own commit)

Kernel: curvature_matrix_via_sparse_operator_from (inversion_imaging_numba_util.py ~:496).

1. Write the oracle first — it has no direct unit test today. New
   test__curvature_matrix_via_sparse_operator_from__matches_dense_psf_precision_operator in
   test_autoarray/inversion/inversion/imaging/test_inversion_imaging_util.py: dense M^T W M vs the
   kernel, parametrised over KERNELS_ODD (asymmetric, non-square). Plus a symmetry / halved-diagonal
   test. Keep the current kernel as curvature_matrix_via_sparse_operator_reference_from.
2. Hoist data_1's <= u1 (pix, weight) pairs into locals once per stored pair (today re-gathered u0x
   from wide-stride 2-D arrays) and take curvature_matrix[pix_0] as a row view. Estimated 10-20%.
   Verify np.array_equal on inversion.curvature_matrix before/after (euclid + hst).

Step 2 — mapper x mapper, two-stage reformulation (own commit, rtol 1e-6)

1. For each data_0: accumulate the dense a[p] = sum_{data_1} W(data_0,data_1) * sum_b w1_b delta(pix1_b, p) into an L1-resident n_source vector, then A[pix_0,:] += w0 * a for each of
   data_0's u0 mappings; keep the A + A^T fold and the halved-diagonal contract.
2. Valid while n_source <~ P * u1 (~3400 at HST; 784/1250 well inside). Keep the reference loop behind
   a runtime branch on pix_pixels; measure the crossover, and if it is dataset-dependent expose it
   as an explicit setting rather than a silent heuristic.
3. Tests: ..._two_stage_matches_reference_kernel parametrised over KERNELS_ODD x pix_pixels on both
   sides of the threshold; inversion-level test__curvature_matrix_mapper_diag__matches_reference_kernel
   in test_curvature_matrix_func_list_blocks.py; the multi-mapper twin
   curvature_matrix_off_diags_via_sparse_operator_from gets its own oracle if touched. Every rtol test
   gets a recorded control run (a deliberately wrong variant must fail).
4. prange: not planned. Zero gain under Nautilus one-process-per-core, and it needs per-thread
   n_source^2 buffers plus a lazily cached prefix-offset array on SparseLinAlgImagingNumba (pickled to
   workers — no new ctor args). Only if step 2 misses target, default off, both arms measured.
5. Rejected: pruning small psf_precision_operator entries (moves the likelihood; fails the pin).

Step 3 — MGE operated matrix (gated by the step-0 split)

1. PyAutoArray, small, bit-identical: OverSampler.binned_array_2d_from
   (autoarray/operators/over_sampling/over_sampler.py ~:202) recomputes the model-independent
   np.bincount(segment_ids) divisor on every one of ~120 calls/eval. Cache it as an immutable
   cached_property (today the code mutates counts in place — the cached array must not be mutated).
   Wider blast radius (every numpy light-profile eval) -> np.array_equal test across repeated calls
   plus the harness Blurred image row as witness.
2. PyAutoGalaxy: shared-geometry fast path in LightProfileLinearObjFuncList
   (autogalaxy/profiles/light/linear/abstract.py, mapping_matrix ~:291,
   operated_mapping_matrix_override ~:320): when all profiles are the same class with identical
   centre/ell_comps, compute the transform and eccentric radii once per grid and evaluate the 60
   Gaussians vectorised over sigma; binning stays a loop initially (binned_array_2d_from is 1-D).
   Fall back to the current loop otherwise. Tests: ..._shared_geometry_matches_per_profile_loop
   (assert bit-identical first; relax to rtol 1e-6 only with a recorded reason) and
   ..._mixed_geometry_falls_back. Estimated MGE 0.22 -> ~0.10-0.13 s.
3. Not pursued: reusing ConvolverState.fft_kernel on the numpy path (scipy already avoids the 3-D FFT;
   single-digit %; a 10-line microbench only), and a narrower memo keyed on centre+ell_comps (weaker
   than the hoist).

Step 4 — measure + ship

1. Paired same-session before/after (host variance 20-30%; alternate B/A/B/A), three paired cells + RTU
   once for currency (unpaired, GPU-only by decision). Discard the first post-rebuild run (numba
   cache: true; documented cold-cache psf_weighted_data_from hazard).
2. Pins at explicit rtol=1e-6 (the harness default is 1e-4) plus a curvature_matrix array comparison.
3. Pool run: cpu_fast_modeling.py with PYAUTO_TEST_MODE=1, PYAUTO_SMALL_DATASETS unset, cores =
   host, plus a linear MGE bulge (the smoke profile runs no pool); compare the parallel speed-up ratio.
4. Note section appended to autolens_profiling/results/notes/numba_curvature_matrix_f_split.md.
5. ship_library PyAutoArray -> ship_library PyAutoGalaxy (if touched) -> ship_workspace
   autolens_profiling. Merge is human.

Key Files

• PyAutoArray: autoarray/inversion/inversion/imaging_numba/inversion_imaging_numba_util.py,
  autoarray/inversion/inversion/imaging_numba/sparse.py,
  autoarray/operators/over_sampling/over_sampler.py,
  test_autoarray/inversion/inversion/imaging/test_inversion_imaging_util.py,
  test_autoarray/inversion/inversion/test_curvature_matrix_func_list_blocks.py
• PyAutoGalaxy: autogalaxy/profiles/light/linear/abstract.py,
  autogalaxy/profiles/light/standard/gaussian.py, autogalaxy/profiles/geometry_profiles.py,
  new tests under test_autogalaxy/profiles/light/linear/
• autolens_profiling: scripts/imaging/likelihood_breakdown/pixelization_numba.py,
  scripts/imaging/likelihood_breakdown/delaunay_numba.py, _profile_cli.py,
  results/notes/numba_curvature_matrix_f_split.md

Verification

pytest test_autoarray (baseline 1296) and pytest test_autogalaxy green; the oracle test passes on the
unmodified kernel before any change; control runs recorded for every rtol test; bit-identical
curvature_matrix for step 1 and the OverSampler change; three pins at rtol 1e-6 (hst rectangular
27661.910133664103, hst-rtu 27180.704715696862, hst-delaunay 29090.52721044813); smoke
cpu_fast_modeling.py; paired before/after table + pool ratio in the note and on this issue.

Marked undetermined (step 0 measures these, never guessed)

The MGE internal split; the actual psf_precision_lengths.sum() / pix_lengths.mean(); whether LLVM
already hoists the scalar product; the two-stage crossover n_source; the numba threading layer on WSL
(only if prange ever runs).

Original Prompt

Click to expand starting prompt

Numba CPU likelihood at HST resolution, phase 2: the mapper×mapper block and the MGE operated matrix

Type: feature
Epic: none (successor to numba_cpu_hst_curvature_matrix_speedup, PyAutoArray#505)
Target: autoarray
Repos:

• @PyAutoArray
• @PyAutoGalaxy
• @autolens_profiling
  Themes:
• numba-cpu
• pixelization
• profiling
  Difficulty: large
  Autonomy: supervised
  Priority: high
  Status: draft
  Filed: 2026-08-28

Original request

yeah do that phase 2

(in reply to: "What's left for a phase 2, if you want one: mapper×mapper (0.28 s) and the MGE
operated matrix (0.22 s) now lead at HST.")

Context

PyAutoArray#505 (branch feature/numba-hst-curvature-matrix-speedup, shipping 2026-08-28) routed
the mapper × linear-func block of the curvature matrix F through the batched FFT Convolver
(Convolver.reversed_kernel), taking F 1.20 → 0.36 s and the HST rectangular numba evaluation
1.56 → 0.60 s (OMP_NUM_THREADS=1, AUTOARRAY_NUMBA_OPERATED_MEMO=0, paired same-session
measurements; note in autolens_profiling/results/notes/numba_curvature_matrix_f_split.md).

Post-#505 HST rectangular split (0.60 s/eval):

step	s	share
F: mapper×mapper [numba sparse-op / w-tilde contraction]	0.277	46%
MGE operated mapping matrix (60 funcs, per-eval convolution of varying profiles)	~0.22	~37%
F: mapper×linear-func [FFT conv + scatter]	0.054	9%
everything else	≤ 0.02 each

Delaunay-1250: mapper×mapper 0.123 s of 0.76 s; F is no longer dominant there — the inversion build
(~0.5 s) is, which is out of scope here.

Plan findings (2026-08-28)

Two read-through findings changed the task's shape at planning time. First, the MGE half of the
cost lives in PyAutoGalaxy, not PyAutoArray: the #505 Convolver batching is already applied in
LightProfileLinearObjFuncList.operated_mapping_matrix_override (the 60 Gaussians are stacked into
one convolution) and scipy already skips the length-60 axis, so the bulk of the ~0.22 s is the 120
per-profile image evaluations, all of which recompute an identical transform and eccentric-radius
grid because the MGE basis shares centre/ell_comps and varies only in sigma — hence
@PyAutoGalaxy added to Repos above. Second, two of the four mapper×mapper candidate levers are
dead on arrival: upper-triangle symmetry is already exploited (the sparse preload stores
ip1 >= ip0 and the kernel folds A + Aᵀ) and unique-mappings compression is already what the
kernel iterates over. The live lever is a two-stage reformulation — a per-data-pixel dense
source-space accumulator followed by contiguous AXPYs.

Goal

Take the HST rectangular numba evaluation from ~0.60 s to ≤ ~0.35 s with the log-likelihood
unchanged to pinned tolerance (pins: hst rectangular 27661.910133664103, hst-rtu
27180.704715696862, hst-delaunay 29090.52721044813; rtol 1e-6 where summation order changes,
bit-identical otherwise).

1. Decompose first, as in perf: numba CPU curvature matrix F at HST resolution (phase 1: split + FFT mapper×func block) #505. Instrument the mapper×mapper kernel
   (inversion_imaging_numba_util.py, sparse-op diag + off-diag kernels) and the MGE operated
   mapping matrix step in the autolens_profiling breakdown harness; record the split at hst +
   euclid rectangular and hst Delaunay-1250. Checkpoint: the split picks the lever.
2. mapper×mapper candidates, cheapest first: unique-mappings compression (is the inner loop over
   data pixels × PSF footprint × source pixels where a per-unique-mapping formulation is smaller);
   upper-triangle-only + mirror; whether the preloaded PSF-precision products are fully exploited
   (nothing mapper-independent recomputed per evaluation); prange measured both at
   OMP_NUM_THREADS=1 and under the Nautilus pool.
3. MGE operated matrix candidates: determine FFT- vs scatter-bound; batch the varying profiles
   through the same Convolver path as perf: numba CPU curvature matrix F at HST resolution (phase 1: split + FFT mapper×func block) #505 if the per-func convolution is the cost; check whether
   any profile subset is invariant across evaluations for a given model (memo, already
   AUTOARRAY_NUMBA_OPERATED_MEMO).
4. Ship behind parity tests; test_autoarray green; smoke
   autolens_workspace/scripts/imaging/features/pixelization/cpu_fast_modeling.py; paired
   before/after on the four breakdown cells + one Nautilus pool run; note in autolens_profiling.

Out of scope: the JAX path; RTU / kernel-CDF meshes (GPU-only by decision); the Delaunay inversion
build; the NNLS solve.

[Jammy2211]

Step 0 checkpoint — instrumented split + measured geometry

Harness commit 3f69b34 (autolens_profiling, feature/numba-hst-curvature-matrix-phase2).
All three cells run with OMP_NUM_THREADS=1 AUTOARRAY_NUMBA_OPERATED_MEMO=0, n_repeats 10,
against the task worktree's PyAutoArray (PYTHONPATH-pinned, autoarray.__file__ asserted).
First post-rebuild run discarded (numba cache: true).

Per-step split (s/eval)

Step	hst rect (bilinear, 784)	euclid rect (bilinear, 784)	hst delaunay (Hilbert, 1250)
Blurred image (FFT convolve)	0.010794	0.005689	0.011387
Inversion build (trace+mesh+mapper)	0.026765	0.012052	0.028456
PSF-weighted data [numba]	0.009156	0.001931	0.009175
MGE: image_2d_from(grid) x60	0.116026	0.030895	0.116086
MGE: image_2d_from(blurring_grid) x60	0.045998	0.024407	0.046184
MGE: batched PSF convolution	0.054078	0.014001	0.064275
MGE residual (dict assembly)	-0.007544	0.001098	-0.003476
(MGE row total, un-split)	0.208558	0.070401	0.223069
Mapper sparse triplets	0.004870	0.001153	0.019793
Data vector D [numba]	0.004511	0.000750	0.005662
F: mapper x mapper [sparse-op]	0.255055	0.049054	0.130403
F: mapper x linear-func [FFT conv + scatter]	0.049128	0.014918	0.054698
F: linear-func x linear-func [BLAS]	0.018941	0.003284	0.019782
F residual (mirror + diag-add)	-0.006710	0.000521	-0.013915
Regularization matrix H	0.001925	0.002463	0.072216
Reconstruction solve (BLAS)	0.009757	0.006513	0.044785
log det (F+H) + log det H	0.016614	0.013514	0.058052
Log evidence (figure of merit)	0.015803	0.012274	0.059074
TOTAL (step-by-step)	0.633439	0.197839	0.734423
Direct log_likelihood_function	0.585031	0.175982	0.663736

The four MGE rows sum to the raw un-split MGE step, as the four F rows do (both
residuals are differences of averaged timings and can come out slightly negative;
recorded as measured, never clipped). mge_split_reproduces_step: true on all three
cells — the three timed pieces reconstruct the step's dict bit-identically, so
the attribution is not an inference.

Geometry constants (model-independent, recorded once into the JSON)

	hst rect	euclid rect	hst delaunay
PSF shape (native)	21x21	21x21	21x21
image pixels (masked)	15361	3841	15361
sparse_operator stored pairs = lengths.sum()	11,080,573	2,318,869	11,080,573
lengths.mean() / .max()	721.34 / 841	603.71 / 841	721.34 / 841
data_to_pix_unique.shape	(15361, 4)	(3841, 4)	(15361, 3)
pix_lengths.mean() / .max() / .sum()	4.0 / 4 / 61444	4.0 / 4 / 15364	1.5523 / 3 / 23845
mapper.params (n_source)	784	784	1250
mapper x mapper inner ops	1.7729e8	3.7102e7	3.2692e7
two-stage stage-1 scatters	4.4322e7	9.2755e6	1.8187e7
two-stage stage-2 AXPY	4.8172e7	1.2045e7	2.9806e7
two-stage accumulator clear	1.2043e7	3.0113e6	1.9201e7
two-stage total ops	1.0454e8	2.4332e7	6.7195e7

mapper x mapper inner ops is the exact count of
curvature_matrix_via_sparse_operator_from's innermost accumulation,
sum over stored (data_0, data_1) pairs of u0(data_0) * u1(data_1) — not an estimate.

Verdict

(a) MGE — profile evaluation, not convolution. At HST the MGE row is 0.2086 s, of
which 0.1620 s (78 %) is the 120 numpy profile evaluations (0.1160 s on the
over-sampled data grid + 0.0460 s on the blurring grid) and 0.0541 s (26 %) is the single
batched real-space convolution; the residual is -0.0075 s. The plan's inference from #505
("~0.17 of the 0.22 s is the 120 per-profile image evaluations") is confirmed within
5 %. Step 3's shared-geometry fast path in PyAutoGalaxy is therefore aimed at the right
0.16 s, and the ConvolverState.fft_kernel idea correctly stayed out of scope. The 2.5x
asymmetry between the grid and blurring-grid evaluations is the 4x4 radial over-sampling
of the data grid, not a per-call overhead.

(b) mapper x mapper op count — confirmed, but the two-stage margin is dataset-dependent.
HST rectangular measures 1.773e8 innermost accumulations against the plan's ~2e8: right
to within 12 %. At 0.2551 s that is 1.44 ns/accumulation, ~4.6 cycles — the signature of an
irregular read-modify-write into a 4.9 MB matrix, exactly the cost the two-stage
reformulation targets. Its op count on the same geometry is 1.045e8, a 1.70x reduction
before any locality or vectorisation gain, so the plan's 1.7-2.5x estimate for HST
rectangular stands.

The Delaunay cell is the caveat, and it is a new finding from this measurement: at
1250 vertices with pix_lengths.mean() = 1.55 (barycentric, vs 4.0 for bilinear) the
reference kernel does only 3.27e7 ops while the two-stage form does 6.72e7 — 2.05x
more. The crossover is not a function of n_source alone: it depends on n_source,
pix_lengths.mean() (u0) and lengths.mean() (P) together, because stage 2 costs
u0 x n_source per data pixel regardless of how sparse the mapping is. Step 2 will
measure whether vectorised dense AXPYs still win at that op ratio, and will report the
measured crossover; the branch will be an explicit, documented threshold, never a silent
heuristic.

Proceeding to step 1 (oracle + bit-identical hoist).

[Jammy2211]

Phase 2 complete — steps 3 and 4, and the result

Continuing from the step-0 checkpoint and the step-1/2
commits. The goal is met: HST rectangular 0.6214 -> 0.3334 s per evaluation
(step total), 0.6184 -> 0.3013 s directly timed — 1.86x / 2.05x, against the
target of ~0.35 s.

Step 3 — the MGE operated matrix

3a, PyAutoArray d8bc3bac — the OverSampler divisor.
binned_array_2d_from recomputed np.bincount(segment_ids) and the
sub_is_uniform np.isclose sweep on every call — 120 calls per evaluation for
a 60-Gaussian MGE, and once more for every ordinary numpy light profile. Both
depend only on the constructor's sub_size; both are now cached_property. On
the HST over sampler (17980 sub-pixels into 15361 pixels) the divisor cost
38.8 us per call and the uniformity check 61.3 us, against 81.5 us for the whole
cached call. Bit-identical. The cached divisor is writeable = False and built
with np.maximum(counts, 1), because the old zero-guard mutated in place and
the array is now shared by every caller. OverSampler is pickled to Nautilus
workers, so a pickle round-trip is asserted.

3b, PyAutoGalaxy 6031aa0a — one eccentric-radius grid per MGE basis. This
is where the MGE time actually was. Each of the 60 Gaussians re-ran the same
reference-frame transform (arctan2 + sin + cos per coordinate) and formed
the same eccentric radii from it; cProfile on the HST cell put the transform
alone at 951 us of the 1.3 ms each profile took.
LightProfileLinearObjFuncList._image_slim_list_from now groups profiles by
(class, centre, ell_comps) and computes each group's transform and radii once.
Bit-identical. Two decisions worth recording:

• Grouping, not an all-or-nothing check. The workspace's canonical MGE
  recipe stacks two sets of 30 Gaussians which share a centre but carry their own
  ell_comps, and GalaxiesToInversion puts both in one func list. An
  all-or-nothing test would have fallen back to 60 independent evaluations on
  exactly the model users are told to write.

• Sersic is deliberately excluded, even in a same-geometry Sersic basis: it
  already avoids the polar transform via _eccentric_radii_grid_from_cartesian,
  a branch taken only when the grid handed to it is not pre-transformed — so
  hoisting a transformed grid into it would be both slower and a different
  floating-point expression. The membership test is the identity of the
  image_2d_from function object against Gaussian.image_2d_from, which also
  correctly excludes GaussianMultipole.

  Follow-up worth filing: giving Gaussian that same Cartesian shortcut would
  remove the polar transform for every Gaussian evaluation, not just an MGE
  basis. It moves values, so it is a separate task.

Binning stays a per-profile call (binned_array_2d_from takes a 1-D array), and
so does the exp — measured, that arithmetic is transcendental-bound rather than
overhead-bound, so the plan's "vectorise over sigma into one block" would buy
nothing and would have to be sliced column-wise to bin.

Step 4 — measured

Three arms, paired B/A/B/A in one session (base = PyAutoArray main @
1b89404b + PyAutoGalaxy main @ 0fbe863d; after steps 1-2; final), three
rounds per cell in rotating arm order, first discarded for the numba recompile,
mean of the last two. OMP_NUM_THREADS=1, AUTOARRAY_NUMBA_OPERATED_MEMO=0,
n_repeats=10.

row	base	step 2	final	step 3	whole phase
hst rect, step total	0.6214	0.4324	0.3334	1.30x	1.86x
hst rect, direct eval	0.6184	0.4207	0.3013	1.40x	2.05x
hst rect, F mapper x mapper	0.2581	0.0819	0.0825	0.99x	3.13x
hst rect, MGE total	0.1970	0.1937	0.0844	2.29x	2.33x
hst rect, MGE 60x image_2d_from(grid)	0.1042	0.0995	0.0233	4.28x	4.48x
hst rect, MGE 60x image_2d_from(blurring)	0.0436	0.0420	0.0115	3.67x	3.81x
hst rect, MGE batched PSF convolution	0.0512	0.0504	0.0523	0.96x	0.98x
hst rect, Blurred image (FFT)	0.0101	0.0097	0.0099	0.98x	1.02x
euclid rect, step total	0.2226	0.1870	0.1417	1.32x	1.57x
euclid rect, F mapper x mapper	0.0531	0.0188	0.0191	0.98x	2.78x
euclid rect, MGE total	0.0802	0.0797	0.0316	2.52x	2.54x
hst delaunay-1250, step total	0.6331	0.6076	0.4884	1.24x	1.30x
hst delaunay, F mapper x mapper	0.1088	0.0697	0.0676	1.03x	1.61x
hst delaunay, MGE total	0.1881	0.1932	0.0839	2.30x	2.24x

The last two hst-rect rows are the control: nothing outside the two terms step 0
identified moves. The Delaunay cell gains least because its evaluation is now
spread across terms this issue did not touch — after phase 2 its largest steps
are the ConstantSplit regularization matrix (0.060 s), the batched MGE
convolution (0.057 s) and the mapper x linear-func block (0.049 s).

hst rectangular RTU was re-run once for currency, unpaired: it is a 7.1 s cell
dominated by its 6.73 s mapper-sparse-triplets step and is GPU-only by the
2026-08-28 decision.

Correctness

Pinned log likelihoods at the explicit rtol=1e-6, PASSED on every arm and
unchanged to every recorded digit: hst bilinear 27661.91013366411 (pinned
27661.910133665442), hst RTU 27180.70471569685, hst Delaunay
29090.527210448134. euclid, unpinned, 6213.306873885871 on all three arms.

Array-level, same instance on each arm:

comparison	MGE operated matrix	inversion.curvature_matrix
step 3 alone (step2 -> final), hst	bit-identical	bit-identical
step 3 alone, euclid	bit-identical	bit-identical
whole phase (base -> final), hst	bit-identical	max rel 1.24e-14, allclose(rtol=1e-12) True
whole phase, euclid	bit-identical	max rel 2.63e-15, allclose(rtol=1e-12) True

Step 3 is bit-identical by construction in both libraries; the only
floating-point movement in the whole phase is step 2's reassociation, twelve
orders inside the pins.

Tests: pytest test_autoarray 1326 passed (1296 at the phase-1 baseline),
pytest test_autogalaxy 1149 passed (1144 before). Control runs recorded for
every rtol test and for every branch of the new step-3b membership test
(dropping the transform, elliptical instead of eccentric radii, transforming the
un-over-sampled grid, removing the image_2d_from-identity guard — each fails
the new tests).

Pool run

cpu_fast_modeling.py's first fit, 8 cores, PYAUTO_TEST_MODE=1,
PYAUTO_SMALL_DATASETS unset (2828 masked pixels), 60-Gaussian linear MGE bulge
added to the lens. Per evaluation:

	pool of 8 (median of 3)	serial
base	0.2000	0.4099
final	0.1697	0.3634

Parallel speed-up ratio 2.05x -> 2.14x — flat to slightly up, which is the
pass condition; it is the number that would fall if any lever had introduced
hidden threads. End-to-end the pool improves 15 % and the single process 11 % on
this small workspace dataset, where F and the MGE are a much smaller share of the
evaluation than in the HST harness cells. Read it as a regression check, not a
speed-up measurement.

One harness bug found and fixed

The MGE: image_2d_from(blurring_grid) x60 row hand-rolled the per-profile loop
that operated_mapping_matrix_override used to contain. After step 3b that loop
no longer exists in the library, so the harness would have kept timing the old
code on the new library and mis-attributed the entire win to the residual row —
the first final run showed exactly that (a -0.030 s residual). It now calls the
same path the override calls, dispatching on the attribute so both libraries stay
measurable, and records which branch it took in the result JSON as
mge_blurring_stack_path.

PRs (do not merge out of order)

• perf(numba): the mapper x mapper block and the MGE operated matrix (#507) #508 — steps 1, 2, 3a
• perf(linear): share one eccentric-radius grid across an MGE basis (PyAutoArray#507) PyAutoGalaxy#590 — step 3b
• profiling: phase-2 instrumentation and the measured result (#507) autolens_profiling#190 — steps 0 and 4 (harness, note, artifacts)

Merge order: #508 -> #590 -> #190. PyAutoGalaxy's CI clones PyAutoArray from
source and checks out the same-named branch when it exists, so #590 is presently
testing against the PyAutoArray branch and falls back to main once #508 merges
and that branch is deleted. No PyAutoGalaxy test depends on the new OverSampler
behaviour, so #590 is green either way.

Heart is RED at ship time for a pre-existing, unrelated reason
("release validation FAILED (stage integrate)"); the human authorisation to
ship over it is quoted verbatim in each PR body.

[Jammy2211]

Shipped

Phase 2 is merged. The issue's goal — HST rectangular numba evaluation ~0.60 s → ≤ ~0.35 s with the log likelihood unchanged to pinned tolerance — is met.

Merged PRs (library-first gate honoured, in this order):

PR	merge commit
#508 — steps 1, 2, 3a	f4b8e758 (head d8bc3bac)
PyAutoLabs/PyAutoGalaxy#590 — step 3b	90317122 (head 6031aa0a)
PyAutoLabs/autolens_profiling#190 — steps 0 and 4	ed1e265a (head 7cc4bbd7)

Result (paired B/A/B/A in one session, OMP_NUM_THREADS=1, AUTOARRAY_NUMBA_OPERATED_MEMO=0, n_repeats=10):

cell	base (1b89404b)	final	whole phase
hst rectangular, step total	0.6214 s	0.3334 s	1.86x
hst rectangular, direct eval	0.6184 s	0.3013 s	2.05x
euclid rectangular	0.2226 s	0.1417 s	1.57x
hst Delaunay-1250	0.6331 s	0.4884 s	1.30x
hst: F mapper × mapper	0.2581 s	0.0825 s	3.13x
hst: MGE total	0.1970 s	0.0844 s	2.33x

The two levers, both chosen by step 0's instrumentation rather than the prompt's candidate list: a two-stage source-space accumulator for the mapper × mapper block (w0 does not depend on data_1, so the sum factorises into a per-data-pixel dense accumulator plus contiguous AXPYs), and one shared eccentric-radius grid per MGE group in PyAutoGalaxy (78% of the MGE row was profile evaluation, not the PSF convolution — the transform and radii were recomputed for all 60 Gaussians that share centre/ell_comps). Plus a cached OverSampler binning divisor and uniformity check, called 120× per evaluation.

Verification: test_autoarray 1326 passed, test_autogalaxy 1149 passed. All three pinned log likelihoods unchanged to every digit at rtol=1e-6 (hst bilinear 27661.91013366411, hst RTU 27180.70471569685, hst Delaunay 29090.527210448134). curvature_matrix bit-identical for steps 1 and 3, allclose(rtol=1e-12) across the phase. Nautilus pool ratio 2.05x → 2.14x on 8 cores — no oversubscription.

CI green on every run and every matrix leg before each merge: PyAutoArray Tests [pull_request] (unittest (3.12), unittest (3.13), unittest-nojax); PyAutoGalaxy Tests [pull_request] (same three legs) + Docs [pull_request] (docs-build); autolens_profiling lint [pull_request].

Heart RED at merge — human override. Merged over the pre-existing Heart RED (pyauto-heart readiness --json 2026-08-28T21:34:42Z, score 45; red_reasons: "release validation FAILED (stage integrate)"). Human authorisation 2026-08-28, verbatim: "prm and then kick off phase 2, I authorize the heart RED thing" — scoped to this task only; release stays human. The RED is pre-existing on main and unrelated: the failing workspace scripts sit on the JAX likelihood path, while this diff touches only the numba imaging-inversion path plus additive-only changes.

Record: PyAutoMind/complete/2026/08/numba-hst-curvature-matrix-phase2.md.