Repository navigation
SO2onDisk: default to the Nyquist band-limit so the DHT is well conditioned - #54
Merged
Merged
Conversation
torch.special.bessel_j0/j1 are only ~5e-7 accurate in float64 on
x in [2.8, 14.5] (torch 2.11, arm64). bessel_jn used them as base
cases for n = 0, 1 and as seeds of the forward recurrence, so every
order, every root and the whole disk-harmonic basis inherited the
error (root j_{0,5} was off by 1e-6; J_n at roots was ~1e-7, not 0).
Evaluate all orders with Miller's backward recurrence, started above
max(n, ceil(max|x|)) so it is stable in both the x < n and the
oscillatory regime. A single pass now returns several orders
(_bessel_jn_multi), which the Newton root finders use to get
J_{n-1}, J_n, J_{n+1} in one evaluation.
Also fix a convergence pathology in the bracketed Newton solvers: once
Newton had converged (3 iterations) the iterate sat on a bracket
endpoint, the next step fell outside by rounding, and the solver fell
back to ~45 bisection iterations. A converged Newton step now pins
the bracket. Net effect: J_n accurate to 1e-15 on [0, 160] for
n <= 60, roots to 1e-14, and compute_all_bessel_roots is 3x faster.
Tests: mpmath-referenced J_n values (atol 1e-14) for n in
{0, 1, 2, 5, 20} on x up to 150, roots to 1e-13, J_0(0) = 1,
root residuals < 1e-13; the two tests that asserted equality with
torch.special are replaced.
…tioned
The default truncation followed fle_2d: ne = floor(pi L^2 / 4) basis
functions, one per disk pixel. That makes the real basis matrix Phi
nearly square and admits radial frequencies ~25% above the grid
Nyquist rate, so cond(Phi) ~ 1e16 at every L. pinv(rcond=1e-10) then
silently discarded directions: DHT roundtrips, inversion and even
exact 90-degree pixel rotations were only accurate to O(1), which the
tests tolerated (rel_err < 1.5, < 0.5, 85% of coefficients within 20%).
Default bandlimit is now pi * L / 2 (new public helper
default_bandlimit). cond(Phi) is ~2 for L = 8..64; DHT roundtrip,
coefficient inversion and bispectrum roundtrip are exact to 1e-14, and
rot90 invariance to 1e-14. The ne-truncation branch and its
conjugate-pair patch are removed; an explicit bandlimit still works.
Coefficient counts change (L=16: 144 DH coefficients, 76 selective
bispectrum entries instead of 105). README table updated.
Tests: cond(Phi) < 10 for L in {16, 32, 64}; exact rank; roundtrip
tolerances tightened to 1e-9..1e-12; new rot90 test; the arbitrary-angle
spatial test now uses band-limited content and bicubic interpolation
(its 0.1 tolerance is interpolation-bound, measured ~2e-2).
This was referenced Oct 14, 2026
johmathe
changed the base branch from
johmathe/fix-bessel-j0-j1-accuracy
to
main
October 14, 2026 16:41
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.
Stacked on #53 (needs the accurate Bessel roots). Merge #53 first, then retarget this to
main.Problem
The default truncation copied fle_2d's coefficient count
ne = ⌊πL²/4⌋, which makes the analysis matrix Φ nearly square on an L×L grid: cond(Φ) ≈ 1e16.pinv(rcond=1e-10)then silently truncated, and the tests toleratedrel_err < 1.5to pass.Fix
default_bandlimit(L) = π·L/2(grid Nyquist) is the new default; cond(Φ) ≈ 2 for L = 8…64.netruncation branch and the conjugate-pair-completion hack from_build_truncated_indices. Explicitbandlimitstill works.Breaking
Default
output_size/ feature layout changes (L=16: 105 → 76). README table updated.Tests
cond(Φ) < 10(L = 16, 32, 64); full-rank basis; DHT roundtrip 0.5 → 1e-12;invertexact coefficients (all, 1e-9) and bispectrum roundtrip 0.5 → 1e-10; rot90 invariance 1e-12; spatial rotation invariance on band-limited inputs; new Nyquist coefficient-count test.Part of the v0.4.0 stack (B).