PaperlayerView exact revision on arXiv

Source-linked AI summary

A higher order lattice-based method for high-dimensional numerical integration using periodization

Felix Bartel, Alexander D. Gilbert, Michael Griebel, Frances Y. Kuo, Ian H. Sloan

arXiv:2610.01113v1math.NAcs.NA

TL;DR

High-dimensional integration poses a challenge: classical deterministic product quadrature slows with dimension, while Monte Carlo converges slowly regardless of smoothness. The paper uses transformed lattice points in a weighted quadrature rule and reports dimension-independent error for p=∞, with at most (1.110721)^s growth otherwise.

  • Problem

    Classical deterministic product quadrature suffers the curse of dimensionality, while Monte Carlo quadrature has a slow rate regardless of smoothness.

    Source: §1.1 · p. 2

  • Method

    The paper constructs unequally weighted quadrature rules using transformed lattice points, optionally subsampled, for non-periodic integrands over the high-dimensional unit cube.

    Source: §1 · p. 1

  • Results

    For p=∞ there is no dependence on s; otherwise growth is at most (1.110721)^s, occurring when p=2.

    Source: Error analysis · p. 6

  • Takeaways & Limitations

    The two error bounds can be interpolated to trade convergence rate against growth in s.

    Source: §1.1 · p. 2

    The random approach requires logarithmic oversampling and accepts a loss in error decay.

    Source: §4.2 · p. 18

Abstract

from arXiv · show

We present a novel numerical integration method for non-periodic functions over the high-dimensional unit cube $[0,1]^s$ by a specially crafted (unequally-)weighted quadrature rule using transformed lattice points, with the additional option of subsampling. The method is designed for integrands whose mixed derivatives up to order $\alpha$ are square-integrable with respect to a product Chebyshev density. The ingredients are (i) periodizing the integrand, (ii) approximating the smooth periodic function by a kernel interpolant at lattice points, (iii) integrating exactly the product of the kernel interpolant and the nonsmooth density, together with subsampling alternatives for (ii) and~(iii). With $|J|$ corresponding to the number of function evaluations, we achieve the convergence rate close to the order $|J|^{-(\alpha-1/2)}$ and $|J|^{-\alpha}\,($1.110721$)^s$, where the implied constants are independent of $s$ under favorable conditions. We can interpolate these two results to trade between the convergence rate and the growth in $s$. We provide numerical experiments to demonstrate our theory.

1 Introduction · 1.1 Background

The paper introduces unequally weighted quadrature using transformed lattice points for non-periodic integration, aiming for higher-order convergence with little or no growth in dimension. Its introduction situates the method among existing approaches and describes its procedure and smoothness assumptions.

1 Introduction

  • The method transforms the integrand into a smooth periodic function times a non-smooth density, approximates it at lattice points with a kernel, and exactly integrates the interpolant times the density; subsampling is optional.The introduction presents the method as a weighted quadrature rule using transformed lattice points.Source: §1 · p. 1
  • The method targets functions with mixed derivatives through order α that behave well at cube boundaries; bounded partial derivatives are one sufficient example.Source: §1 · pp. 1–2

1.1 Background

  • Classical deterministic product quadrature has error O(N^-α/s), while Monte Carlo has root-mean-square error O(N^-1/2), independent of s but slow regardless of α.Source: §1.1 · p. 2
  • Conventional QMC and higher-order methods can approach O(N^-α), but logarithmic factors and implied constants that grow exponentially with s limit their practical use to moderate dimensions.The passage also notes logarithmic factors that may grow exponentially in α.Source: §1.1 · p. 2
  • Weighted mixed Sobolev spaces with summable weights allow lattice and higher-order QMC rates O(N^-α+δ) without log N factors and with dimension-independent implied constants.The passage identifies this dimension-independent behavior as strong tractability and notes similar results for anisotropically adapted sparse grids.Source: §1.1 · p. 2
  • Error-bound constants for methods with the same convergence rate can depend very differently on α, with some growing super-exponentially and requiring many points to resolve the preasymptotic regime.Source: §1.1 · p. 2
  • The proposed unequally weighted rules aim for higher-order convergence with little or no growth in dimension; their bounds can be independent of s under favorable conditions, with an interpolatable rate–dimension tradeoff.The convergence rate is expressed in terms of |J|, the number of function evaluations.Source: §1.1 · p. 2

1.2 Our approach

The approach periodizes the integrand, kernel-interpolates the transformed function at lattice points, and integrates the interpolant against the transformed density; it also allows random subsampling of the lattice.

Step 1. Periodizing the integrand

  • The coordinate map traverses [0, 1] twice, making the transformed function smooth and 1-periodic in each coordinate; the transformed univariate density is half the Jacobian.Source: Step 1. Periodizing the integrand · p. 3

Step 2. Kernel approximation of the function F

Step 3. Exact integration of the approximate integrand

  • The method exactly integrates the kernel interpolant against the density using weights w = K^-1κ, which are independent of F and optimal for weighted integration in the corresponding space.With product kernels, the final step requires accurate computation of only s one-dimensional integrals.Source: Step 3. Exact integration of the… · pp. 4–5

Step 2′ and Step 3′. The subsampling alternative

Subsampling lattice points improves the convergence rate of the method, while the error analysis bounds approximation across Lp norms and exposes a trade-off between convergence rate and dimension dependence. The analysis also extends subsampled-lattice approximation results to L∞ and interpolates them to all p in [2, ∞].

Step 2′ and Step 3′. The subsampling alternative

Error analysis

  • For subsampled lattices, the paper extends prior L2 approximation results to L∞ and interpolates L2 and L∞ bounds to obtain Lp results for every p ∈ [2, ∞].A single lattice generating vector is suited to all p, and the bounds depend on the weighted Korobov norm of F.Source: Error analysis · p. 5
  • The error bounds trade polynomial convergence in N or |J| against dimension growth: p=∞ has no s-dependence, while p=2 grows at most (1.110721)^s.The method does not require choosing p in advance and automatically yields the smallest bound over p ∈ [2, ∞].Source: Error analysis · p. 6
  • As ϵ → 0, the convergence rates approach their stated limits but the constants diverge; under appropriate weight conditions, the constants and norm of f can be bounded independently of s.The remaining error dependence on s is through ξ_s^p.Source: Error analysis · p. 6
  • If f is twice as smooth as the reproducing kernel with smoothness α, a generalization of recent work yields double the rates.The resulting norm has smoothness 2α and squared and rescaled weights η².Source: Error analysis · p. 6

1.3 Structure of this paper

The paper moves from defining weighted function spaces and constructing quadrature weights to analyzing approximation and method errors, then presents numerical results and comparisons with related strategies.

Sections 2–3

  • Section 2 defines weighted Korobov and Sobolev spaces and relates the Korobov norm of F to the Sobolev norm of f; Section 3 details kernel interpolation, its cost, and exact integration for quadrature weights.Source: §1.3 · p. 6

Section 4

  • Section 4 analyzes full and subsampled lattice interpolation under L2 and L∞ settings, interpolates results for Lp approximation, and constructs one lattice generating vector.Source: §1.3 · pp. 6–7

Section 5

  • Section 5 uses the approximation bounds to establish the method’s final error bound in Theorem 5.1 and doubled rates in Corollary 5.3 when f is twice as smooth as the kernel.Source: §1.3 · p. 7

Sections 6–7

  • Section 6 presents numerical results, and Section 7 compares the method with related strategies and gives concluding remarks.Source: §1.3 · p. 7

2 Function spaces

The section defines weighted periodic Korobov and non-periodic mixed Sobolev spaces, then relates the norm of the periodized integrand to the original integrand’s Sobolev norm.

2.1 Weighted Korobov space of periodic functions

  • The weighted Korobov space comprises periodic functions with absolutely convergent Fourier series, with smoothness α > 1/2 and weights γ_u controlling the importance of variable subsets.For integer α, the space has a reproducing kernel.Source: §2.1 · p. 7

2.2 Weighted Sobolev space of non-periodic functions

  • For integer α ≥ 1, the weighted Sobolev space contains functions with mixed square-integrable derivatives up to order α; β weights encode variable importance, while the Chebyshev density weight controls boundary behavior.With β_u = 1 and constant density weight ω ≡ 1, its norm is the classical mixed Sobolev norm.Source: §2.2 · p. 8

2.3 Relating the Korobov norm of F to the Sobolev norm of f

  • The periodization F = f ◦ ψ has Korobov norm bounded by f’s weighted Sobolev norm, and the norm inequality is sharp for constant functions.The derivative relationship is established using the multivariate Faà di Bruno chain rule.Source: §2.3 · pp. 8, 10

3 Kernel interpolation and integration

The section turns kernel interpolation of the transformed function into a weighted quadrature rule, then describes FFT-based full-lattice construction and methods for computing subsampled weights.

