Skip to content

Commit 105f6ea

Browse files
Jammy2211claude
andcommitted
perf: warm-start the positive-only NNLS solve from the previous evaluation's passive set (#498)
Phase 3a of the numba CPU likelihood speed restoration. fnnls_cholesky gains solver diagnostics (stats dict), accepts a mask or index warm start, reuses one Cholesky factorisation for the warm start instead of a dense solve plus a from-scratch refactorisation, and repairs an infeasible warm passive set before the outer loop (an all-True seed previously bypassed the loop and returned the clipped vector). A process-local FIFO memo (nnls_memo.py) seeds each solve from the previous evaluation's final passive set, keyed on the index space (mesh/data shapes + edge-zeroed subset); Settings.nnls_warm_start_memo, default on, AUTOARRAY_NNLS_WARM_START=0 disables; a failed seeded solve retries once from the dense-sign start. Delaunay-1250 fiducial, random-walk sequence: median active-set iterations 69.5 -> 7.0 (euclid), 32 -> 8 (hst); solve 0.54 -> 0.05 s / 0.25 -> 0.07 s; log-likelihood parity 3e-14. Solve time is never worse than the dense-sign start, including on uncorrelated i.i.d. jumps. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01XySQZP6Npj4adZMkYZ96Bo
1 parent 86e2944 commit 105f6ea

10 files changed

Lines changed: 763 additions & 23 deletions

File tree

‎autoarray/config/general.yaml‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@ inversion:
88
use_border_relocator: false # If True, by default a pixelization's border is used to relocate all pixels outside its border to the border.
99
nnls_jacobi_preconditioning: true # If True (default), the curvature matrix passed to jaxnnls.solve_nnls_primal is Jacobi-preconditioned (D Q D y = D q, x = D y). Fixes NaN backward-pass gradients on ill-conditioned Q and roughly halves forward solve time. Set False to restore the raw unpreconditioned solve.
1010
nnls_target_kappa: 1.0e-11 # Central-path relaxation parameter passed to jaxnnls.solve_nnls_primal. Larger values smooth the relaxed-KKT backward pass and prevent NaN gradients on ill-conditioned Q; smaller values tighten the primal solve. Verified finite gradients across all MGE/rectangular/delaunay pipelines (imaging + interferometer) with scale invariance over 5 orders of magnitude in noise. jaxnnls's own default (1e-3) is too aggressive for the backward pass.
11+
nnls_warm_start_memo: true # If True (default), the NumPy/numba positive-only (fnnls) solve warm-starts its active set from the previous likelihood evaluation's passive set, cutting active-set iterations on successive sampler evaluations. The NNLS optimum is unique so the reconstruction is unchanged. On by default as of PyAutoArray#498, measured on the euclid+hst Delaunay-1250 fiducial (9.9x / 4.0x fewer active-set iterations on successive evaluations, reconstruction unchanged). Set false, or AUTOARRAY_NNLS_WARM_START=0, to disable. JAX path unaffected.
1112
reconstruction_vmax_factor: 0.5 # Plots of an Inversion's reconstruction use the reconstructed data's bright value multiplied by this factor.
1213
log_det_method: cholesky # How the Bayesian-evidence log-determinant terms are computed. "cholesky" (default) is the historical 2*sum(log(diag(cholesky(M)))); "slogdet" uses logabsdet of slogdet(M), which is identical where M is positive-definite but finite (not NaN) where the Cholesky fails, for gradient-based searches (opt-in, non-default; does not change the default evidence). Under "slogdet" the kernel regularization schemes (Matern/Gaussian/Exponential) also compute the regularization log-det analytically from a Cholesky of their covariance instead of factorizing the formed inverse. See PyAutoArray#391.
1314
regularization_term_method: matmul # How the Bayesian-evidence regularization term s^T H s is computed. "matmul" (default) is the historical s @ (H @ s) against the explicitly formed regularization matrix; "cho_solve" evaluates coefficient * s^T C^-1 s for the kernel schemes (Matern/Gaussian/Exponential/MaternAdapt) via one Cholesky solve of their covariance C, avoiding the explicit inverse whose round-off is amplified by cond(C) (~1e9 on clustered traced mesh vertices). Opt-in, non-default; does not change the default evidence. Schemes with no such factorization fall back to the formed matrix.

