fix(inversion): JAX PDIP positive-only solve fails to converge on SLaM MGE systems

[Jammy2211]

Overview

The JAX positive-only inversion solve (jaxnnls primal-dual interior point, max_iter=50) fails to converge on the real SLaM source_lp[1] MGE model (2 lens bases x 20 Gaussians with sigma_min = pixel_scale/10 + 20 source Gaussians = 60 linear columns) and returns the unconverged iterate silently. At 14/48 parameter vectors within ±0.005 of truth the log-likelihood is wrong (-2e5 down to -1e138 vs 15k-22k from NumPy fnnls_cholesky); with cap 200, 9 converge (72-195 it) and 5 never do. Under Nautilus vmap (n_batch 50) one bad lane also drives every batch to the full cap (NNLS ≈ 28% of a GPU batch). Found in the 2026-09-24 MGE likelihood speed audit; the autolens_profiling imaging/mge cell (20+20 Gaussians, mass fixed) hides it (1/32).

Plan

• Capture the failing (Q, q) linear systems from the SLaM model into a small fixture via a capture script in autolens_profiling.
• Add a PyAutoArray regression test that is red on main: JAX PDIP must report converged and match the fnnls_cholesky objective on every fixture system, single and under vmap.
• Diagnose the mechanism on the fixture: unreachable absolute KKT tolerance (6.7e-11 at cond ~1e11) vs genuine divergence in the KKT Cholesky; effect of Jacobi preconditioning; whether the certified active-set solver certifies these pure-MGE systems.
• Checkpoint with the human on mechanism + fix choice.
• Fix so the solve never silently returns an unconverged iterate; converged systems stay within fp tolerance of today; NumPy path untouched.
• Ship library PR (library-first, Heart RED needs the development override); file the follow-on SLaM timing prompt and re-check downstream MGE JAX likelihood pins.

Detailed implementation plan

Affected Repositories

• PyAutoArray (primary)
• autolens_profiling (capture script + hazard JSON only; parallel-claim, disjoint files)

Branch Survey

Repository	Current Branch	Dirty?
array/PyAutoArray	main (7fa8d27)	clean, unclaimed
lens/autolens_profiling	main	clean (untracked dataset/abell_1201 pre-existing); claimed by certified-solver-phase-c1-lane-rate, hst-gpu-residue-p4, point-source-cpu-p3 → parallel-claim with disjoint file set

Suggested branch: feature/mge-pdip-nnls-convergence
Worktree: ~/Code/PyAutoLabs-wt/mge-pdip-nnls-convergence/

Code facts (main 7fa8d27)

• Entry reconstruction_positive_only_from — autoarray/inversion/inversion/inversion_util.py:291-486; JAX branch :390; config nnls_jacobi_preconditioning (True), nnls_target_kappa (1e-11); settings.nnls_solver_tol / nnls_max_iter (None → 50); Jacobi D = 1/sqrt(diag Q) with no tiny-diagonal guard (:429-438); returns only x, stats gets solver="pdip" only.
• Loop solve_nnls — autoarray/util/jax_nnls.py:33-76: lax.while_loop(pdip_iter < max_iter and converged == 0, pdip_pc_step); stop rule is an absolute KKT inf-norm < min(n*EPSILON, 1e-2), EPSILON = eps*5e3 (jaxnnls pdip.py:32,206) → 6.7e-11 for n=60; returns (x, s, z, converged, pdip_iter) but both callers discard converged (:91, :94); custom_vjp :79-108.
• NumPy oracle: fnnls_cholesky (autoarray/util/fnnls.py:27) via inversion_util.py:516-570.
• Certified active-set solver autoarray/util/jax_active_set.py (solve_certified :265, solve_certified_with_fallback :302); dispatch refuses any inversion containing AbstractLinearObjFuncList (inversion/inversion/abstract.py:584-596) because a joint 60-MGE + Delaunay-1500 system failed to certify in 40 passes; pure-MGE systems never tested.
• Settings autoarray/settings.py (:18 nnls_solver_tol, :19 nnls_max_iter, :25-28 positive_only_solver/certified_*); yaml autoarray/config/general.yaml:5-16.
• Tests: test_autoarray/util/test_jax_nnls.py, test_jax_active_set.py, test_autoarray/inversion/inversion/test_inversion_util.py (:263+, :390+), test_positive_only_dispatch.py; JAX guard requires_jax = pytest.mark.skipif(find_spec("jax") is None); x64 enabled inline. No test compares JAX PDIP with fnnls or checks the cap.
• Related: PyAutoArray#377 (pixelized knife-edge flips, ΔLL ~1e-3) — different symptom, cross-reference only.

