Sparse Probabilistic Richardson Extrapolation

1. Problem Statement

The paper addresses the computational challenge of accelerating the convergence of complex numerical simulators (e.g., whole-heart models, fluid dynamics, agent-based systems) to their exact continuum limit $f(0)$.

• Context: Most numerical tasks involve discretization parameters $x$ (e.g., time steps, mesh widths). The simulator output $f(x)$ approximates the true value $f(0)$, with error decreasing as $x \to 0$. However, high-fidelity simulations are computationally expensive (often requiring hundreds of CPU hours).
• The Bottleneck: Traditional extrapolation methods (like Richardson Extrapolation) require fitting polynomials to data. While effective for single parameters, they struggle with multiple parameters ($d > 1$).
• Limitation of Previous Work: A recent probabilistic approach, Gauss–Richardson Extrapolation (GRE), modeled the simulator using a Gaussian Process (GP) to provide uncertainty quantification and convergence acceleration. However, GRE suffers from a super-exponential data requirement in the dimension $d$ of the parameters. Specifically, the minimum number of simulations required scales as $n_{min} \ge (2^d d!)^d$. For $d=4$, this requires $\sim 2 \times 10^{10}$ simulations, rendering it impractical for complex simulators.

2. Methodology: Sparse Probabilistic Richardson Extrapolation (SPRE)

The authors propose SPRE, a novel method that combines probabilistic modeling with the concept of extrapolation sparsity.

A. Core Concept: Extrapolation Dimension

The paper introduces the notion of Extrapolation Dimension, denoted $d_{ext}(f)$.

• Definition: It is the minimal cardinality of an index set $A \subset \mathbb{N}_0^d$ such that the simulator output can be expanded as:
  $$f(x) = f(0) + \sum_{\alpha \in A} \beta_\alpha x^\alpha + R(x)$$
  where $R(x)$ is a residual term that vanishes faster than the expansion terms as $x \to 0$.
• Sparsity Hypothesis: The authors argue that for most complex simulators, the effective number of terms required to describe the error (the size of $A$) is small ($d_{ext}(f) \ll \text{total possible terms}$). This is analogous to sparsity in linear regression.

B. The SPRE Algorithm

SPRE models the simulator $f$ as a Gaussian Process with a specific structure:

1. Mean Function: The mean is a linear combination of monomials from the sparse index set $A$:
   $$E[f(x)] = \sum_{\alpha \in A} \beta_\alpha x^\alpha$$
   The coefficients $\beta_\alpha$ are assigned an improper flat prior and integrated out.
2. Covariance Function: A symmetric, positive-definite kernel $k(x, x')$ models the residual $R(x)$.
3. Posterior Estimation: The posterior mean at $x=0$ (the estimate of $f(0)$) is computed analytically. Crucially, the method requires only $n \ge d_{ext}(f)$ data points to be well-defined, rather than the super-exponential number required by GRE.

C. Implementation Details

• Data-Driven Selection: Since the index set $A$ is unknown, SPRE uses Leave-One-Out Cross-Validation (LOOCV) to select the optimal set $A$ and kernel hyperparameters. It employs a stepwise search, adding monomials that reduce the LOOCV error.
• Experimental Design: The method supports cost-constrained experimental design to select the most informative simulation configurations $\{x_i\}$ within a fixed computational budget.

3. Key Contributions

1. Theoretical Breakthrough: The paper establishes that the data requirement for probabilistic extrapolation can be polynomial in the dimension $d$ (specifically dependent on $d_{ext}(f)$) rather than super-exponential.
2. Convergence Guarantees:
   • Theorem 3: Proves that SPRE achieves convergence acceleration comparable to classical Multivariate Richardson Extrapolation (MRE) but with probabilistic uncertainty quantification.
   • Theorem 4: Demonstrates that if the residual $R(x)$ possesses additional smoothness, SPRE can achieve further acceleration beyond MRE by leveraging the kernel's properties, without needing prior knowledge of the smoothness order.
3. Sparsity Discovery: The authors provide empirical evidence that "extrapolation sparsity" is a generic phenomenon in complex simulators (e.g., contact mechanics, flocking models), meaning $d_{ext}(f)$ is typically small (e.g., 2–4) even when the total number of parameters $d$ is larger.
4. Practical Framework: The method unifies numerical analysis (polynomial expansion) with statistical learning (GPs), offering uncertainty quantification, model criticism, and experimental design capabilities absent in classical methods.

4. Results and Empirical Assessment

The authors validated SPRE on three increasingly complex case studies:

1. Two Spheres 3D Model: A simulation of colliding spheres with contact constraints ($d=3$).
   • Result: SPRE achieved convergence rates similar to MRE and significantly better than raw estimates. It correctly identified a sparse index set ($d_{ext} \approx 2-4$). GRE failed to accelerate convergence due to insufficient data ($n=64$ was too small for $d=3$).
2. Five Shapes 3D Model: A more complex collision scenario with five interacting bodies ($d=3$).
   • Result: SPRE again outperformed raw estimates and GRE. The method remained robust even as the simulation complexity increased.
3. Agent-Based Flocking Model: A simulation of 60 interacting agents with nonlinear dynamics ($d=3$).
   • Result: Despite the highly nonlinear and emergent behavior, SPRE successfully accelerated convergence. The estimated index set suggested a first-order expansion in time step and repulsion parameters, confirming sparsity.

Key Observations:

• Accuracy: SPRE estimates were often two orders of magnitude more accurate than raw low-fidelity simulations.
• Uncertainty Quantification: The predictive intervals provided by SPRE were well-calibrated (the true value fell within the error bars).
• Kernel Choice: The "white noise" kernel and Matérn-1/2 kernel performed best, likely because the complex simulators did not satisfy the high smoothness assumptions required by Gaussian kernels.

5. Significance

• Feasibility of High-Dimensional Extrapolation: SPRE removes the "curse of dimensionality" barrier that previously made probabilistic extrapolation impractical for multi-parameter simulators.
• Cost Reduction: By drastically reducing the number of required simulations (from super-exponential to polynomial), SPRE enables high-accuracy inference in fields like cardiac modeling, climate science, and engineering where simulations are prohibitively expensive.
• Unified Framework: It bridges the gap between deterministic numerical analysis (Richardson extrapolation) and probabilistic numerics, allowing for rigorous uncertainty quantification in computational science.
• Generalizability: The discovery that complex, nonlinear simulators often exhibit low extrapolation dimension suggests this method is widely applicable across scientific domains.

In summary, Sparse Probabilistic Richardson Extrapolation represents a paradigm shift in how numerical errors are handled, transforming extrapolation from a data-hungry, deterministic procedure into a data-efficient, probabilistic tool capable of accelerating complex simulations with rigorous uncertainty bounds.