Work Classification

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

Affected Repositories

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

Branch Survey

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

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.