[NAThompson]

Leja rules reuse earlier evaluations; the sparse grid refines where contributions are largest.

leja_legendre and leja_hermite are nested weighted Leja rules (one new point per rule, weights from a √w-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, over [−1, 1]^d or against exp(−|x|²) on ℝ^d. 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), Monte Carlo and randomized Sobol' QMC (medians of 32):

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.
• Bookkeeping 2–4 µs per point at d = 8–16, up to ~40 µs at high levels in d = 2–4 (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), 378 claimed convergence; 4 of those had a larger error, all for the kinked exp(−Σc|x−u|) in d = 2, by at most 3.8×. Without the extra round it was 38, up to 13×.

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.7 / 14.5 ms
32-point rule, cos 3x	0.28 µs
adaptive, Runge to 1e-6 (82 calls)	4.1 µs
sparse grid, exp(Σ xᵢ/2ⁱ), d = 8, 1e-8 (6111 calls)	6.2 ms
sparse grid, Hermite cos(Σ xᵢ/(i+1)), d = 8, 1e-6 (15697 calls)	20 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