  Gaussian process modulated Poisson processes provide a flexible framework for
modelling spatiotemporal point patterns. So far this had been restricted to one
dimension, binning to a pre-determined grid, or small data sets of up to a few
thousand data points. Here we introduce Cox process inference based on Fourier
features. This sparse representation induces global rather than local
constraints on the function space and is computationally efficient. This allows
us to formulate a grid-free approximation that scales well with the number of
data points and the size of the domain. We demonstrate that this allows MCMC
approximations to the non-Gaussian posterior. We also find that, in practice,
Fourier features have more consistent optimization behavior than previous
approaches. Our approximate Bayesian method can fit over 100,000 events with
complex spatiotemporal patterns in three dimensions on a single GPU.