3.1 Kernel interpolation

  • For a Korobov kernel on a full lattice, the circulant kernel matrix enables three-FFT interpolation in O(N log N) arithmetic operations.Kernel-evaluation cost scales linearly in s for product weights and quadratically in s for POD and SPOD weights [25, Table 1].Source: §3.1 · p. 11

3.2 The numerical integration step

  • The analytically computable kernel integrals yield an exactly integrated interpolant as a weighted quadrature rule with optimal weights for μ-weighted integration in the Korobov space.Here, optimal means minimizing the quadrature rule’s worst-case integration error for the specified point set.Source: §3.2 · pp. 11–12
  • After applying ψ, lattice points are related to a Chebyshev lattice and cluster near the boundary, with some points duplicated.Figure 2 shows 73 distinct transformed points from a lattice with N = 144.Source: §3.2 · pp. 11–12, Figure 2 · p. 12

3.3 Implementation for full lattice

  • Full-lattice weights are obtained by solving Kw = κ with FFTs in O(N log N), while total quadrature-rule construction cost is dominated by the CBC construction.The fast CBC algorithm constructs the generating vector; its cost depends on the weights used.Source: §3.3 · pp. 12–13

3.4 Implementation for subsampled lattice

  • Subsampled weights can be computed directly in O(|J|^3), iteratively in O(|J|^2) per iteration, or using FFT-based full-matrix products in O(N log N) per iteration.The iterative methods’ iteration counts depend on the matrix condition number and can be high in practice.Source: §3.4 · p. 14

4 Theory for lattice-based kernel interpolant and subsampling

The section analyzes full-lattice kernel interpolation and subsampled least-squares alternatives, then specializes their error bounds to weighted Korobov spaces using one lattice generating vector.

4.1 Classical lattice algorithm and kernel interpolant

  • The kernel interpolant is optimal among algorithms using the same evaluations, but the analysis avoids a potentially tighter α > 1 bound to retain L2 and L∞ bounds for one lattice vector.The alternative bound requires a different criterion and would produce a different generating vector.Source: §4.1 · pp. 16–17

4.2 Subsampled least squares approximation and subsampled kernel interpolant

  • The random result uses |J| := ⌈12|B|(log |B| + t)⌉, with the stated event holding with probability exceeding 1 − 3 exp(−t).The random techniques enforce logarithmic oversampling; mentioned constructive approaches rely on BSS and have arithmetic complexity O(|B|3N).Source: §4.2 · p. 18

4.3 Specializing to weighted Korobov space

  • The unified weighted Korobov-space theorem recovers prior full-lattice L2 and L∞ rates and subsampled L2 rates, while adding subsampled L∞ and intermediate-Lp results [1].All bounds use the same generating vector; constants may differ from the cited results for that reason.Source: §4.3 · p. 19

A , HKor

This passage specifies the size and sampling scheme for J when t ≥ 4.

A , HKor

  • For t ≥ 4, J is a multiset of uniformly i.i.d. draws from {0, …, N − 1}, with |J| = ⌈12 |B| (log |B| + t)⌉.Source: A , HKor · p. 20

B , HKor

The section derives bounds for full and subsampled lattices by treating p=2 and p=∞ separately, then interpolating for intermediate p with λ=1/2(α−ε), where ε>0 is arbitrarily small.

B , HKor

  • The proof treats p=2 and p=∞ separately for both lattice variants, then uses interpolation to obtain bounds for p∈[2,∞].It chooses λ=1/2(α−ε) with arbitrarily small ε>0.Source: B , HKor · pp. 20–21
  • The connection between |J| and N is expressed using φ(N); when N is a prime power, the expression for φ(N) simplifies that connection.The section also gives bounds relating |J| to |B| and φ(N).Source: B , HKor · p. 21
  • The constants depend on s, α, γ, and ε; assuming the weights γ satisfy (4.24) ensures they are bounded independently of s.The constants for intermediate p also depend trivially on p.Source: B , HKor · p. 21
  • For sufficiently large |J|, (log|J|+t)≤c_δ|J|^δ for arbitrarily small δ>0, simplifying (4.23); the stated result is for prime-power N, and its constants diverge as ε→0.The estimate applies for |J| sufficiently large.Source: B , HKor · pp. 20, 22

5 Error analysis

The error analysis derives lattice-based integration bounds from Lp approximation errors of the transformed integrand, then extends the results to smoother integrands using lower-smoothness kernels.

5 Error analysis

  • The error bound uses Hölder’s inequality and a worst-case Lp approximation bound for F, then relates the transformed Korobov norm to the original integrand’s Sobolev norm.An analogous bound holds with a subsampled lattice.Source: §5 · p. 22

