Fast and Provably Accurate Sequential Designs using Hilbert Space Gaussian Processes

1. Problem Statement

Sequential design with Gaussian Process (GP) surrogates is a standard approach for optimizing expensive simulation experiments. A critical component of this process is the acquisition function, which determines where to sample next. The Integrated Mean Squared Error (IMSE) is a highly effective acquisition function for global surrogate modeling, as it aims to minimize the average prediction error over the entire design domain $\Omega$.

However, implementing IMSE-based sequential design faces a significant computational bottleneck:

• Intractable Integrals: The IMSE criterion requires evaluating integrals of the form $\int_\Omega k(x, t)k(x, x') dP(x)$, where $k$ is the kernel function.
• Lack of Closed Forms: For most kernel functions (e.g., Matérn, Generalized Wendland), these integrals lack closed-form expressions.
• Current Limitations: Existing efficient implementations (e.g., in the hetGP package) rely on pre-computed analytic integrals, restricting them to a narrow set of kernels (Gaussian, Matérn-3/2, Matérn-5/2) with uniform measures.
• Computational Cost: Numerical integration (Monte Carlo or quadrature) within an optimization loop is prohibitively expensive, especially as the number of samples $N$ increases, because the integrand changes at every iteration.

2. Methodology

The authors propose a novel framework using Hilbert Space Gaussian Processes (HSGP) to approximate the IMSE acquisition function, enabling closed-form evaluation for a broad class of kernels.

A. HSGP Approximation

The method approximates the stationary kernel $k(x, x')$ using a truncated eigenbasis expansion of the Laplace operator on a bounded domain $(-L, L)^d$ (where $L > B$ and the design domain is $\Omega = (-B, B)^d$).

• Basis Functions: The eigenfunctions are sine functions: $\phi_j(x) = L^{-d/2} \prod \sin(\frac{\pi j_k}{2L}(x_k + L))$.
• Spectral Density: The kernel is approximated by truncating the series:
  $$\hat{k}_m(x, x') = \sum_{j \in [m]^d} S\left(\frac{\pi j}{2L}\right) \phi_j(x) \phi_j(x')$$
  where $S(\omega)$ is the spectral density of the kernel (available in closed form for Gaussian, Matérn, and Generalized Wendland kernels).
• Closed-Form Integration: Because the basis functions are trigonometric, the integral of the product of kernels required for IMSE can be computed analytically. This transforms the IMSE calculation from a numerical integration problem into a matrix-vector multiplication problem.

B. The $\gamma$-Stabilizing Strategy

To ensure numerical stability and prevent the algorithm from selecting points too close to existing data (which causes the denominator in the IMSE formula to approach zero), the authors introduce a $\gamma$-stabilizing constraint:

• The next candidate point $t$ is selected from a restricted domain $\Omega_\gamma = \{t \in \Omega : \text{dist}(t, X_N) \geq \gamma h_N\}$, where $h_N$ is the fill distance.
• This ensures the resulting design is quasi-uniform, which is crucial for theoretical convergence guarantees.

C. Algorithm

The proposed sequential design algorithm (Algorithm 2) iteratively:

1. Fits a GP model.
2. Computes the approximate IMSE, $\widehat{\text{IMSE}}_m(t)$, using the closed-form expression derived from the HSGP expansion.
3. Selects the next point by maximizing $\widehat{\text{IMSE}}_m(t)$ over the constrained domain $\Omega_\gamma$.
4. Updates the dataset and repeats.

3. Key Contributions

Methodological Innovation

• General Applicability: Unlike previous methods restricted to specific kernels, the HSGP-IMSE approach applies to any kernel with a closed-form spectral density (including Matérn kernels with arbitrary smoothness $\nu$ and Generalized Wendland kernels).
• Computational Efficiency: The method reduces the complexity of evaluating the acquisition function. While standard approaches require $O(N^2)$ or higher operations for matrix updates involving pre-computed integrals, the HSGP approach achieves a cost of $o(N^2)$ (specifically $O(NM)$ where $M$ is the number of basis functions, and $M \ll N$).

Theoretical Guarantees

• Sharp Error Bounds: The authors derive non-asymptotic, explicit error bounds for both the kernel approximation and the resulting IMSE acquisition function.
  • For Gaussian kernels, the error decays exponentially with respect to the domain parameter $L$ and the truncation level $m$.
  • For Matérn kernels, the error decays polynomially with $m$ at a rate determined by the smoothness $\nu$.
• Convergence Rates: The paper establishes that the approximation error of the acquisition function converges at rates comparable to the optimal convergence rates of the underlying GP regression, provided the HSGP parameters ($L, m$) are tuned correctly relative to the sample size $N$.
• Stability: The $\gamma$-stabilizing strategy is proven to guarantee quasi-uniformity of the design points, ensuring the denominator in the IMSE formula remains bounded away from zero.

4. Results

Numerical Experiments

The authors conducted extensive experiments in 1D and 2D domains using various kernels (Gaussian, Matérn-3/2, Matérn-5/2, Matérn-2, and Generalized Wendland).

1. Approximation Accuracy:

   • Visual comparisons (Figures 1 & 2) show that the HSGP approximation $\widehat{\text{IMSE}}_m(t)$ almost perfectly overlaps with the numerically integrated reference IMSE(t).
   • The geometric features (modes, ridges, basins) of the acquisition landscape are preserved, ensuring that the optimization of the surrogate yields the same maximizers as the exact criterion.

2. Sequential Design Performance:

   • Predictive Accuracy: The HSGP-IMSE method achieved lower Root Mean Square Error (RMSE) and faster reduction in posterior variance compared to both the standard hetGP (IMSPE) baseline and Latin Hypercube Sampling (LHS) baselines.
   • Computational Time: The runtime of HSGP-IMSE was comparable to the highly optimized hetGP implementation for supported kernels and significantly faster than numerical integration approaches. For unsupported kernels (e.g., Matérn-2, Generalized Wendland), it outperformed LHS in accuracy while maintaining similar computational costs.

3. Robustness: The method performed well across different noise levels (nugget effects) and kernel types, demonstrating its versatility.

5. Significance

This work addresses a fundamental limitation in the field of sequential design and surrogate modeling:

• Breaking the Kernel Barrier: It liberates IMSE-based sequential design from the restriction to a few specific kernel functions, allowing practitioners to use kernels better suited to their specific physical or statistical problems (e.g., anisotropic or less smooth kernels).
• Scalability: By providing a closed-form solution for the acquisition function, it makes global sequential design computationally feasible for larger datasets and more complex kernels.
• Theoretical Rigor: The paper provides the first sharp, non-asymptotic error bounds for HSGP approximations in the context of sequential design, offering a solid theoretical foundation for parameter selection ($L$ and $m$).

In summary, the proposed HSGP-IMSE framework offers a fast, accurate, and theoretically grounded alternative to existing methods, enabling efficient global emulation for a much broader class of Gaussian Process models.