Scalable Gaussian Process Flows

Outside of the linear-Gaussian regime, conditional sampling from Gaussian processes (GPs) is challenging. Recent sampling algorithms provide machinery allowing conditioning on arbitrary statements, including non-linear and non-Gaussian conditions. However, important bottlenecks remain: typically the additional flexibility comes at an increased computational cost. This is seen prominently in (Moss et al. (2026)), which vastly increases the scope of GP modelling capabilities, but at the cost of executing an expensive iterative diffusion model. In this paper, we alleviate two significant drawbacks of (Moss et al. (2026)) by (1) introducing kernel approximations that improve scalability to high-resolution domains and (2) allowing for kernel hyperparameter optimisation. We demonstrate the utility on two important downstream tasks: probabilistic downsampling using areal summary statistics, and inference of PDE solutions on irregular domains from noisy observations.