Skip to content

feat(inversion): read log det(F + λH) off the NNLS Cholesky factor; cache curvature_reg_matrix (#267 lever 3) - #555

Merged
Jammy2211 merged 1 commit into
mainfrom
feature/fixed-light-numba-levers-l3
Sep 16, 2026
Merged

Jammy2211 merged 1 commit into
mainfrom
feature/fixed-light-numba-levers-l3

Conversation

@Jammy2211

Copy link
Copy Markdown
Collaborator

Summary

Lever 3 of PyAutoLabs/autolens_profiling#267 (phase 3 of the fixed-lens-light-numba-cpu epic), stacked on lever 2: log det(F + λH) is read off the Cholesky factor the NNLS solver already built, and curvature_reg_matrix is cached. Resolves the folded-in PyAutoMind draft curvature_reg_matrix_rebuilt_every_access.md.

After levers 1 and 2, one S3 likelihood evaluation on the production numba CPU path (HST Delaunay N=1500, AdaptSplit, InversionImagingSparseNumba, NNLS memo ON, 1 thread) still factorised the same 1500×1500 F + λH twice: once inside fnnls_cholesky on its ~1485-column passive set (the factor lived in a local U_buffer and was never returned), and once more as a fresh full potrf in log_det_curvature_reg_matrix_term (~40 ms of a ~270 ms call). curvature_reg_matrix was a plain property reached twice per evaluation (~4 ms) whose docstring described an in-place add removed in 0766edd4 (2025-06-26) and a cache dropped in the blanket e819fa12 sweep (2025-11-05).

What changed

  • fnnls_cholesky(..., factor: Optional[dict] = None) publishes U_buffer, k_active, passive_set and matrix_shape into the out-dict. No arithmetic change; the solve is byte-identical.
  • reconstruction_positive_only_from(..., factor=None) clears the dict at entry and threads it into both the memo-seeded and the dense-seed/fallback fnnls_cholesky calls, so a warm_start_fallback cannot leave a stale factor.
  • New inversion_util.log_det_from_passive_cholesky_from(matrix, U_buffer, k_active, passive_set): log det(M) = log det(M_PP) + log det(M_AA − M_APᵀ M_PP⁻¹ M_PA); the first term is 2·Σ log diag(U), the second one triangular solve with |active| right-hand sides plus an |active|×|active| Cholesky (|active| ≈ 15 on the production system). Raises LinAlgError rather than returning a non-finite value.
  • AbstractInversion.reconstruction stores the factor on the instance (one per inversion, never a module global). log_det_curvature_reg_matrix_term evaluates reconstruction first, then takes the fast path only when all guards hold: numpy backend, log_det_method == "cholesky", all linear objects regularised, factor present with k_active > 0, shape match, and the factor's index set covers the whole reduced system. Otherwise, and on any LinAlgError from the fast path, the unchanged dense route runs (which carries the test-mode guard). "slogdet" mode is untouched.
  • curvature_reg_matrix → cached_property, docstring rewritten. The two stale docstrings on curvature_matrix in imaging/sparse.py and imaging/mapping.py (which claimed it was not cached and was overwritten in place — both already cached_property, both out-of-place since 0766edd4) are corrected, docstring only. curvature_reg_matrix_reduced is left alone (reached once, ~0.01 ms).
  • JAX branch untouched: the JAX path never runs fnnls, so no factor exists there and the guard is numpy-only by construction. The GPU-side question (does XLA collapse the duplicate adds under jit) is answered separately by the hst-gpu-non-solver-residue HLO census on build: raise caps on scipy, astropy, scikit-image, scikit-learn #267.

API Changes

  • New optional keyword factor on fnnls_cholesky and reconstruction_positive_only_from (default None, no behaviour change when omitted).
  • curvature_reg_matrix returns the same object on repeat access (was a fresh array each time); no caller in PyAutoArray, PyAutoGalaxy or PyAutoLens mutates it in place. Holds one extra (n, n) array alive per inversion (18 MB at n=1500).
  • log_det_curvature_reg_matrix_term now evaluates reconstruction first. The likelihood path already did so (regularization_term runs before the log-dets), so production order is unchanged; a diagnostic asking only for this term now pays for the solve.
  • Values: on the production S3 system the log-det term differs from the dense potrf by +1.8e-12 nats and the fit's log evidence is bit-identical (19705.717582263595 on both routes). log_det_method="slogdet" is untouched by construction.

Measurements

Local (real S3 system, 1 thread, medians of 20; laptop under load so absolute numbers are ~3x RAL): log_det_curvature_reg_matrix_term dense 149.5 ms → fast path 13.1 ms (11.5x); curvature_reg_matrix computations per evaluation 2 → 1; factor on the fiducial instance k_active = 1485 of 1500, |active| = 15, solve_ids_to_keep = arange(1500) (the coverage guard accepts the identity subset, which is the production configuration).

RAL A/B (control = lever 2 0b17c292, feature = lever 3 b4322c3e, gpu partition CPUs-only, 1 thread, phase-2 Leg D cell, route b memo ON):

arm PyAutoArray whole call (ms) log_det_curvature_reg_matrix_term (ms) curvature_reg_matrix (ms / n_calls) cell gates job
control (lever 2) 0b17c292 273.194 40.849 4.536 / 2 4/4 PASS 343356
feature (lever 3) b4322c3e 230.031 6.318 2.281 / 1 4/4 PASS 343356

1.188x on the whole call (−43.163 ms), 6.47x on the site, over 32 counterbalanced ABBA blocks at --n-repeats 64; log_likelihood bit-identical at 18665.014330695605 in both arms, witness W1–W5 PASS with the fast and dense-forced log dets agreeing at Δ exactly 0.0 nats and the figure_of_merit and reconstruction bit/byte-identical. (Job 343355 ran the same A/B at 16 repeats; its feature arm measured 224.330 ms and then tripped the harness's 1.03 instrumentation-overhead cap over 8 blocks, so it is kept only as the bridge row — its control, 268.681 ms, reproduces lever 2's recorded 267.448 ms to +0.46 %.)

A100 row (fixed_light.py Delaunay cell, fp64, job 343357, COMPLETED): identity, as constructed. s3/rows["Log det Cholesky (F+λH reduced)"] 1.125 → 1.139 ms and (H reduced) 1.150 → 1.147 ms; all six log_evidence_terms bit-identical in both systems (S3 log_evidence 29140.295897816348, S0 29140.295897816537), and the eager xp=np family — which does take the fast path — is bit-identical too, log_evidence_eager 29140.295876199474 in both arms, where lever 2's moved −3.4466e-06 nats. Both arms' S0 log-det pins PASSED, mapper_block_log_det_status ASSERTED, pinned_drift empty. The A100's own dense potrf of the same F + λH is 1.139 ms against 40.849 ms on one CPU core and 5.344 ms on this PR's route, so there is nothing to carry to the JAX path — and structurally nothing to carry it from, since the JAX solver never builds a factor. Cumulative levers 1+2+3 against the merge base: 413.301 → 230.031 ms, 1.797x (1.80x).

Testing

  • New test_autoarray/inversion/inversion/test_log_det_shared_factor.py (35 tests): the Schur identity vs 2·Σ log diag(cholesky(M)) at rtol 1e-12 over sizes 64/256/512 and active-set sizes 0/1/5/40 with unordered passive sets; the factor contract (UᵀU reproduces M[P][:,P], solution byte-identical with and without factor, 4 seeds); end-to-end on a real aa.Inversion with seeded signed data so the active set is non-empty, asserting the fast path fired and the value; the arange coverage trap; fallbacks (route c takes the dense route and does not evaluate the reconstruction, slogdet unchanged, forced LinAlgError falls back to the dense value, warm_start_fallback still publishes the retry's factor).
  • pytest test_autoarray/inversion -q — 605 passed; pytest test_autoarray/util -q — 127 passed (one test double for fnnls_cholesky in test_jax_nnls.py gained the new keyword and forwards it; no assertion changed).

Not done here (recorded for follow-up on #267)

  • Option A′ (permute active-last, one potrf, folds away both fancy-index copies): ~15–20 ms more but perturbs the passive factor at 5e-13 and so the knife-edge active-set decisions; candidate lever 4.
  • Edge-zeroed systems (zeroed_pixels > 0) fall back to the dense route; generalising the index mapping is a follow-up.
  • S0 systems with linear funcs: a different identity (log det(M_mm) = log det(M) − log det((M⁻¹)_ff)); out of scope.
  • reconstruction_covariance_matrix factorises F + λH a third time for diagnostics; could reuse the same factor.

Tracked by PyAutoLabs/autolens_profiling#267 (lever 3 of 3). Merge order: #553, then lever 2, then this.

Heart gate

pyauto-heart readiness was RED (score 20) at ship time with two release-chain reasons unrelated to this branch — "install verification FAILED (testpypi; checks F)" and "release validation FAILED (stage integrate)" — acknowledged by the human on 2026-09-16 for PR-open only (levers 2 and 3). Merge stays a human act via /prm.

🤖 Generated with Claude Code

…ache curvature_reg_matrix

log_det_curvature_reg_matrix_term re-factorised the 1500x1500 F + λH that
fnnls_cholesky had just factorised on its passive set (~40 ms of a ~270 ms HST
Delaunay N=1500 numba CPU call). fnnls_cholesky now publishes its factor through
an optional out-dict (no arithmetic change; the solve is byte-identical), and on
the numpy path with the default "cholesky" method the term reads
log det(M) = log det(M_PP) + log det(M_AA - M_AP^T M_PP^-1 M_PA) off that factor,
one triangular solve with |active| right-hand sides plus an |active|x|active|
Cholesky; agrees with the dense potrf to <= 2e-12 nats. Guards: numpy, cholesky
mode, all-regularised, factor present, shape and index coverage match; anything
else, and any LinAlgError on the fast path, takes the unchanged dense route.

curvature_reg_matrix becomes a cached_property (it was reached twice per
evaluation); its docstring described an in-place add removed in 0766edd and a
cache dropped in the blanket e819fa1 sweep. Resolves the folded-in
PyAutoMind draft curvature_reg_matrix_rebuilt_every_access.md.

Lever 3 of PyAutoLabs/autolens_profiling#267 (stacked on lever 2).

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
@Jammy2211

Copy link
Copy Markdown
Collaborator Author

Paired measurement PR: PyAutoLabs/autolens_profiling#272 — RAL jobs 343356 (CPU A/B, 1.19x, control = lever 2; 343355 = 16-repeat bridge row) and 343357 (A100, identity; CPU-only lever).

Base automatically changed from feature/fixed-light-numba-levers-l2 to main September 16, 2026 17:35
@Jammy2211
Jammy2211 merged commit 91240e4 into main Sep 16, 2026
3 checks passed
@Jammy2211
Jammy2211 deleted the feature/fixed-light-numba-levers-l3 branch September 16, 2026 21:50
@Jammy2211 Jammy2211 removed the pending-release PR queued for the next release build label Sep 26, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant