Add nested Leja quadrature and a dimension-adaptive sparse grid

doc/math.qbk

@@ -774,6 +774,7 @@ and as a CD ISBN 0-9504833-2-X  978-0-9504833-2-0, Classification 519.2-dc22.
 [include quadrature/gauss.qbk]
 [include quadrature/gauss_kronrod.qbk]
 [include quadrature/gauss_hermite.qbk]
+[include quadrature/leja.qbk]
 [include quadrature/double_exponential.qbk]
 [include quadrature/ooura_fourier_integrals.qbk]
 [include quadrature/naive_monte_carlo.qbk]

doc/quadrature/leja.qbk

@@ -0,0 +1,184 @@
+[/
+Copyright (c) 2026 Nick Thompson
+Use, modification and distribution are subject to the
+Boost Software License, Version 1.0. (See accompanying file
+LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+]
+
+[section:leja Leja quadrature]
+
+[heading Synopsis]
+
+`#include <boost/math/quadrature/leja.hpp>`
+
+   namespace boost{ namespace math{ namespace quadrature{
+
+   template <class Real, class ``__Policy`` = boost::math::policies::policy<> >
+   class leja_legendre  // the integral of f over [-1, 1]
+   {
+   public:
+      leja_legendre();
+
+      const std::vector<Real>& abscissa() const;
+      const std::vector<Real>& weights(unsigned n) const;
+      static constexpr unsigned max_points();
+
+      template <class F>
+      auto integrate(F f, Real tolerance = tools::root_epsilon<Real>(),
+                     Real* error = nullptr, Real* L1 = nullptr) const -> decltype(f(Real(0)));
+
+      // The integral over [a, b]; leja_legendre and leja_jacobi only.
+      template <class F>
+      auto integrate(F f, Real a, Real b, Real tolerance = tools::root_epsilon<Real>(),
+                     Real* error = nullptr, Real* L1 = nullptr) const -> decltype(f(Real(0)));
+   };
+
+   template <class Real, class ``__Policy`` = boost::math::policies::policy<> >
+   class leja_hermite;  // the integral of f(x) exp(-x[super 2]) over the real line, with the same members except integrate(f, a, b)
+
+   // The integral of f(x) x[super alpha] exp(-x) over \[0, [infin]), alpha > -1; as leja_hermite, with the constructor
+   //    explicit leja_laguerre(Real alpha = 0);
+   template <class Real, class ``__Policy`` = boost::math::policies::policy<> >
+   class leja_laguerre;
+
+   // The integral of f(x) (1 - x)[super alpha] (1 + x)[super beta] over \[-1, 1\], alpha, beta > -1, with the constructor
+   //    leja_jacobi(Real alpha, Real beta);
+   template <class Real, class ``__Policy`` = boost::math::policies::policy<> >
+   class leja_jacobi;
+
+   }}} // namespaces
+
+[heading Description]
+
+A weighted Leja sequence is built greedily: each new point maximizes [radic]/w/ times the product of its distances to the points
+already chosen (Leja; Narayan and Jakeman). `leja_legendre` has the Legendre weight /w/(/x/) = 1 on \[-1, 1\], which gives
+the classical Leja points, here in the order 0, -1, 1; `leja_hermite` has /w/(/x/) = exp(-/x/[super 2]), `leja_laguerre`
+/w/(/x/) = /x/[super [alpha]] exp(-/x/) on \[0, [infin]), and `leja_jacobi` /w/(/x/) = (1 - /x/)[super [alpha]] (1 + /x/)[super [beta]] on \[-1, 1\].
+All start from the mean of /w/, so that the one-point rule is the one-point Gauss rule, exact to degree 1, the natural first level
+of a sparse grid. Where an exponent is negative, [radic]/w/ has no maximum, and the points are chosen as if it were 0; the weights
+are those of /w/. `leja_jacobi` with [alpha] = [beta] = 0 has the points of `leja_legendre`.
+The constructors raise a __domain_error unless the exponents are > -1.
+The /n/-point rule uses the first /n/ points and has the weights that make it exact for polynomials of degree /n/ - 1.
+
+The rules are nested: going from /n/ to /n/ + 1 points costs one more evaluation of /f/, where a Gauss rule must start again.
+This makes them the usual building block of sparse grids, and lets the adaptive `integrate` add one point at a time.
+The price is degree: a Leja rule with /n/ points is exact to degree /n/ - 1, a Gauss rule to 2/n/ - 1.
+
+The points and weights are computed when they are first needed, up to `max_points()`: 1023 for the built-in floating point
+types and 127 for others. `weights(n)` returns the /n/ weights of the /n/-point rule, computing it and any smaller ones first,
+and raises a __domain_error unless 1 [le] /n/ [le] `max_points()`. `abscissa()` returns the points computed so far, in
+the order they were chosen. All members are `const` and safe to call from several threads, which share what any of them
+computes, as `tanh_sinh` does; points that another thread is computing read as NaN until it is done.
+`example/leja_parallel_assembly.cpp` assembles the Hamiltonian of the anharmonic oscillator /x/[super 4] in 60 Hermite functions
+on all cores with one `leja_legendre`, in agreement with the exact entries to 2[times]10[super -15] times the largest: 0.21 s on one core
+of an Apple M-series and 0.026 s on ten, once the rule has its points.
+`leja_hermite` suits expectations over a normal distribution, such as the mean of the logistic function of /X/ ~ N(1, 1):
+
+   leja_hermite<double> rule;
+   double error;
+   double mean = rule.integrate([](double x) { return 1 / (1 + std::exp(-1 - std::sqrt(2.0) * x)); }, 1e-10, &error)
+               / std::sqrt(boost::math::constants::pi<double>());
+
+takes 105 points for 0.69673467014367874, within 4.6[times]10[super -15] of the true value; the estimate of its error, `error` divided by [radic][pi], is 5.6[times]10[super -12].
+
+When several points maximize the product to within rounding, the leftmost is taken; in the tests float, double, long double and
+`cpp_bin_float_50` choose the same points. Some weights are negative, but not by much: up to 100 points the sum of |/w[sub i]/|
+is at most 1.32 times the total mass for the exponents tried, but 2.44 times for [alpha] = -0.9, [beta] = 1/2.
+
+`integrate(f, tolerance, error, L1)` evaluates /f/ once at each point in turn, until /n/ [ge] 4 and
+|/Q[sub m]/ - /Q[sub m - max(1, m/4)]/| [le] `tolerance` * /L1/ for every /m/ from /n/ - /n//4 to /n/, where /Q[sub n]/ is the /n/-point rule
+and /L1/ = [sum] |/w[sub i]/| |/f/(/x[sub i]/)|. If `error` or `L1` are not null they are set to the largest of those differences,
+a heuristic estimate of the error, and to /L1/.
+If `max_points()` is reached first, its result is returned with its error estimate, as `tanh_sinh` does; no error is raised.
+/f/ may return a real or complex value, or a vector of them: any type with `operator[]` and `size()`, such as `std::array` or
+`std::vector`, whose values all have the same size. A vector is integrated componentwise, and |/y/| above is its max norm.
+`leja_legendre` evaluates /f/ at the endpoints [plusminus]1, and `leja_laguerre` and `leja_jacobi` at an end of the support where the exponent is [le] 0.
+
+`integrate(f, a, b, tolerance, error, L1)` integrates /f/ over \[/a/, /b/\] with `leja_legendre`, and /f/(/x/) (/b/ - /x/)[super [alpha]] (/x/ - /a/)[super [beta]]
+with `leja_jacobi`, by the affine map from \[-1, 1\], which also scales the error and /L1/. If /a/ > /b/ the result is minus the
+integral over \[/b/, /a/\], with [alpha] still at /b/. Infinite or NaN bounds raise a __domain_error.
+
+[heading Accuracy and cost]
+
+The rules up to /N/ points take O(/N/[super 2]) operations: about 0.3 ms for 64 points, 1 ms for 128,
+16 ms for 512 and 60 ms for 1023 in double precision, and 0.4 s for 127 in `cpp_bin_float_50`.
+Each rule follows from the last in Newton form, with the integral of the nodal polynomial from a recurrence carried to twice the
+working precision: the rules integrate /x/[super /k/] to within 20 /n/ [epsilon] [sum] |/w[sub i]/ /x[sub i]/[super /k/]|, and in double
+precision the weights are within an ulp of the largest weight of the same rules in 50 digits, up to 1023 points.
+A bounded sum of |/w[sub i]/| does not help integrands that grow like exp(/x/[super 2]) against `leja_hermite`, such as products of Hermite polynomials
+of high degree: /L1/ is 10[super 19] times the integral of /h/[sub 40][super 2]. Integrate such products with their Gaussian over
+a finite interval instead, as the example does.
+
+[heading References]
+
+* F. Leja, ['Sur certaines suites liées aux ensembles plans et leur application à la représentation conforme], Annales Polonici Mathematici 4 (1957), 8-13.
+* L. Reichel, ['Newton interpolation at Leja points], BIT 30 (1990), 332-346, [@https://doi.org/10.1007/BF02017352 doi:10.1007/BF02017352].
+* A. Narayan and J. D. Jakeman, ['Adaptive Leja sparse grid constructions for stochastic collocation and high-dimensional approximation], SIAM Journal on Scientific Computing 36 (2014), no. 6, A2952-A2983, [@https://doi.org/10.1137/140966368 doi:10.1137/140966368].
+
+[endsect] [/section:leja Leja quadrature]
+
+[section:sparse_grid Sparse grids]
+
+[heading Synopsis]
+
+`#include <boost/math/quadrature/sparse_grid.hpp>`
+
+   namespace boost{ namespace math{ namespace quadrature{
+
+   // Rule: leja_legendre<Real>, leja_hermite<Real>, leja_laguerre<Real> or leja_jacobi<Real>, Real = typename Rule::value_type; f(const std::vector<Real>& x), or f(x) with x a std::array<Real, N>, Eigen::VectorXd or other container with operator[] and size().
+   template <class F, class Rule>
+   auto sparse_grid(F f, const Rule& rule, unsigned dimension, std::size_t max_calls = 100000,
+                    Real tolerance = tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr)
+      -> decltype(f(x));
+
+   // Over the box [a[0], b[0]] x ... x [a[d-1], b[d-1]], d = a.size(), for leja_legendre and leja_jacobi; Box has operator[] and size().
+   template <class F, class Rule, class Box>
+   auto sparse_grid(F f, const Rule& rule, const Box& a, const Box& b, std::size_t max_calls = 100000,
+                    Real tolerance = tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr)
+      -> decltype(f(x));
+
+   }}} // namespaces
+
+[heading Description]
+
+`sparse_grid` integrates /f/ over \[-1, 1\][super /d/] with `leja_legendre`, or /f/(/x/) exp(-|/x/|[super 2]) over the whole space
+with `leja_hermite`, or against the product of the weights of `leja_laguerre` or `leja_jacobi`, by the dimension-adaptive sparse grid of Gerstner and Griebel. Coordinate /i/ at level /k[sub i]/ uses the
+first 2/k[sub i]/ - 1 points of the rule, and the grid is the sum of the tensor products of differences of consecutive levels
+over a set of multi-indices /k/. The set grows where these contributions are largest, so the grid spends its points on the
+coordinates, and the interactions between them, that matter. Exact integrals of products of polynomials of degree 2/k[sub i]/ - 2
+need only the multi-index /k/ and those below it, and a function of a few coordinates costs little more than in those dimensions.
+
+The levels add two points at a time, not one: for an even /f/ such as /x/[super 2] the rules with one and two points agree, and a
+zero difference would stop the refinement of that coordinate. For the same reason, when the error estimate (the sum of the
+norms of the contributions not yet refined) falls below `tolerance` * /L1/, every such contribution is refined once more,
+and the result is accepted only if the estimate is still below the tolerance. Nothing is accepted while /L1/ is 0, as it is while /f/ has been
+zero at every point: /x/[sub 1][super 2] /x/[sub 2][super 2], with `leja_legendre`, is zero at the first point, the mean of the weight,
+and on the axes through it. After /d/ + 1 such rounds of refinement with /f/ still zero everywhere, the result is 0 with an error estimate of 0. An interaction of three or more variables
+that vanishes on the planes through the mean can still be missed if /f/ is not zero there, such as 1 + /x/[sub 1][super 2] /x/[sub 2][super 2] /x/[sub 3][super 2]. /L1/ is the grid's estimate of the integral of |/f/| times the
+weight, the scale of the tolerance; before the grid has converged it can be far off, even negative. It is not a condition
+number: the sparse grid's weights have mixed signs, and the sum of |/w[sub i] f/(/x[sub i]/)| can be much larger.
+If `error` or `L1` are not null they are set to these. The grid computes the points of `rule` it needs, and keeps them there for the next call.
+
+`max_calls` bounds the number of evaluations of /f/: a refinement that would exceed it is not made, and the result has an
+error estimate above the tolerance. A contribution that would need more points than that in a coordinate is not refined
+further, and its norm stays in the error estimate;
+if the calls run out during the final refinements the error estimate is NaN. No error is raised. A __domain_error is
+raised if `dimension` or `max_calls` is zero, or if /f/ takes a fixed-size container whose size is not `dimension`. /f/ may return a real or complex value or a vector of them, as for the
+one-dimensional rules; a vector is integrated componentwise, and its error estimate and /L1/ use the max norm.
+With bounds `a` and `b`, each coordinate is mapped as `integrate(f, a[i], b[i])` maps it, and the result, error and /L1/ are scaled by the
+product of the factors; a __domain_error is raised if a bound is infinite or NaN or `b.size()` is not `a.size()`.
+
+The bookkeeping costs under 2 microseconds per evaluation of /f/ in double precision on Genz's test integrands, more where
+the levels are high, so a sparse grid pays off when /f/ costs more than that. On Genz's test integrands, with 10[super 4] calls or more, it beats
+quasi-Monte Carlo (Boost.Random's Sobol' sequence) by orders of magnitude up to four dimensions; in eight it still wins for
+entire functions such as cos(0.7 + [sum] /c[sub i] x[sub i]/), but not for peaked ones such as Genz's product peak and Gaussian,
+and in sixteen quasi-Monte Carlo is better unless /f/ depends mostly on a few coordinates. For integrands with kinks, such as
+exp(-[sum] /c[sub i]/|/x[sub i]/ - /u[sub i]/|), quasi-Monte Carlo is better in any dimension, and the error estimate, a heuristic
+as for the one-dimensional rule, can be too small.
+
+[heading References]
+
+* T. Gerstner and M. Griebel, ['Dimension-adaptive tensor-product quadrature], Computing 71 (2003), 65-87, [@https://doi.org/10.1007/s00607-003-0015-5 doi:10.1007/s00607-003-0015-5].
+* A. Narayan and J. D. Jakeman, as above.
+
+[endsect] [/section:sparse_grid Sparse grids]

example/Jamfile.v2

@@ -151,6 +151,7 @@ test-suite examples :
    [ run reciprocal_fibonacci_constant.cpp : : : [ check-target-builds ../config//has_mpfr "MPFR Support" : <linkflags>"-lmpfr -lgmp" : <build>no ] ]
    [ run reverse_mode_linear_regression_example.cpp ] 
    [ run autodiff_reverse_black_scholes.cpp ]
+   [ compile leja_parallel_assembly.cpp : [ requires cxx17_if_constexpr cxx11_hdr_thread ] [ check-target-builds ../../multiprecision/config//has_eigen : : <build>no ] ]
 ;
 
 run root_elliptic_finding.cpp /boost/timer : : : release <link>static  [ requires cxx11_unified_initialization_syntax cxx11_defaulted_functions ] <target-os>freebsd:<linkflags>"-lrt" <target-os>linux:<linkflags>"-lrt -lpthread" ;

example/leja_parallel_assembly.cpp

@@ -0,0 +1,111 @@
+// Copyright Nick Thompson 2026.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Parallel assembly of a Hamiltonian with one shared Leja rule.
+//
+// The anharmonic oscillator H = -d^2/dx^2 + x^4, in the basis of Hermite functions psi_n(x / s), n < N. With y = x / s,
+// H_ij = int psi_i(y) ((2j + 1 - y^2) / s^2 + s^4 y^4) psi_j(y) dy, since -psi_j'' = (2j + 1 - y^2) psi_j. The integrands
+// are entire and negligible beyond |y| = L, and the threads integrate them over [-L, L] with one leja_legendre: whichever
+// first needs more points computes them, and the others use them. Each entry takes as many points as it needs, and gets
+// an error estimate from the points it has already used.
+//
+// The matrix elements are also known exactly, from x = (a + a^+) / sqrt(2), and the lowest eigenvalues agree with those
+// of Hioe and Montroll (1975).
+//
+// A rule for exp(-y^2), leja_hermite, would be the wrong choice: h_i h_j grows like exp(y^2), the Leja points spread
+// further than the Gauss points, and sum |w_k h_i h_j(x_k)| is 1e19 times H_ij for i = j = 40.
+
+#include <boost/math/quadrature/leja.hpp>
+#include <Eigen/Dense>
+#include <algorithm>
+#include <atomic>
+#include <chrono>
+#include <cmath>
+#include <iomanip>
+#include <iostream>
+#include <thread>
+#include <vector>
+
+// psi_i(y) psi_j(y), from h_{n+1} = y sqrt(2 / (n + 1)) h_n - sqrt(n / (n + 1)) h_{n-1} times exp(-y^2 / 2).
+double psi_product(unsigned i, unsigned j, double y)
+{
+   double previous = 0;
+   double psi = std::exp(-y * y / 2) / std::sqrt(std::sqrt(boost::math::constants::pi<double>()));
+   double psi_i = i == 0 ? psi : 0;
+   double psi_j = j == 0 ? psi : 0;
+   for (unsigned n = 0; n < (std::max)(i, j); ++n)
+   {
+      const double next = y * std::sqrt(2.0 / (n + 1)) * psi - std::sqrt(double(n) / (n + 1)) * previous;
+      previous = psi;
+      psi = next;
+      if (n + 1 == i)
+         psi_i = psi;
+      if (n + 1 == j)
+         psi_j = psi;
+   }
+   return psi_i * psi_j;
+}
+
+// H from the matrix X of y, tridiagonal with X_{n,n+1} = sqrt((n + 1) / 2): X^2 and X^4 are exact in the first N rows
+// when X has N + 4.
+Eigen::MatrixXd exact(unsigned N, double s)
+{
+   Eigen::MatrixXd X = Eigen::MatrixXd::Zero(N + 4, N + 4);
+   for (unsigned n = 0; n < N + 3; ++n)
+      X(n, n + 1) = X(n + 1, n) = std::sqrt((n + 1) / 2.0);
+   const Eigen::MatrixXd X2 = X * X;
+   const Eigen::MatrixXd D = Eigen::VectorXd::LinSpaced(N, 1, 2 * N - 1).asDiagonal();
+   return (D - X2.topLeftCorner(N, N)) / (s * s) + s * s * s * s * (X2 * X2).topLeftCorner(N, N);
+}
+
+int main()
+{
+   const unsigned N = 60;
+   const double s = 0.75;
+   const double L = 14;
+   const unsigned threads = (std::max)(1u, std::thread::hardware_concurrency());
+   const boost::math::quadrature::leja_legendre<double> rule;
+
+   Eigen::MatrixXd H(N, N);
+   Eigen::MatrixXi calls(N, N);
+   // The entries i <= j, taken in turn by the threads.
+   auto assemble = [&] {
+      std::atomic<unsigned> next{0};
+      auto work = [&] {
+         for (unsigned e = next++; e < N * N; e = next++)
+         {
+            const unsigned i = e / N, j = e % N;
+            if (i > j)
+               continue;
+            unsigned count = 0;
+            H(i, j) = H(j, i) = rule.integrate([&](double y) {
+               ++count;
+               return psi_product(i, j, y) * ((2 * j + 1 - y * y) / (s * s) + s * s * s * s * y * y * y * y);
+            }, -L, L, 1e-13);
+            calls(i, j) = calls(j, i) = count;
+         }
+      };
+      const auto start = std::chrono::steady_clock::now();
+      std::vector<std::thread> pool;
+      for (unsigned t = 0; t < threads; ++t)
+         pool.emplace_back(work);
+      for (std::thread& t : pool)
+         t.join();
+      return std::chrono::duration<double>(std::chrono::steady_clock::now() - start).count();
+   };
+   const double cold = assemble();
+   const double warm = assemble();
+
+   const Eigen::MatrixXd H_exact = exact(N, s);
+   const double worst = (H - H_exact).cwiseAbs().maxCoeff(), largest = H_exact.cwiseAbs().maxCoeff();
+   std::cout << std::setprecision(17) << N * (N + 1) / 2 << " entries on " << threads << " threads, with "
+             << calls.minCoeff() << " to " << calls.maxCoeff()
+             << " points: " << cold << " s computing the rule as it goes, " << warm << " s with it computed.\n"
+             << "max |H - H_exact| = " << worst << ", " << worst / largest << " of the largest entry.\n";
+   const Eigen::VectorXd E = Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd>(H, Eigen::EigenvaluesOnly).eigenvalues();
+   const double reference[] = {1.0603620904841829, 3.7996730298013941, 7.4556979379867383, 11.644745511378162};
+   for (unsigned n = 0; n < 4; ++n)
+      std::cout << "E_" << n << " = " << E(n) << ", Hioe and Montroll " << reference[n] << '\n';
+}