  Gaussian processes (GPs) offer a flexible class of priors for nonparametric
Bayesian regression, but popular GP posterior inference methods are typically
prohibitively slow or lack desirable finite-data guarantees on quality. We
develop an approach to scalable approximate GP regression with finite-data
guarantees on the accuracy of pointwise posterior mean and variance estimates.
Our main contribution is a novel objective for approximate inference in the
nonparametric setting: the preconditioned Fisher (pF) divergence. We show that
unlike the Kullback--Leibler divergence (used in variational inference), the pF
divergence bounds the 2-Wasserstein distance, which in turn provides tight
bounds the pointwise difference of the mean and variance functions. We
demonstrate that, for sparse GP likelihood approximations, we can minimize the
pF divergence efficiently. Our experiments show that optimizing the pF
divergence has the same computational requirements as variational sparse GPs
while providing comparable empirical performance--in addition to our novel
finite-data quality guarantees.
