Parametric Numerical Integration

• Parametric numerical integration is a technique to approximate families of integrals dependent on variable parameters, ensuring rigorous error control and reduced computational cost.
• The Magic Point Empirical Interpolation method uses an offline-online decomposition to construct efficient quadrature schemes with exponential convergence for analytic integrands.
• Recent advances using machine learning and neural operators enhance accuracy for high-dimensional, complex integrals, offering significant speed-ups over traditional quadrature methods.

Parametric numerical integration refers to the rigorous approximation and efficient numerical evaluation of families of integrals that depend on one or more external parameters, typically of the form

$$I(\theta) = \int_{\Omega} h(x;\theta)\,dx, \quad \theta \in \mathcal{P} \subset \mathbb{R}^d,$$

where $h$ is often available in closed or computable form, and the goal is to produce accurate, often fast, approximations for a large set of parameters $\theta$. Such integrals arise ubiquitously in scientific computing, engineering design, uncertainty quantification, financial modeling (e.g., option pricing), and machine learning. The core challenge is to address the computational complexity of evaluating many integrals for variable parameters, while rigorously controlling error and computational cost.

1. Problem Classes and Fundamental Challenges

Parametric numerical integration spans both finite- and infinite-dimensional parameter spaces and a variety of application domains. The main classes include:

• Low- to moderate-dimensional parameterized quadrature (e.g., Fourier-based option pricing, statistical functionals): Here, $h(x;\theta)$ is smooth or analytic, and $|\mathcal{P}|$ is small to moderate. The focus is on minimal-cost evaluation across $\theta$ (Gaß et al., 2015, Gaß et al., 2015).
• High- and infinite-dimensional parametric integrals: Arising in uncertainty quantification for PDEs with random inputs, these feature $h(x;\theta)$ defined on a space with $d \gg 1$ or $d = \infty$ (Guth et al., 2022, Dũng, 2019, Dick et al., 2014).
• Geometric parametric integration over domains specified by parameter-dependent or high-order parametric surfaces, such as in trimmed isogeometric geometries (Antolin et al., 2021, Chin et al., 2020).
• Data-driven and operator learning approaches, where the mapping $\theta \mapsto I(\theta)$ is learned via machine learning surrogates (Leitao et al., 12 Dec 2025, Maître et al., 2022).

The critical challenge is reduction of computational complexity as a function of the number of queries in $\theta$, dimension $d$, and target error. Other key issues include analytic regularity, adaptivity, the role of randomized versus deterministic methods, and efficient treatment of singular, oscillatory, or near-singular integrands.

2. Magic Point Empirical Interpolation: Offline-Online Decomposition

The Magic Point Empirical Interpolation method provides an explicit, constructive approach to parametric integration for low- to moderate-dimensional, analytic problems (Gaß et al., 2015, Gaß et al., 2015). The method consists of two stages:

• Offline phase: A greedy algorithm selects "magic" points $x_1^*, \ldots, x_N^*$ in $\Omega$ and builds a corresponding basis $\{q_m\}_{m=1}^N$ for the space of possible integrands $U = \{ h(\cdot;\theta): \theta \in \mathcal{P} \}$. At each step, the function in $U$ with maximal residual error is identified, and its maximal residue over $\Omega$ sets the next interpolation point. Basis construction ensures lower-triangular interpolation matrices, and the associated quadrature weights $\omega_m(\theta)$ are computed by integrating (parameter-independent) linear combinations of the basis.
• Online phase: For any new $\theta$, evaluate $h(x_m^*;\theta)$ at the magic points and compute $$I_N(\theta) = \sum_{m=1}^N \omega_m(\theta) h(x_m^*;\theta)$$ with cost $\mathcal{O}(N)$.

Rigorous exponential convergence holds if $x \mapsto h(x;\theta)$ is analytic in a complex neighborhood of $\Omega$, with error $$|I(\theta) - I_N(\theta)| \leq C N e^{-(\alpha-\ln 4)N}$$ uniformly over $\theta$, where $\alpha$ depends on the width of the strip of analyticity. Numerical experiments confirm that in practice, $N \asymp 20$–$50$ suffices for $10^{-8}$ accuracy in prototypical Fourier-based finance applications, and the approach is markedly more efficient than conventional quadrature or COS methods (Gaß et al., 2015).

Application domain	Offline cost (per model)	Online cost (per $\theta$)	Typical $N$ for $10^{-8}$ error
Option pricing	$\sim$ few hours	$\mathcal{O}(N)$ ops, closed-form $h$	$20$–$50$
Fourier inversion	(Identical process)	$\mathcal{O}(N)$	$20$–$40$

3. High-Dimensional and Infinite-Dimensional Approaches

For parametric integrals over high- or infinite-dimensional parameter domains, two principal approaches provide tractable complexity and rigorous error control:

• Dimension truncation and Taylor analysis: If $g(y)$ (often $g(y) = G(u(y))$ for PDE solution $u$ and quantity-of-interest $G$) is analytic and satisfies weighted $\ell^p$-summability of derivatives, truncating to the first $s$ coordinates gives error $\mathcal{O}(s^{-2/p+1})$, with $p \in (0,1)$ determined by decay of coefficient norms (Guth et al., 2022).
• Sparse-grid and QMC methods: For isotropic or anisotropic regularity, sparse-grid Hermite interpolation (Dũng, 2019) or higher-order QMC (Dick et al., 2014) are effective. Under suitable summability and holomorphy conditions, QMC achieves dimension-independent rates of $\mathcal{O}(N^{-1/p})$ for sample size $N$, where $p$ quantifies regularity, using SPOD (smoothness-driven product and order-dependent) weights.
• Hybrid error balancing (truncation, cubature, discretization): In uncertainty quantification workflows, the overall accuracy is determined by balancing dimension truncation, spatial discretization (e.g., finite element mesh size $h$), and cubature error, with