5.1 Main theorem

  • Theorem 5.1’s bounds apply under its weight relation and sufficiently large N, using a component-by-component constructed z and arbitrarily small ϵ > 0.The construction requires S_N(z) ≤ 1/4; the theorem covers f in the associated weighted Sobolev space.Source: §5.1 · p. 22
  • At the subsampled-lattice extremes, the error rate is close to |J|^-(α-1/2) for p = ∞ and |J|^-α (1.110721)^s for p = 2.Taking the infimum over p ∈ [2, ∞] trades the polynomial rate in |J| against exponential dependence on s.Source: §5.1 · p. 23

5.2 Doubling the rate

  • The doubling-rate result extends to intermediate smoothness σ ∈ [1, 2] and general Lp approximation for p ∈ [1, ∞]; analogous results hold for subsampled lattices.Changing weights changes the size of the weighted Korobov unit ball, not membership in the space.Source: §5.2 · p. 23
  • For integrands with twice the kernel smoothness, Corollary 5.3 replaces Theorem 5.1’s full- and subsampled-lattice error bounds; lower-order kernels are often more stable and easier to compute.The corollary states the result for double smoothness, σ = 2.Source: §5.2 · pp. 23–24

6 Numerical experiments

Experiments compare full and subsampled lattice rules for α = 1, 2, 3 in dimensions up to 64. Results show how convergence depends on smoothness and on whether the integrand belongs to the expected function space.

6.1 Example 1: reciprocal function

  • For q = 6, the full lattice achieves the best rate N^-α for all α, while full and subsampled rules behave similarly in dimension 64.Convergence is even faster for α = 1 and for α = 2, 3 in lower dimensions.Source: §6.1 · p. 26, Figure 3 · p. 26
  • At q = 2.5, the best rate N^-α appears only for α = 1 or in lower dimensions; for α = 3, the observed slower rate is 1.25.For α = 2 in dimension 64, convergence approaches rate 2 only at higher N; full and subsampled rules behave comparably with comparable point counts.Source: §6.1 · pp. 26–27, Figure 4 · p. 27

6.2 Example 2: product function

  • For the product function at α = 3 and q = 6, convergence has an intermediate rate rather than the full doubling-the-rate effect.The integrand is not in HSob_s,α+1,β,ω, so the full effect is not expected; results for q = 2.5 were similar and omitted.Source: §6.2 · p. 27, Figure 5 · p. 28

7 Other related strategies

The section situates the paper’s transformed-lattice method among periodization, product-integration, tent-transform, Chebyshev-lattice, and higher-order digital-net approaches. It highlights differences in transformation and weighting, as well as limitations and comparison criteria for these alternatives.

7.1 A common periodization method

  • Common monotone periodizations use smooth increasing endpoint-flat maps such as Korobov or Sidi, whereas this paper’s ψ is nonmonotone and its transformed density has a first-derivative discontinuity.The traditional transformed integrand can be treated as a somewhat smooth periodic function and integrated using the trapezoidal rule.Source: §7.1 · pp. 28–29

7.2 The method of product integration

  • This method generalizes product integration by replacing polynomial interpolation and a weakly singular ν with translated kernels and weight μ; asymptotic positivity of its weights is not known.In product integration, polynomial interpolation approximates the smooth factor while weighted basis integrals account for ν.Source: §7.2 · p. 29

7.3 Integrating the tent-transformed integrand

  • The tent transform traverses [0, 1] twice and yields density 1, but its nonsmooth endpoint behavior precludes arbitrary-order convergence; it can attain the order for smoothness α = 2.It is also described as an extreme periodization and product-integration case, using a periodizing function far from smooth [9].Source: §7.3 · pp. 29–30

7.4 Chebyshev lattices

  • Chebyshev lattices relate to the paper’s transformed lattice points; tent transformation duplicates points, leaving ⌊N/2 + 1⌋ distinct points under the stated coprimality condition.The condition is gcd(z_j, N) = 1 for at least one component z_j; the alternative tent-then-cosine interpretation permits odd N [33].Source: §7.4 · p. 30

7.5 Higher-order digital nets

  • Higher-order digital nets are equal-weight quasi-Monte Carlo competitors that can achieve full convergence order in H_mix^α and admit component-by-component construction.The paper notes that comparisons should also consider implied constants and their growth with dimension and α, as well as reported high-precision arithmetic difficulties; detailed comparison is left for future work.Source: §7.5 · p. 30
Loading 2610.01113v1…