Repository navigation
Fix J_0/J_1 accuracy: Miller recurrence for every Bessel order - #53
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.
This was referenced Oct 14, 2026
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #53 +/- ##
==========================================
+ Coverage 95.46% 95.56% +0.09%
==========================================
Files 11 11
Lines 2227 2231 +4
==========================================
+ Hits 2126 2132 +6
+ Misses 101 99 -2 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
johmathe
added a commit
that referenced
this pull request
Oct 14, 2026
* Release v0.4.0 Bump version to 0.4.0 and record the feature-layout changes in the SO3onS2, DnonDn and SO2onDisk class docstrings (versionchanged notes). Breaking: output_size / layout changes for SO3onS2 (lmax >= 3), DnonDn (odd n >= 5) and the default SO2onDisk truncation. * Raise coverage to ~98.6%: drop dead code, test dense-CG cache and graph wiring Dead code removed (never reachable from the library): * _bessel: scalar _jn_scalar/_djn_scalar/_bisect_newton (only tests used them; the batch solver is the single root path) and the unreachable num_brackets <= 0 branch in compute_all_bessel_roots. * dn_on_dn: n2d == 0 / n3 == 0 branches (n >= 3 implies n2d >= 1). * so3_on_s2: empty-slice guard in the power loop of _ensure_precomputed_on_device (empty entries are never flagged power). Bug found by the new tests: _load_or_compute_reduced_cg did not catch pickle.UnpicklingError, so a corrupt reduced-CG cache file crashed module construction instead of being recomputed. New tests: package version fallback; _bessel_jn_multi(()) and order-0 batch roots; _save_cache write failure; triangle-violating sparse CG entries (empty tables, exact zeros in forward); SO3onS2 dense tables vs sparse forward (incl. CG power); reduced-CG disk cache compute/load/ corrupt/write-failure paths; lazy sparse gather-table construction; CUDA graph capture wiring with torch.cuda stream/graph APIs stubbed on CPU. Remaining uncovered lines are defensive numerical fallbacks in the D_n / octahedral CG constructions and the CUDA dispatch in SO3onS2.forward.
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.
Problem
bessel_jnseeded orders 0/1 fromtorch.special.bessel_j0/j1, which are only ~5e-7 accurate on x ∈ [2.8, 14.5] (torch 2.11, arm64). Forward recurrence inherited the error, and so did Bessel roots and theSO2onDiskDHT basis.Fix
nandx(_bessel_jn_multi), start order frommax(n, ⌈max|x|⌉), normalised withJ_0 + 2ΣJ_{2k} = 1, periodic rescale against overflow;J_0(0) = 1._bisect_newton_batchpins the bracket on a tiny Newton step (the old loop converged in 3 steps but then bisected ~45 times because the bracket never collapsed)._bessel_jn_forward/_bessel_jn_miller.Result
|J_n| error ≲ 1e-15 on [0, 160] for n ≤ 60; roots to 1.4e-14; root finding 3× faster than before (0.27 s vs 0.83 s).
Tests
mpmath-referenced
J_nvalues (30 points, incl. x = 3.5, 8, 12 where torch is worst; atol 1e-14) and roots (1e-13); float32 input; batched; high-order residuals. Replaces the oldtest_j0/j1_matches_torch, which asserted agreement with the inaccurate reference.Part of the v0.4.0 stack (A). Stacked on top: #54 (disk conditioning).
Made with Cursor