Embarrassingly parallel inference for Gaussian processes

Training Gaussian process-based models typically involves an $ O(N^3)$ computational bottleneck. Popular methods for overcoming this matrix inversion problem cannot adequately model all types of latent functions, and are often not parallelizable. We present an embarrassingly parallel method that takes advantage of inverting block diagonal approximations, while maintaining much of the expressivity of a full covariance matrix. By using importance sampling to average over different realizations of low-rank GP approximations, we ensure our algorithm is both asymptotically unbiased and embarrassingly parallel. We show comparable or improved performance over competing methods, on a range of synthetic and real datasets.