Equispaced Fourier Representations for Efficient Gaussian Process Regression from a Billion Data Points
Philip Greengard, Manas Rachh, Alex H. Barnett · SIAM/ASA Journal on Uncertainty Quantification · 2025
Abstract. We introduce a Fourier-based fast algorithm for Gaussian process regression in low dimensions. It approximates a translationally invariant covariance kernel by complex exponentials on an equispaced Cartesian frequency grid of [Formula: see text] nodes. This results in a weight-space [Formula: see text] system matrix with Toeplitz structure, which can thus be applied to a vector in [Formula: see text] operations via the fast Fourier transform (FFT), independent of the number of data points [Formula: see text]. The linear system can be set up in [Formula: see text] operations using nonuniform FFTs. This enables efficient massive-scale regression via an iterative solver, even for kernels with fat-tailed spectral densities (large [Formula: see text]). We provide bounds on both kernel approximation and posterior mean errors. Numerical experiments for squared-exponential and Matérn kernels in one, two, and three dimensions often show 1–2 orders of magnitude acceleration over state-of-the-art rank-structured solvers at comparable accuracy. Our method allows two-dimensional Matérn-[Formula: see text] regression from [Formula: see text] data points to be performed in two minutes on a standard desktop, with posterior mean accuracy [Formula: see text]. This opens up spatial statistics applications 100 times larger than previously possible.