  Gaussian processes are flexible probabilistic regression models which are
widely used in statistics and machine learning. However, a drawback is their
limited scalability to large data sets. To alleviate this, we consider
full-scale approximations (FSAs) that combine predictive process methods and
covariance tapering, thus approximating both global and local structures. We
show how iterative methods can be used to reduce the computational costs for
calculating likelihoods, gradients, and predictive distributions with FSAs. We
introduce a novel preconditioner and show that it accelerates the conjugate
gradient method's convergence speed and mitigates its sensitivity with respect
to the FSA parameters and the eigenvalue structure of the original covariance
matrix, and we demonstrate empirically that it outperforms a state-of-the-art
pivoted Cholesky preconditioner. Further, we present a novel, accurate, and
fast way to calculate predictive variances relying on stochastic estimations
and iterative methods. In both simulated and real-world data experiments, we
find that our proposed methodology achieves the same accuracy as Cholesky-based
computations with a substantial reduction in computational time. Finally, we
also compare different approaches for determining inducing points in predictive
process and FSA models. All methods are implemented in a free C++ software
library with high-level Python and R packages.