Implementation Steps

1. Fixture capture (autolens_profiling) — new scripts/imaging/hazards/mge_nnls_capture.py: build the SLaM source_lp[1] model as in autolens_workspace/scripts/imaging/features/pixelization/slam.py:66-116 on the profiling HST dataset; sample a seeded ±0.005 near-truth box; wrap reconstruction_positive_only_from to record the (curvature_reg_matrix, data_vector) passed (post-regularisation, pre-Jacobi), PDIP converged/pdip_iter at cap 50 and 200, fnnls objective and logL. Write results/hazards/component/mge/nnls_capture_slam_hst_<version>.json + a compressed .npz of 8 systems (5 never-converging, 2 cap-hit-converge-by-200, 1 healthy; ~29 KB raw each). Copy the .npz to test_autoarray/inversion/inversion/files/mge_slam_nnls_systems.npz with a README line naming the generator and versions.
2. Red regression test (PyAutoArray) — test_autoarray/inversion/inversion/test_nnls_mge_convergence.py (JAX-guarded, x64): per system x_np = fnnls_cholesky(Q, q); x, s, z, converged, it = solve_nnls(Q_pc, q_pc) mirroring the Jacobi path; assert converged == 1, it < 50, obj(x*D) <= obj(x_np)*(1+1e-8)+1e-8 with obj = 0.5 xᵀQx − qᵀx; end-to-end reconstruction_positive_only_from(..., xp=jnp, stats=stats) matches the NumPy call within 1e-8 relative in objective and (post-fix) stats["converged"]; vmap variant over the 8 systems. Run red on main first; record counts here.
3. Diagnosis — per system cond(Q), cond(Q_pc), min diag, column collinearity, fnnls active-set size; unroll pdip_pc_step 200 steps logging KKT inf-norm and max|x| → plateau above tolerance vs blow-up; toggles: Jacobi off, relative solver_tol, max_iter 200; solve_certified(Q, q) directly on each system. Post the mechanism note here; checkpoint with the human before step 4.
4. Fix — F1 (always): keep converged/pdip_iter through _solve_nnls_primal_with; reconstruction_positive_only_from writes stats["converged"], stats["iterations"]; non-converged solves never pass through silently. F2 (if plateau): problem-scaled tolerance (max(n·EPSILON, rel·‖q_pc‖∞) or relative KKT residual), verified bit-identical on converging systems. F3 (if blow-up): lax.cond(converged, x, solve_certified(...)[0]) fallback for linear-object-only inversions, gated on step 3 showing certification on pure-MGE systems (extend abstract.py:584-596 only for not has(Mapper); the joint refusal stays); else return the best-KKT iterate tracked in the loop. Never lower the cap, add a ridge to Q, or alter the NumPy path. Docstrings/yaml (try/except KeyError) for any new knob.
5. Verification — python -m pytest test_autoarray/ green; existing PDIP-vs-scipy and certified tests unchanged; random well-conditioned systems bit-identical before/after (asserted in the new test); autolens_workspace_test JAX MGE likelihood pins against the worktree library (pins may only move at previously failing vectors); re-run step 1 capture with the fixed library: 0/48 unconverged, logL within 1 nat of fnnls.
6. Ship + follow-ons — ship_library (Heart RED today: integrate:fail + 3 unrelated autolens cluster/weak script failures → human development override at ship time); file draft/research/autolens_profiling/mge_nnls_fix_slam_timing.md (SLaM 60-column single + vmap16, CPU/RTX/A100, NNLS share of a batch).

