Fast Adaptive Fourier Integration for Spectral Densities of Gaussian Processes

The specification of a covariance function is of paramount importance when employing Gaussian process models, but the requirement of positive definiteness severely limits those used in practice. Designing flexible stationary covariance functions is, however, straightforward in the spectral domain, where one needs only to supply a positive and symmetric spectral density. In this work, we introduce an adaptive integration framework for efficiently and accurately evaluating covariance functions and their derivatives at irregular locations directly from any continuous, integrable spectral density. In order to make this approach computationally tractable, we employ high-order panel quadrature, the nonuniform fast Fourier transform, and a Nyquist-informed panel selection heuristic, and derive novel algebraic truncation error bounds which are used to monitor convergence. As a result, we demonstrate several orders of magnitude speedup compared to naive uniform quadrature approaches, allowing us to evaluate covariance functions from slowly decaying, singular spectral densities at millions of locations to a user-specified tolerance in seconds on a laptop. We then apply our methodology to perform gradient-based maximum likelihood estimation using a previously numerically infeasible long-memory spectral model for wind velocities below the atmospheric boundary layer.