Gaussian process regression is widely used in geostatistics, time-series analysis, and machine learning. It infers an unknown continuous function in a principled fashion from noisy measurements at $N$ scattered data points. The prior on the function is Gaussian, with covariance given by some user-chosen translationally invariant kernel. Yet $N$ has been limited to about $10^6$, even with modern low-rank methods. Focusing on low spatial dimension (1--3), we present a GP regression method using kernel approximation by an equispaced quadrature grid in the Fourier domain. This enables the iterative solution of a smaller Toeplitz linear system, exploiting both the FFT and the nonuniform FFT to give ${\mathcal O}(N)$ cost. The result is often one to two orders of magnitude faster than state of the art methods, and enables cheap massive-scale regressions. For example, for a 2D Matérn-3/2 kernel and $N = 10^9$ points, the posterior mean function is found to 3-digit accuracy in two minutes on a desktop.
Joint work with Philip Greengard (Columbia) and Manas Rachh (Flatiron Institute)