  One of the most compelling features of Gaussian process (GP) regression is
its ability to provide well-calibrated posterior distributions. Recent advances
in inducing point methods have sped up GP marginal likelihood and posterior
mean computations, leaving posterior covariance estimation and sampling as the
remaining computational bottlenecks. In this paper we address these
shortcomings by using the Lanczos algorithm to rapidly approximate the
predictive covariance matrix. Our approach, which we refer to as LOVE (LanczOs
Variance Estimates), substantially improves time and space complexity. In our
experiments, LOVE computes covariances up to 2,000 times faster and draws
samples 18,000 times faster than existing methods, all without sacrificing
accuracy.
