  We generalize fast Gaussian process leave-one-out formulae to multiple-fold
cross-validation, highlighting in turn in broad settings the covariance
structure of cross-validation residuals. The employed approach, that relies on
block matrix inversion via Schur complements, is applied to both Simple and
Universal Kriging frameworks. We illustrate how resulting covariances affect
model diagnostics and how to properly transform residuals in the first place.
Beyond that, we examine how accounting for dependency between such residuals
affect cross-validation-based estimation of the scale parameter. It is found in
two distinct cases, namely in scale estimation and in broader covariance
parameter estimation via pseudo-likelihood, that correcting for covariances
between cross-validation residuals leads back to maximum likelihood estimation
or to an original variation thereof. The proposed fast calculation of Gaussian
Process multiple-fold cross-validation residuals is implemented and benchmarked
against a naive implementation, all in R language. Numerical experiments
highlight the accuracy of our approach as well as the substantial speed-ups
that it enables. It is noticeable however, as supported by a discussion on the
main drivers of computational costs and by a dedicated numerical benchmark,
that speed-ups steeply decline as the number of folds (say, all sharing the
same size) decreases. Overall, our results enable fast multiple-fold
cross-validation, have direct consequences in GP model diagnostics, and pave
the way to future work on hyperparameter fitting as well as on the promising
field of goal-oriented fold design.