$s^{-2/p+1} \sim h^r \sim N^{-1/p} \sim \varepsilon$

for target tolerance $\varepsilon$ (Guth et al., 2022, Dick et al., 2014).

4. Machine Learning-Based and Neural-Operator Methods

Parametric numerical integration via neural networks, especially with differential information, has recently demonstrated significant empirical advantages in sample efficiency and scalability, particularly for integrals with high-dimensional parameter spaces or vector-valued outputs (Leitao et al., 12 Dec 2025, Maître et al., 2022). Approaches include:

• Surrogate regression: Standard feedforward ANNs are trained to learn the map $\theta \mapsto I(\theta)$, using single-sample Monte Carlo targets for supervised regression.
• Differential machine learning (DML): Loss functions are augmented to include both value and gradient (with respect to $\theta$) information. Sampling of MC gradients and analytic differentiation enables unbiased simultaneous estimation of $I(\theta)$ and $\nabla_\theta I(\theta)$, substantially reducing variance and accelerating convergence.
• Automatic antiderivative fitting: For smooth $f(x;\theta)$, neural networks are trained to approximate an $n$-fold primitive, with the loss enforcing the correct mixed partials. The integral $I(\theta)$ is then exactly recovered via the generalized fundamental theorem of calculus from boundary evaluations (Maître et al., 2022).

These methods yield uniformly lower mean squared error and improved sample efficiency relative to value-only surrogates, and post-training inference cost is reduced to milliseconds per query, with speed-ups of $40$–$125\times$ over standard numerical integration software (Leitao et al., 12 Dec 2025, Maître et al., 2022).

5. Geometric and High-Order Parametric Integration on Curved Domains

Parametric quadrature for domains bounded by parametric curves or surfaces (planar curves, polyhedral patches, trimmed NURBS, etc.) demands flexible domain decomposition and rigorous treatment of boundary or singularity effects:

• Scaled Boundary Cubature (SBC) (Chin et al., 2020): Parameterizes star-convex or nonconvex regions via a scaling map from a center $x_0$, enabling tensor-product quadrature on a reference domain. Special transforms handle homogeneous integrands and weak or near singularities, with analytic error control for polynomial and non-polynomial cases.
• Folded decomposition (Antolin et al., 2021): Curved polyhedra are decomposed into generally pyramidal cells via arbitrary "seed" vertices, resulting in integration subdomains where the volume map Jacobian may change sign ("folded cells"). Rigorous analysis shows that the accuracy and convergence of tensor-product Gauss quadrature is unaffected by sign changes, provided the mapped integrand remains smooth and is suitably extended outside the principal domain if necessary.

These methods enable robust, high-order quadrature for trimmed geometric domains, directly accommodating real-world complex geometries and sharp features.

6. Complexity Theory and Adaptivity in Parametric Integration

The information-based complexity of parametric numerical integration, including the minimal errors of deterministic and randomized algorithms, is characterized by the smoothness of the input space and the power of adaptation (Heinrich, 2023). Main findings include:

• Randomized $n$-th minimal error rates: For mean computation over $L_p$ spaces, deterministic (and non-adaptive randomized) algorithms achieve error rates $\mathcal{O}(n^{-(1-1/p)})$ (up to logarithmic factors), while adaptive randomized methods achieve $\mathcal{O}(n^{-1/2})$ whenever $p > 2$ and $p < q$.
• Benefit of adaptation: There exists a strict, unbounded polynomial factor by which adaptively sampling in the randomized setting outperforms non-adaptive approaches for $p > 2$. This resolves an open problem regarding the necessity of adaptation for optimal convergence even in linear parametric integration settings.
• Extension to Sobolev-space and continuous problems: These results extend to infinite-dimensional Sobolev models via suitable discretizations, ensuring transfer of optimal rates and adaptation advantages.

7. Specialized High-Frequency and Singular Kernel Quadrature

For integrals involving sharply peaked, oscillatory, or singular kernels, such as $I_\sigma = \int_a^b f(x)\,G_\alpha(x)\,dx$ with $G_\alpha(x)$ a sharply localized Gaussian ($\alpha \gg 1$), specialized mesh and interpolation strategies are required (Ma et al., 2018). Techniques include:

• Graded meshes: Subdivide the integration domain with subintervals increasingly clustered near singularities or localization centers via geometric or adaptive meshes, balancing quadrature errors.
• Panelwise Chebyshev interpolation: On each subinterval, interpolate $f$ at Chebyshev points and compute moments of kernel-weighted polynomials exactly, enabling polynomial or exponential convergence in $\alpha$.
• Extensions: The method generalizes to oscillatory or strongly decaying weights, yielding black-box quadrature rules whose work scales polylogarithmically in the kernel parameter or singularity strength.

---

Through these analytic, algorithmic, and computational developments, parametric numerical integration provides a unified theoretical and practical framework for efficient, high-accuracy quadrature in the context of parameterized physical, financial, and statistical modeling, with extensions spanning machine learning surrogates and high-dimensional probabilistic computation. The continued refinement of hybrid, data-driven, and problem-specific algorithms, and the exploitation of regularity and analytics, remains a focal direction for further advances.