[NAThompson]

Legendre, Hermite, Laguerre and Jacobi Leja rules reuse earlier evaluations; the sparse grid refines where contributions are largest.

leja_legendre, leja_hermite, leja_laguerre(max_points, alpha = 0) (x^α e^{−x} on [0, ∞)) and leja_jacobi(max_points, alpha, beta) ((1−x)^α (1+x)^β on [−1, 1]) are nested weighted Leja rules (one new point per rule, starting at the mean of the weight, weights from a scaled moment system), with an adaptive integrate. sparse_grid(f, rule, dimension, max_calls, tolerance) builds a dimension-adaptive sparse grid (Gerstner–Griebel) from them, against a product of any of these weights. Both are in the nanobind layer.

One dimension

Error vs function calls: Leja (nested, one new call per rule), its free estimate, one-shot Gauss, and Gauss re-run at each doubling.

Where Leja loses, Hermite with a pole at ±i:

Sparse grid

Genz integrands, relative error vs calls, against Tasmanian 8.2 (Leja rules, surplus and anisotropic refinement), one Monte Carlo run and Boost.Random's unshifted Sobol' sequence, as a user would run them:

Where QMC wins, a peaked integrand in eight dimensions:

• Comparable to Tasmanian up to d = 4; at d = 8 ahead for the oscillatory integrand and comparable for the Gaussian and peaked ones; neither reaches 1e-4 at d = 16. Peak memory 44 MB vs 150 MB–1.6 GB.
• Our curves sweep max_calls (tolerance 1e-14) so every point is a budget-limited result; the honesty numbers below come from a tolerance sweep.
• Time to reach 1e-4 and 1e-8, single-threaded, against Tasmanian's faster mode (it through its Python bindings): faster in all 30 cases both reach, median 16×, range 1.3× (the kinked integrand in d = 2, where we need 9× more calls) to 53×; only we reach 1e-4 for the anisotropic product peak in d = 8. Bookkeeping is under 2 µs per call on these integrands (Tasmanian 6–200 µs).
• Stopping: when the estimate passes, every active index is refined once more. Over 828 runs (d = 2–16, tol 1e-1 to 1e-12), 398 claimed convergence; 12 of those had a larger error, all for the kinked exp(−Σc|x−u|) in d = 2, by at most 6.7×. Without the extra round, smooth integrands were wrong too (38 claims, up to 13×, with the earlier 61-point rule).
• A coordinate that needs more points than the rule has gets a larger copy of the rule (2^m − 1 points, up to 1023 for built-in types, 127 otherwise). The weights are accurate to beyond 2048 points in double (moment residual 3e-14); the cap is the O(N⁴) build, about 30 s at 1023.

Laguerre and Jacobi, exponents from −0.9 to 10: stable to 1024 points (cos x, 1/(2+x) within 1.2e-14), points within 8 ulps and weights within 160 ulps of 50-digit builds up to 100 points (α = −0.9: 112 / 1340). Laguerre moments lose n²ε·L1 in the LU solve rather than nε.

Evidence (H5Web, VS Code extension): leja_correctness.h5.zip (exactness for every n, weights vs 100 digits, greedy choices, Lebesgue constants), sparse_grid_evidence.h5.zip (all convergence plots, honesty of the estimate, chosen index sets).

Benchmark (reporting/performance/leja_quadrature_performance.cpp, Apple M-series, -O3, CPU time)

build 64 / 128 points	1.5 / 13 ms (all four rules)
32-point rule, cos 3x	0.25 µs
adaptive, Runge to 1e-6 (82 calls)	3.5 µs
sparse grid, exp(Σ xᵢ/2ⁱ), d = 8, 1e-8 (6111 calls)	2.5 ms
sparse grid, Hermite cos(Σ xᵢ/(i+1)), d = 8, 1e-6 (15697 calls)	6.2 ms

🤖 Generated with Claude Code

[NAThompson]

@JacobHass8 : Check it out at your convenience.

[JacobHass8]

I wonder if there is a way to generalize this for any weight function if you don't know the recurrence relation. I think there is a way if you know the moments of the weight function but I'm not exactly sure.

[NAThompson]

Yes. The two halves of the construction need different things. The points only use √w through log √w and its derivatives (the greedy argmax), so any weight we can evaluate works as is. The orthonormal recurrence (beta(k), mu0) is only used for the weights, to keep the moment system well conditioned.

Recovering the recurrence from raw moments ∫xᵏw is badly conditioned: the Hankel matrix's condition number grows exponentially. The stable routes are Gautschi's:

• the modified Chebyshev algorithm, from modified moments ∫pₖ(x)w(x)dx against a known basis such as Legendre or Chebyshev, which is stable on a bounded interval;
• discretized Stieltjes or Lanczos, which only needs w, discretized by a fine Gauss or tanh-sinh rule.

Either would give a general leja<Real>(w, support). I kept this PR to the two closed-form weights, but a Jacobi weight (also closed form) or a general weight via discretized Stieltjes would be a natural follow-up.

-claude

[NAThompson]

@JacobHass8 : I believe that performance-wise, API wise, feature-wise, and otherwise this is good.

Please try these commands against the wheel at your convenience:

import boost.math as bm, numpy as np

# The rule is nested: adding a point keeps all the earlier ones,
# so the adaptive loop never evaluates f at the same point twice.
r = bm.leja_legendre(max_points=64)
r.abscissa[:5]                       # [0, -1, 1, -0.577, 0.659]
r.integrate_adaptive(np.cos, tolerance=1e-12)
# (1.682941969615793, 6e-13, ...)      exact 2 sin 1 = 1.682941969615793

# Chebyshev weight 1/sqrt(1-x^2) via Jacobi(-1/2, -1/2)
cheb = bm.leja_jacobi(max_points=64, alpha=-0.5, beta=-0.5)
cheb.integrate_adaptive(np.cos, tolerance=1e-12)
# (2.403939430634413, 1e-12, ...)      exact π J0(1) = 2.403939430634413

# Singular weight on a half-line: ∫₀^∞ cos(x) e^{-x} / sqrt(x) dx
bm.leja_laguerre(max_points=128, alpha=-0.5).integrate_adaptive(np.cos, tolerance=1e-10)
# (1.3769963318531544, 6e-11, ...)     exact 1.376996331853153

# An 8-dimensional Gaussian integral on a sparse grid
c = 1 / np.arange(1, 9)
bm.sparse_grid(lambda x: np.cos(c @ x), rule=bm.leja_hermite(max_points=64),
               dimension=8, max_calls=200_000, tolerance=1e-8)
# (66.49082912155218, 3.3e-7, ...)     exact π⁴ exp(-Σcᵢ²/4) = 66.49082913179191

# A product of exponential weights on [0, ∞)³
bm.sparse_grid(lambda x: np.exp(-x[0]/2 - x[1]/3 - x[2]/4),
               rule=bm.leja_laguerre(max_points=64), dimension=3,
               max_calls=100_000, tolerance=1e-12)
# (0.3999999999999982, 1e-13, ...)     exact 2/5

[JacobHass8]

Wow, this is excellent! I had tried to implement sparse grids but it was quite slow and inaccurate. I'm impressed it's so accurate! Thank you!