Activity

1. Jammy2211 commented on Sep 24, 2026

   

   Jammy2211

   on Sep 24, 2026

   CollaboratorAuthor

   More actions

   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.
2. Jammy2211 commented on Sep 24, 2026

   

   Jammy2211

   on Sep 24, 2026

   CollaboratorAuthor

   More actions

   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.
3. Jammy2211 commented on Sep 24, 2026

   

   Jammy2211

   on Sep 24, 2026

   CollaboratorAuthor

   More actions

   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.
4. Jammy2211

   added a commit that references this issue on Sep 24, 2026

   Merge pull request #572 from PyAutoLabs/feature/mge-pdip-nnls-converg…

   ...

   3de624b
5. Jammy2211

   closed this as completedin #572on Sep 24, 2026
6. Jammy2211 commented on Sep 24, 2026

   

   Jammy2211

   on Sep 24, 2026

   CollaboratorAuthor

   More actions

   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.