Linear cost and exponentially convergent approximation of Gaussian Matérn processes on intervals

TL;DR

This work addresses the large-$n$ challenge in Gaussian processes with Matérn covariance by introducing a generally applicable, linear-cost approximation on bounded intervals. It develops an optimal rational approximation of the spectral density that represents the target process as a sum of independent Gaussian Markov processes, enabling a sparse, banded precision structure and straightforward inclusion in Bayesian software. Theoretical guarantees show exponential convergence of the covariance approximation in the order $m$, with exactness when $ u$ is integer, and practical efficiency as a Markov-based representation. Empirical results demonstrate superior accuracy-at-cost compared with state-of-the-art methods across covariance, prediction, and process-approximation tasks, and the method is implemented in the rSPDErSPDE package to integrate with tools like R-INLA. This provides a scalable, accurate alternative for Matérn GP modeling on intervals, with potential extensions to more general domains via a covariance-operator rational framework.

Abstract

The computational cost for inference and prediction of statistical models based on Gaussian processes with Matérn covariance functions scales cubicly with the number of observations, limiting their applicability to large data sets. The cost can be reduced in certain special cases, but there are currently no generally applicable exact methods with linear cost. Several approximate methods have been introduced to reduce the cost, but most of these lack theoretical guarantees for the accuracy. We consider Gaussian processes on bounded intervals with Matérn covariance functions and for the first time develop a generally applicable method with linear cost and with a covariance error that decreases exponentially fast in the order $m$ of the proposed approximation. The method is based on an optimal rational approximation of the spectral density and results in an approximation that can be represented as a sum of $m$ independent Gaussian Markov processes, which facilitates easy usage in general software for statistical inference, enabling its efficient implementation in general statistical inference software packages. Besides the theoretical justifications, we demonstrate the accuracy empirically through carefully designed simulation studies which show that the method outperforms all state-of-the-art alternatives in terms of accuracy for a fixed computational cost in statistical tasks such as Gaussian process regression.

Figures (8)

• Figure 1: True and approximate spectral densities for $\nu=0.4$, $1.4$ and $2.4$, and for different rational approximation orders (top row), and corresponding absolute errors for the rational approximations (bottom row).
• Figure 2: Theoretical error bound from \ref{['total_error_bound']} in Theorem \ref{['ra_bound']} using the approximation for $C_{\{\alpha\}}$ given in Remark \ref{['rem:constant_C_alpha']}.
• Figure 3: Covariance errors for different values of $\nu$ and different orders of approximation.
• Figure 4: Prediction errors, measured as mean squared errors of the posterior means and posterior standard deviations, evaluated on the interval $I = [0, 50]$ for $\sigma_e = 0.1$, $\rho = 2$, various values of $\nu$, and different orders of approximation.
• Figure 5: Prediction errors, measured as mean squared errors of the posterior means and posterior standard deviations, evaluated on the interval $I = [0, 50]$ for $\sigma_e = \sqrt{0.1}$, $\rho = 2$, various values of $\nu$, and different orders of approximation.

Theorems & Definitions (14)

• Theorem 1
• Remark 2
• Proposition 3
• Remark 4
• Proposition 5
• Proposition 6
• Remark 7
• Proposition 8
• Proposition 9
• Proposition 10