‎autoarray/inversion/inversion/abstract.py‎

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
import copy
2+
import hashlib
23
import warnings
34

45
import numpy as np
@@ -534,6 +535,37 @@ def solve_ids_to_keep(self) -> Optional[np.ndarray]:
534535

535536
return self.zeroed_ids_to_keep
536537

538+
def _nnls_warm_start_fingerprint(self, ids_to_keep=None) -> Optional[str]:
539+
"""
540+
Identify the index space this inversion's positive-only solve works in, so the
541+
cross-evaluation passive-set memo (`autoarray.inversion.inversion.nnls_memo`) only
542+
reuses a seed across evaluations whose reconstruction entries mean the same thing.
543+
544+
Deliberately built from shapes, NOT from the curvature matrix or data vector: the
545+
memo exists to carry a passive set between *neighbouring* parameter points, whose
546+
matrices differ. What must not change is the mapping from index to source pixel --
547+
the mesh sizes and the parameter-vector length. Under edge zeroing the passive set
548+
lives in the subset index space, so `ids_to_keep` is part of the identity too.
549+
550+
Returns `None` on the JAX path, which does not run fnnls at all.
551+
"""
552+
if self._xp.__name__.startswith("jax"):
553+
return None
554+
555+
parts = [str(tuple(np.shape(self.data_vector)))]
556+
557+
for mapper in self.cls_list_from(cls=Mapper):
558+
parts.append(str(tuple(np.shape(mapper.source_plane_mesh_grid))))
559+
560+
if ids_to_keep is not None:
561+
parts.append(
562+
hashlib.sha256(
563+
np.ascontiguousarray(ids_to_keep, dtype=np.int64).tobytes()
564+
).hexdigest()
565+
)
566+
567+
return "|".join(parts)
568+
537569
@cached_property
538570
def reconstruction(self) -> np.ndarray:
539571
"""
@@ -576,6 +608,9 @@ def reconstruction(self) -> np.ndarray:
576608
curvature_reg_matrix=curvature_reg_matrix,
577609
settings=self.settings,
578610
xp=self._xp,
611+
fingerprint=self._nnls_warm_start_fingerprint(
612+
ids_to_keep=ids_to_keep
613+
),
579614
)
580615
)
581616

@@ -599,6 +634,7 @@ def reconstruction(self) -> np.ndarray:
599634
curvature_reg_matrix=self.curvature_reg_matrix,
600635
settings=self.settings,
601636
xp=self._xp,
637+
fingerprint=self._nnls_warm_start_fingerprint(),
602638
)
603639

604640
return inversion_util.reconstruction_positive_negative_from(

‎autoarray/inversion/inversion/inversion_util.py‎

Lines changed: 55 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -258,6 +258,7 @@ def reconstruction_positive_only_from(
258258
curvature_reg_matrix: np.ndarray,
259259
settings: Settings = None,
260260
xp=np,
261+
fingerprint=None,
261262
):
262263
"""
263264
Solve the linear system Eq.(2) (in terms of minimizing the quadratic value) of
@@ -296,6 +297,9 @@ def reconstruction_positive_only_from(
296297
settings
297298
Controls the settings of the inversion, for this function where the solution is checked to not be all
298299
the same values.\
300+
fingerprint
301+
Identifies the index space this solve's passive set lives in, enabling the cross-evaluation warm-start
302+
memo (`Settings.nnls_warm_start_memo`) on the NumPy path. `None` disables the memo for this call.
299303
300304
Returns
301305
-------
@@ -371,13 +375,60 @@ def reconstruction_positive_only_from(
371375
try:
372376

373377
from autoarray.util.fnnls import fnnls_cholesky
378+
from autoarray.inversion.inversion import nnls_memo
374379

375-
return fnnls_cholesky(
376-
curvature_reg_matrix,
377-
(data_vector).T,
378-
P_initial=np.linalg.solve(curvature_reg_matrix, data_vector) > 0,
380+
use_memo = (
381+
settings is not None
382+
and settings.nnls_warm_start_memo
383+
and nnls_memo.memo_enabled()
384+
and fingerprint is not None
385+
)
386+
387+
stats = {}
388+
389+
key = (
390+
nnls_memo.memo_key(n=data_vector.shape[0], fingerprint=fingerprint)
391+
if use_memo
392+
else None
393+
)
394+
395+
warm_start = (
396+
nnls_memo.passive_set_get(key=key, n=data_vector.shape[0])
397+
if use_memo
398+
else None
379399
)
380400

401+
if warm_start is not None:
402+
try:
403+
reconstruction = fnnls_cholesky(
404+
curvature_reg_matrix,
405+
(data_vector).T,
406+
P_initial=warm_start,
407+
stats=stats,
408+
)
409+
except (RuntimeError, np.linalg.LinAlgError, ValueError):
410+
# A seed from a previous evaluation is a guess about a
411+
# different matrix, so it can factorise badly where the
412+
# dense-sign start would not. That must cost one retry, never a
413+
# resample: fall back to exactly the un-memoized computation
414+
# before the InversionException path below is reached.
415+
reconstruction = None
416+
else:
417+
reconstruction = None
418+
419+
if reconstruction is None:
420+
reconstruction = fnnls_cholesky(
421+
curvature_reg_matrix,
422+
(data_vector).T,
423+
P_initial=np.linalg.solve(curvature_reg_matrix, data_vector) > 0,
424+
stats=stats,
425+
)
426+
427+
if use_memo:
428+
nnls_memo.passive_set_put(key=key, passive_set=stats["passive_set"])
429+
430+
return reconstruction
431+
381432
except (RuntimeError, np.linalg.LinAlgError, ValueError) as e:
382433
if is_test_mode():
383434
# See reconstruction_positive_negative_from: benign dummy in test mode,
Lines changed: 81 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,81 @@
1+
import os
2+
3+
import numpy as np
4+
from typing import Dict, Optional
5+
6+
# Cross-evaluation memo for the positive-only (fnnls) solve's passive set.
7+
#
8+
# The Bro & De Jong active set iteration is warm-started from a guess at which
9+
# reconstruction entries are non-zero. The production guess -- the sign of the
10+
# unconstrained dense solve -- gets ~150 of ~1560 entries wrong on a Delaunay
11+
# euclid fit, and each wrong entry costs an active-set iteration in a solve
12+
# that is ~70% of the whole numba CPU likelihood evaluation. Successive
13+
# sampler evaluations sit close together in parameter space and share nearly
14+
# all of their passive set, so seeding from the previous evaluation's FINAL
15+
# passive set is a much better guess.
16+
#
17+
# Entries are keyed by the index space the passive set lives in (mesh/data
18+
# shapes, plus the edge-zeroed subset when one is in use) -- never by the
19+
# matrix values, since the whole point is to hit across nearby parameter
20+
# points. A wrong seed cannot corrupt the answer: the NNLS optimum is unique
21+
# and the solver adds and removes entries until the KKT conditions hold, so a
22+
# stale entry costs iterations, not correctness. Forked pool workers inherit a
23+
# copy at fork and diverge from there, which is fine for the same reason.
24+
#
25+
# Disable with AUTOARRAY_NNLS_WARM_START=0.
26+
_nnls_passive_set_memo: Dict[str, np.ndarray] = {}
27+
28+
_NNLS_PASSIVE_SET_MEMO_MAX_ENTRIES = 8
29+
30+
31+
def memo_enabled() -> bool:
32+
"""
33+
Whether the passive-set memo is active in this process.
34+
"""
35+
return os.environ.get("AUTOARRAY_NNLS_WARM_START", "1") != "0"
36+
37+
38+
def memo_key(n: int, fingerprint) -> str:
39+
"""
40+
The memo key for a solve of size `n` whose index space is described by
41+
`fingerprint` (see `AbstractInversion._nnls_warm_start_fingerprint`).
42+
"""
43+
return f"{n}:{fingerprint}"
44+
45+
46+
def passive_set_get(key: str, n: int) -> Optional[np.ndarray]:
47+
"""
48+
The memoized passive set for `key`, or None on a miss.
49+
50+
An entry whose indices do not all fit a size-`n` solve is a miss, not a
51+
hit: `n` is already part of the key, so this only fires if a caller
52+
fingerprints two different index spaces identically, and a miss is always
53+
a safe outcome.
54+
"""
55+
passive_set = _nnls_passive_set_memo.get(key)
56+
57+
if passive_set is None:
58+
return None
59+
60+
if passive_set.size and passive_set.max() >= n:
61+
return None
62+
63+
return passive_set
64+
65+
66+
def passive_set_put(key: str, passive_set: np.ndarray) -> None:
67+
"""
68+
Store a solve's final passive set, evicting the oldest entry once the memo
69+
is full (FIFO; the memo tracks one inversion's recent history, not a
70+
working set worth ranking).
71+
"""
72+
stored = np.asarray(passive_set, dtype=int).copy()
73+
stored.setflags(write=False)
74+
75+
if (
76+
key not in _nnls_passive_set_memo
77+
and len(_nnls_passive_set_memo) >= _NNLS_PASSIVE_SET_MEMO_MAX_ENTRIES
78+
):
79+
_nnls_passive_set_memo.pop(next(iter(_nnls_passive_set_memo)))
80+
81+
_nnls_passive_set_memo[key] = stored

‎autoarray/settings.py‎

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@ def __init__(
1717
no_regularization_add_to_curvature_diag_value: float = None,
1818
nnls_solver_tol: Optional[float] = None,
1919
nnls_max_iter: Optional[int] = None,
20+
nnls_warm_start_memo: Optional[bool] = None,
2021
log_det_method: Optional[str] = None,
2122
regularization_term_method: Optional[str] = None,
2223
):
@@ -99,6 +100,16 @@ def __init__(
99100
jaxnnls's own hard-coded cap of 50 (production HST pixelization+MGE systems converge in ~19-21
100101
iterations). Under `vmap` the solve runs until the slowest lane in the batch converges, so this
101102
also caps the worst-case batched cost. Only the JAX (`xp=jnp`) path honors this.
103+
nnls_warm_start_memo
104+
Whether the NumPy / numba positive-only (NNLS) solve warm-starts its active-set iteration from the
105+
passive set of the *previous* likelihood evaluation, held in a small process-local memo
106+
(:mod:`autoarray.inversion.inversion.nnls_memo`). Successive sampler evaluations sit close together in
107+
parameter space and share most of their passive set, so the seeded start removes active-set iterations
108+
from the solve, which is ~70% of a numba CPU likelihood evaluation on production Delaunay meshes. The
109+
NNLS optimum is unique, so the reconstruction is unchanged. `None` (default) reads the packaged value
110+
(`true`); setting ``AUTOARRAY_NNLS_WARM_START=0`` is the process-wide kill-switch that disables the
111+
memo whatever this setting says. Only the NumPy (`xp=np`) fnnls path honors this; the JAX path
112+
ignores it.
102113
log_det_method
103114
Which computation is used for the two Bayesian-evidence log-determinant terms
104115
(``log_det_curvature_reg_matrix_term`` and ``log_det_regularization_matrix_term``). `None`
@@ -151,6 +162,7 @@ def __init__(
151162
self.use_mixed_precision = use_mixed_precision
152163
self.nnls_solver_tol = nnls_solver_tol
153164
self.nnls_max_iter = nnls_max_iter
165+
self._nnls_warm_start_memo = nnls_warm_start_memo
154166
self._use_positive_only_solver = use_positive_only_solver
155167
self._use_edge_zeroed_pixels = use_edge_zeroed_pixels
156168
self._use_border_relocator = use_border_relocator
@@ -201,6 +213,22 @@ def no_regularization_add_to_curvature_diag_value(self):
201213

202214
return self._no_regularization_add_to_curvature_diag_value
203215

216+
@property
217+
def nnls_warm_start_memo(self) -> bool:
218+
"""
219+
Whether the NumPy fnnls solve is warm-started from the previous evaluation's passive set.
220+
"""
221+
if self._nnls_warm_start_memo is None:
222+
try:
223+
return conf.instance["general"]["inversion"]["nnls_warm_start_memo"]
224+
except KeyError:
225+
# Workspaces ship their own general.yaml that shadows autoarray's
226+
# and will not carry this key, so in practice this fallback *is*
227+
# the production default: on, matching the packaged value.
228+
return True
229+
230+
return self._nnls_warm_start_memo
231+
204232
@property
205233
def log_det_method(self):
206234
if self._log_det_method is None:

0 commit comments

Comments
 (0)