fix: adaptive Richardson step in LensCalc NumPy Hessian (too coarse near compact deflectors)

[Jammy2211]

Overview

LensCalc._hessian_via_richardson differences the deflections at a hardcoded 0.01 arcsec step and discards the h/(h/2) pair's truncation estimate. Where a ray passes ~3.5e-4" from a compact isothermal centre (the PyAutoLens#480 fixture) the magnification comes out 102–122 % wrong with a flipped sign, while the same route agrees to 1e-8 at a smooth plane. The NumPy path is the default for point-source flux fits and PointSolver's magnification threshold. Reproduced digit-for-digit by PyAutoLens#715's strict xfail and printed live in autolens_workspace#517.

Plan

• Make the Richardson step adaptive per point: use the error estimate the pair already gives (|H(h/2) − H(h)|/3), halve the step for unconverged points only (reusing the previous half-step evaluation), stop at rtol=1e-6, atol=1e-8 or after 12 halvings, and warn loudly (never silently) for points that do not converge.
• Public hessian_from(grid, xp) signature unchanged; JAX path untouched; smooth fields exit after the first pair at today's cost.
• Regression tests: compact SIS near its centre vs the profile's analytic shear/convergence (fails on main first — the control); the existing diagonal-grid literals unchanged to 1e-10; the unconverged-point warning.
• Follow-up in PyAutoLens (separate PR, merged after this one): convert the strict xfail in test_multi_plane_cross_validation.py into a positive assertion.
• Workspace follow-up (filed, not done here): the guide's live "not fixed" warning section needs rewording once the fix ships in a release.

Detailed implementation plan

Work Classification

Library (PyAutoGalaxy primary; PyAutoLens test follow-up).

Affected Repositories

• PyAutoGalaxy (primary)
• PyAutoLens (test-only follow-up PR)

Branch Survey

Repository	Current Branch	Dirty?
./PyAutoGalaxy	main	clean
./PyAutoLens	main	clean

Suggested branch: feature/lenscalc-adaptive-hessian-step
Worktree root: ~/Code/PyAutoLabs-wt/lenscalc-adaptive-hessian-step/

Implementation Steps

1. autogalaxy/operate/lens_calc.py _hessian_via_richardson(grid, buffer=0.01, rtol=1e-6, atol=1e-8, max_halvings=12): compute H(h), H(h/2) via _hessian_via_finite_difference; R = (4H(h/2) − H(h))/3; E = |H(h/2) − H(h)|/3; converged where E <= atol + rtol·|R| on all four components; loop: for the unconverged subset h ← h/2, previous H(h/2) becomes H(h), one new FD evaluation on the subset; after max_halvings keep the last R and warnings.warn with the unconverged count and the largest relative estimate.
2. Update the hessian_from docstring (drop the "matches JAX to float64 precision" claim; describe the adaptive step and tolerances).
3. Tests in test_autogalaxy/operate/test_deflections.py: test__hessian_from__adaptive_step__compact_sis_near_centre (IsothermalSph R_E=0.2, points 3e-4"–1e-3" from centre, LensCalc.from_mass_obj, shear/convergence via Hessian vs analytic at rtol 1e-4; confirm it FAILS on main before the fix), test__hessian_from__adaptive_step__smooth_field_unchanged (existing literals at :100/:116 unchanged to 1e-10), test__hessian_from__unconverged_points_warn (point exactly on the SIS centre: warns, finite, no raise).
4. PyAutoLens (own branch, second PR after the Galaxy merge): test_autolens/lens/test_multi_plane_cross_validation.py:1126-1158 — replace the strict xfail with a positive assertion at rtol 1e-3 vs the ray-traced Jacobian.

Key Files

• autogalaxy/operate/lens_calc.py (:382 hessian_from, :417 _hessian_via_richardson, :460 _hessian_via_finite_difference)
• test_autogalaxy/operate/test_deflections.py
• PyAutoLens/test_autolens/lens/test_multi_plane_cross_validation.py (follow-up)
• callers (read-only): PyAutoLens/autolens/point/fit/abstract.py:131, autolens/point/solver/shape_solver.py

Original Prompt

Click to expand starting prompt

LensCalc NumPy Hessian step is too coarse for multi-plane tracers

Type: bug
Target: PyAutoGalaxy
Repos:

• PyAutoGalaxy
• PyAutoLens
  Themes:
• point-source
• jax-gradient
  Difficulty: large
  Autonomy: supervised
  Priority: high
  Status: formalised

Filed: 2026-08-27

LensCalc's NumPy Hessian uses a hardcoded finite-difference step that is too coarse for
multi-plane configurations, returning magnifications that are wrong by >100% with flipped signs.

