  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 $M$
nodes. This results in a weight-space $M\times M$ system matrix with Toeplitz
structure, which can thus be applied to a vector in ${\mathcal O}(M \log{M})$
operations via the fast Fourier transform (FFT), independent of the number of
data points $N$. The linear system can be set up in ${\mathcal O}(N + M
\log{M})$ operations using nonuniform FFTs. This enables efficient
massive-scale regression via an iterative solver, even for kernels with
fat-tailed spectral densities (large $M$). We provide bounds on both kernel
approximation and posterior mean errors. Numerical experiments for
squared-exponential and Mat\'ern 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 2D
Mat\'ern-$\mbox{$\frac{3}{2}$}$ regression from $N=10^9$ data points to be
performed in 2 minutes on a standard desktop, with posterior mean accuracy
$10^{-3}$. This opens up spatial statistics applications 100 times larger than
previously possible.