Key Files

• autoarray/util/jax_nnls.py — PDIP loop, custom_vjp, discarded converged
• autoarray/inversion/inversion/inversion_util.py — reconstruction_positive_only_from, Jacobi, stats
• autoarray/inversion/inversion/abstract.py — solver dispatch/refusal
• autoarray/util/jax_active_set.py — certified solver (candidate fallback)
• autoarray/settings.py, autoarray/config/general.yaml — knobs
• test_autoarray/inversion/inversion/test_nnls_mge_convergence.py (new), test_autoarray/inversion/inversion/files/mge_slam_nnls_systems.npz (new)
• autolens_profiling/scripts/imaging/hazards/mge_nnls_capture.py (new)

Original Prompt

Click to expand starting prompt

JAX positive-only (PDIP NNLS) solve does not converge on the SLaM source_lp[1] MGE model and returns wrong log-likelihoods

Type: bug
Priority: high
Repos:

• PyAutoArray
• autolens_profiling
  Witness: A PyAutoArray regression test builds the 2-basis SLaM MGE model (2 x 20 lens Gaussians with sigma_min = pixel_scale/10 plus 20 source Gaussians) and, at every one of a fixed set of >= 48 near-truth parameter vectors, the JAX positive-only reconstruction reports converged=True and its log-likelihood agrees with NumPy fnnls within 1 nat on CPU fp64; the test fails red on current main (14/48 unconverged).

Target: PyAutoArray

Witness (2026-09-24 audit, autolens 2aaa1c1a8 / autoarray 681938a, jax 0.10.2, CPU fp64)

The SLaM source_lp[1] imaging model (autolens_workspace scripts/imaging/features/pixelization/slam.py:66-116: lens light = 2 bases x 20 Gaussians, sigma_min = pixel_scale/10; source = 20 Gaussians; free Isothermal + ExternalShear; 60 linear columns, 17 free parameters) was evaluated at 48 parameter vectors within +/-0.005 of the simulated truth on the HST-like simple__no_lens_light-style dataset used by autolens_profiling.