Found on 2026-08-27 while cross-checking the fix for PyAutoLens#480. It is a separate, pre-existing
bug: #480's fix is accurate, and this was found by the control arm of that check.

LensCalc._hessian_via_richardson (autogalaxy/operate/lens_calc.py) evaluates the Hessian by
central finite differences at a hardcoded buffer=0.01 arcsec, Richardson-extrapolated at h and
h/2. That step is fixed regardless of the scale the deflection field actually varies on.

Measured, on the tracer from PyAutoLens#480 (lens z=0.5 Isothermal R_E=1.6; source z=1.0 with its
own Isothermal R_E=0.2; source z=2.0), magnification at the four image positions of the z=1.0
source, computed three ways:

last plane (z=2.0)
numpy Richardson FD -0.00694 -0.00221 0.00139 0.00246
jax exact autodiff 0.04508 0.01099 -0.08602 -0.01118
ray-traced Jacobian 0.04508 0.01099 -0.08602 -0.01101

JAX autodiff (float64) and a Jacobian derived independently from traced_grid_2d_list_from agree
with each other; the NumPy path disagrees by 122% and has the WRONG SIGN on three of four points.
The ray-traced values are stable across step sizes h=1e-4 to 1e-7, so this is not noise in the
cross-check.

The same three-way comparison at the intermediate plane (z=1.0) agrees to 1.7e-08. So the failure
is configuration-dependent, not general: the map to z=1.0 involves only the smooth main lens, while
the map to z=2.0 additionally passes the compact z=1.0 deflector (R_E=0.2), whose deflection field
varies on scales where a 0.01 arcsec step is far too coarse.

Why it matters: the NumPy path is the default. AbstractFitPoint.magnifications_at_positions uses
it, so point-source flux fits and source-plane chi-squareds on multi-plane models with a compact
intermediate deflector are exposed, as is PointSolver's magnification threshold. A sign flip on a
magnification is not a small error.

Scope to consider: scale the step to the local deflection scale rather than hardcoding it; or
error-estimate from the Richardson pair (the h vs h/2 difference already bounds the truncation
error and is currently discarded) and warn or refine when it is large; or make the JAX path
reachable from NumPy callers. A regression test should pin the multi-plane configuration above
against the ray-traced Jacobian, which is the independent oracle used to find this.

[Jammy2211]

Library PR Created: #593

PyAutoLens follow-up (xfail → regression assert, theta_for_beta tol 1e-12) is committed on the same-named branch and will be opened after this merges.

[Jammy2211]

PyAutoLens follow-up PR: PyAutoLabs/PyAutoLens#717

[Jammy2211]

Shipped — 2026-08-29

PRs

• fix: adaptive Richardson step in LensCalc NumPy Hessian (#591) #593 (MERGED) — the fix
• test: Richardson-step regression asserts; theta_for_beta tol 1e-12 (PyAutoGalaxy#591) PyAutoLens#717 (MERGED) — test follow-up

Summary

LensCalc._hessian_via_richardson is now adaptive per point: successive Richardson estimates from (h, h/2) and (h/2, h/4), halving on the unconverged subset only, converged at |ΔR| ≤ 1e-8 + 1e-7·|R|, a roundoff-floor stop that returns the best extrapolant when the change starts growing, at most 20 halvings, and exactly one UserWarning per call for genuinely singular points. Compact-SIS control: κ 1.1e-9 / shear 7.6e-9 vs analytic (was 54–100 % off on main); the #480 multi-plane fixture agrees with the ray-traced Jacobian to 3.4e-4 with no warning. Smooth fields exit after 3 FD evaluations (was 2). PyAutoLens: the strict xfail became a regression assert; theta_for_beta root tol relaxed 1e-14 → 1e-12 (scipy status-3 artefact at residual 1e-16). Galaxy 1152 passed, Lens 570 passed; CI green on all legs.

Session Notes

• rtol=1e-8 was tried and rejected: the multi-plane path's FD roundoff floor is ~3e-8 (at 14 halvings), so a tighter gate warned permanently and over-refined. The floor rule makes the result tolerance-independent at the floor.
• Remaining warnings in the Lens suite (1 of N points … 5.00e-01) come from fixtures with a point exactly on a profile centre — the designed loud case, not a regression.
• Workspace follow-up (gated on the release carrying fix: adaptive Richardson step in LensCalc NumPy Hessian (#591) #593): PyAutoMind/draft/docs/autolens_workspace/multi_plane_guide_richardson_warning_update.md — the guide's live "not fixed" warning section must be reworded once users can see the fix.