Fast Adaptive Fourier Integration for Spectral Densities of Gaussian Processes
Paul G. Beckman, Christopher J. Geoga
TL;DR
This work addresses the challenge of constructing flexible, positive-definite Gaussian process covariances from arbitrary spectral densities by transforming the problem to the spectral domain via Bochner's theorem. It introduces an adaptive panel Fourier integration framework, enhanced by Gauss-Legendre quadrature and a Nyquist-guided panel strategy, with acceleration from a type-3 nonuniform FFT to achieve high-accuracy kernel evaluations and derivatives for irregularly sampled data. Key contributions include rigorous quadrature and truncation error bounds, robust handling of origin singularities, and automatic differentiation-based derivative computation, enabling efficient gradient-based maximum likelihood estimation for slow-decaying or long-memory spectra (e.g., singular Matérn) and real-world applications such as Doppler LiDAR wind data. The methodology delivers orders-of-magnitude speedups over naive approaches, supports modeling with broad spectral densities, and comes with open-source Julia code, making advanced spectral GP modeling practical for irregular data and complex covariance structures.
Abstract
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 \textit{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.