• 14/48 points hit the max_iter=50 cap of solve_nnls_primal (jaxnnls primal-dual interior point, lax.while_loop, autoarray/util/jax_nnls.py:33-78, called from reconstruction_positive_only_from in autoarray/inversion/inversion/inversion_util.py:391-466 after Jacobi preconditioning) and return log-likelihoods from -2e5 down to -1e138. NumPy fnnls on the same curvature/data vector gives sane values (15k-22k).
• Raising the cap to 200: 9 of the 14 converge (72-195 iterations), 5 never converge and return NaN.
• Converged results are not always right either: one GPU point was 1 nat off fnnls.
• Results differ between CPU, GPU, single and vmap evaluation and under 1e-16 perturbations of the images. The Jacobi-preconditioned curvature matrix has condition number ~1.5e11. Suspect: near-collinear narrow Gaussians across the two lens bases.
• The autolens_profiling imaging/mge cell (20+20 Gaussians, mass FIXED) hides this: 1/32 bad points. Any regression test must use the 2-basis SLaM model.
• Under Nautilus use_jax_vmap=True (default, autofit/non_linear/search/nest/nautilus/search.py:198, n_batch=50) one failing lane makes every batch run all 50 iterations: NNLS is ~28% of a GPU vmap batch and ~44% of a single A100 evaluation, so convergence is also the largest remaining GPU speed lever (~15-17% of a batch).
• Lowering the cap is not safe: 15 iterations already moves the objective by 2e-8 at converging points.
• The certified active-set solver (autoarray/util/jax_active_set.py, PyAutoArray#566/feat(inversion): certified active-set positive solver on the JAX path, opt-in, mapper-only dispatch (#566) #567) is refused whenever the inversion contains AbstractLinearObjFuncList (autoarray/inversion/inversion/abstract.py:578-598), so the MGE path cannot opt into it today.

Ask

1. Reproduce with a regression test on the 2-basis SLaM MGE model (compare the JAX positive-only reconstruction / log-likelihood against fnnls at a set of near-truth vectors; assert agreement within a nats tolerance and that converged is true).
2. Isolate the root cause (conditioning of the collinear 2-basis block vs PDIP tolerance/nnls_target_kappa vs cap) and fix convergence: candidates are better preconditioning, a tolerance/adaptive-cap policy, a fallback when the loop hits the cap (never return the unconverged iterate as a likelihood), or extending the certified active-set solver dispatch to linear-object-function inversions.
3. Numerics: converged points must stay within fp tolerance of current results; only the failing points may change (towards fnnls).
4. Report the CPU and GPU timing impact on the SLaM model single and vmap16 (autolens_profiling imaging/mge runtime cell is NOT representative; measure the 60-column model).

Audit scripts and JSON artefacts (nnls_fail.py, nnls_fail_slam_cpu.json, lanes.py) exist in the 2026-09-24 session scratch dir; regenerate rather than rely on them.

🤖 Generated with Claude Code

[Jammy2211]

Step 3 diagnosis

Steps 1-2 done (worktree feature/mge-pdip-nnls-convergence, not pushed): the capture (autolens_profiling d6926af, scripts/imaging/hazards/mge_nnls_capture.py) reproduces the audit exactly. 14/48 vectors hit cap 50 and 5/48 never converge in 200. The fixture test_autoarray/inversion/inversion/files/mge_slam_nnls_systems.npz is 158 KB and holds 8 systems (60 columns each): keys 0-4 never converge, 5-6 converge by 200, 7 is healthy. The red test test_nnls_mge_convergence.py (PyAutoArray a417cb28) gives 15 failed / 6 passed: 7 of 7 fail the single-solve test, 7 of 7 fail the end-to-end objective test, and the vmap test fails (lane iterations [50]*7 + [19]). The fixture check, the healthy system x2 and the 3 well-conditioned controls (bit-identical to upstream jaxnnls, within 1e-8 of scipy) all pass. The existing NNLS/certified/dispatch tests still pass (52).

Mechanism: BLOW-UP from degenerate empty columns. It is not a tolerance plateau and not collinearity. All numbers are CPU fp64 (EPSILON=1.1e-12, so tol = 6.66e-11).

key	cond(Q) / cond(Q_pc)	max col cos (19,39)	empty cols (diag = 1e-3 floor)	fnnls #x>0	KKT at fnnls sol	PDIP KKT it0 → it10	outcome
0-4	9.7e10 / 1.5e11	0.99995-0.99999	40-43(44)	23-26	3-6e-13	6e3 → 1e86…1e212	frozen, then NaN at 24-145 (key 3 decays to 6e14 at 200)
5, 6	9.7e10 / 1.5e11	0.99997	40-43(44)	22	3e-13	6e3 → 1e127, 1e151	decays about 1e-2 per step, converges at 195 / 183
7 healthy	9.7e10 / 1.5e11	0.99996	40 only (q=3.6e-5)	27	5e-13	6e3 → 1.7e-4	converges at 19

• The narrowest source Gaussians (columns 40-44) carry raw curvature of 1e-13 to 1e-5, against a max diag of 2.1e7. Their diagonal is therefore just the no_regularization_add_to_curvature_diag_value floor (1e-3). Their q is 1e-10 down to 1e-258, and the off-diagonals are about 3x smaller than q. After Jacobi scaling each is an almost decoupled coordinate with q_pc ≈ 0, so the optimum sits at the degenerate corner y* = z* = 0.
• On the first PDIP step, the slack for column 40 collapses to 1e-12…1e-51 and its dual z jumps to 1e85…1e217. The step length then collapses, x freezes at about 1e2-1e4, and r2 = s·z dominates the KKT residual. That is divergence (BLOW-UP), not a plateau.
• Neither conditioning nor collinearity separates the failures: the healthy system has the same cond(Q_pc) of 1.47e11 and the same 0.99996 cosine. The tolerance is reachable, because the KKT residual at the exact fnnls solution is 3-6e-13, below 6.7e-11.
• Counterfactuals on the Jacobi system: dropping the floored columns makes all 8 converge in 14-16 iterations. Dropping collinear column 39 still leaves key 1 at 200 and the others at 88-187.
• Hard-deleting the "empty" columns is not exact. On key 4, fnnls puts x = 9.1 on a column with raw curvature 4e-6, and deleting it shifts the objective by 4e-8 relative.

Toggles (Δobj is relative to fnnls):

toggle	keys 0-4	keys 5-6	key 7
Jacobi off (raw Q), cap 50/200	never "converged", but Δobj ≤ 6e-16 (a plateau: right answer, absolute tol unreachable at raw scale)	same	cap 50 unconverged; cap 200 NaN
tol = max(1,‖q_pc‖∞)·n·EPSILON (=6.9e-8)	unchanged: Δobj 4.5…3e132 or NaN	193 / 180 it	16 it
Settings(nnls_solver_tol=1e-6 or 1e-8)	unchanged (garbage / NaN)	garbage at cap 50; with nnls_max_iter=200, Δobj ≤ 2.5e-13	Δobj ≤ 4e-14
solve_certified, raw or Jacobi, pass_budget 16/40	certified only key 4 (Jacobi, 9 passes); the rest are uncertified, Δobj 1.6e-6…2.8e-5	certified (10, 13 passes), Δobj ≤ 4e-16	not certified (budget 16 and 40)

Fix options, ranked by what the evidence supports:

1. F1 (surface converged / pdip_iter, never pass an unconverged iterate silently) is needed in every case and fully supported. It is not sufficient alone, because 14/48 would still need a correct answer.
2. F3 (fallback on non-convergence) is only partly supported as specified. The certified solver fixes 3 of the 7 cap-hit systems and fails to certify the healthy one. Its uncertified iterates are off by up to 2.8e-5 relative, about 27 in the objective. A fallback target that the data does support is the un-preconditioned PDIP solve: every failing system gets Δobj ≤ 6e-16. Its absolute-tol flag is never met, though, and it NaNs on key 7 at 200 iterations, so it would need its own acceptance test, for example the objective or a relative KKT check.
3. F2 (problem-scaled tolerance) is contradicted. The residual diverges by 90+ orders of magnitude rather than stalling, and a scaled or looser tol changes nothing on keys 0-4.

Not in the plan, but the evidence points at it: the trigger is the Jacobi scaling of floor-dominated, signal-free columns. Candidates are a guard in the Jacobi/initialisation step for columns where diag(Q) − floor is negligible (for example, pin D for those columns, or start them at the boundary), or a Jacobi-off retry. Any of these needs care, because a hard column drop is off by 4e-8 on key 4. The fix choice is left to the human checkpoint.

[Jammy2211]

Step 4 implementation

4a probe (8 fixture systems, 3 random controls and a pixelized-like n=400 SPD sentinel with cond about 1e8; cap 50, CPU fp64; Δobj is relative to fnnls):

variant	fixture 0-6	fixture 7 (healthy)	controls + sentinel	worst Δobj
Jacobi (today)	7/7 unconverged, garbage or NaN	19 it	9-22 it	1e132
(i) Q/σ, q/σ with σ = max diag Q, absolute tol	4/7 unconverged	30 it	ok	2e197
(ii) raw Q, tol = n·EPS·max(1,‖q‖∞)	7/7 converged, 14-17 it	16 it	9-17 it	3.9e-11
(ii) with that tol × 1e-2	7/7 converged, 16-19 it	18 it	10-19 it	4.1e-13
(ii) × 1e-3 / × 1e-4	18-21 it	20-21 it	10-20 it	3.6e-14 / 4.2e-15
(iii) σ = ‖q‖∞	converged, 16-17 it	16 it	sentinel 7e-12	3e-11
(iii) σ = ‖Q‖∞	5/7 unconverged	29 it	ok	5e143

With tol × 1, variant (ii) converges everywhere but stops about 1e-11 relative short of fnnls. That is about 4e-5 in the objective, which misses the 1e-6-nat logL target. Chosen: (ii) with factor 1e-2, the fewest iterations that meet ≤1e-12. No max|x| blow-up anywhere (max|x| ≤ 9e2), and the healthy key 7 stops at 18 iterations.

Gradient check changed the design. A relaxed-KKT backward pass on the raw Q gives NaN gradients on fixture keys 2, 3 and 4, which is why Jacobi scaling was introduced in the first place. The implemented "raw" mode therefore does raw PDIP forward plus today's Jacobi-space backward pass. The raw iterate is mapped to scaled coordinates (y = x/D, s/D, z·D), then solve_relaxed_nnls and diff_nnls run on D Q D. With that:

• gradients are finite and non-zero on all 8 keys (the Jacobi-mode gradient is NaN on key 1);
• on healthy key 7 the gradient is within 2.8e-3 of today's Jacobi gradient;
• the custom_vjp is exact for x, and (Q, q, D) get zero cotangents, since y depends only on (Q_pc, q_pc).

Files (PyAutoArray 9d223748; autolens_profiling a284862; not pushed)

• autoarray/util/jax_nnls.py: adds data_scaled_solver_tol (DATA_SCALED_TOL_FACTOR = 1e-2) and solve_nnls_primal_with_status, which returns (x, converged, pdip_iter) with no cotangents for the integer flags. solve_nnls_primal keeps its drop-in x-only signature, because about 20 downstream autolens_profiling / workspace_developer scripts import it. Also adds solve_nnls_primal_raw_forward.
• inversion_util.reconstruction_positive_only_from: new preconditioning="jacobi"|"raw", default "jacobi". The PDIP paths now record stats["converged"], stats["iterations"] and stats["preconditioning"] as traced scalars, the same out-dict channel as the certified path's certified/passes, which nothing in AbstractInversion consumes either. "raw" combined with solver="certified" raises.
• AbstractInversion.positive_only_preconditioning_used: inversions with a Mapper get "jacobi"; mapper-less ones get Settings.nnls_preconditioning_no_mapper, backed by config key nnls_preconditioning_no_mapper: raw (read with the KeyError fallback). The NumPy path is untouched.
• The Jacobi path is bit-identical to main: x, jitted x, and dQ/dq on 11 systems.

Tests. test_nnls_mge_convergence.py and test_positive_only_dispatch.py give 42 passed, including:

• raw solve converges in under 50 iterations with Δobj ≤ 1e-12 on all 8 systems;
• end-to-end stats["converged"] == 1;
• the Jacobi mode reports converged == 0 with 50 iterations on keys 0-6;
• the vmap+jit test: all 8 lanes converge in under 50 iterations;
• gradients finite and non-zero;
• controls: Jacobi bit-identical to upstream jaxnnls, raw within 1e-8 of scipy;
• 5 dispatch/Settings tests.

Full suite: 1681 passed. black is clean on the touched files.

Capture re-run against the worktree library, results/hazards/component/mge/nnls_capture_slam_hst_v2026.8.17.1_postfix.json: library preconditioning is raw, 0/48 unconverged (15-19 iterations), and max |logL_jax − logL_numpy| = 5.6e-7 on all 48 vectors (pre-fix: 14/48 wrong, down to −1e140).

Timing (SLaM 60-column model, CPU fp64, 3 interleaved rounds, minimum per side; this laptop is shared, so noise is a few ms):

	main	fix
single, healthy vector	27.3 ms	25.8 ms
single, formerly capped vector	26.2 ms	26.6 ms
vmap16, per evaluation	31.3 ms	30.3 ms

Timing is neutral on CPU. The GPU vmap gain (no lane stuck at the cap) is not measured here.

Downstream pins (14 mapper-less jax_likelihood scripts in autolens_workspace_test, run with main vs the fix):

• The largest change is 2.5e-11 relative, in imaging/mge_group (1185.53571327 → 1185.53571324).
• imaging/mge and multi_dataset/mge move about 1e-8 absolute; the rest move ≤ 8e-13 or not at all. Every pin moves far inside its rtol 1e-4.
• 4 scripts fail identically on main and on the fix, with the existing "JAX vmap likelihood mismatch" at rtol 1e-4: interferometer/lp, interferometer/mge, multi_dataset/dataset_model, multi_dataset/lp. This predates the change.

Not verified: GPU timing and GPU parity.

[Jammy2211]

Library PR Created

• PyAutoArray: fix(inversion): raw-forward PDIP for mapper-less positive-only JAX solves, surface convergence (#571) #572 (label pending-release) — branch feature/mge-pdip-nnls-convergence @ 07145ea
• autolens_profiling (draft, depends on fix(inversion): raw-forward PDIP for mapper-less positive-only JAX solves, surface convergence (#571) #572): hazards: SLaM MGE NNLS capture for PyAutoArray#571 autolens_profiling#307 — branch feature/mge-pdip-nnls-convergence @ 9b19d1e

Test counts: python -m pytest test_autoarray/ — 1681 passed (worktree, JAX fp64, re-run at ship time). Regression fixture (8 SLaM systems): 15 failed / 6 passed on main, green after the fix. Capture re-run: 0/48 unconverged, max |ΔlogL| 5.6e-7 vs NumPy. 14 mapper-less autolens_workspace_test JAX pins: ≤ 2.5e-11 relative movement. autolens_profiling: ruff check, ruff format --check, build_readme.py --check pass.

Heart RED development override

Human authorisation, verbatim, given live in the session on 2026-09-24 in reply to the override request naming #571:

I authorise you to continue

Heart RED reasons at ship time (pyauto-heart readiness, 2026-09-24), unrelated to this branch:

• release validation FAILED (stage integrate)
• workspace validation not passing (4 failed, cloud#35579888156: autolens notebooks/cluster/modeling.ipynb, autolens notebooks/weak/a2744.ipynb, autolens scripts/cluster/modeling.py, +1 more)
• manifest drift: hub organism blurb (organs present) — 7 mismatch(es) vs PyAutoMind/repos.yaml

Branch gates passed: full suite 1681 passed at 07145ea; regression fixture red→green; capture 0/48 unconverged; workspace_test pins ≤ 2.5e-11; profiling lint clean.

Scope of the override: commit / push / PR-open only. No merge, no release — merging requires a separate /prm with green checks.

[Jammy2211]

Shipped

PRs (both merged 2026-09-24):

• fix(inversion): raw-forward PDIP for mapper-less positive-only JAX solves, surface convergence (#571) #572 (merge 3de624b) — library fix, pending release
• hazards: SLaM MGE NNLS capture for PyAutoArray#571 autolens_profiling#307 (merge 00211d36) — capture script + results

Summary

• Bug: the jaxnnls PDIP positive-only solve (max_iter 50) diverged on signal-free, floor-only MGE source-Gaussian columns after Jacobi scaling, returning garbage log-likelihoods on 14/48 near-truth SLaM source_lp[1] vectors.
• Fix: mapper-less JAX inversions now run PDIP on the raw (unscaled) forward system with a data-scaled tolerance 1e-2·n·EPSILON·max(1,‖q‖∞); the Jacobi-space backward pass is kept; Mapper inversions are bit-identical. converged / iterations are surfaced in the solver stats, and Settings.nnls_preconditioning_no_mapper (default raw) selects the mode.
• Evidence: red regression fixture (8 real systems, 158 KB) went from 15 failed to green; PyAutoArray suite 1681 passed; capture re-run 0/48 unconverged, max |ΔlogL| 5.6e-7; 14 workspace_test mapper-less JAX pins agree to ≤2.5e-11 (4 fail identically before and after — a pre-existing vmap rtol 1e-4 mismatch, not filed).
• Governance: shipped under the human's Heart RED development override ("I authorise you to continue", 2026-09-24, recorded above); merged on the human's "ok merge it" via /prm.

Not done / follow-ups

• GPU timing and parity — PyAutoMind draft/research/autoarray/mge_nnls_fix_pyautoarray_571_slam_60.md (blocked on release).
• Already filed: draft/bug/workspaces/mge_likelihood_breakdown_steps_are_cumulative_an.md.
• The certified active-set solver certifies only 3/7 of these systems.