diff --git a/complete/2026/08/numerical-inversion-failures.md b/complete/2026/08/numerical-inversion-failures.md new file mode 100644 index 00000000..81dd4131 --- /dev/null +++ b/complete/2026/08/numerical-inversion-failures.md @@ -0,0 +1,238 @@ +# numerical-inversion-failures — release-run non-positive-definite inversion failures + +**Date:** 2026-08-22 +**Issue:** [PyAutoArray#467](https://github.com/PyAutoLabs/PyAutoArray/issues/467) (closed) +**PRs:** none — **no code was changed in any repo** +**Outcome:** investigated to a definitive verdict; no defect exists to fix. + +## What this task was + +Release run `28784914443` (PyAutoHeart#27, the first real release-profile validation) reported 42 +script failures, split across seven prompts in `PyAutoMind/draft/bug/health_fixes/`. This task owned +two, both alleging `LinAlgError`-class failures from non-positive-definite matrices in inversion +paths: + +- `autolens_workspace_test/scripts/interferometer/model_fit.py` +- `autogalaxy_workspace/scripts/interferometer/features/pixelization/galaxy_reconstruction.py` + +The prompt asserted the autolens leg "reproduces on current `main`" and prescribed a five-step +repair: capture curvature/regularization matrix properties at failure, localise the defect between +sampled parameters, regularization construction, numerical stabilization or an underdetermined +script model, then fix the owning library. + +Brain sized it `too-large` (score 16), fix-locus "library source", strategy "split into phases". + +## What was actually done + +Not started as the prescribed repair. The premise was gated first, for the same reason its sibling +`autofit_sampler_database` was: the claim was six weeks old, and the cluster around it had already +produced two independent findings that it had aged out. + +Both scripts were re-run on current `main` from a **cleared** `output/`, under each workspace's +`config/build/profile_release.yaml`, env resolved by `autohands.env_config.build_env_for_script` at +workspace CWD, 1800s `mode=release` cap. Libraries at `main`: PyAutoFit `248ca971f`, PyAutoArray +`b808a9b1`, PyAutoGalaxy `7e3856dd`, PyAutoLens `d8f6bb3df`, PyAutoNerves `f6d6d52`. Three workspace +checkouts were behind `origin/main` and were synced first. + +## Result: 0 / 2 reproduce + +| Script | Result | Secs | +|---|---|--:| +| `autolens_workspace_test/scripts/interferometer/model_fit.py` | PASS | 78 | +| `autogalaxy_workspace/scripts/interferometer/features/pixelization/galaxy_reconstruction.py` | PASS | 70 | + +No non-positive-definite failure and no `LinAlgError` in either. The prompt's specific claim that the +autolens leg reproduces on `main` is false. + +## Why this verdict is strong: two prior independent refutations of the same hypothesis + +This is not a lone green run. The same hypothesis — that PyAuto's inversion path produces +non-positive-definite curvature/regularization matrices — had already been tested and refuted twice: + +1. **`complete/2026/07/pix-inversion-not-positive-definite.md`** (2026-07-21) investigated a + six-marker `LinAlgError: matrix not positive definite / singular` cluster across + autogalaxy_workspace + HowToGalaxy. Outcome **inverted**: all six markers were stale and **no code + was fixed**. The pix `LinAlgError` had been cured on 2026-04-10 — the same day the markers were + filed — by PyAutoArray's `GaussianKernel` PD-guarantee `f1817af0` (symmetrise + trace-scaled + diagonal jitter). Evidence was a 40-draw numpy inversion A/B across the full `GaussianKernel` + LogUniform prior (coeff/scale `1e-6`..`~5e5`): **0 raises / 0 non-finite** on both `cholesky` and + `slogdet`. +2. **`complete/2026/08/autofit-sampler-database.md`** (PyAutoFit#1508, 2026-08-21) — the sibling from + this same release run, **0/9 reproduce**, closed with no code change. + +A PD-guarantee landed in the owning library four months ago, was independently verified by direct +matrix probing, and both of this prompt's scripts now pass. The defect described here does not exist +on `main`. + +## Why the issue was closed rather than parked open + +The same test the sibling applied, and it passes here: **neither script is parked.** Checked against +`main` on 2026-08-22 — + +- `autogalaxy_workspace/config/build/no_run.yaml` — 8 entries, none matching + `interferometer/features/pixelization/galaxy_reconstruction` (GUI scripts, fits/png_make, + search-viz, and one SLOW shapelets entry). +- `autolens_workspace_test/config/build/no_run.yaml` — no `interferometer/model_fit` entry. It + appears only as a *consumer* on line 62, where `interferometer/simulator/with_lens_light.py` is + marked `BOOTSTRAP-TARGET` because it produces `model_fit`'s dataset — which confirms `model_fit` + itself runs. + +So both scripts re-execute under exactly this profile in **every** `mode=release` pass. Re-validation +is automatic; a surviving defect fails the next release run loudly and earns a fresh issue with fresh +evidence. There is no human reminder to lose — the same reasoning that closed PyAutoFit#1508, and the +reason this prompt closes while `samples_parameter_paths` (#1327) stays parked. + +This is now the **fourth** independent finding that this cluster aged out. The consistent explanation +across all of them is the one #1327 reached: stale cached `output/` in the 2026-07 release run, +against libraries that have since absorbed dozens of fixes. + +## What this does NOT establish + +1. **The autolens leg ran on numpy, not JAX.** `autolens_workspace_test`'s release profile *defaults* + `PYAUTO_DISABLE_JAX="1"`; scripts opt back in with an in-file `ENV: jax` declaration, and + `model_fit.py` carries none. That is release-faithful — it is what the release run itself executes + — but JAX-on and JAX-off are different numerical code paths, and this refutation covers only the + numpy one. A JAX-only inversion conditioning defect would not appear here. +2. **These were source-tree runs**, not the TestPyPI wheels the release run installed. A wheel-only + packaging defect would not show. + +## Incidental finding — filed as its own prompt + +`galaxy_reconstruction.py` passes while emitting 4x `RuntimeWarning: invalid value encountered in +sqrt` from `PyAutoArray/autoarray/inversion/inversion/abstract.py:859`. **This looks exactly like +evidence for the prompt's hypothesis and is not** — it is a separate, unconditional defect: + +```python +@property +def reconstruction_noise_map_with_covariance(self) -> np.ndarray: + return np.sqrt(np.linalg.inv(self.curvature_reg_matrix)) +``` + +`sqrt` is applied **elementwise to the whole inverse matrix**, whose off-diagonal entries are +covariances and are generally negative — so those entries are NaN *by construction*, for any input +matrix, however well-conditioned. It is not a conditioning symptom and does not rescue the prompt. + +Confirmed still present on PyAutoArray `main` @ `a6b07cd` (2026-08-22). The 1D +`reconstruction_noise_map` is unaffected — it takes the diagonal, and elementwise-sqrt commutes with +taking the diagonal — so the science path is correct; only the covariance-aware consumer and the +warning spam are hit. + +Filed as **`draft/bug/autoarray/reconstruction_noise_map_covariance_sqrt.md`**, not fixed here: the +correct off-diagonal semantics are an API/science decision, and the fix has a real trap (the 1D +science path is derived from this property's diagonal, so a naive change silently converts a +standard-deviation into a variance). + + +**2026-08-22 follow-up — the incidental finding got bigger.** Research into *why* source-reconstruction +noise maps have been unreliable found the sqrt bug is **not** the cause: `np.sqrt` is elementwise, so it +commutes with taking the diagonal and provably cannot reach the 1D noise map. Three deeper defects were +found in the same property, and split into a second prompt, +`draft/bug/autoarray/reconstruction_noise_map_solver_mismatch.md`: + +1. **Estimator mismatch (the big one).** `inv(curvature_reg_matrix)` is the posterior covariance of the + *unconstrained* Warren & Dye solve. But `use_positive_only_solver: true` is the shipped default, so the + reconstruction is an NNLS active-set solve. Imposing `s >= 0` truncates the posterior, so the reported + noise is overstated near the boundary and meaningless for pinned pixels — and a compact lensed source + pins a large fraction of the mesh at zero, so the formula is worst exactly where it is most used. +2. **Edge-zeroed pixels ignored.** `use_edge_zeroed_pixels: true` is also default; the reconstruction + solves on `zeroed_ids_to_keep` and scatters back zeros, while the noise map inverts the *full* matrix — + re-admitting the poorly-constrained boundary vertices the zeroing exists to remove. +3. **`use_edge_zeroed_pixels` is nested inside the positive-only branch**, so turning the positive-only + solver off silently disables edge-zeroing with no warning. Two orthogonal settings, silently coupled. + +Corroboration for the numerics half: `abstract.py:805` already documents `~1e-6` evidence round-off from +"factorizing the explicitly formed inverse" at `cond(C) ~ 1e9` on clustered traced mesh vertices — applied +to the log-det, never to the noise map. And `inversion_plots.py:395` already wraps the noise map in +`except np.linalg.LinAlgError`, writing NaN to the CSV: a guard that exists because this fails in practice. + +## Follow-on state of the cluster + +`draft/bug/health_fixes/README.md` row struck through. Of the original seven prompts: two shipped or +closed with code (`aggregator_output_contracts`), two closed as refuted with no code change +(`autofit_sampler_database`, this one), one parked not-reproducing (`samples_parameter_paths`, +#1327), and three remain in `draft/` with dated gate annotations +(`jax_runtime_and_parity`, `jit_visualization_outputs`, `release_timeout_policy`) — those are *not* +complete, because their SLOW/NEEDS_FIX parkings describe *intermittent* failures that a single green +run cannot clear. + +## Original prompt + +# Fix release-profile numerical inversion failures + +Type: bug +Target: health_fixes +Difficulty: too-large +Autonomy: supervised +Priority: high +Status: formalised + +## Context + +Two interferometer scripts fail in inversion paths with non-positive-definite matrices. +The Autolens test failure reproduces on current `main`; the Autogalaxy script passed in a +stateful local checkout and needs a clean confirmation. + +Owners: @PyAutoArray, @PyAutoGalaxy, @PyAutoLens, @autogalaxy_workspace, and +@autolens_workspace_test. + +## Scripts + +- `autogalaxy_workspace/scripts/interferometer/features/pixelization/galaxy_reconstruction.py` +- `autolens_workspace_test/scripts/interferometer/model_fit.py` + +## Required work + +1. Reproduce in clean output/worktrees with deterministic seeds and release settings. +2. Capture the curvature and regularization matrix properties at failure: symmetry, + conditioning, eigenvalue range, dtype, backend, and mapper configuration. +3. Identify whether the defect is invalid sampled parameters, regularization construction, + numerical stabilization, or a script model that permits an undefined inversion. +4. Fix the owning library for valid inputs. Do not catch `LinAlgError` or alter the script + to hide a genuine inversion failure. +5. Add numerical regression tests and rerun both scripts repeatedly under the profile. + + + +## 2026-08-21 — REPRODUCTION GATE RUN: **2/2 PASS — prompt refuted** + +Method (identical to the gate that closed the sibling `autofit_sampler_database`, PyAutoFit#1508): +every script run from a **cleared** `output/`, under its workspace's +`config/build/profile_release.yaml`, env resolved by `autohands.env_config.build_env_for_script` +at workspace CWD, 1800s `mode=release` cap. Libraries at `main`: PyAutoFit `248ca971f`, +PyAutoArray `b808a9b1`, PyAutoGalaxy `7e3856dd`, PyAutoLens `d8f6bb3df`, PyAutoNerves `f6d6d52`. +Three workspace checkouts were **behind `origin/main`** and were synced first. + +| Script | Result | Secs | +|---|---|--:| +| `autolens_workspace_test/scripts/interferometer/model_fit.py` | PASS | 78 | +| `autogalaxy_workspace/scripts/interferometer/features/pixelization/galaxy_reconstruction.py` | PASS | 70 | + +The prompt states the autolens leg "reproduces on current `main`". It does not. No +non-positive-definite failure, no `LinAlgError`, in either. + +**Note which numerical path each took.** `autolens_workspace_test`'s release profile *defaults* +`PYAUTO_DISABLE_JAX="1"`, and scripts opt back in with an in-file `ENV: jax` declaration. +`model_fit.py` has no such declaration, so it ran on **numpy** — release-faithful, but worth +knowing for a claim about inversion numerics, since JAX-on and JAX-off are different code paths. + +### Incidental finding — a real defect, but NOT this prompt's + +`galaxy_reconstruction.py` passes while emitting 4x +`RuntimeWarning: invalid value encountered in sqrt` from +`PyAutoArray/autoarray/inversion/inversion/abstract.py:859`: + +```python +def reconstruction_noise_map_with_covariance(self): + return np.sqrt(np.linalg.inv(self.curvature_reg_matrix)) +``` + +`sqrt` is applied **elementwise to the whole inverse matrix**, whose off-diagonal entries are +covariances and are generally negative — so those entries are NaN *by construction*, for any +matrix, however well-conditioned. + +**This is not evidence of a non-positive-definite matrix** and does not rescue the prompt's +hypothesis, despite looking exactly like it would. It is a separate defect: a property whose +docstring promises a matrix that "accounts for the covariance of the noise between pixels" returns +NaN wherever that covariance is negative. The 1D `reconstruction_noise_map` is unaffected — it +takes the diagonal, and `diag(sqrt(M)) == sqrt(diag(M))` — so the science path is correct; only +the covariance-aware consumer and the warning spam are hit. Worth its own PyAutoArray prompt. diff --git a/complete/2026/08/reconstruction-noise-map-covariance-sqrt.md b/complete/2026/08/reconstruction-noise-map-covariance-sqrt.md new file mode 100644 index 00000000..a3e02330 --- /dev/null +++ b/complete/2026/08/reconstruction-noise-map-covariance-sqrt.md @@ -0,0 +1,320 @@ +# reconstruction-noise-map-covariance-sqrt — covariance NaNs and the Cholesky rewrite + +**Date:** 2026-08-22 +**Issue:** [PyAutoArray#468](https://github.com/PyAutoLabs/PyAutoArray/issues/468) (closed) +**PRs:** [PyAutoArray#469](https://github.com/PyAutoLabs/PyAutoArray/pull/469) (MERGED, `2784056`) +**Outcome:** shipped — phase 1 of a two-phase cluster; the estimator half is still open. + +## What this task was + +`AbstractInversion.reconstruction_noise_map_with_covariance` was one line: + +```python +return np.sqrt(np.linalg.inv(self.curvature_reg_matrix)) +``` + +`np.sqrt` applied **elementwise to the whole inverse**. The off-diagonals of a covariance matrix are +covariances and are routinely negative, so every one was `NaN` by construction — for any matrix, +however well-conditioned — with a `RuntimeWarning` on every call. The docstring promised a matrix +that "accounts for the covariance of the noise between pixels"; the entries carrying that covariance +were precisely the broken ones. + +Found as an incidental finding while reproduction-gating +`complete/2026/08/numerical-inversion-failures.md`, whose own hypothesis was refuted. + +## The trap this cluster kept setting + +**The elementwise sqrt does not affect the 1D `reconstruction_noise_map`.** `np.sqrt` is elementwise, +so it commutes with taking the diagonal: `diagonal(sqrt(C))[i] == sqrt(C[i,i])`. This was mistaken +for the cause of unreliable source noise maps twice — once by the original gate, once during +planning. It is not. It looks exactly like evidence for a non-positive-definite matrix and is not +that either. + +## The A/B refuted two of the arguments for the fix + +Run **before** writing any code, which is why the reasoning in the shipped docstring is narrower than +the reasoning in the original prompt: + +| Claim | Verdict | +|---|---| +| `inv` gives negative diagonals on well-formed SPD | **REFUTED** — 0 across cond 1e3–1e15, n=400, 20 trials each | +| `inv` is materially less accurate on the diagonal | **REFUTED** — matches `cho_solve`; at cond 1e15 `inv` was marginally *better* | +| Near-coincident mesh vertices degrade the inverse | **REFUTED** — with regularization the matrix stays PD (cond ~6.8e7 even at exactly duplicated columns) | +| `inv` returns asymmetric output | **CONFIRMED** — 5.2e-7 at cond 1e12 vs 2.6e-16 | +| `inv` silently succeeds on indefinite matrices | **CONFIRMED** | + +**The case for Cholesky is detection, not accuracy.** `cho_factor` raises `LinAlgError` on a negative +eigenvalue; `np.linalg.inv` raises only on an *exactly* singular matrix. At eigenvalue `-1e-8` all +300 diagonals came back negative (whole noise map NaN); at `-1.0`, **zero** did — no NaN, no warning, +no error, and wrong numbers. That silent case is the failure mode worth fixing. + +The `abstract.py:805` note about "the round-off of factorizing the explicitly formed inverse" at +`cond ~ 1e9` was cited as corroboration during planning. The A/B did not support it as an *accuracy* +argument for the noise map; keep it as history, not as evidence. + +## What shipped + +- **new `reconstruction_covariance_matrix`** — `cho_solve(cho_factor(...))`, input *and* output + symmetrized, explicit finiteness guard raising `LinAlgError`. +- **`reconstruction_noise_map`** — decoupled to `sqrt(diag(C))`. It was previously correct only + *incidentally*, via the elementwise sqrt; the invariant is now stated so it cannot silently become + a variance. +- **`reconstruction_noise_map_with_covariance`** — deprecated alias, warning states the value change. +- 5 new regression tests plus 2 from review; 2 plotter monkeypatch sites repointed. + +## The review caught a regression the first commit introduced + +Worth recording, because it is a trap anyone swapping `inv` for a scipy factorization will hit: + +**scipy's `cho_factor`/`cho_solve` default to `check_finite=True` and raise `ValueError`, not +`LinAlgError`.** Both call sites (`inversion_plots.py:169`, `:397`) catch only `LinAlgError`, and the +CSV writer's docstring explicitly promises a failure there may not abort the enclosing model-fit. So +a NaN-contaminated curvature matrix would have killed a fit that previously wrote a `nan` column and +continued. Measured: old code returned `[nan, nan]`; first fix raised `builtins.ValueError` past the +guards. + +Fixed *inside* the property with an explicit finiteness check rather than by broadening the callers' +`except` clauses — the contract "raises `LinAlgError` for any input that has no covariance" then +holds for downstream callers too. `check_finite=False` is passed onward, so the check costs nothing. + +Second review find: **`cho_factor` reads only the upper triangle**, so an asymmetric input was +silently inverted as though its lower triangle matched (`[[2.0, 0.5], [0.1, 2.0]]` → diag `0.5333` +vs the true `0.5063`). Output symmetrization does not fix that; the input is symmetrized now too. + +Also corrected: the "no value change" claim for `reconstruction_noise_map` was overstated — +algebraically identical, only *numerically* equivalent (~7e-15 at cond 1e3, ~4e-5 at cond 1e13). +And the symmetry test was tautological, since `0.5 * (C + C.T)` is bitwise symmetric for any `C`; it +now also asserts accuracy against an exactly-constructed ground truth. + +## Downstream API risk — resolved by grep, not assumption + +The deprecated alias changes values under an unchanged name, and `DeprecationWarning` is invisible by +default when raised from library code. So it was checked: + +| Repo | `with_covariance` | `reconstruction_noise_map` | +|---|---|---| +| PyAutoGalaxy `3ca31bf` | none | none | +| PyAutoLens `87e5827` | none | none | +| autolens_workspace | none | **4 scripts + notebooks** | + +Control greps confirm the checkouts were real (394 / 237 / 465 `.py` files), so the nulls are genuine. +**Not checked:** `autogalaxy_workspace`, the HowTo repos, external user code. + +`reconstruction_noise_map` **is** used, and this is the finding that matters for the sibling prompt: +`autolens_workspace/scripts/{imaging,interferometer,group,multi_galaxy}/features/pixelization/source_science.py` +compute `signal_to_noise_map = reconstruction / reconstruction_noise_map`. The noise map feeds +published S/N maps on source reconstructions. + +## Still open — the larger half + +`draft/bug/autoarray/reconstruction_noise_map_solver_mismatch.md` (`Priority: high`, +`Autonomy: human-required`) holds the estimator-level defects, deliberately excluded from this PR: + +1. The covariance is that of the **unconstrained** Warren & Dye solve, but + `use_positive_only_solver: true` is the shipped default, so the reconstruction is an NNLS + active-set solve. A compact source pins a large fraction of the mesh at zero. +2. The noise map ignores `zeroed_ids_to_keep` under `use_edge_zeroed_pixels: true`, re-admitting the + poorly-constrained boundary vertices the zeroing exists to remove. +3. `use_edge_zeroed_pixels` is nested inside the positive-only branch, so turning that solver off + silently disables edge-zeroing. + +First job there is to instrument a real fit: **it was never established that a real +`curvature_reg_matrix` is indefinite in a converged fit**, nor that the NNLS pinned fraction is +actually large. Both were reasoned, not measured, and the whole estimator argument rests on the +second. + +## Environment note + +Python 3.11 is too old for this repo (`requires-python >= 3.12`); `python3.12` was present and a venv +there installed cleanly. Three `test_transformer.py` pynufft failures reproduce on clean `a6b07cd` +in that sandbox and were **green in CI** — a local dependency artefact, not a repo problem. + +## Original prompt + +# `reconstruction_noise_map_with_covariance` — form the covariance properly, fix the sqrt + +Type: bug +Target: autoarray +Repos: +- @PyAutoArray +Difficulty: small +Autonomy: supervised +Priority: medium +Status: in-progress + +## Why this exists + +Found during the 2026-08-21 reproduction gate for +`draft/bug/health_fixes/numerical_inversion_failures.md` (record: +`complete/2026/08/numerical-inversion-failures.md`). That prompt alleged +non-positive-definite inversion matrices and was **refuted**; this was the real +defect the gate turned up. + +**Scope note (2026-08-22).** Deeper research found the noise map is wrong in +*four* distinct ways. This prompt owns the two that need no science decision — +the numerics and the semantics. The estimator-level defects (the noise map +describes a different estimator than the default solver, and ignores edge-zeroed +pixels) are **`draft/bug/autoarray/reconstruction_noise_map_solver_mismatch.md`**, +which is the larger and more consequential of the two. Do that one second; this +one first, because it is small, safe and unblocks the other. + +## Defect A — elementwise sqrt NaNs every off-diagonal + +`autoarray/inversion/inversion/abstract.py:839-859`, verified on `main` @ `a6b07cd`: + +```python +@property +def reconstruction_noise_map_with_covariance(self) -> np.ndarray: + """... a two dimension matrix which accounts for the covariance of the noise between pixels.""" + return np.sqrt(np.linalg.inv(self.curvature_reg_matrix)) +``` + +`np.sqrt` is applied **elementwise to the entire inverse**. That inverse is the +covariance matrix `C`, whose off-diagonals are covariances and are generally +negative — so they become `NaN`, and each call emits +`RuntimeWarning: invalid value encountered in sqrt`. Unconditional: any matrix +with an anti-correlated pixel pair NaNs, however well-conditioned. + +**This does NOT affect the 1D noise map.** `np.sqrt` is elementwise, so it +commutes with taking the diagonal — `np.diagonal(np.sqrt(C))[i] == sqrt(C[i,i])`. +The off-diagonal NaNs never reach `reconstruction_noise_map`. Do not cite this +defect as the cause of unreliable 1D noise maps; that is the sibling prompt. + +## Defect B — `np.linalg.inv` is the wrong routine, and the repo already says so + +The same file, 50 lines up, documents exactly this hazard for the log-det path +(`abstract.py:805-806`): + +> the analytically exact `pixels * log(coeff) - log det C` from a single Cholesky +> of their covariance `C`, **avoiding the round-off of factorizing the explicitly +> formed inverse** (which reaches ~1e-6 absolute in the evidence at +> **cond(C) ~ 1e9 on clustered traced mesh vertices**) + +That reasoning was applied to the log-det and never to the noise map, which still +forms the explicit inverse — of the same matrix, at the same conditioning, on the +same clustered-mesh geometry. + +Three consequences, all pointing the same way: + +1. `np.linalg.inv` is LU-based. It exploits neither symmetry nor + positive-definiteness, both of which this matrix has (when it is well-posed). +2. It **raises only on exactly-singular input.** Near-singular passes through with + amplified error, so a diagonal entry can come back negative — impossible for a + true PD inverse — and `sqrt` turns it into NaN, or leaves it barely positive and + yields a wildly wrong RMS. Silently. +3. The reconstruction path never inverts: `reconstruction_positive_negative_from` + uses `xp.linalg.solve` and `fnnls_cholesky` uses + `slg.solve(..., assume_a="pos")`. The noise map is the only place in the + inversion that forms an explicit inverse. + +**This is already biting users.** `inversion_plots.py:395` wraps the noise map in +`except np.linalg.LinAlgError` and writes the CSV column as NaN with a warning — +a guard that exists because this fails in practice. + +## The fix + +Option 1 from the original draft, chosen 2026-08-22: the property should return +the actual covariance matrix, computed via Cholesky. `scipy` is already a hard +dependency (`pyproject.toml`). + +```python +from scipy.linalg import cho_factor, cho_solve + +@property +def reconstruction_covariance_matrix(self) -> np.ndarray: + """The covariance matrix C = [F + λH]^-1 of the reconstruction.""" + matrix = np.asarray(self.curvature_reg_matrix) + covariance = cho_solve(cho_factor(matrix), np.eye(matrix.shape[0])) + return 0.5 * (covariance + covariance.T) # remove rounding asymmetry + +@property +def reconstruction_noise_map(self) -> np.ndarray: + """1D RMS noise: sqrt of the diagonal of the covariance matrix.""" + return np.sqrt(np.diag(self.reconstruction_covariance_matrix)) +``` + +Why this shape: + +- **`cho_factor` raises `LinAlgError` on a non-PD matrix**, so the noise map now + fails loudly on exactly the matrices the reconstruction already rejects. Today + the two disagree: `solve` raises and resamples, `inv` returns garbage. +- **`reconstruction_noise_map` is decoupled** and computes `sqrt(diag(C))` + directly. Today it is correct only *incidentally*, because sqrt happens to be + elementwise — change the matrix and it silently becomes a variance. Decoupling + removes that trap permanently. +- **Off-diagonals become real covariances**, so the docstring's promise holds. + +**Naming.** `..._with_covariance` returning a covariance matrix should be +`reconstruction_covariance_matrix`. Keep the old name as a `DeprecationWarning` +alias returning the new matrix. Its values *do* change — diagonal from std-dev to +variance, off-diagonals from NaN to covariances — but every off-diagonal consumer +was reading NaN, so nothing correct can break. Note the change in the release +notes regardless. + +Optional, only if profiling asks for it: if just the diagonal is needed, +`diag(C)` is available from the Cholesky factor as the squared row-norms of +`L^-1`, avoiding the full `n x n` product. Not worth the complexity up front — +this is computed once per fit, not per-likelihood. + +## Verification + +- **Off-diagonals finite** for a well-conditioned matrix with an anti-correlated + pixel pair. **No such test exists today** — the only assertion on this property + (`test_autoarray/inversion/inversion/test_abstract.py:684`) checks `[0, 0]`, a + *diagonal* element. That gap is why this shipped. +- **No `RuntimeWarning`.** Run the regression test under `-W error::RuntimeWarning` + so a regression fails rather than warns. +- **`reconstruction_noise_map` still returns `sqrt(diag(C))`** — assert against a + hand-computed value, and assert the invariant explicitly, not just the numbers. +- **The `inv`-vs-`cho_solve` A/B was run on 2026-08-22 and refuted two of the + claims above.** Recorded so nobody re-derives the wrong reasoning: + + | Claim | Verdict | + |---|---| + | `inv` gives negative diagonals on well-formed SPD | **REFUTED** — 0 across cond 1e3–1e15, n=400, 20 trials each | + | `inv` is materially less accurate on the diagonal | **REFUTED** — matches `cho_solve`; at cond 1e15 `inv` was marginally *better* | + | Near-coincident mesh vertices degrade the inverse | **REFUTED** — with regularization the matrix stays PD (cond ~6.8e7 even at exactly duplicated columns) | + | `inv` returns asymmetric output | **CONFIRMED** — 5.2e-7 at cond 1e12 vs 2.6e-16 | + | `inv` silently succeeds on indefinite matrices | **CONFIRMED** | + + The surviving argument is **detection, not accuracy**. `cho_factor` raises + `LinAlgError` on a negative eigenvalue; `inv` raises only on an *exactly* + singular matrix and otherwise returns a plausible-looking covariance. At + eigenvalue `-1e-8` all 300 diagonals came back negative (whole noise map NaN); + at `-1.0`, **zero** did — no NaN, no warning, no error, wrong numbers. + + **Not established:** that a real `curvature_reg_matrix` *is* indefinite in a + converged fit. `Settings.no_regularization_add_to_curvature_diag_value` and the + `curvature_matrix_with_added_to_diag_from` docstring ("it is common for the + `curvature_matrix` computed to not be positive-definite") say it happens, but no + fit was instrumented to confirm it. Worth doing under the sibling prompt. +- `test_autoarray/inversion/plot/test_inversion_plotters.py:82,110` monkeypatch + this property to force a `LinAlgError` and check plots/CSV degrade gracefully. + Confirm the same exception still escapes — `cho_factor` also raises + `LinAlgError`, so this should hold, but assert it. +- **Downstream:** this sweep covered PyAutoArray only, where the sole in-repo + consumer is `reconstruction_noise_map`. Grep @PyAutoGalaxy and @PyAutoLens for + `reconstruction_noise_map_with_covariance` before assuming containment. + +## Also fold in + +`reconstruction_noise_map`'s docstring claims it "is computed as the square root +of the diagonal of the `reconstruction_noise_map_with_covariance` matrix". The +code takes the diagonal of an already-square-rooted matrix — no second sqrt. The +two agree today only because sqrt is elementwise, which is precisely the bug. +Rewrite the sentence to match whatever ships. + +## Note on the JAX path + +This property uses bare `np`, not `self._xp`, so it is already numpy-only even +under a JAX fit — a JAX `curvature_reg_matrix` is coerced via `__array__`, forcing +a device→host sync. The scipy fix does not regress that (there was no JAX support +to lose) but it does make it explicit. Add `np.asarray` at the boundary, as above, +and note the limitation in the docstring rather than leaving it implicit. + +## Provenance + +- Found during: `complete/2026/08/numerical-inversion-failures.md` (2026-08-22) +- Sibling: `draft/bug/autoarray/reconstruction_noise_map_solver_mismatch.md` +- **Not** a symptom of the refuted non-positive-definite hypothesis in that record, + nor of `complete/2026/07/pix-inversion-not-positive-definite.md` (also refuted). diff --git a/complete/index.md b/complete/index.md index 396191e3..3ef4dc23 100644 --- a/complete/index.md +++ b/complete/index.md @@ -6,7 +6,7 @@ Token-light navigation over the finished-work records (schema: only then grep a dated bucket. Curators: edit the band between the CURATED markers; everything below GENERATED is rebuilt. -1074 records across 7 buckets. +1076 records across 7 buckets. ## Highlights @@ -123,6 +123,7 @@ _(curate hard-won records here — survives regeneration.)_ - [numba-cpu-likelihood-profiling](2026/08/numba-cpu-likelihood-profiling.md) - [numba-first-call-garbage-psf-weighted-data](2026/08/numba-first-call-garbage-psf-weighted-data.md) - [numerical-hazard-profiling](2026/08/numerical-hazard-profiling.md) — Phase 1 established the reusable numerical-hazard profiling framework, typed risk measurements, reachability a… +- [numerical-inversion-failures](2026/08/numerical-inversion-failures.md) - [one-tap-dashboard-rollout](2026/08/one-tap-dashboard-rollout.md) — the pattern (a generated page + one-tap 📋 copy-for-Claude payloads, phone - [optional-dependency-skip-guards](2026/08/optional-dependency-skip-guards.md) — `test_nautilus.py::test__single_core_builds_no_pool` had no skip guard - [organism-board-final-readmes](2026/08/organism-board-final-readmes.md) — close on record @@ -156,6 +157,7 @@ _(curate hard-won records here — survives regeneration.)_ - [queue-intake-claude](2026/08/queue-intake-claude.md) — auto-closed on merge - [readme-scientist-buttons](2026/08/readme-scientist-buttons.md) - [reconcile-upstream-repo-mode](2026/08/reconcile-upstream-repo-mode.md) +- [reconstruction-noise-map-covariance-sqrt](2026/08/reconstruction-noise-map-covariance-sqrt.md) - [rectangular-spline-adapt-image-jit-fix](2026/08/rectangular-spline-adapt-image-jit-fix.md) - [registry-integrity-check](2026/08/registry-integrity-check.md) - [regularization-jax-gradient-gaps](2026/08/regularization-jax-gradient-gaps.md) diff --git a/dashboard.html b/dashboard.html index 3194c422..64ac3eea 100644 --- a/dashboard.html +++ b/dashboard.html @@ -40,7 +40,7 @@

