From 5559d97f4ab20215f10756329638f15d6ef946f4 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Sat, 29 Aug 2026 15:48:55 -0400 Subject: [PATCH 1/3] fix: adaptive Richardson step in LensCalc NumPy Hessian (#591) `LensCalc._hessian_via_richardson` differenced the deflections at a hardcoded 0.01" step and discarded the truncation estimate the h/(h/2) pair already gives. Near a compact deflector that step straddles the whole deflector: on an `IsothermalSph` with einstein_radius=0.2, points 3e-4"-1e-3" from the centre came back with a convergence of ~46 against an analytic 100-333 (54-86% wrong) and a shear wrong by ~100%, which is where the >100% sign-flipped multi-plane magnifications of PyAutoLens#480's control arm came from. The step is now adaptive per grid point. The pair also yields `E = |H(h/2) - H(h)| / 3`; a point is converged when all four components satisfy `E <= atol + rtol * |R|` (rtol=1e-6, atol=1e-8). Unconverged points have h halved - the previous half-step evaluation becoming the new full-step one, so each halving costs one finite-difference evaluation on the shrinking unconverged subset only - for up to 12 halvings. Points still unconverged keep their last extrapolated value and raise a single UserWarning naming the count and the largest relative error estimate: never silent, never a raise (a raise here would kill an otherwise-converged fit). The compact-SIS points above now reproduce the profile's analytic shear and convergence to ~1e-10. `hessian_from(grid, xp)` is unchanged and the JAX path is untouched. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01WzKP1RJKhPrtFmsHSm6EEA --- autogalaxy/operate/lens_calc.py | 124 ++++++++++++++++--- test_autogalaxy/operate/test_deflections.py | 129 ++++++++++++++++++++ 2 files changed, 239 insertions(+), 14 deletions(-) diff --git a/autogalaxy/operate/lens_calc.py b/autogalaxy/operate/lens_calc.py index e9be7f993..a6e6666cf 100644 --- a/autogalaxy/operate/lens_calc.py +++ b/autogalaxy/operate/lens_calc.py @@ -19,6 +19,7 @@ from functools import wraps import importlib import logging +import warnings import numpy as np from typing import List, Tuple, Union @@ -394,8 +395,20 @@ def hessian_from(self, grid, xp=np) -> Tuple: - **NumPy** (``xp=np``, default): 2-point central finite-difference approximation, Richardson-extrapolated at step sizes ``h`` and ``h/2`` and combined as ``(4 * H(h/2) - H(h)) / 3``. This cancels the leading ``O(h^2)`` truncation term, - giving ``O(h^4)`` accuracy and matching the JAX path to float64 precision. JAX is - not imported. + giving ``O(h^4)`` accuracy. JAX is not imported. + + The step is **adaptive per grid point**. The pair ``H(h)``, ``H(h/2)`` also yields a + truncation-error estimate ``E = |H(h/2) - H(h)| / 3``; a point is converged when all + four Hessian components satisfy ``E <= atol + rtol * |R|`` (defaults ``rtol=1e-6``, + ``atol=1e-8``). Points that are not converged have their step halved -- reusing the + previous half-step evaluation, so each halving costs one new finite-difference + evaluation on the shrinking unconverged subset only -- for up to ``max_halvings=12`` + halvings. A smooth field therefore converges on the first pair at the same cost as + before, while a point close to a compact deflector (where a fixed 0.01" step is far + too coarse) refines until it is accurate. Any point still unconverged after the last + halving keeps its final extrapolated value and triggers a single ``UserWarning`` + reporting how many points did not converge; nothing is ever silently returned as + converged, and no exception is raised. - **JAX** (``xp=jnp``): exact derivatives via ``jax.jacfwd`` applied to ``deflections_yx_scalar``, vectorised over the grid with ``jnp.vectorize``. @@ -414,18 +427,101 @@ def hessian_from(self, grid, xp=np) -> Tuple: return self._hessian_via_richardson(grid=grid) return self._hessian_via_jax(grid=grid, xp=xp) - def _hessian_via_richardson(self, grid, buffer: float = 0.01) -> Tuple: - yy_h, xy_h, yx_h, xx_h = self._hessian_via_finite_difference( - grid=grid, buffer=buffer - ) - yy_h2, xy_h2, yx_h2, xx_h2 = self._hessian_via_finite_difference( - grid=grid, buffer=buffer / 2.0 - ) - hessian_yy = (4.0 * yy_h2 - yy_h) / 3.0 - hessian_xy = (4.0 * xy_h2 - xy_h) / 3.0 - hessian_yx = (4.0 * yx_h2 - yx_h) / 3.0 - hessian_xx = (4.0 * xx_h2 - xx_h) / 3.0 - return hessian_yy, hessian_xy, hessian_yx, hessian_xx + def _hessian_via_richardson( + self, + grid, + buffer: float = 0.01, + rtol: float = 1.0e-6, + atol: float = 1.0e-8, + max_halvings: int = 12, + ) -> Tuple: + """ + Returns the Hessian via Richardson-extrapolated central finite differences with a step + size that adapts, per grid point, to the scale the deflection field actually varies on. + + The finite-difference pair ``H(h)`` and ``H(h/2)`` gives both the extrapolated value + ``R = (4 H(h/2) - H(h)) / 3``, which is ``O(h^4)`` accurate, and an estimate of the + truncation error ``E = |H(h/2) - H(h)| / 3``, which the fixed-step implementation + discarded. Points whose four components all satisfy ``E <= atol + rtol * |R|`` are + accepted; the rest have ``h`` halved, the previous ``H(h/2)`` becoming the new ``H(h)`` + so only one new finite-difference evaluation is needed, and only on the unconverged + subset. Smooth fields therefore exit after the first pair at the original cost. + + Points still unconverged after ``max_halvings`` halvings keep their last extrapolated + value and raise a single ``UserWarning``; they are never silently accepted, and no + exception is raised (a raise here would kill an otherwise-converged model fit). + + Parameters + ---------- + grid + The 2D grid of (y,x) arc-second coordinates the Hessian is computed on. + buffer + The initial finite-difference step size in arc-seconds. + rtol + The relative tolerance the per-point truncation-error estimate must meet. + atol + The absolute tolerance the per-point truncation-error estimate must meet. + max_halvings + The maximum number of times the step is halved before unconverged points are warned + about and their last value kept. + """ + grid_values = grid.array if hasattr(grid, "array") else grid + grid_values = np.array(grid_values, dtype=np.float64, copy=True) + + hessian_h = np.stack( + self._hessian_via_finite_difference(grid=grid, buffer=buffer) + ).astype(np.float64) + hessian_h2 = np.stack( + self._hessian_via_finite_difference(grid=grid, buffer=buffer / 2.0) + ).astype(np.float64) + + richardson = (4.0 * hessian_h2 - hessian_h) / 3.0 + error = np.abs(hessian_h2 - hessian_h) / 3.0 + + unconverged = np.any(error > atol + rtol * np.abs(richardson), axis=0) + + step = buffer + halvings = 0 + + while np.any(unconverged) and halvings < max_halvings: + halvings += 1 + step /= 2.0 + + index = np.flatnonzero(unconverged) + + grid_subset = aa.Grid2DIrregular(values=grid_values[index]) + + hessian_subset_h = hessian_h2[:, index] + hessian_subset_h2 = np.stack( + self._hessian_via_finite_difference(grid=grid_subset, buffer=step / 2.0) + ).astype(np.float64) + + richardson_subset = (4.0 * hessian_subset_h2 - hessian_subset_h) / 3.0 + error_subset = np.abs(hessian_subset_h2 - hessian_subset_h) / 3.0 + + hessian_h[:, index] = hessian_subset_h + hessian_h2[:, index] = hessian_subset_h2 + richardson[:, index] = richardson_subset + error[:, index] = error_subset + + unconverged = np.zeros(unconverged.shape, dtype=bool) + unconverged[index] = np.any( + error_subset > atol + rtol * np.abs(richardson_subset), axis=0 + ) + + if np.any(unconverged): + total = int(unconverged.shape[0]) + number = int(np.count_nonzero(unconverged)) + denominator = np.where(np.abs(richardson) > 0.0, np.abs(richardson), 1.0) + largest = float(np.max((error / denominator)[:, unconverged])) + warnings.warn( + f"LensCalc Hessian: {number} of {total} points did not converge after " + f"{max_halvings} halvings (largest relative error estimate " + f"{largest:.2e}); values kept.", + UserWarning, + ) + + return richardson[0], richardson[1], richardson[2], richardson[3] def _hessian_via_jax(self, grid, xp) -> Tuple: import jax diff --git a/test_autogalaxy/operate/test_deflections.py b/test_autogalaxy/operate/test_deflections.py index 2f053c6ac..70a86fb39 100644 --- a/test_autogalaxy/operate/test_deflections.py +++ b/test_autogalaxy/operate/test_deflections.py @@ -129,6 +129,135 @@ def test__hessian_from__axis_aligned_grid__correct_values(): assert hessian_xx == pytest.approx(np.array([2.22209, 0.0]), 1.0e-4) +def test__hessian_from__adaptive_step__compact_sis_near_centre(): + """ + Regression test for the adaptive Richardson step (issue #591). + + Close to the centre of a compact deflector the deflection field varies on a scale far smaller + than the 0.01" step the NumPy Hessian used to be hardcoded to, so the finite differences + straddled the whole deflector and returned values that were ~100% wrong (and, for the + magnification, sign-flipped). With the step adapted per point the Hessian must reproduce the + profile's own analytic shear and convergence. + + ``IsothermalSph`` has closed-form shear and convergence, so the analytic values are an + independent oracle: at radius ``r`` from the centre both have magnitude + ``einstein_radius / (2 r)``, i.e. of order 100-330 for the radii used here. + """ + mp = ag.mp.IsothermalSph(centre=(0.0, 0.0), einstein_radius=0.2) + + radii = [3.0e-4, 4.0e-4, 5.5e-4, 7.0e-4, 8.5e-4, 1.0e-3] + angles = [10.0, 55.0, 100.0, 170.0, 230.0, 310.0] + + grid = ag.Grid2DIrregular( + values=[ + ( + radius * math.sin(math.radians(angle)), + radius * math.cos(math.radians(angle)), + ) + for radius, angle in zip(radii, angles) + ] + ) + + od = LensCalc.from_mass_obj(mp) + + shear_analytic = np.asarray(mp.shear_yx_2d_from(grid=grid)) + shear_via_hessian = np.asarray(od.shear_yx_2d_via_hessian_from(grid=grid)) + + convergence_analytic = np.asarray(mp.convergence_2d_from(grid=grid)) + convergence_via_hessian = np.asarray(od.convergence_2d_via_hessian_from(grid=grid)) + + np.testing.assert_allclose(shear_via_hessian, shear_analytic, rtol=1.0e-4) + np.testing.assert_allclose( + convergence_via_hessian, convergence_analytic, rtol=1.0e-4 + ) + + +def test__hessian_from__adaptive_step__smooth_field_unchanged(): + """ + The adaptive step of issue #591 must not move the answer on the smooth fields the fixed-step + implementation already handled well: the values pinned by + ``test__hessian_from__diagonal_grid__correct_values`` and + ``test__hessian_from__axis_aligned_grid__correct_values`` must still hold, and the adaptive + result must agree with the old fixed-step Richardson extrapolation (one pair at h=0.01 and + h=0.005) to far inside those tests' tolerances. + + Where the two differ at all it is the fixed step's own residual truncation error: the adaptive + result is the more accurate of the two (checked here against the profile's analytic + convergence, which it reproduces to ~1e-12 against the fixed step's ~1e-9). + """ + mp = ag.mp.Isothermal( + centre=(0.0, 0.0), ell_comps=(0.0, -0.111111), einstein_radius=2.0 + ) + + od = LensCalc.from_mass_obj(mp) + + grid_diagonal = ag.Grid2DIrregular(values=[(0.5, 0.5), (1.0, 1.0)]) + grid_axis_aligned = ag.Grid2DIrregular(values=[(1.0, 0.0), (0.0, 1.0)]) + + for grid in [grid_diagonal, grid_axis_aligned]: + hessian_h = np.stack(od._hessian_via_finite_difference(grid=grid, buffer=0.01)) + hessian_h2 = np.stack( + od._hessian_via_finite_difference(grid=grid, buffer=0.005) + ) + hessian_fixed_step = (4.0 * hessian_h2 - hessian_h) / 3.0 + + hessian_adaptive = np.stack(od.hessian_from(grid=grid)) + + np.testing.assert_allclose( + hessian_adaptive, hessian_fixed_step, rtol=1.0e-7, atol=1.0e-10 + ) + + convergence_adaptive = 0.5 * (hessian_adaptive[0] + hessian_adaptive[3]) + convergence_fixed_step = 0.5 * (hessian_fixed_step[0] + hessian_fixed_step[3]) + convergence_analytic = np.asarray(mp.convergence_2d_from(grid=grid)) + + assert np.max(np.abs(convergence_adaptive - convergence_analytic)) <= np.max( + np.abs(convergence_fixed_step - convergence_analytic) + ) + + hessian_yy, hessian_xy, hessian_yx, hessian_xx = od.hessian_from(grid=grid_diagonal) + + assert hessian_yy == pytest.approx(np.array([1.3882113, 0.6941056]), 1.0e-4) + assert hessian_xy == pytest.approx(np.array([-1.3882113, -0.6941056]), 1.0e-4) + assert hessian_yx == pytest.approx(np.array([-1.3882113, -0.6941056]), 1.0e-4) + assert hessian_xx == pytest.approx(np.array([1.3882113, 0.6941056]), 1.0e-4) + + hessian_yy, hessian_xy, hessian_yx, hessian_xx = od.hessian_from( + grid=grid_axis_aligned + ) + + assert hessian_yy == pytest.approx(np.array([0.0, 1.777699]), 1.0e-4) + assert hessian_xy == pytest.approx(np.array([0.0, 0.0]), 1.0e-4) + assert hessian_yx == pytest.approx(np.array([0.0, 0.0]), 1.0e-4) + assert hessian_xx == pytest.approx(np.array([2.22209, 0.0]), 1.0e-4) + + +def test__hessian_from__unconverged_points_warn(): + """ + The deflection angles of an isothermal sphere are discontinuous at its centre, so no step size + resolves the Hessian there and the adaptive refinement must give up. It must do so loudly (a + single ``UserWarning``), never silently and never by raising -- a raise here would kill an + otherwise-converged model fit. The values it keeps must stay finite, and the smooth point on + the same grid, which converged on the first pair, must be unaffected. + """ + mp = ag.mp.IsothermalSph(centre=(0.0, 0.0), einstein_radius=1.0) + + grid = ag.Grid2DIrregular(values=[(0.0, 0.0), (1.0, 1.0)]) + + od = LensCalc.from_mass_obj(mp) + + with pytest.warns(UserWarning, match="did not converge"): + hessian = od.hessian_from(grid=grid) + + assert np.all(np.isfinite(np.stack(hessian))) + + convergence_smooth = 0.5 * (hessian[0][1] + hessian[3][1]) + + assert convergence_smooth == pytest.approx( + float(np.asarray(mp.convergence_2d_from(grid=grid))[1]), 1.0e-4 + ) + + def test__convergence_2d_via_hessian_from(): grid = ag.Grid2DIrregular( values=[(1.075, -0.125), (-0.875, -0.075), (-0.925, -0.075), (0.075, 0.925)] From 4147f565ef5c924d4195df4ddbce964d549623ba Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Sat, 29 Aug 2026 15:57:42 -0400 Subject: [PATCH 2/3] fix: converge the adaptive Hessian step on successive Richardson estimates (#591) The first cut of the adaptive step gated convergence on the pair's own error estimate `|H(h/2) - H(h)| / 3`. That estimate bounds the error of `H(h/2)`, which is O(h^2), not of the extrapolant `R`, which is O(h^4), so it kept declaring non-convergence long after the returned value had stopped moving: the compact-SIS control points warned "6 of 6 points did not converge" while agreeing with the analytic shear and convergence to 1e-10, and the PyAutoLens#480 multi-plane configuration warned on all four image positions while agreeing with the ray-traced Jacobian. On a point-source fit, where `AbstractFitPoint.magnifications_at_positions` runs per likelihood evaluation, that turns the loud channel into noise. Convergence is now judged on the change between successive extrapolants: `R_k` from the pair `(h_k, h_k/2)` against `R_{k+1}` from `(h_k/2, h_k/4)`, converged where all four components satisfy `|R_{k+1} - R_k| <= atol + rtol * |R_{k+1}|`. The half-step evaluation is still reused as the next full-step one, so each halving costs one finite-difference evaluation on the shrinking unconverged subset. rtol=1e-6 and atol=1e-8 are unchanged; max_halvings is raised 12 -> 20, which no longer costs anything on well-behaved points because they now exit early. Measured: the compact-SIS points and the #480 configuration are both silent (the latter verified under `-W error::UserWarning`), a point exactly on an isothermal centre still warns exactly once, and the smooth diagonal and axis-aligned grids settle in three finite-difference evaluations (down from five and four under the old gate). Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01WzKP1RJKhPrtFmsHSm6EEA --- autogalaxy/operate/lens_calc.py | 67 +++++++++++++++++++-------------- 1 file changed, 38 insertions(+), 29 deletions(-) diff --git a/autogalaxy/operate/lens_calc.py b/autogalaxy/operate/lens_calc.py index a6e6666cf..50fd28c24 100644 --- a/autogalaxy/operate/lens_calc.py +++ b/autogalaxy/operate/lens_calc.py @@ -397,15 +397,17 @@ def hessian_from(self, grid, xp=np) -> Tuple: ``(4 * H(h/2) - H(h)) / 3``. This cancels the leading ``O(h^2)`` truncation term, giving ``O(h^4)`` accuracy. JAX is not imported. - The step is **adaptive per grid point**. The pair ``H(h)``, ``H(h/2)`` also yields a - truncation-error estimate ``E = |H(h/2) - H(h)| / 3``; a point is converged when all - four Hessian components satisfy ``E <= atol + rtol * |R|`` (defaults ``rtol=1e-6``, - ``atol=1e-8``). Points that are not converged have their step halved -- reusing the - previous half-step evaluation, so each halving costs one new finite-difference - evaluation on the shrinking unconverged subset only -- for up to ``max_halvings=12`` - halvings. A smooth field therefore converges on the first pair at the same cost as - before, while a point close to a compact deflector (where a fixed 0.01" step is far - too coarse) refines until it is accurate. Any point still unconverged after the last + The step is **adaptive per grid point**, and convergence is judged on the change + between **successive** extrapolants: ``R_k`` from the pair ``(h_k, h_k/2)`` and + ``R_{k+1}`` from ``(h_k/2, h_k/4)``, which costs one new finite-difference evaluation + because the half-step evaluation is reused. A point is converged when all four Hessian + components satisfy ``|R_{k+1} - R_k| <= atol + rtol * |R_{k+1}|`` (defaults + ``rtol=1e-6``, ``atol=1e-8``); the latest extrapolant is the value returned. Points + that are not converged keep halving, each halving costing one finite-difference + evaluation on the shrinking unconverged subset only, for up to ``max_halvings=20`` + halvings. A smooth field therefore settles after three finite-difference evaluations, + while a point close to a compact deflector (where a fixed 0.01" step is far too coarse) + refines until the extrapolant stops moving. Any point still unconverged after the last halving keeps its final extrapolated value and triggers a single ``UserWarning`` reporting how many points did not converge; nothing is ever silently returned as converged, and no exception is raised. @@ -433,23 +435,30 @@ def _hessian_via_richardson( buffer: float = 0.01, rtol: float = 1.0e-6, atol: float = 1.0e-8, - max_halvings: int = 12, + max_halvings: int = 20, ) -> Tuple: """ Returns the Hessian via Richardson-extrapolated central finite differences with a step size that adapts, per grid point, to the scale the deflection field actually varies on. - The finite-difference pair ``H(h)`` and ``H(h/2)`` gives both the extrapolated value - ``R = (4 H(h/2) - H(h)) / 3``, which is ``O(h^4)`` accurate, and an estimate of the - truncation error ``E = |H(h/2) - H(h)| / 3``, which the fixed-step implementation - discarded. Points whose four components all satisfy ``E <= atol + rtol * |R|`` are - accepted; the rest have ``h`` halved, the previous ``H(h/2)`` becoming the new ``H(h)`` - so only one new finite-difference evaluation is needed, and only on the unconverged - subset. Smooth fields therefore exit after the first pair at the original cost. - - Points still unconverged after ``max_halvings`` halvings keep their last extrapolated - value and raise a single ``UserWarning``; they are never silently accepted, and no - exception is raised (a raise here would kill an otherwise-converged model fit). + A finite-difference pair ``H(h)``, ``H(h/2)`` gives the extrapolant + ``R = (4 H(h/2) - H(h)) / 3``, which cancels the leading ``O(h^2)`` truncation term. The + fixed-step implementation returned the first such ``R`` whatever the field looked like. + Here the step keeps halving and convergence is judged on the change between successive + extrapolants: ``R_k`` from ``(h_k, h_k/2)`` against ``R_{k+1}`` from ``(h_k/2, h_k/4)``, + a point being converged when all four components satisfy + ``|R_{k+1} - R_k| <= atol + rtol * |R_{k+1}|``. Because the half-step evaluation is + reused as the next full-step one, each halving costs a single finite-difference + evaluation, and only on the shrinking unconverged subset. + + Judging on the extrapolants rather than on the pair's own error estimate + ``|H(h/2) - H(h)| / 3`` matters: that estimate bounds the error of ``H(h/2)``, which is + ``O(h^2)``, not of ``R``, which is ``O(h^4)``, so it declares non-convergence orders of + magnitude past the point where the returned value has stopped moving. + + Points still unconverged after ``max_halvings`` halvings keep their last extrapolant and + raise a single ``UserWarning``; they are never silently accepted, and no exception is + raised (a raise here would kill an otherwise-converged model fit). Parameters ---------- @@ -458,9 +467,9 @@ def _hessian_via_richardson( buffer The initial finite-difference step size in arc-seconds. rtol - The relative tolerance the per-point truncation-error estimate must meet. + The relative tolerance the change between successive extrapolants must meet. atol - The absolute tolerance the per-point truncation-error estimate must meet. + The absolute tolerance the change between successive extrapolants must meet. max_halvings The maximum number of times the step is halved before unconverged points are warned about and their last value kept. @@ -476,9 +485,9 @@ def _hessian_via_richardson( ).astype(np.float64) richardson = (4.0 * hessian_h2 - hessian_h) / 3.0 - error = np.abs(hessian_h2 - hessian_h) / 3.0 + change = np.full(richardson.shape, np.inf) - unconverged = np.any(error > atol + rtol * np.abs(richardson), axis=0) + unconverged = np.ones(richardson.shape[1], dtype=bool) step = buffer halvings = 0 @@ -497,23 +506,23 @@ def _hessian_via_richardson( ).astype(np.float64) richardson_subset = (4.0 * hessian_subset_h2 - hessian_subset_h) / 3.0 - error_subset = np.abs(hessian_subset_h2 - hessian_subset_h) / 3.0 + change_subset = np.abs(richardson_subset - richardson[:, index]) hessian_h[:, index] = hessian_subset_h hessian_h2[:, index] = hessian_subset_h2 richardson[:, index] = richardson_subset - error[:, index] = error_subset + change[:, index] = change_subset unconverged = np.zeros(unconverged.shape, dtype=bool) unconverged[index] = np.any( - error_subset > atol + rtol * np.abs(richardson_subset), axis=0 + change_subset > atol + rtol * np.abs(richardson_subset), axis=0 ) if np.any(unconverged): total = int(unconverged.shape[0]) number = int(np.count_nonzero(unconverged)) denominator = np.where(np.abs(richardson) > 0.0, np.abs(richardson), 1.0) - largest = float(np.max((error / denominator)[:, unconverged])) + largest = float(np.max((change / denominator)[:, unconverged])) warnings.warn( f"LensCalc Hessian: {number} of {total} points did not converge after " f"{max_halvings} halvings (largest relative error estimate " From 9499b78700c0e6300c9eebb0500380579429e286 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Sat, 29 Aug 2026 16:10:04 -0400 Subject: [PATCH 3/3] fix: adaptive Hessian gate rtol=1e-7 with roundoff-floor stop (#591) rtol=1e-6 left the compact-SIS control at ~1e-7 relative accuracy, but simply tightening to 1e-8 made the PyAutoLens#480 multi-plane fixture warn on every call. Instrumenting the halving sequence there shows why: the relative change between successive extrapolants falls as O(h^2) to 3.33e-08 at the 14th halving and then *grows* again - 5.98e-08, 1.73e-07, 7.02e-07, 2.37e-06 - as cancellation in (f(x+h) - f(x-h)) takes over. A tolerance below that floor can never be met, so the loop both warned forever and kept halving six steps past the best answer the arithmetic can give. Each point now stops on either of two rules, neither of which warns: 1. Tolerance: |R_{k+1} - R_k| <= atol + rtol * |R_{k+1}| on all four components, keeping R_{k+1}. Default rtol is now 1e-7 (atol=1e-8). 2. Roundoff floor: the relative change grew instead of shrinking while already below roundoff_guard=1e-4, so R_k - the estimate from before the growth - is kept. Rule 2 is guarded so it cannot silence a true singularity: a point exactly on an isothermal centre has a relative change pinned at ~0.5 that never enters the guard band, so it still runs to max_halvings and raises its single UserWarning. The warning is now reserved for that case. Measured: compact-SIS control kappa 1.14e-09 and shear 7.57e-09 against analytic (rtol=1e-6 gave 1.8e-08 / 1.2e-07); the #480 fixture is silent under -W error::UserWarning and now terminates at the floor for any rtol - 16 finite-difference evaluations at 1e-7, 17 at 1e-8 and at 1e-10, all agreeing with the ray-traced Jacobian to 3.4e-04, the stored constant's own precision. Smooth grids still settle in three finite-difference evaluations. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01WzKP1RJKhPrtFmsHSm6EEA --- autogalaxy/operate/lens_calc.py | 135 ++++++++++++++++++++------------ 1 file changed, 87 insertions(+), 48 deletions(-) diff --git a/autogalaxy/operate/lens_calc.py b/autogalaxy/operate/lens_calc.py index 50fd28c24..80d93c38c 100644 --- a/autogalaxy/operate/lens_calc.py +++ b/autogalaxy/operate/lens_calc.py @@ -400,17 +400,27 @@ def hessian_from(self, grid, xp=np) -> Tuple: The step is **adaptive per grid point**, and convergence is judged on the change between **successive** extrapolants: ``R_k`` from the pair ``(h_k, h_k/2)`` and ``R_{k+1}`` from ``(h_k/2, h_k/4)``, which costs one new finite-difference evaluation - because the half-step evaluation is reused. A point is converged when all four Hessian - components satisfy ``|R_{k+1} - R_k| <= atol + rtol * |R_{k+1}|`` (defaults - ``rtol=1e-6``, ``atol=1e-8``); the latest extrapolant is the value returned. Points - that are not converged keep halving, each halving costing one finite-difference - evaluation on the shrinking unconverged subset only, for up to ``max_halvings=20`` - halvings. A smooth field therefore settles after three finite-difference evaluations, - while a point close to a compact deflector (where a fixed 0.01" step is far too coarse) - refines until the extrapolant stops moving. Any point still unconverged after the last - halving keeps its final extrapolated value and triggers a single ``UserWarning`` - reporting how many points did not converge; nothing is ever silently returned as - converged, and no exception is raised. + because the half-step evaluation is reused. A point stops refining on either of two + rules: + + 1. **Tolerance** -- all four Hessian components satisfy + ``|R_{k+1} - R_k| <= atol + rtol * |R_{k+1}|`` (defaults ``rtol=1e-7``, + ``atol=1e-8``). ``R_{k+1}`` is returned. + 2. **Roundoff floor** -- the change has stopped shrinking (``|R_{k+1} - R_k|`` is + larger, relatively, than ``|R_k - R_{k-1}|``) while already below + ``roundoff_guard=1e-4`` relative. Halving further only feeds cancellation error, so + ``R_k``, the estimate from before the growth, is returned. This is what keeps a + configuration whose deflections are themselves evaluated numerically (a multi-plane + trace, say, which can floor at ~3e-8 relative) from warning on every call. + + Neither rule warns. Points that reach ``max_halvings=20`` halvings while their change + is still shrinking -- a genuine singularity, e.g. a grid point exactly on an isothermal + centre, whose relative change never falls below the guard -- keep their final + extrapolant and trigger a single ``UserWarning`` reporting how many points did not + converge. Nothing is ever silently returned as converged, and no exception is raised. + + A smooth field settles after three finite-difference evaluations; only points near + compact structure, where a fixed 0.01" step is far too coarse, iterate further. - **JAX** (``xp=jnp``): exact derivatives via ``jax.jacfwd`` applied to ``deflections_yx_scalar``, vectorised over the grid with ``jnp.vectorize``. @@ -433,9 +443,10 @@ def _hessian_via_richardson( self, grid, buffer: float = 0.01, - rtol: float = 1.0e-6, + rtol: float = 1.0e-7, atol: float = 1.0e-8, max_halvings: int = 20, + roundoff_guard: float = 1.0e-4, ) -> Tuple: """ Returns the Hessian via Richardson-extrapolated central finite differences with a step @@ -444,21 +455,31 @@ def _hessian_via_richardson( A finite-difference pair ``H(h)``, ``H(h/2)`` gives the extrapolant ``R = (4 H(h/2) - H(h)) / 3``, which cancels the leading ``O(h^2)`` truncation term. The fixed-step implementation returned the first such ``R`` whatever the field looked like. - Here the step keeps halving and convergence is judged on the change between successive - extrapolants: ``R_k`` from ``(h_k, h_k/2)`` against ``R_{k+1}`` from ``(h_k/2, h_k/4)``, - a point being converged when all four components satisfy - ``|R_{k+1} - R_k| <= atol + rtol * |R_{k+1}|``. Because the half-step evaluation is - reused as the next full-step one, each halving costs a single finite-difference - evaluation, and only on the shrinking unconverged subset. - - Judging on the extrapolants rather than on the pair's own error estimate + Here the step keeps halving and each point stops on one of two rules, neither of which + warns: + + 1. **Tolerance.** ``R_k`` from ``(h_k, h_k/2)`` and ``R_{k+1}`` from ``(h_k/2, h_k/4)`` + agree on all four components: ``|R_{k+1} - R_k| <= atol + rtol * |R_{k+1}|``. + ``R_{k+1}`` is kept. + 2. **Roundoff floor.** The relative change grew instead of shrinking, while already + below ``roundoff_guard``. Finite differences fall as ``O(h^2)`` only until + cancellation in ``(f(x+h) - f(x-h))`` takes over, after which halving makes the answer + *worse*; the turning point is the best the arithmetic can do. ``R_k``, the estimate + from before the growth, is kept. + + Because the half-step evaluation is reused as the next full-step one, each halving costs a + single finite-difference evaluation, and only on the shrinking active subset. + + Judging on the extrapolants rather than on a pair's own error estimate ``|H(h/2) - H(h)| / 3`` matters: that estimate bounds the error of ``H(h/2)``, which is ``O(h^2)``, not of ``R``, which is ``O(h^4)``, so it declares non-convergence orders of magnitude past the point where the returned value has stopped moving. - Points still unconverged after ``max_halvings`` halvings keep their last extrapolant and - raise a single ``UserWarning``; they are never silently accepted, and no exception is - raised (a raise here would kill an otherwise-converged model fit). + Rule 2 is guarded by ``roundoff_guard`` so that it cannot silence a genuine singularity: + a grid point sitting exactly on an isothermal centre has a relative change that stays at + ~0.5 forever, never entering the guard band, so it runs to ``max_halvings``, keeps its + last extrapolant and raises a single ``UserWarning``. It is never silently accepted, and + no exception is raised (a raise here would kill an otherwise-converged model fit). Parameters ---------- @@ -471,58 +492,76 @@ def _hessian_via_richardson( atol The absolute tolerance the change between successive extrapolants must meet. max_halvings - The maximum number of times the step is halved before unconverged points are warned - about and their last value kept. + The maximum number of times the step is halved before points that are still refining + are warned about and their last value kept. + roundoff_guard + The relative change below which a growing change is read as the roundoff floor rather + than as a field the step has not resolved yet. """ grid_values = grid.array if hasattr(grid, "array") else grid grid_values = np.array(grid_values, dtype=np.float64, copy=True) - hessian_h = np.stack( + hessian_full = np.stack( self._hessian_via_finite_difference(grid=grid, buffer=buffer) ).astype(np.float64) - hessian_h2 = np.stack( + hessian_half = np.stack( self._hessian_via_finite_difference(grid=grid, buffer=buffer / 2.0) ).astype(np.float64) - richardson = (4.0 * hessian_h2 - hessian_h) / 3.0 - change = np.full(richardson.shape, np.inf) + richardson = (4.0 * hessian_half - hessian_full) / 3.0 - unconverged = np.ones(richardson.shape[1], dtype=bool) + total = richardson.shape[1] + relative_change = np.full(total, np.inf) + refining = np.ones(total, dtype=bool) step = buffer halvings = 0 - while np.any(unconverged) and halvings < max_halvings: + while np.any(refining) and halvings < max_halvings: halvings += 1 step /= 2.0 - index = np.flatnonzero(unconverged) + index = np.flatnonzero(refining) grid_subset = aa.Grid2DIrregular(values=grid_values[index]) - hessian_subset_h = hessian_h2[:, index] - hessian_subset_h2 = np.stack( + full_subset = hessian_half[:, index] + half_subset = np.stack( self._hessian_via_finite_difference(grid=grid_subset, buffer=step / 2.0) ).astype(np.float64) - richardson_subset = (4.0 * hessian_subset_h2 - hessian_subset_h) / 3.0 - change_subset = np.abs(richardson_subset - richardson[:, index]) + richardson_subset = (4.0 * half_subset - full_subset) / 3.0 - hessian_h[:, index] = hessian_subset_h - hessian_h2[:, index] = hessian_subset_h2 - richardson[:, index] = richardson_subset - change[:, index] = change_subset + change = np.abs(richardson_subset - richardson[:, index]) - unconverged = np.zeros(unconverged.shape, dtype=bool) - unconverged[index] = np.any( - change_subset > atol + rtol * np.abs(richardson_subset), axis=0 + within_tolerance = np.all( + change <= atol + rtol * np.abs(richardson_subset), axis=0 ) - if np.any(unconverged): - total = int(unconverged.shape[0]) - number = int(np.count_nonzero(unconverged)) - denominator = np.where(np.abs(richardson) > 0.0, np.abs(richardson), 1.0) - largest = float(np.max((change / denominator)[:, unconverged])) + denominator = np.where( + np.abs(richardson_subset) > 0.0, np.abs(richardson_subset), 1.0 + ) + change_relative = np.max(change / denominator, axis=0) + + at_roundoff_floor = (change_relative > relative_change[index]) & ( + relative_change[index] <= roundoff_guard + ) + + # Points at the roundoff floor keep R_k, the extrapolant from before the change + # started growing, so their columns are left untouched; every other point adopts + # the new extrapolant and the half-step evaluation that will seed its next pair. + adopted = index[~at_roundoff_floor] + + richardson[:, adopted] = richardson_subset[:, ~at_roundoff_floor] + hessian_half[:, adopted] = half_subset[:, ~at_roundoff_floor] + relative_change[adopted] = change_relative[~at_roundoff_floor] + + refining = np.zeros(total, dtype=bool) + refining[index] = ~(within_tolerance | at_roundoff_floor) + + if np.any(refining): + number = int(np.count_nonzero(refining)) + largest = float(np.max(relative_change[refining])) warnings.warn( f"LensCalc Hessian: {number} of {total} points did not converge after " f"{max_halvings} halvings (largest relative error estimate "