fix: adaptive Richardson step in LensCalc NumPy Hessian (#591) - #593
Merged
Conversation
`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 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01WzKP1RJKhPrtFmsHSm6EEA
…mates (#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 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01WzKP1RJKhPrtFmsHSm6EEA
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 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WzKP1RJKhPrtFmsHSm6EEA
Collaborator
Author
|
Downstream PyAutoLens test PR: PyAutoLabs/PyAutoLens#717 |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
LensCalc._hessian_via_richardsondifferenced the deflections at a hardcoded 0.01" step and discarded the h/(h/2) pair's error information, so near compact structure the NumPy Hessian was badly wrong: on an isothermal sphere with R_E = 0.2" evaluated 3e-4"–1e-3" from its centre, κ came out ~46.5 where the analytic value is 100–333 (54–86 % off) and the shear was ~100 % off; on the PyAutoLens#480 multi-plane fixture the last-plane magnification was 102–122 % off with the sign flipped on all four image positions. The NumPy path is the default for point-source flux fits (magnifications_at_positions) andPointSolver's magnification threshold.This makes the step adaptive per point: Richardson estimates R_k are formed from (h, h/2) and (h/2, h/4), the step is halved on the unconverged subset only (reusing the previous half-step evaluation), convergence is
|R_{k+1} − R_k| ≤ atol + rtol·|R_{k+1}|withrtol=1e-7,atol=1e-8plus the roundoff-floor stop described below, at most 20 halvings, and any point that still has not converged (a genuinely singular point, e.g. exactly on an isothermal centre) keeps its last value and raises exactly oneUserWarningper call — never silent, never an exception. Smooth fields exit after three finite-difference evaluations (was two; a 50×50 grid costs 0.007 s); only points near compact structure iterate.The roundoff-floor stop is the second termination rule. Finite differences improve as O(h²) only until cancellation in
(f(x+h) − f(x−h))takes over; past that turning point halving makes the answer worse. On the #480 fixture the relative change between successive extrapolants bottoms out at 3.3e-08 at the 14th halving and then grows again (5.98e-08, 1.73e-07, 7.02e-07, …). So a point whose change grows instead of shrinking, while already belowroundoff_guard=1e-4relative, stops there and keeps R_k — the estimate from before the growth — with no warning. The guard means a true singularity (relative change pinned at ~0.5, never entering the guard band) still runs tomax_halvingsand still warns. No configuration can warn permanently: the #480 fixture is silent at rtol 1e-7, 1e-8 and 1e-10 alike.Result on the compact-SIS control: κ relative error 1.14e-09, shear 7.57e-09 vs analytic (was 0.54–1.0). On the #480 fixture the NumPy Hessian now agrees with the ray-traced Jacobian and JAX autodiff to 3e-4 (the stored constant's own precision), with no warning.
API Changes
No public signature changes.
LensCalc.hessian_from(grid, xp=np)and everything built on it (convergence_2d_via_hessian_from,shear_yx_2d_via_hessian_from,magnification_2d_via_hessian_from,jacobian_from, critical curves / caustics) return more accurate values near compact structure; smooth-field results move by ≲ 6e-9 relative (the adaptive result is the more accurate one — verified against analytic convergence). AUserWarningis now emitted for points that cannot converge. The JAX path is unchanged.See full details below.
Test Plan
main(κ 54–86 % off, shear ~100 % off) and passes with the fixtest_autogalaxy/operate+test_isothermal.pygreen; fulltest_autogalaxy1152 passedtest_autolens570 passed with that branch)Full API Changes (for automation & release notes)
Removed
Added
LensCalc._hessian_via_richardson(grid, buffer=0.01, rtol=1e-7, atol=1e-8, max_halvings=20, roundoff_guard=1e-4)— private; adaptive per-point Richardson extrapolation with a roundoff-floor stop (was a fixedbuffer=0.01)UserWarningfrom the NumPy Hessian path when points do not converge aftermax_halvingsMigration
Closes #591.
Generated by the PyAutoLabs agent workflow.
🤖 Generated with Claude Code
https://claude.ai/code/session_01WzKP1RJKhPrtFmsHSm6EEA