📋 PyAutoMind Dashboard

Every task the Mind is holding. Tap a task's 📋 and its /start_dev command is on your clipboard — paste it into a Claude Code chat to route Claude straight to that task.

In flight 2 · Parked 3 · Planned 6 · Backlog 146 · markdown version

Start here

-

Highest priority (filed as high) — showing 12 of 19

+

Highest priority (filed as high) — showing 12 of 18

TRIAGE: needs manual review before routingmedium · safe · high

Numba CPU likelihood phase 1: batched MGE convolution + operated-matrixautoarray · medium · supervised · high

Rectangular mesh split: Bilinear (fast CPU default) vs RTU (advanced/GPU)autoarray · medium · supervised · high

@@ -114,7 +114,6 @@

Backlog

Stamp the small-datasets regime at the FITS writer funnelpyautonerves · large · supervised · high

Fix release JAX runtime compatibility and likelihood parityhealth_fixes · too-large · supervised · high

Fix JIT quick-update visualization output regressionshealth_fixes · too-large · supervised · high

-

Fix release-profile numerical inversion failureshealth_fixes · too-large · supervised · high

Fix release result/sample parameter-path regressionshealth_fixes · too-large · supervised · high

pixel_scales given as an int (or np.float64) is never widenedautoarray · small · supervised · medium

jax 0.11 breaks beta/gamma message log_partition under jit ('tuple' objectautofit · small · supervised · medium

@@ -143,6 +142,7 @@

Backlog

Priors & Messages cleanup — trackerpriors · too-large · supervised · normal

autolens_workspace_test jax_likelihood pins: 4 scripts fail smoke on mainautolens · low · supervised · normal

status.sh --repos sources a file that no longer existspyautomind · small · supervised · low

+

The reconstruction noise map describes a different estimator than theautoarray · medium · human-required · low

interferometer/start_here.py OOM in nightly release-validation integrate legautolens

diff --git a/dashboard.md b/dashboard.md index df646139..2778aa1a 100644 --- a/dashboard.md +++ b/dashboard.md @@ -15,7 +15,7 @@ Every task the Mind is holding, on one page: what is in flight, what is parked, ## Start here -**Highest priority** (filed as `high`) — showing 12 of 19 +**Highest priority** (filed as `high`) — showing 12 of 18
📋 TRIAGE: needs manual review before routing — medium · safe · high @@ -487,14 +487,6 @@ Scoped but not started; some are not yet prompt files. Full detail in [`planned.
-
📋 Fix release-profile numerical inversion failures — health_fixes · too-large · supervised · high - -``` -/start_dev draft/bug/health_fixes/numerical_inversion_failures.md -``` - -
-
📋 Fix release result/sample parameter-path regressions — health_fixes · too-large · supervised · high ``` @@ -719,6 +711,14 @@ Scoped but not started; some are not yet prompt files. Full detail in [`planned.
+
📋 The reconstruction noise map describes a different estimator than the — autoarray · medium · human-required · low + +``` +/start_dev draft/bug/autoarray/reconstruction_noise_map_solver_mismatch.md +``` + +
+
📋 interferometer/start_here.py OOM in nightly release-validation integrate leg — autolens ``` diff --git a/draft/bug/autoarray/reconstruction_noise_map_solver_mismatch.md b/draft/bug/autoarray/reconstruction_noise_map_solver_mismatch.md new file mode 100644 index 00000000..7c214246 --- /dev/null +++ b/draft/bug/autoarray/reconstruction_noise_map_solver_mismatch.md @@ -0,0 +1,429 @@ +# The reconstruction noise map describes a different estimator than the default solver + +Type: bug +Target: autoarray +Repos: +- @PyAutoArray +Difficulty: medium +Autonomy: human-required +Priority: low +Status: draft + +## Why this exists + +Found 2026-08-22 while researching +`draft/bug/autoarray/reconstruction_noise_map_covariance_sqrt.md` (the numerics +and semantics half). The question that started it: *why has the noise map on +source reconstructions not always been reliable?* + +The answer is not the elementwise-sqrt bug — that only touches off-diagonals and +provably cannot reach the 1D noise map. It is this: **the noise map is computed +from a formula that describes an estimator PyAutoArray no longer uses by +default.** `Autonomy: human-required` because the primary fix is a statistical +decision about what the reported uncertainty *means*, not a code repair. + +## DECISION MADE 2026-08-22 — posterior, not diagnostic + +The human was asked whether this noise map is a **posterior** (a calibrated uncertainty on the +reconstructed flux, used for error bars and S/N) or a **diagnostic** (a relative readout of how +well-constrained each pixel is). Answer: **posterior.** + +That settles the direction and closes the diagnostic-only escape hatch offered under Defect 1 +("rewrite the docstring rather than the maths"). Defects 1 and 2 are real defects, not documentation +drift. What remains open is the *design* of the fix, not whether one is needed. + +## Defect 1 — the covariance formula assumes an unconstrained solve; the default is NNLS + +`abstract.py:859` computes `C = inv(curvature_reg_matrix)` = `[F + λH]^-1`. That +is the posterior covariance of the **classical semi-linear inversion** — the +unconstrained linear-Gaussian solution of Warren & Dye (2003) eq. 12, which is +what `reconstruction_positive_negative_from` computes. + +But `config/general.yaml` ships: + +```yaml +use_positive_only_solver: true # DEFAULT +``` + +So the default reconstruction is `fnnls_cholesky` — a **non-negative least +squares** solve, minimising `||Zs - x||²` subject to `s >= 0`. That is a +*constrained* estimator with an active set: `fnnls` maintains a passive set `P` +and solves `slg.solve(ZTZ[P][:,P], ..., assume_a="pos")` on the free pixels only, +pinning the rest at exactly zero. + +Imposing `s >= 0` is equivalent to truncating the Gaussian prior to the +non-negative orthant. The posterior is then a **truncated** multivariate Gaussian, +whose covariance is *not* `[F + λH]^-1`: + +- For pixels well inside the positive region, the constraint is inactive and + `[F + λH]^-1` is a good approximation. +- For pixels near the boundary, truncation **reduces** the variance, so the + reported noise is systematically **overstated**. +- For pixels pinned at exactly zero, the marginal posterior is not Gaussian at + all — it piles up at the boundary. The reported number is meaningless. + +**Why this bites source reconstructions hardest.** A lensed source is compact; most +of the mesh is empty sky. So NNLS pins a *large fraction* of pixels at zero, the +active set is large, and the unconstrained formula is worst exactly where it is +most used. That matches the reported symptom precisely. + +The docstring makes an uncertainty claim, not a diagnostic one — "the RMS standard +deviation of the noise in every pixel ... should be used for any scientific +analysis (e.g. source reconstructions of strong lenses)" — so this matters. + +**The honest counter-argument, which a human should weigh:** if the noise map is +meant only as a *diagnostic* of how well each pixel is constrained by data plus +regularization, `[F + λH]^-1` is defensible for any solver, and the fix is to +rewrite the docstring rather than the maths. Decide which of the two it is before +writing code. That decision is the point of this prompt. + +## Defect 2 — the noise map ignores edge-zeroed pixels the reconstruction excluded + +`config/general.yaml` also ships `use_edge_zeroed_pixels: true`. Under it, the +reconstruction (`abstract.py:511-539`) subsets the system: + +```python +curvature_reg_matrix = self.curvature_reg_matrix[self.zeroed_ids_to_keep][:, self.zeroed_ids_to_keep] +``` + +solves the reduced problem, then scatters back with **exact zeros** at the zeroed +pixels. `reconstruction_noise_map_with_covariance` inverts `self.curvature_reg_matrix` +— the **full** matrix, respecting neither this reduction nor the separate +`mapper_indices` reduction that `curvature_reg_matrix_reduced` applies for the +log-det. + +Those excluded rows are, by the Delaunay mesh's own docstring, the +"poorly constrained boundary vertices" whose zeroing exists to "stabilize the +linear inversion" and "prevent poorly constrained boundary vertices from absorbing +flux". **The noise map re-admits into an explicit inverse precisely the degenerate +rows the reconstruction deliberately dropped to stay stable.** + +The user-visible result: a pixel whose reconstruction reads exactly `0` (meaning +"not solved for") gets a noise value computed as though it had been solved. +Reconstruction and noise map disagree about what those pixels mean. + +Scope: only bites when `zeroed_pixels > 0`. `Delaunay.__init__` defaults it to `0`, +so this is opt-in per mesh — check `rectangular_rtu_adapt_density` and any +workspace configs before sizing the blast radius. + +## Defect 3 — `use_edge_zeroed_pixels` is silently ignored unless the positive-only solver is on + +This is the "does the edge-pixel handling make sense next to positive-only?" +question, and the answer is no. The control flow (`abstract.py:509-554`): + +```python +if self.settings.use_positive_only_solver: # default True + if self.settings.use_edge_zeroed_pixels and self.has(cls=Mapper): + ...subset, fnnls, scatter back... + else: + return reconstruction_positive_only_from(FULL matrix) +return reconstruction_positive_negative_from(FULL matrix) # edge-zeroing never consulted +``` + +`use_edge_zeroed_pixels` is nested **inside** the positive-only branch. Setting +`use_positive_only_solver: false` — a reasonable thing to do, for speed or to +permit negative values — **silently disables edge-zeroing too**, with no warning. +The poorly-constrained boundary vertices come straight back into the solve and +results change for a reason the config does not express. + +These are orthogonal concerns. Which parameters are *solvable* (edge-zeroing) is a +statement about the mesh; which solver walks them is a separate choice. Edge-zeroing +should apply to both branches, or the coupling should be made explicit and +documented. + +## Suggested direction (interpretation now settled; design still open) + +If the noise map is to describe the estimator actually used, the covariance should +be formed on the **same index set the reconstruction solved**, and scattered back: + +1. Determine the kept set exactly as `reconstruction` does — `zeroed_ids_to_keep` + under edge-zeroing, and, for the NNLS answer to Defect 1, further restricted to + the free set (pixels with `reconstruction > 0`). +2. Cholesky-invert that submatrix (per the sibling prompt's `cho_factor` / + `cho_solve` fix). +3. Scatter back into full shape. **Decide what the excluded pixels report** — `0` + matches the reconstruction's own convention and keeps plots working; `NaN` is + more honest ("never estimated") but breaks colourbars and the CSV. Recommend + `0` with an explicit docstring statement, since the reconstruction already + reports `0` there and consumers handle it. + +Restricting to the NNLS free set gives the covariance *conditional on the active +set* — standard practice for constrained least squares, and a defensible, +documentable choice. It is still an approximation: it ignores the uncertainty in +the active set itself. Say so in the docstring rather than implying exactness. + +## Downstream evidence for why this matters (added 2026-08-22) + +`autolens_workspace` uses the 1D noise map directly in user-facing science scripts — +`scripts/{imaging,interferometer,group,multi_galaxy}/features/pixelization/source_science.py` +compute: + +```python +reconstruction_noise_map = inversion.reconstruction_noise_map +signal_to_noise_map = reconstruction / reconstruction_noise_map +``` + +So the quantity this prompt argues is computed for the wrong estimator is divided into the +reconstruction to produce a **signal-to-noise map on a source reconstruction** — the number +that ends up in papers. If the NNLS/unconstrained mismatch is real, it propagates straight +into published S/N. That raises the stakes on the Defect 1 decision and is the concrete +reason to instrument a real fit rather than reason about it further. + +## SUPERSEDED — synthetic proxy, badly underestimated the effect (kept as a caution) + +This prompt's load-bearing claim was that NNLS pins a large fraction of a compact source's mesh, +flagged as reasoned-not-measured. Now measured with PyAutoArray's **real solver** +(`autoarray.util.fnnls.fnnls_cholesky`) on a problem shaped like a source-plane inversion — compact +Gaussian source, 4-neighbour gradient regularization, mesh mostly covering empty sky. + +**CONFIRMED — the pinned fraction is large, and scales with source compactness:** + +| mesh | pixels | source area | pinned at 0 | +|---|--:|--:|--:| +| 20x20 | 400 | 2% | **46.0%** | +| 20x20 | 400 | 5% | **24.8%** | +| 20x20 | 400 | 15% | 1.2% | +| 30x30 | 900 | 2% | **47.8%** | +| 30x30 | 900 | 5% | **26.4%** | +| 30x30 | 900 | 15% | 0.4% | +| 40x40 | 1600 | 2% | **50.0%** | +| 40x40 | 1600 | 5% | **26.8%** | +| 40x40 | 1600 | 15% | 0.5% | + +A quarter to a half of the mesh is pinned for a genuinely compact source; the effect vanishes for an +extended one. The premise holds. + +**But the numerical consequence is far smaller than this prompt implied.** Shipped full-matrix noise +map vs the active-set-conditional one (covariance restricted to the free set), on the 900-pixel / +5%-source case (226 pinned): + +| quantity | value | +|---|---| +| median noise, shipped (full matrix) | 0.002226 | +| median noise, active-set-conditional | 0.002159 | +| ratio shipped / conditional, free pixels | **median 1.025, min 1.002, max 1.123** | + +The shipped map **overstates** uncertainty on the free pixels by ~2.5% median, up to ~12% worst case. +The bias is systematic and one-directional (`min 1.002` — it never understates), exactly as the +truncation argument predicts. But it is a few percent, **not** the order-of-magnitude error the +`Priority: high` grading assumed. + +The pinned pixels are less alarming than feared: their reconstruction is exactly `0.0`, so +`signal_to_noise_map = reconstruction / reconstruction_noise_map` yields `0 / 0.00228 = 0`. Zero S/N +for an unlit pixel is defensible — though the reported noise value itself is still meaningless, the +posterior there being a spike at the boundary rather than a Gaussian. + +**Re-graded `Priority: high` -> `medium`.** A systematic, always-one-direction ~2.5% (up to 12%) +overstatement of source-plane error bars is worth correcting for a quantity now confirmed to be a +**posterior** and which feeds published S/N maps. It is not an emergency. + +**Caveat, and a real one.** This is a *structural proxy*, not a lens fit: the mapping matrix is random +rather than produced by ray tracing, so neighbouring image pixels do not map to neighbouring source +pixels as they do in reality. That structure affects conditioning and could move the magnitude either +way. The pinned *fraction* is robust to it (it follows from source compactness, not mapping +structure); the *2.5% / 12%* figures are indicative only. **Re-measure on a real Delaunay fit before +quoting them anywhere.** + +Defect 3 (the config coupling) is untouched by all of this — it is an unambiguous bug at any priority, +but it sits on the reconstruction path and so changes fit results. It needs its own sign-off. + +## MEASURED ON A REAL LENS FIT 2026-08-22 — the proxy was wrong; re-graded back to high + +The synthetic measurement above was re-run on a **real ray-traced fit**, reproducing +`autolens_workspace/scripts/imaging/features/pixelization/source_science.py`: Isothermal +`einstein_radius=1.6` + shear, `RectangularBilinearAdaptDensity(28, 28)`, `Constant` regularization, +`r=3.0"` mask, `over_sample_size_pixelization=4`, compact Sersic source, PSF and Poisson noise. +Defaults confirmed live: `use_positive_only_solver=True`, `use_edge_zeroed_pixels=True`, +`mesh.zeroed_pixels = 108` of 784. + +**The proxy understated the effect by an order of magnitude, and its caveat was the reason.** A random +mapping matrix spreads data support across the whole mesh; real ray tracing concentrates it in the +arc, so far more of the mesh is unconstrained. + +Sensitivity sweep (`noise x med` = median shipped/active-set-conditional on free pixels; `flux %` = +change in source flux passing the workspace's `S/N >= 5` cut): + +| r_eff | reg coeff | mesh | params | pinned % | noise x med | max | flux % | +|--:|--:|--:|--:|--:|--:|--:|--:| +| 0.05 | 1.0 | 28 | 784 | **97.1%** | **2.837** | 9.34 | **-24.1%** | +| 0.10 | 1.0 | 28 | 784 | **88.6%** | **1.692** | 10.60 | -11.1% | +| 0.30 | 1.0 | 28 | 784 | 53.2% | 1.208 | 8.02 | -3.9% | +| 0.60 | 1.0 | 28 | 784 | 21.0% | 1.025 | 4.06 | -1.1% | +| 0.10 | 0.1 | 28 | 784 | 87.1% | 1.713 | 10.39 | **-29.8%** | +| 0.10 | 10.0 | 28 | 784 | 82.3% | 1.117 | 1.90 | -2.4% | +| 0.10 | 100.0 | 28 | 784 | 83.2% | 1.009 | 1.57 | -0.2% | +| 0.10 | 1.0 | 20 | 400 | 91.0% | 1.481 | 4.93 | 0.0% | +| 0.10 | 1.0 | 40 | 1600 | 90.4% | 1.421 | 10.36 | **-48.9%** | + +Two clean trends: the gap **grows with source compactness** and **shrinks with regularization +strength**. At `coeff=100` it essentially vanishes (1.009); at realistic `coeff~1` with a compact +source it is a factor of **1.7-2.8**, with individual pixels up to **10x**. + +### The framing that matters — these are two bounds, not right-vs-wrong + +Do not read "1.8x" as "the shipped noise map is wrong by 1.8x". The two quantities **bracket** the +truth: + +- **Shipped (full-matrix)** ignores the `s >= 0` constraint entirely, so it **overstates** — an upper + bound. +- **Active-set-conditional** treats the active set as *known*, ignoring uncertainty about which + pixels are pinned, so it **understates** — a lower bound. +- The true truncated-Gaussian posterior lies between them. + +So the finding is: **the source-plane noise map is ambiguous at the factor-of-2 level for compact +sources**, and the shipped value sits at the pessimistic end of that range. That is still a serious +problem for a published error bar — but "switch to the conditional covariance" is **not** the +correct fix, it just swaps one bound for the other. Designing the real fix means either computing the +truncated posterior properly, or documenting the bracket honestly. + +### Downstream consequence is a threshold effect, which amplifies it + +`source_science.py` does not merely display the noise map: + +```python +signal_to_noise_map = reconstruction / reconstruction_noise_map +mesh_pixel_mask = signal_to_noise_map < 5.0 +reconstruction_masked[mesh_pixel_mask] = 0.0 +``` + +Source flux and magnification are computed from the surviving pixels. A hard cut turns a smooth noise +bias into a discrete one: pixels near `S/N = 5` flip in or out. Measured flux shifts of **-24%**, +**-30%** and **-49%** in the table above come from a handful of pixels crossing that line. + +**Re-graded `Priority: medium` -> `high`.** The earlier downgrade was made on the synthetic proxy and +was wrong. + +### What is still not established + +- One lens configuration, one mass model, one noise realization per row. No error bars on these numbers. +- `reg coeff` was fixed by hand. In a real model-fit λ is a **free parameter the sampler optimises**, + and the `coeff=100` row shows the effect nearly vanishes when λ is large. Where fitted λ actually + lands for these datasets decides whether this bites in practice — **measure that before acting**. +- `RectangularBilinearAdaptDensity`, not Delaunay. Delaunay is the other common source mesh and was + not tested. + +Scripts: `scratchpad/real_fit_measure.py`, `scratchpad/sensitivity.py` (session artefacts). + +## SETTLED 2026-08-22 — at the FITTED lambda the effect nearly vanishes. Re-graded to low. + +The open question from the sweep above was: those numbers used a hand-set regularization +coefficient, but in a real fit `Constant.coefficient` is a free parameter under +`LogUniform(1e-6, 1e6)` (workspace `config/priors/regularization/constant.yaml`), and the figure of +merit for a pixelized fit is the **Bayesian log evidence**. So `argmax_lambda log_evidence` is what +the sampler converges on — computed deterministically over a 17-point grid rather than by running +Nautilus. + +| r_eff | lambda* | log evidence | pinned % | noise x med | max | S/N>=5 ship | cond | flux % | +|--:|--:|--:|--:|--:|--:|--:|--:|--:| +| 0.05 | **10** | 27605.4 | 96.6% | **1.263** | 1.91 | 12 | 12 | **0.0%** | +| 0.10 | **10** | 26932.8 | 87.1% | **1.055** | 2.45 | 29 | 30 | **-1.3%** | +| 0.30 | **10** | 24635.2 | 42.3% | **1.007** | 2.10 | 138 | 138 | **0.0%** | + +`lambda* = 10` in all three cases, comfortably inside the scanned grid (`1e-3 .. 1e5`) — not an edge +artefact. + +**The Bayesian evidence self-selects away from the problematic regime.** The factor-of-2.8 gap found +above occurs at `lambda ~ 0.1-1`, i.e. under-regularized solutions the evidence *penalises*. At the +lambda a real fit chooses, the median gap is **1.007-1.263** and the downstream science outputs — +source flux and magnification through the `S/N >= 5` cut — move by **0.0%, -1.3%, 0.0%**. The pixel +counts either side of the cut are all but identical (12/12, 29/30, 138/138). + +The pinned fraction stays large (42-97%), so the *mechanism* in Defect 1 is real and confirmed. It +simply does not have a large numerical consequence at the operating point. + +**Re-graded `Priority: high` -> `low`.** Honest accounting: this prompt has been graded high -> +medium -> high -> low across three measurements. The swings came from measuring progressively less +wrong things — synthetic proxy, then a real fit at hand-set lambda, then a real fit at *fitted* +lambda. Only the last is the operating point, and it is the one that governs. + +### What remains true and worth doing + +- For a **very compact source** (`r_eff = 0.05`) the shipped noise map still overstates by **~26% + median**, up to ~2x on individual pixels. That is a real bias on a published error bar even though + it moves no flux. Anyone quoting per-pixel source uncertainties on a compact source should know. +- The **bracket framing stands**: full-matrix overstates, active-set-conditional understates, and the + truncated-Gaussian posterior lies between. A fix that simply swaps to the conditional covariance + would be wrong at any lambda. +- **Recommended immediate action is documentation, not code**: state in + `reconstruction_noise_map`'s docstring that the covariance is that of the unconstrained solve, + that the default solver is NNLS, and that for compact sources this overstates per-pixel noise by a + few tens of percent at most. Cheap, honest, no API change. +- A proper truncated-posterior implementation is only worth it if someone needs calibrated per-pixel + error bars on very compact sources. Not now. + +### Caveats on this result + +- **The lens mass was fixed at truth.** In a real model-fit the mass is free too, and a poor mass + model may need a lower lambda to absorb residuals — which is the regime where the gap opens. This + is the most likely way the conclusion could be wrong. +- One noise realization per row; no error bars on lambda*. +- `RectangularBilinearAdaptDensity` only; Delaunay untested. +- Evidence was maximised on a grid, not sampled — Nautilus explores a posterior over lambda, so some + posterior mass sits at lower lambda where the gap is larger. + +Script: `scratchpad/fitted_lambda.py` (session artefact). + +## PARTIALLY ADDRESSED 2026-08-22 — the docstring caveat shipped + +The recommended immediate action was taken: **[PyAutoArray#472](https://github.com/PyAutoLabs/PyAutoArray/pull/472)** +documents on `reconstruction_noise_map` that the covariance is that of the *unconstrained* solve +while the default solver is NNLS, quantifies the overstatement at the evidence-optimal coefficient +(x1.01 to x1.26 median, up to ~2x per pixel) and below it (~2.8x median, ~10x per pixel), and records +that restricting to the free set is **not** the correction because the two bracket the truth. + +Documentation only — no behaviour change, no API change. + +### What this prompt still owns + +1. **Defect 1's actual maths.** Computing the truncated-Gaussian posterior properly. Graded `low`: + worth it only if someone needs calibrated per-pixel error bars on a very compact source. +2. **Defect 2** — the covariance ignores `zeroed_ids_to_keep` while the reconstruction subsets by it. + Never measured in isolation; the real fits here had 108 zeroed pixels of 784 and they are folded + into the "pinned" counts throughout, so its separate contribution is unknown. +3. **Defect 3** — `use_edge_zeroed_pixels` nested inside the `use_positive_only_solver` branch, so + turning the positive-only solver off silently disables edge-zeroing. **Untouched by any of the + measurement above and unambiguous at any priority.** It sits on the reconstruction path, so it + changes fit results and needs its own sign-off. This is the most likely next piece of real work + here. +4. **The open measurement**: re-run the evidence-optimal lambda with the lens mass **free** rather + than fixed at truth. That is the single result most likely to overturn the `low` grading. + +## Verification + +- **Reproduce the symptom first.** Take a real Delaunay source fit, compute the + noise map under the current code and under the free-set-restricted covariance, + and compare. Quantify how many mesh pixels are pinned at zero by NNLS — the + claim that this fraction is large for compact sources is **reasoned, not + measured**, and the whole prompt rests on it. If the fraction turns out small, + Defect 1 is a much smaller problem than stated here and should be re-graded. +- Confirm reconstruction and noise map agree on which pixels were solved: every + pixel the reconstruction reports as an exact structural zero should be + identifiable in the noise map by the documented convention. +- With `zeroed_pixels > 0`, assert the covariance is formed on the reduced matrix + — regression-test the shape and the scatter-back, not just values. +- For Defect 3, assert `use_edge_zeroed_pixels: true` + `use_positive_only_solver: + false` either applies edge-zeroing or raises/warns. It must not silently ignore + the setting. +- Check whether `curvature_reg_matrix_reduced`'s `mapper_indices` reduction should + apply to the noise map too. The log-det uses it; the noise map does not. Decide + deliberately — this is a third, separate index set and the inconsistency between + all three is itself a finding. + +## Prior art — read before starting + +- `complete/2026/08/numerical-inversion-failures.md` — this cluster's refutation. +- `complete/2026/07/pix-inversion-not-positive-definite.md` — an earlier + non-positive-definite hypothesis, also refuted; documents the `GaussianKernel` + PD-guarantee `f1817af0`. +- `autoarray/util/cholesky_funcs.py:50-80` — near-coincident mesh vertices make the + Schur pivot's sign depend on BLAS thread count. The degeneracy is real and + documented; this prompt is about not feeding it into an explicit inverse. +- `abstract.py:805` — the repo already documents `~1e-6` evidence round-off from + "factorizing the explicitly formed inverse" at `cond(C) ~ 1e9` on clustered + traced mesh vertices. + +## Provenance + +- Found during: research for the sibling prompt, 2026-08-22. +- Do the sibling first — it is small, needs no science decision, and its Cholesky + covariance helper is the building block this prompt reuses. diff --git a/draft/bug/health_fixes/README.md b/draft/bug/health_fixes/README.md index 7464ca8c..0acd0f92 100644 --- a/draft/bug/health_fixes/README.md +++ b/draft/bug/health_fixes/README.md @@ -18,7 +18,7 @@ Each failing script is assigned to exactly one prompt: | ~~aggregator_output_contracts.md~~ — ✅ **SHIPPED 2026-07-07**, record `complete/2026/07/aggregator-output-contracts.md` (PyAutoFit#1324; autogalaxy_workspace#122, autolens_workspace#229, autolens_workspace_test#146 all merged) | 7 | Result/aggregator prerequisites and generated paths | | [jax_runtime_and_parity.md](jax_runtime_and_parity.md) — ⚠️ **6/6 pass 2026-08-21**, defect refuted; parkings NOT cleared (intermittent) | 6 | JAX/TFP compatibility and likelihood parity | | [jit_visualization_outputs.md](jit_visualization_outputs.md) — ⚠️ **4/4 pass 2026-08-21**, refuted; point_source parking is stale | 4 | Quick-update visualizations not producing images | -| [numerical_inversion_failures.md](numerical_inversion_failures.md) — ⚠️ **2/2 pass 2026-08-21**, refuted; incidental PyAutoArray sqrt-NaN found | 2 | Non-positive-definite inversion matrices | +| ~~numerical_inversion_failures.md~~ — ✅ **CLOSED 2026-08-22**, record `complete/2026/08/numerical-inversion-failures.md` ([PyAutoArray#467](https://github.com/PyAutoLabs/PyAutoArray/issues/467); **0/2 reproduce** on current `main` — no defect, no code changed in any repo; neither script is parked, so re-validation is automatic in every `mode=release` pass). Incidental PyAutoArray sqrt-NaN filed as `draft/bug/autoarray/reconstruction_noise_map_covariance_sqrt.md` | 2 | Non-positive-definite inversion matrices | | [release_timeout_policy.md](release_timeout_policy.md) — ⚠️ **4/4 measured pass far under cap 2026-08-21**; start_here not measured | 5 | 300-second release-surface decisions | Total: **42 scripts**. Scripts that pass on current `main` remain listed because they @@ -69,3 +69,34 @@ consistent explanation across all of them is the one #1327 reached: stale cached candidate. 2. **These were source-tree runs**, not the TestPyPI wheels the release run installed. A wheel-only defect would not show here. + +## 2026-08-22 — `numerical_inversion_failures` closed (the "four" above are now three) + +Struck through in the table. The 2026-08-21 sweep left four prompts in `draft/`; one of them has +now closed, on the same ground that closed `autofit_sampler_database` rather than on its green run +alone: + +**Neither of its two scripts is parked.** Verified against `main`: +`autogalaxy_workspace/config/build/no_run.yaml` has no +`interferometer/features/pixelization/galaxy_reconstruction` entry, and +`autolens_workspace_test/config/build/no_run.yaml` has no `interferometer/model_fit` entry — it +appears only as a *consumer*, where its simulator is marked `BOOTSTRAP-TARGET` for producing +`model_fit`'s dataset. Both scripts therefore re-execute in **every** `mode=release` pass, so +re-validation is automatic and there is no human reminder to lose. + +That is the distinction that decides this whole folder: caveat 1 of the 2026-08-21 sweep blocks the +*parked* prompts from closing, and it simply does not apply here. The three that remain +(`jax_runtime_and_parity`, `jit_visualization_outputs`, `release_timeout_policy`) each still carry +SLOW/NEEDS_FIX parkings describing *intermittent* failures a single green run cannot clear. + +It also makes **four** independent refutations of this cluster, not three — the fourth predates the +release run: `complete/2026/07/pix-inversion-not-positive-definite.md` (2026-07-21) tested the same +non-positive-definite hypothesis across six markers, found all six stale, and changed **no code**. +The `LinAlgError` had been cured on 2026-04-10 by PyAutoArray's `GaussianKernel` PD-guarantee +`f1817af0`, confirmed by a 40-draw inversion A/B across the full prior range (0 raises, 0 non-finite). + +One real defect came out of the gate and is filed separately: +`draft/bug/autoarray/reconstruction_noise_map_covariance_sqrt.md`. **It looks exactly like evidence +for the non-positive-definite hypothesis and is not** — `abstract.py:859` applies `np.sqrt` +elementwise to a whole covariance matrix, so negative off-diagonals are NaN unconditionally, for any +matrix, however well-conditioned. It has been mistaken for a conditioning symptom once already. diff --git a/draft/bug/health_fixes/numerical_inversion_failures.md b/draft/bug/health_fixes/numerical_inversion_failures.md deleted file mode 100644 index b0e3de4d..00000000 --- a/draft/bug/health_fixes/numerical_inversion_failures.md +++ /dev/null @@ -1,79 +0,0 @@ -# Fix release-profile numerical inversion failures - -Type: bug -Target: health_fixes -Difficulty: too-large -Autonomy: supervised -Priority: high -Status: formalised - -## Context - -Two interferometer scripts fail in inversion paths with non-positive-definite matrices. -The Autolens test failure reproduces on current `main`; the Autogalaxy script passed in a -stateful local checkout and needs a clean confirmation. - -Owners: @PyAutoArray, @PyAutoGalaxy, @PyAutoLens, @autogalaxy_workspace, and -@autolens_workspace_test. - -## Scripts - -- `autogalaxy_workspace/scripts/interferometer/features/pixelization/galaxy_reconstruction.py` -- `autolens_workspace_test/scripts/interferometer/model_fit.py` - -## Required work - -1. Reproduce in clean output/worktrees with deterministic seeds and release settings. -2. Capture the curvature and regularization matrix properties at failure: symmetry, - conditioning, eigenvalue range, dtype, backend, and mapper configuration. -3. Identify whether the defect is invalid sampled parameters, regularization construction, - numerical stabilization, or a script model that permits an undefined inversion. -4. Fix the owning library for valid inputs. Do not catch `LinAlgError` or alter the script - to hide a genuine inversion failure. -5. Add numerical regression tests and rerun both scripts repeatedly under the profile. - - - -## 2026-08-21 — REPRODUCTION GATE RUN: **2/2 PASS — prompt refuted** - -Method (identical to the gate that closed the sibling `autofit_sampler_database`, PyAutoFit#1508): -every script run from a **cleared** `output/`, under its workspace's -`config/build/profile_release.yaml`, env resolved by `autohands.env_config.build_env_for_script` -at workspace CWD, 1800s `mode=release` cap. Libraries at `main`: PyAutoFit `248ca971f`, -PyAutoArray `b808a9b1`, PyAutoGalaxy `7e3856dd`, PyAutoLens `d8f6bb3df`, PyAutoNerves `f6d6d52`. -Three workspace checkouts were **behind `origin/main`** and were synced first. - -| Script | Result | Secs | -|---|---|--:| -| `autolens_workspace_test/scripts/interferometer/model_fit.py` | PASS | 78 | -| `autogalaxy_workspace/scripts/interferometer/features/pixelization/galaxy_reconstruction.py` | PASS | 70 | - -The prompt states the autolens leg "reproduces on current `main`". It does not. No -non-positive-definite failure, no `LinAlgError`, in either. - -**Note which numerical path each took.** `autolens_workspace_test`'s release profile *defaults* -`PYAUTO_DISABLE_JAX="1"`, and scripts opt back in with an in-file `ENV: jax` declaration. -`model_fit.py` has no such declaration, so it ran on **numpy** — release-faithful, but worth -knowing for a claim about inversion numerics, since JAX-on and JAX-off are different code paths. - -### Incidental finding — a real defect, but NOT this prompt's - -`galaxy_reconstruction.py` passes while emitting 4x -`RuntimeWarning: invalid value encountered in sqrt` from -`PyAutoArray/autoarray/inversion/inversion/abstract.py:859`: - -```python -def reconstruction_noise_map_with_covariance(self): - return np.sqrt(np.linalg.inv(self.curvature_reg_matrix)) -``` - -`sqrt` is applied **elementwise to the whole inverse matrix**, whose off-diagonal entries are -covariances and are generally negative — so those entries are NaN *by construction*, for any -matrix, however well-conditioned. - -**This is not evidence of a non-positive-definite matrix** and does not rescue the prompt's -hypothesis, despite looking exactly like it would. It is a separate defect: a property whose -docstring promises a matrix that "accounts for the covariance of the noise between pixels" returns -NaN wherever that covariance is negative. The 1D `reconstruction_noise_map` is unaffected — it -takes the diagonal, and `diag(sqrt(M)) == sqrt(diag(M))` — so the science path is correct; only -the covariance-aware consumer and the warning spam are hit. Worth its own PyAutoArray prompt.