Fast Calculation of Gaussian Process Multiple-Fold Cross-Validation Residuals and their Covariances

Abstract

We generalize fast Gaussian process leave-one-out formulas to multiple-fold cross-validation, highlighting in turn the covariance structure of cross-validation residuals in simple and universal kriging frameworks. We illustrate how resulting covariances affect model diagnostics. We further establish in the case of noiseless observations that correcting for covariances between residuals in cross-validation-based estimation of the scale parameter leads back to maximum likelihood estimation. Also, we highlight in broader settings how differences between pseudo-likelihood and likelihood methods boil down to accounting or not for residual covariances. The proposed fast calculation of cross-validation residuals is implemented and benchmarked against a naive implementation, all in R. Numerical experiments highlight the substantial speed-ups that our approach enables. However, as supported by a discussion on main drivers of computational costs and by a numerical benchmark, speed-ups steeply decline as the number of folds (say, all sharing the same size) decreases. An application to a contaminant localization test case illustrates that the way of grouping observations in folds may affect model assessment and parameter fitting compared to leave-one-out. Overall, our results enable fast multiple-fold cross-validation, have consequences in model diagnostics, and pave the way to future work on hyperparameter fitting as well as on goal-oriented fold design. Supplementary materials for this article are available online.

Keywords:

• Cross-validation
• Diagnostics
• Gaussian process
• Hyperparameter estimation
• Universal kriging
• Woodbury formula

1 Introduction and Notation

Gaussian process (GP) models are at the heart of a number of prominent methods in spatial statistics, machine learning, and beyond. Properly validating such models and fitting underlying hyperparameters thus appear as crucial and impactful endeavors. It is often the case that, for reasons of data scarcity or other, the latter need to be conducted based on a unique dataset. Cross-validation is commonly used in such a context, not only in GP modeling but also as a general approach to assess statistical models and tune parameters in a vast class of prediction algorithms. Fast leave-one-out cross-validation formulas are known for some models and approaches including GP prediction. However, leave-one-out is known to suffer from some pitfalls, and multiple-fold cross-validation is increasingly preferred over it in broader contexts, notably in the realm of machine learning.

We focus here on cross-validation in the context of GP modeling, and more specifically on adapting fast formulas for leave-one-out (LOO) to multiple-fold cross-validation and exploring how these could be efficiently exploited in topics such as model diagnostics and covariance hyperparameter fitting. One of the essential bricks of the proposed approaches turns out to be closed-form formulas for the covariance structure of cross-validation residuals. Our main results are presented both in the simple kriging and universal kriging frameworks, also with links to ridge regression. They contribute to a lively stream of research on cross-validation that pertains to GP modeling, and also to further related classes of statistical models.

Leave-one-out (LOO) and multiple-fold cross-validation in the GP framework were tackled as early as in Dubrule (1983), and some further references as well as a more recent account of their use in machine learning is given in chap. 5 of Rasmussen and Williams (2006). On a different note, cross-validation has led to a number of investigations in the context of spline modeling, for instance in Golub, Heath, and Wahba (1979) where generalized cross-validation was used for the selection of a ridge (regularization) parameter. In the domain of GP modeling for computer experiments, cross-validation was tackled as early as in Currin et al. (1988) (with a pioneering discussion on links to MLE), and leave-one-out was investigated in Bachoc (2013) and found to deliver a valuable alternative to MLE in misspecified cases. A general survey of cross-validation procedures for model selection is presented in Arlot and Celisse (2010). Cross-validation for model selection is further tackled in Zhang and Yang (2015). Le Gratiet, Cannamela, and Iooss (2015) proposed cokriging-based sequential design strategies using fast cross-validation for multi-fidelity computer codes. Fischer (2016) studied model selection for GP regression by Approximation Set Coding. Martino, Laparra, and Camps-Valls (2017) suggested new probabilistic cross-validation Estimators for Gaussian process regression. Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure were investigated in Roberts et al. (2017). New prediction error estimators were defined in Rabinowicz and Rosset (2020a), where novel model selection criteria in the same spirit as AIC and Mallow’s C_p were suggested. LOO cross-validation for Bayesian model comparison in large data was recently tackled in Magnusson et al. (2020). Kaminsky et al. (2021) introduced an efficient batch multiple-fold cross-validation Voronoi adaptive sampling technique for global surrogate modeling. Rabinowicz and Rosset (2020b) introduced a bias-corrected cross-validation estimator for correlated data. Bates, Hastie, and Tibshirani (2023) tackles the question of what one actually does estimate in cross-validation. In contrast with the random design settings considered in the latter as in many theoretical works pertaining to cross-validation, our baseline settings here will assume a fixed design.

LOO cross-validation consists in predicting at each of the observation locations when (temporarily) removing the corresponding observation from the dataset, and then comparing obtained predictions to left-out observations. As a first example, we consider a case of GP prediction for a one-dimensional test function (which definition is recalled in Section A.2 of the supplementary material) observed at a regular design, here a 10-point subdivision of $[0,1]$, and a GP model assuming a Matérn 5/2 stationary covariance kernel. The test function is represented in black in Figure 1, along with the GP predictor (dashed blue line) and LOO predictions at the design locations (red stars). Modeling is performed with the DiceKriging package Roustant, Ginsbourger, and Deville (2012) and cross-validation relies on a fast implementation using the cv function newly available in the DiceKriging package (version 1.6.0). This simple example in a misspecified case illustrates how LOO residuals may be more representative of actual prediction errors than the built-in GP prediction standard deviation, the latter being in this case not depending on the actual observations and calculated under a questionable stationarity assumption.

Fig. 1 On the upper panel, GP mean predictor (dashed blue line) of the test function (black line) defined by (A.4) based on 10 evaluations at a regular grid, LOO cross-validation predictions (red stars). Lower panel: absolute prediction errors associated with GP (black line) and LOO (red stars) predictions, and GP prediction standard deviation (dashed blue line).

It has in fact been established in the seminal paper Dubrule (1983) that the LOO prediction errors and variance could be calculated quite elegantly and efficiently by relying on the inverse of a certain matrix (the main focus in Dubrule (1983) being on a variogram matrix extended by monomial basis functions) and quantities derived thereof. We work here in terms of covariances and related matrices as in Bachoc (2013), where fast LOO formulas were obtained for simple kriging in the square-integrable case.

An elegant by-product formula for the square norm of LOO residulas has been leveraged in Bachoc (2013) for hyperparameter estimation in the stationary case, namely by minimizing this quantity as a function of covariance parameters (excluding the “variance” parameter). In particular, numerical experiments conducted in Bachoc (2013) suggested that leave-one-out error minimization is more robust than maximum likelihood estimation (MLE) in the case of model misspecification, a common situation in practice. Yet, depending on the experimental design, the errors obtained by leaving one point at a time may be far from representative of generalization errors, as illustrated in Section 6 and in Section A.2.1 of the supplementary material. While remaining with a prescribed design, multiple-fold cross-validation allows to mitigate shortcomings of leave-one-out by considering situations where larger parts of the design (i.e., several points) are removed at a time.

Our primary aim in the present paper is to generalize fast-leave-one-out formulas to multiple-fold cross-validation (MFCV), hence, facilitating their use in both contexts of model diagnostics and hyperparameter optimization. We derive such formulas and highlight in turn the covariance structure of cross-validation residuals, a crucial ingredient in investigated diagnostics and in approaches pertaining to parameter fitting. Furthermore, following ideas already sketched in Dubrule (1983), we present extensions of our results to the case of universal kriging. The obtained results clarify some links between cross-validation and MLE, and also open new perspectives regarding the design of folds and ultimately also of the underlying experimental design points.

Section 2 presents the considered class of statistical models and recalls the basics of simple and universal kriging under Gaussian process assumptions, with some brief detours in Section 2.1 on how Gaussian linear and ridge regression may fit into this framework as well as on Bayesian interpretations of universal kriging and ridge regression. Prediction equations are revisited in Section 2.2 in transductive settings thanks to Schur complement approaches, hence, setting the decor for fast cross-validation. Section 3 is dedicated to fast multiple-fold cross-validation, starting in Section 3.1 with generic results that apply in Gaussian vector conditioning. We then present in Section 3.2 our main results on the fast calculation and the joint probability distribution of cross-validation residuals in the simple kriging case. The latter are then extended to universal kriging settings in Section 3.3, notably retrieving as a particular case fast cross-validation results for linear regression presented in Zhang (1993). In Section 4 we present consequences of the main results in the context of model fitting. Section 4.1 focuses on graphical diagnostics via pivotal statistics for cross-validation residuals, while the two subsequent sections deal with cross-validation-based covariance parameter estimation. More specifically, Section 4.2 is devoted to the estimation of the scale parameter; Section 4.3 then pertains to the estimation of further hyperparameters building upon either the norm or the (pseudo-)likelihood of multiple-fold cross-validation residuals. Section 5 mostly consists in numerical experiments illustrating the accuracy and speed-ups offered by the closed-form formulas (Section 5.1), and some considerations regarding associated computational complexities (Section 5.2). Section 6 presents an application to a contaminant localization test case and illustrate that grouping observations in folds may lead to improved model assessment and parameter fitting compared to LOO. The ultimate Section 7 is a discussion opening perspectives on how the established results could be further exploited to improve diagnostics and parameter estimation procedures in GP modeling.

2 Models and Transductive Prediction Settings

2.1 Of Universal Kriging and Ridge Regression

Throughout the article we consider statistical models of the form (1) $$Z_i=\xi ({\mathbf{x}}_i)+{\epsilon }_i({\mathbf{x}}_i\in D,i\in \{1,\dots ,n\}),$$(1) with $\xi (\mathbf{x})=\mu (\mathbf{x})+\eta (\mathbf{x})$ ($\mathbf{x}\in D$) where $\eta$ is a centered GP with covariance kernel k, μ is a trend function, and $\mathit{\epsilon }={({\epsilon }_1,\dots ,{\epsilon }_n)}^{\top }\sim \mathcal{N}(\mathbf{0},{\Sigma }_{\epsilon })$ independently of $\eta$. In some cases, μ is known up to some coefficients: unknown constant (ordinary kriging settings) or linear combination of basis functions with unknown coefficients (universal kriging settings). When μ is known, one speaks of simple kriging settings. While several variants of Kriging can be defined without distributional assumptions (and without even assuming square-integrability of $\eta$), we restrict the exposition here to the broadly used Gaussian case. However, we do not require k to be stationary. Let us remark that throughout the following exposition on kriging/GP prediction, the kernel (including hyperparameters) is considered as given. In practice, the equations are also employed with plugged in estimates of hyperparameters. Full Bayesian approaches where the uncertainty on kernel hyperparameters is further propagated to the predictions are beyond the considered scope. However, in Sections 4 and 6 we tackle the issue of model validation and hyperparameter fitting via the established cross-validation results.

We consider a linear trend case that encompasses the above settings, and also has some interesting connections to linear (ridge) regression as we will see below. We hence assume μ to be defined by (2) $$\mu (\mathbf{x})=\sum _{j=1}^p{\beta }_j{f}_j(\mathbf{x})(\mathbf{x}\in D),$$(2) where the f_j ($1\le j\le p$) are prescribed basis functions and β_j ($1\le j\le p$) are real-valued coefficients. Simple kriging can be retrieved by taking p = 1 with ${\beta }_1=1$ and f₁ arbitrary (possibly null), and universal kriging with $p\ge 1$ and $\mathit{\beta }=({\beta }_1,\dots ,{\beta }_p)$ unknown (ordinary kriging is a special case with p = 1 and $f_1\equiv 1$). In simple kriging, a predictor of $\xi$ of the form $\widehat{\xi }(\mathbf{x})=\mu (\mathbf{x})+\mathit{\lambda }{(\mathbf{x})}^{\top }\mathbf{Z}(\mathbf{x}\in D)$ is sought, where $\mathbf{Z}={(Z_1,\dots ,Z_n)}^{\top }$. $\mathit{\lambda }(\mathbf{x})$ is determined by minimizing the residual variance (3) $$\begin{array}{c}\text{var}[\xi (\mathbf{x})-\widehat{\xi }(\mathbf{x})]=k(\mathbf{x},\mathbf{x})+\mathit{\lambda }{(\mathbf{x})}^{\top }(K+{\Sigma }_{\epsilon })\mathit{\lambda }(\mathbf{x})\\ -2\mathbf{k}{(\mathbf{x})}^{\top }K\mathit{\lambda }(\mathbf{x}),\end{array}$$(3)

where $K={(k({\mathbf{x}}_i,{\mathbf{x}}_j))}_{i,j\in \{1,\dots ,n\}}$ and $\mathbf{k}(\mathbf{x})={(k(\mathbf{x},{\mathbf{x}}_j))}_{j\in \{1,\dots ,n\}}$. In universal kriging, the predictor is of the form $$\widehat{\xi }(\mathbf{x})=\mathit{\lambda }{(\mathbf{x})}^{\top }\mathbf{Z}(\mathbf{x}\in D),$$ and p linear constraints $f_j(\mathbf{x})={\sum }_{i=1}^n{\lambda }_i(\mathbf{x})f_j({\mathbf{x}}_j)$ for $j\in \{1,\dots ,p\}$ are added to the minimization of (3) to ensure unbiasedness of $\widehat{\xi }(\mathbf{x})$. Minimizing (3) under these unbiasedness constraints can be addressed by Lagrange duality, leading to a linear problem featuring for each $\mathbf{x}\in D$a p-dimensional Lagrange multiplier $@@\mathcal{l}(\mathbf{x})$: (4) $$\left(\begin{array}{cc}K+{\Sigma }_{\epsilon }& F\\ F^{\top }& 0\end{array}\right)\left(\begin{array}{c}\mathit{\lambda }(\mathbf{x})\\ \mathbf{\ell }(\mathbf{x})\end{array}\right)=\left(\begin{array}{c}\mathbf{k}(\mathbf{x})\\ \mathbf{f}(\mathbf{x})\end{array}\right),$$(4) where $F={(f_j({\mathbf{x}}_i))}_{1\le i\le n,1\le j\le p}\in {\mathbb{R}}^{n\times p}$ and $\mathbf{f}(\mathbf{x})=$ ${(f_1(\mathbf{x}),\dots ,f_p(\mathbf{x}))}^{\top }\in {\mathbb{R}}^p$. Assuming that $\Sigma =K+{\Sigma }_{\epsilon }$ and $\mathfrak{I}=F^{\top }{\Sigma }^{-1}F$ are invertible, solving for $\mathit{\lambda }(\mathbf{x})$ delivers the universal kriging predictor, that can ultimately be written, as further detailed in Section A.3 of the supplementary material, as (5) $$\widehat{\xi }(\mathbf{x})=\mathbf{f}{(\mathbf{x})}^{\top }\widehat{\mathit{\beta }}+\mathbf{k}{(\mathbf{x})}^{\top }{\Sigma }^{-1}(\mathbf{Z}-F\widehat{\mathit{\beta }})(\mathbf{x}\in D)$$(5) where $\widehat{\mathit{\beta }}={\mathfrak{I}}^{-1}{F}^{\top }{\Sigma }^{-1}\mathbf{Z}$. In simple kriging with $\mu (\mathbf{x})=\mathbf{f}{(\mathbf{x})}^{\top }\mathit{\beta }$, the predictor has the same form as (5) yet with $\mathit{\beta }$ instead of $\widehat{\mathit{\beta }}$. Furthermore, the universal kriging residual covariance writes, for arbitrary $\mathbf{x},\mathbf{x}\prime \in D$: (6) $$\begin{array}{cc}\hfill & \text{cov}[\xi (\mathbf{x})-\widehat{\xi }(\mathbf{x}),\xi (\mathbf{x}\prime )-\widehat{\xi }(\mathbf{x}\prime )]=k(\mathbf{x},\mathbf{x}\prime )-\mathbf{k}{(\mathbf{x})}^{\top }{\Sigma }^{-1}\mathbf{k}(\mathbf{x}\prime )\hfill \\ \hfill & +{(\mathbf{f}(\mathbf{x})-F^{\top }{\Sigma }^{-1}\mathbf{k}(\mathbf{x}))}^{\top }{\mathfrak{I}}^{-1}(\mathbf{f}(\mathbf{x}\prime )-F^{\top }{\Sigma }^{-1}\mathbf{k}(\mathbf{x}\prime )).\hfill \end{array}$$(6)

The simple kriging residual covariance boils down to the right-hand side of the first line of (6). On a different note, it is worth noting that when setting $\eta \equiv 0$ (and hence $\xi (\mathbf{x})={\sum }_{j=1}^p{\beta }_j{f}_j(\mathbf{x})$) one retrieves thereby the equations of (generalized) least squares.

We now review the Bayesian approach to universal kriging Omre (1987), Handcock and Stein (1993), and Helbert, Dupuy, and Carraro (2009), that will in turn give us an occasion to incorporate ridge regression into the discussion and will prove practical in further developments throughout the article. Let us assume to this end that $\mathit{\beta }$ is endowed with a Gaussian prior distribution $\mathcal{N}(\mathbf{0},{\Sigma }_{\beta })$ where ${\Sigma }_{\beta }$ is an invertible covariance matrix. Let us first consider ${\Sigma }_{\beta }$ as fixed. The posterior distribution of $\mathit{\beta }$ knowing Z is Gaussian with $\mathbb{E}[\mathit{\beta }|\mathbf{Z}]={\widehat{\mathit{\beta }}}_{\text{MAP}}:={({\Sigma }_{\beta }^{-1}+F^{\top }{\Sigma }^{-1}F)}^{-1}{F}^{\top }{\Sigma }^{-1}\mathbf{Z}$ and (7) $$\text{cov}\left[\mathit{\beta }\right|\mathbf{Z}]={({\Sigma }_{\beta }^{-1}+F^{\top }{\Sigma }^{-1}F)}^{-1}$$(7)

whereby, for $\mathbf{x},\mathbf{x}\prime \in D,\mathbb{E}[\xi (\mathbf{x})|\mathbf{Z}]=\mathbf{f}{(\mathbf{x})}^{\top }{\widehat{\mathit{\beta }}}_{\text{MAP}}+\mathbf{k}{(\mathbf{x})}^{\top }{\Sigma }^{-1}(\mathbf{Z}-F{\widehat{\mathit{\beta }}}_{\text{MAP}})$ and (8) $$\begin{array}{c}\text{cov}[\xi (\mathbf{x}),\xi (\mathbf{x}\prime )|\mathbf{Z}]=k(\mathbf{x},\mathbf{x})-\mathbf{k}{(\mathbf{x})}^{\top }{\Sigma }^{-1}\mathbf{k}(\mathbf{x}\prime )\\ +{(\mathbf{f}(\mathbf{x})-F^{\top }{\Sigma }^{-1}\mathbf{k}(\mathbf{x}))}^{\top }\text{cov}\left[\mathit{\beta }\right|\mathbf{Z}](\mathbf{f}(\mathbf{x}\prime )-F^{\top }{\Sigma }^{-1}\mathbf{k}(\mathbf{x}\prime )).\end{array}$$(8)

With ${\Sigma }_{\beta }^{-1}\to \mathbf{0}$, we hence have that $\text{cov}[\mathit{\beta }|\mathbf{Z}]\to {\mathfrak{I}}^{-1}$ and so $\mathbb{E}[\xi (\mathbf{x})|\mathbf{Z}]\to \widehat{\xi }(\mathbf{x})$ and $\text{cov}[\xi (\mathbf{x}),\xi (\mathbf{x}\prime )|\mathbf{Z}]\to \text{cov}[\xi (\mathbf{x})-\widehat{\xi }(\mathbf{x}),\xi (\mathbf{x}\prime )-\widehat{\xi }(\mathbf{x}\prime )]$.

We now turn to linear (ridge) regression and present how the resulting equations can be obtained in similar settings. In linear ridge regression, say as classically with ${\Sigma }_{\epsilon }={\sigma }^2{I}_n,$ one estimates $\mathit{\beta }$ via ${\widehat{\mathit{\beta }}}_{\lambda }={(F^{\top }F+\lambda I_n)}^{-1}{F}^{\top }\mathbf{Z}$, where $\lambda >0$ is a regularization parameter and I_n stands for the identity matrix ${({\delta }_{\mathit{\text{ij}}})}_{i,j\in \{1,\dots ,n\}}\in {\mathbb{R}}^{n\times n}$. The canonical way to introduce this estimator is by replacing the ordinary least squares minimization by the minimization of its penalized counterpart $$\mathit{\beta }\mapsto \left|\right|\mathbf{Z}-F\mathit{\beta }|{|}_2^2+\lambda |\left|\mathit{\beta }\right|{|}_2^2.$$

Yet it is well-known that ${\widehat{\mathit{\beta }}}_{\lambda }$ can also be seen as maximum a posterior estimator in a Bayesian framework, under the prior distribution $\mathit{\beta }\sim \mathcal{N}(\mathbf{0},{\gamma }^2{I}_p)$, where ${\gamma }^2=\frac{{\sigma }^2}{\lambda }$. In such a framework, one has indeed $\mathbf{Z}\sim \mathcal{N}(\mathbf{0},{\gamma }^{2}F{F}^{\top }+{\sigma }^2{I}_n)$ and $\mathit{\beta }|\mathbf{Z}=\mathbf{z}\sim \mathcal{N}(F^{\top }{(F{F}^{\top }+\lambda I_p)}^{-1}\mathbf{z},{\gamma }^2{I}_p-{\gamma }^2{F}^{\top }{(F{F}^{\top }+\lambda I_n)}^{-1}F)$. From there one can obtain, for instance by using the Woodbury formula, that $$F^{\top }{(F{F}^{\top }+\lambda I_p)}^{-1}={(F^{\top }F+\lambda I_p)}^{-1}{F}^{\top },$$ whereby $\mathbb{E}[\mathit{\beta }|\mathbf{Z}]={(F^{\top }F+\lambda I_p)}^{-1}{F}^{\top }\mathbf{Z}={\widehat{\mathit{\beta }}}_{\lambda }$, which coincides in the considered Gaussian case with the MAP estimator. For F full column-ranked and $\gamma \to \infty$, the posterior distribution is proper and one recovers as posterior expectation and MAP the ordinary least squares estimator $\widehat{\mathit{\beta }}={(F^{\top }F)}^{-1}{F}^{\top }\mathbf{Z}$.

2.2 Transductive GP Prediction via Schur Complements

One speaks of transductive settings in cases where the point(s) at which one wishes to predict are fixed in advance. We consider in this section the case where $\xi ({\mathbf{x}}_i)$ (resp. Z_i ) are to be predicted for $m\le n-1$ distinct values of i ($m\le n-p$ in the case of universal kriging) from the remaining n–m observations Z_j . Without loss of generality, we will assume that the indices at which the predictions are to be made are ${\mathbf{i}}_o=(1,\dots ,m)$ and denote by ${\mathbf{j}}_o=(m+1,\dots ,n)$ the remaining ones. For any vector $\mathbf{v}\in {\mathbb{R}}^n$, matrix $A\in {\mathbb{R}}^{n\times n}$, and arbitrary vectors of ordered indices $\mathbf{i},\mathbf{j}$, we denote here and in the following by $\mathbf{v}[\mathbf{i}]$ the subvector of v with indices in i and by $A[\mathbf{i},\mathbf{j}]$ the block extracted from A with corresponding indices, using in turn the convention $A[\mathbf{i}]=A[\mathbf{i},\mathbf{i}]$. Similarly, we use the convention $k(\mathbf{x})=k(\mathbf{x},\mathbf{x})$ for $\mathbf{x}\in D$.

In simple kriging settings (here $\mu \equiv 0$), it is common knowledge that $\widehat{\xi }({\mathbf{x}}_{{\mathbf{i}}_o})$ and the associated residual covariance matrix $\text{cov}(\xi ({\mathbf{x}}_{{\mathbf{i}}_o})-\widehat{\xi }({\mathbf{x}}_{{\mathbf{i}}_o}))$ can be obtained elegantly based on manipulating (9) $${\Sigma }_{\star }=\left(\begin{array}{cc}k({\mathbf{x}}_{{\mathbf{i}}_o})& k({\mathbf{x}}_{{\mathbf{i}}_o},{\mathbf{x}}_{{\mathbf{j}}_o})\\ k({\mathbf{x}}_{{\mathbf{j}}_o},{\mathbf{x}}_{{\mathbf{i}}_o})& k({\mathbf{x}}_{{\mathbf{j}}_o})+{\Sigma }_{\epsilon }\left[{\mathbf{j}}_o\right]\end{array}\right).$$(9)

Assuming indeed that ${\Sigma }_{\star }$ is invertible, we obtain by block inversion formula (See Section A.1 in supplementary material, in particular Proposition A.1) that ${\Sigma }_{\star }^{-1}[{\mathbf{i}}_o]={(k({\mathbf{x}}_{{\mathbf{i}}_o})-k({\mathbf{x}}_{{\mathbf{i}}_o},{\mathbf{x}}_{{\mathbf{j}}_o}){\Sigma }_{\star }{[{\mathbf{j}}_o]}^{-1}k({\mathbf{x}}_{{\mathbf{j}}_o},{\mathbf{x}}_{{\mathbf{i}}_o}))}^{-1}$, whereby (10) $$\text{cov}(\xi ({\mathbf{x}}_{{\mathbf{i}}_o})-\widehat{\xi }({\mathbf{x}}_{{\mathbf{i}}_o}))={({\Sigma }_{\star }^{-1}[{\mathbf{i}}_o])}^{-1}$$(10) where ${\Sigma }_{\star }^{-1}[{\mathbf{i}}_o]$ stands for the m × m upper left block of ${\Sigma }_{\star }^{-1}$. Using now the upper right block ${\Sigma }_{\star }^{-1}[{\mathbf{i}}_o,{\mathbf{j}}_o]$ of ${\Sigma }_{\star }^{-1}$ and (10), we obtain similarly that ${\Sigma }_{\star }^{-1}[{\mathbf{i}}_o,{\mathbf{j}}_o]=-{\Sigma }_{\star }^{-1}[{\mathbf{i}}_o]k({\mathbf{x}}_{{\mathbf{i}}_o},{\mathbf{x}}_{{\mathbf{j}}_o}){\Sigma }_{\star }{[{\mathbf{j}}_o]}^{-1}$, so (11) $$\widehat{\xi }({\mathbf{x}}_{{\mathbf{i}}_o})=-{({\Sigma }_{\star }^{-1}[{\mathbf{i}}_o])}^{-1}{\Sigma }_{\star }^{-1}[{\mathbf{i}}_o,{\mathbf{j}}_o]\mathbf{Z}\left[{\mathbf{j}}_o\right].$$(11)

Applying now under invertibility assumption a similar block inversion to (12) $$\begin{array}{c}\Sigma =\left(\begin{array}{cc}\Sigma \left[{\mathbf{i}}_o\right]& \Sigma [{\mathbf{i}}_o,{\mathbf{j}}_o]\\ \Sigma [{\mathbf{j}}_o,{\mathbf{i}}_o]& \Sigma \left[{\mathbf{j}}_o\right]\end{array}\right)\\ =\mathrm{}\left(\mathrm{}\begin{array}{cc}k({\mathbf{x}}_{{\mathbf{i}}_o})+{\Sigma }_{\epsilon }\left[{\mathbf{i}}_o\right]& k({\mathbf{x}}_{{\mathbf{i}}_o},{\mathbf{x}}_{{\mathbf{j}}_o})+{\Sigma }_{\epsilon }[{\mathbf{i}}_o,{\mathbf{j}}_o]\\ k({\mathbf{x}}_{{\mathbf{j}}_o},{\mathbf{x}}_{{\mathbf{i}}_o})+{\Sigma }_{\epsilon }[{\mathbf{j}}_o,{\mathbf{i}}_o]& k({\mathbf{x}}_{{\mathbf{j}}_o})+{\Sigma }_{\epsilon }\left[{\mathbf{j}}_o\right]\end{array}\right)\mathrm{},\end{array}$$(12) we get ${\Sigma }^{-1}[{\mathbf{i}}_o]={(\Sigma [{\mathbf{i}}_o]-\Sigma [{\mathbf{i}}_o,{\mathbf{j}}_o]\Sigma {[{\mathbf{j}}_o]}^{-1}\Sigma [{\mathbf{j}}_o,{\mathbf{i}}_o])}^{-1}$ and ${\Sigma }^{-1}[{\mathbf{i}}_o,{\mathbf{j}}_o]=-{\Sigma }^{-1}[{\mathbf{i}}_o]\Sigma [{\mathbf{i}}_o,{\mathbf{j}}_o]\Sigma {[{\mathbf{j}}_o]}^{-1}$ whereby the BLUP of ${\mathbf{Z}}_{{\mathbf{i}}_o}$ given ${\mathbf{Z}}_{{\mathbf{j}}_o}$, denoted here and in the following by ${\widehat{\mathbf{Z}}}^{(-{\mathbf{i}}_o)}[{\mathbf{i}}_o]$, is given by (13) $$\begin{array}{c}{\widehat{\mathbf{Z}}}^{(-{\mathbf{i}}_o)}\left[{\mathbf{i}}_o\right]=-{({\Sigma }^{-1}[{\mathbf{i}}_o])}^{-1}{\Sigma }^{-1}[{\mathbf{i}}_o,{\mathbf{j}}_o]\mathbf{Z}\left[{\mathbf{j}}_o\right].\hfill \end{array}$$(13)

Considering now the prediction residual ${\mathbf{E}}_{{\mathbf{i}}_o}:=\mathbf{Z}[{\mathbf{i}}_o]-{\widehat{\mathbf{Z}}}^{(-{\mathbf{i}}_o)}[{\mathbf{i}}_o]$, we get from the latter (14) $$\begin{array}{cc}{\mathbf{E}}_{{\mathbf{i}}_o}\hfill & ={({\Sigma }^{-1}[{\mathbf{i}}_o])}^{-1}({\Sigma }^{-1}\mathbf{Z})\left[{\mathbf{i}}_o\right],\hfill \\ \hfill & ={(Q[{\mathbf{i}}_o])}^{-1}(Q\mathbf{Z})\left[{\mathbf{i}}_o\right],\hfill \end{array}$$(14) with $Q={\Sigma }^{-1}$. Let us stress that ${\widehat{\mathbf{Z}}}^{(-{\mathbf{i}}_o)}[{\mathbf{i}}_o]$ departs from the most common simple kriging predictor via the treatment of noise, in the sense that in the considered transductive settings it can be relevant to predict the noisy $\mathbf{Z}[{\mathbf{i}}_o]$ instead of $\xi ({\mathbf{x}}_{{\mathbf{i}}_o})$, and since a general noise covariance matrix is assumed for the sake of generality, this implies a straightforward adaptation of the blocks entering into play (i.e., blocks of Σ instead of those of ${\Sigma }_{\star }$). Of course, in case $\Sigma ={\Sigma }_{\star }$ (e.g., under a null observation noise, as often encountered in computer experiments), the two approaches above coincide.

We will now see that under invertibility conditions, the universal kriging predictor and residual covariance too can be retrieved based on the inverse of a matrix, namely of (15) $$M_{\star }=\left(\begin{array}{ccc}k({\mathbf{x}}_{{\mathbf{i}}_o})& k({\mathbf{x}}_{{\mathbf{i}}_o},{\mathbf{x}}_{{\mathbf{j}}_o})& F[{\mathbf{i}}_o,]\\ k({\mathbf{x}}_{{\mathbf{j}}_o},{\mathbf{x}}_{{\mathbf{i}}_o})& k({\mathbf{x}}_{{\mathbf{j}}_o})+{\Sigma }_{\epsilon }\left[{\mathbf{j}}_o\right]& F[{\mathbf{j}}_o,]\\ F{[{\mathbf{i}}_o,]}^{\top }& F{[{\mathbf{j}}_o,]}^{\top }& \mathbf{0}\end{array}\right),$$(15) where $F[{\mathbf{i}}_o,]$ is the m × p matrix of basis functions evaluated at the prediction points. As we prove below, assuming invertibility of $M_{\star }$ and of $M_{\star }[-{\mathbf{i}}_o]=M_{\star }[m+1:n+p]$, the universal kriging residual covariance matrix can then be plainly obtained as follows: (16) $$\text{cov}(\xi ({\mathbf{x}}_{{\mathbf{i}}_o})-\widehat{\xi }({\mathbf{x}}_{{\mathbf{i}}_o}))={(M_{\star }^{-1}[{\mathbf{i}}_o])}^{-1}.$$(16)

Furthermore, the universal kriging predictor too writes similarly as in the simple kriging case: (17) $$\widehat{\xi }({\mathbf{x}}_{{\mathbf{i}}_o})=-{(M_{\star }^{-1}[{\mathbf{i}}_o])}^{-1}{M}_{\star }^{-1}[{\mathbf{i}}_o,{\mathbf{j}}_o]\mathbf{Z}\left[{\mathbf{j}}_o\right].$$(17)

In fact, in the two last equations, what changed between simple and universal kriging settings is simply the replacement of ${\Sigma }_{\star }$ by $M_{\star }$. Let us now present how this works. Assuming indeed that $M_{\star }$ and $M_{\star }[-{\mathbf{i}}_o]$ are invertible, we first apply Proposition A.1 and obtain (18) $$\begin{array}{c}{(M_{\star }^{-1}[{\mathbf{i}}_o])}^{-1}=k({\mathbf{x}}_{{\mathbf{i}}_o})-[k({\mathbf{x}}_{{\mathbf{i}}_o},{\mathbf{x}}_{{\mathbf{j}}_o}),F[{\mathbf{i}}_o,\left]\right]M_{\star }{[-{\mathbf{i}}_o]}^{-1}\\ {[k({\mathbf{x}}_{{\mathbf{i}}_o},{\mathbf{x}}_{{\mathbf{j}}_o}),F[{\mathbf{i}}_o,\left]\right]}^{\top }.\end{array}$$(18)

Now, provided $\Sigma [{\mathbf{j}}_o]=k({\mathbf{x}}_{{\mathbf{j}}_o})+{\Sigma }_{\epsilon }[{\mathbf{j}}_o]$ and ${\mathfrak{I}}_{-{\mathbf{i}}_o}=F{[{\mathbf{j}}_o,]}^{\top }\Sigma {[{\mathbf{j}}_o]}^{-1}F[{\mathbf{j}}_o,]$ are invertible, we get using a variation of the block inversion formula (detailed in Remark A.3): (19) $$M_{\star }{[-{\mathbf{i}}_o]}^{-1}=\left(\begin{array}{cc}\Sigma {\left[{\mathbf{j}}_o\right]}^{-1}-\Sigma {\left[{\mathbf{j}}_o\right]}^{-1}F[{\mathbf{j}}_o,]{\mathfrak{I}}_{-{\mathbf{i}}_o}^{-1}F{[{\mathbf{j}}_o,]}^{\top }\Sigma {\left[{\mathbf{j}}_o\right]}^{-1}& \Sigma {\left[{\mathbf{j}}_o\right]}^{-1}F[{\mathbf{j}}_o,]{\mathfrak{I}}_{-{\mathbf{i}}_o}^{-1}\\ {\mathfrak{I}}_{-{\mathbf{i}}_o}^{-1}F{[{\mathbf{j}}_o,]}^{\top }\Sigma {\left[{\mathbf{j}}_o\right]}^{-1}& -{\mathfrak{I}}_{-{\mathbf{i}}_o}^{-1}\end{array}\right).$$(19)

Expanding the product on the right-hand side of (18) using (19) then delivers (20) $$\begin{array}{c}(k({\mathbf{x}}_{{\mathbf{i}}_o},{\mathbf{x}}_{{\mathbf{j}}_o}),F[{\mathbf{i}}_o,])M_{\star }^{-1}{\left[{\mathbf{i}}_o\right]}^{-1}{(k({\mathbf{x}}_{{\mathbf{i}}_o},{\mathbf{x}}_{{\mathbf{j}}_o}),F[{\mathbf{i}}_o,])}^{\top }\\ =k({\mathbf{x}}_{{\mathbf{i}}_o},{\mathbf{x}}_{{\mathbf{j}}_o})\Sigma {\left[{\mathbf{j}}_o\right]}^{-1}k({\mathbf{x}}_{{\mathbf{j}}_o},{\mathbf{x}}_{{\mathbf{i}}_o})\\ -(F[{\mathbf{i}}_o,]-k({\mathbf{x}}_{{\mathbf{i}}_o},{\mathbf{x}}_{{\mathbf{j}}_o})\Sigma {\left[{\mathbf{j}}_o\right]}^{-1}F[{\mathbf{j}}_o,]){\mathfrak{I}}_{-{\mathbf{i}}_o}^{-1}(F[{\mathbf{i}}_o,]\\ -k({\mathbf{x}}_{{\mathbf{i}}_o},{\mathbf{x}}_{{\mathbf{j}}_o})\Sigma {\left[{\mathbf{j}}_o\right]}^{-1}F[{\mathbf{j}}_o,]{)}^{\top },\end{array}$$(20) which establishes (16). As for (17), it then follows from Proposition A.2 and (19) that $$\begin{array}{cc}\hfill & M_{\star }^{-1}[{\mathbf{i}}_o,{\mathbf{j}}_o]\mathbf{Z}\left[{\mathbf{j}}_o\right]\hfill \\ \hfill & =-M_{\star }^{-1}\left[{\mathbf{i}}_o\right](k({\mathbf{x}}_{{\mathbf{i}}_o},{\mathbf{x}}_{{\mathbf{j}}_o}),F[{\mathbf{i}}_o,])M_{\star }{[-{\mathbf{i}}_o]}^{-1}\tilde{\mathbf{Z}}[-{\mathbf{i}}_o]\hfill \\ \text{where by }\hfill & -{(M_{\star }^{-1}[{\mathbf{i}}_o])}^{-1}{M}_{\star }^{-1}[{\mathbf{i}}_o,{\mathbf{j}}_o]\mathbf{Z}\left[{\mathbf{j}}_o\right]\hfill \\ \hfill & =k({\mathbf{x}}_{{\mathbf{i}}_o},{\mathbf{x}}_{{\mathbf{j}}_o})\Sigma {\left[{\mathbf{j}}_o\right]}^{-1}(\mathbf{Z}[{\mathbf{j}}_o]-F[{\mathbf{j}}_o,]\hfill \\ \hfill & \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}{\widehat{\mathit{\beta }}}^{(-{\mathbf{i}}_o)})+F[{\mathbf{i}}_o,]{\widehat{\mathit{\beta }}}^{(-{\mathbf{i}}_o)}\hfill \end{array}$$ with $\tilde{\mathbf{Z}}={({\mathbf{Z}}^{\top }{\mathbf{0}}^{\top })}^{\top }\in {\mathbb{R}}^{n+p}$ and ${\widehat{\mathit{\beta }}}^{(-{\mathbf{i}}_o)}$ the GLS estimator of $\mathit{\beta }$ based on $\mathbf{Z}[{\mathbf{j}}_o]$.

Applying now under suitable invertibility assumptions a similar block inversion to (21) $$M=\left(\begin{array}{cc}\Sigma & F\\ F^{\top }& \mathbf{0}\end{array}\right),$$(21) we get the BLUP ${\widehat{\mathbf{Z}}}^{(-{\mathbf{i}}_o)}[{\mathbf{i}}_o]$ of $\mathbf{Z}[{\mathbf{i}}_o]$ given $\mathbf{Z}[{\mathbf{j}}_o]$ under universal kriging settings as (22) $${\widehat{\mathbf{Z}}}^{(-{\mathbf{i}}_o)}\left[{\mathbf{i}}_o\right]=-{(M^{-1}[{\mathbf{i}}_o])}^{-1}{M}^{-1}[{\mathbf{i}}_o,{\mathbf{j}}_o]\mathbf{Z}\left[{\mathbf{j}}_o\right].$$(22)

Considering further the associated prediction residual ${\mathbf{E}}_{{\mathbf{i}}_o}:=\mathbf{Z}[{\mathbf{i}}_o]-{\widehat{\mathbf{Z}}}^{(-{\mathbf{i}}_o)}[{\mathbf{i}}_o]$, we finally get (23) $$\begin{array}{cc}{\mathbf{E}}_{{\mathbf{i}}_o}\hfill & ={(M^{-1}[{\mathbf{i}}_o])}^{-1}(M^{-1}[1:n]\mathbf{Z})\left[{\mathbf{i}}_o\right],\hfill \\ \hfill & ={(\tilde{Q}[{\mathbf{i}}_o])}^{-1}(\tilde{Q}\mathbf{Z})\left[{\mathbf{i}}_o\right],\hfill \end{array}$$(23) with $\tilde{Q}=Q-\mathit{\text{QF}}{(F^{\top }\mathit{\text{QF}})}^{-1}{F}^{\top }Q$.

3 Fast Multiple-Fold Cross-Validation

In this section we leverage transductive GP equations in settings when subsets of the n observations are left out, delivering in turn the joint distribution of multiple-fold cross-validation residuals together with a fast computation approach. The key is that the previous results of (14) and (23) carry over to arbitrary ways of separating the dataset into left out and remaining observations (i.e., between test and learning sets). Throughout the following we consider vectors of strictly ordered indices from $\{1,\dots ,n\}$, and denote by $\mathcal{S}$ the set of all such index vectors. Since these index vectors are completely characterized by non-empty subsets of $\{1,\dots ,n\}$, there are $2^n-1$ elements in $\mathcal{S}$. Elements $\mathbf{i}\in \mathcal{S}$, that may also be thought of as subsets of $\{1,\dots ,n\}$, are called folds.

3.1 The Generic Gaussian Vector Case

We now consider generic Gaussian cross-validation settings and deterministic functions arising in this context. $\mathbf{v}\in {\mathbb{R}}^n$ and $A\in {\mathbb{R}}^{n\times n}$ stand for a vector and a symmetric positive definite matrix, respectively. For any $\mathbf{i}\in \mathcal{S}$, we define the following cross-validation residual mapping: (24) $$\begin{array}{c}{e}_{\mathbf{i}}(·;\mathbf{m},A):\mathbf{v}\in {\mathbb{R}}^n\to e_{\mathbf{i}}(\mathbf{v};\mathbf{m},A)=\mathbf{v}\left[\mathbf{i}\right]-\mathbf{m}\left[\mathbf{i}\right]\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}-A[\mathbf{i},-\mathbf{i}]A{[-\mathbf{i}]}^{-1}(\mathbf{v}[-\mathbf{i}]-\mathbf{m}[-\mathbf{i}]),\end{array}$$(24) delivering for $\mathbf{v}\in {\mathbb{R}}^n$ the deterministic counterpart of the cross-validation residual obtained when predicting $\mathbf{V}[\mathbf{i}]$ by ${\widehat{\mathbf{V}}}^{(-\mathbf{i})}[\mathbf{i}]=\mathbf{m}[\mathbf{i}]+A[\mathbf{i},-\mathbf{i}]A{[-\mathbf{i}]}^{-1}(\mathbf{V}[-\mathbf{i}]-\mathbf{m}[-\mathbf{i}])$, where $\mathbf{V}\sim \mathcal{N}(\mathbf{m},A)$. The notation $-\mathbf{i}$ refers again to the extraction of sub-vectors/matrices by removing components with indices in i.

Proposition 1

(Fast calculation of cross-validation residuals). Let $\mathbf{i}\in \mathcal{S}$ and $P=A^{-1}$. Then, (25) $$e_{\mathbf{i}}(\mathbf{v};\mathbf{m},A)=P{\left[\mathbf{i}\right]}^{-1}(P(\mathbf{v}-\mathbf{m}))\left[\mathbf{i}\right](\mathbf{v}\in {\mathbb{R}}^n).$$(25)

Consequently, for $\mathbf{V}\sim \mathcal{N}(\mathbf{m},A)$, the ${\mathbf{E}}_{\mathbf{i}}=e_{\mathbf{i}}(\mathbf{V};\mathbf{m},A)$ $(\mathbf{i}\in \mathcal{S})$ are jointly Gaussian, centered, and (26) $$\text{cov}({\mathbf{E}}_{\mathbf{i}},{\mathbf{E}}_{\mathbf{j}})=P{\left[\mathbf{i}\right]}^{-1}P[\mathbf{i},\mathbf{j}]P{\left[\mathbf{j}\right]}^{-1}(\mathbf{i},\mathbf{j}\in \mathcal{S}).$$(26)

For the case of folds ${\mathbf{i}}_1,\dots ,{\mathbf{i}}_q$ for which concatenation gives $(1,\dots ,n)$, ${\mathbf{E}}_{{\mathbf{i}}_1,\dots ,{\mathbf{i}}_q}:={[{\mathbf{E}}_{{\mathbf{i}}_1}^{\top },\dots ,{\mathbf{E}}_{{\mathbf{i}}_q}^{\top }]}^{\top }$ has centered n-dimensional Gaussian distribution with covariance matrix (27) $$\text{cov}({\mathbf{E}}_{{\mathbf{i}}_1,\dots ,{\mathbf{i}}_q})=B^{-1}P{B}^{-1}$$(27) where $B=\text{blockdiag}(P[{\mathbf{i}}_1],\dots ,P[{\mathbf{i}}_q])$.

Proof.

Equation (25)(25) $$e_{\mathbf{i}}(\mathbf{v};\mathbf{m},A)=P{\left[\mathbf{i}\right]}^{-1}(P(\mathbf{v}-\mathbf{m}))\left[\mathbf{i}\right](\mathbf{v}\in {\mathbb{R}}^n).$$(25) follows from suitably expressing blocks of A’s inverse using the Schur complement approach of Equations (A.1) and (A.2), delivering, respectively $$\begin{array}{c}{A}^{-1}\left[\mathbf{i}\right]={(A[\mathbf{i}]-A[\mathbf{i},-\mathbf{i}\left]A{[-\mathbf{i}]}^{-1}A\right[-\mathbf{i},\mathbf{i}])}^{-1}\text{, and}\\ A^{-1}[\mathbf{i},-\mathbf{i}]=-{(A[\mathbf{i}]-A[\mathbf{i},-\mathbf{i}\left]A{[-\mathbf{i}]}^{-1}A\right[-\mathbf{i},\mathbf{i}])}^{-1}\\ A[\mathbf{i},-\mathbf{i}]A{[-\mathbf{i}]}^{-1}.\end{array}$$

Considering the rows indexed by i of the product $P(\mathbf{v}-\mathbf{m})=A^{-1}(\mathbf{v}-\mathbf{m})$, we hence get that $$\begin{array}{cc}(P(\mathbf{v}-\mathbf{m}))\left[\mathbf{i}\right]=\hfill & P\left[\mathbf{i}\right](\mathbf{v}[\mathbf{i}]-\mathbf{m}[\mathbf{i}])\hfill \\ \hfill & -P\left[\mathbf{i}\right]A[\mathbf{i},-\mathbf{i}]A{[-\mathbf{i}]}^{-1}(\mathbf{v}[-\mathbf{i}]-\mathbf{m}[-\mathbf{i}])\hfill \end{array}$$ whereby $e_{\mathbf{i}}(\mathbf{v};\mathbf{m},A)={(P[\mathbf{i}])}^{-1}(P\mathbf{v})[\mathbf{i}]=P{[\mathbf{i}]}^{-1}{I}_n[\mathbf{i},]P(\mathbf{v}-\mathbf{m})$, where I_n denotes the n × n identity matrix. The joint Gaussianity of the ${\mathbf{E}}_{\mathbf{i}}$’s and the covariance structure follow from the affine dependence of ${\mathbf{E}}_{\mathbf{i}}={(P[\mathbf{i}])}^{-1}{I}_n[\mathbf{i},]P(\mathbf{V}-\mathbf{m})$ on V, so that concatenating the $q\ge 1$ random vectors ${\mathbf{E}}_{{\mathbf{i}}_1},\dots ,{\mathbf{E}}_{{\mathbf{i}}_q}$ leads to a Gaussian vector by left multiplication of $\mathbf{V}-\mathbf{m}$ by a deterministic matrix. As for the covariance matrix of ${\mathbf{E}}_{{\mathbf{i}}_1,\dots ,{\mathbf{i}}_q}$, it follows from $$\begin{array}{cc}{[{\mathbf{E}}_{{\mathbf{i}}_1}^{\top },\dots ,{\mathbf{E}}_{{\mathbf{i}}_q}^{\top }]}^{\top }\hfill & =\left(\begin{array}{c}P{\left[{\mathbf{i}}_1\right]}^{-1}{I}_n[{\mathbf{i}}_1,]P\mathbf{V}\hfill \\ \dots \hfill \\ P{\left[{\mathbf{i}}_q\right]}^{-1}{I}_n[{\mathbf{i}}_q,]P\mathbf{V}\hfill \end{array}\right)\hfill \\ \hfill & =\text{blockdiag}\left(P{\left[{\mathbf{i}}_1\right]}^{-1},\dots ,P{\left[{\mathbf{i}}_q\right]}^{-1}\mathrm{}\right)\mathrm{}\left(\mathrm{}\begin{array}{c}{I}_n[{\mathbf{i}}_1,]\hfill \\ \dots \hfill \\ I_n[{\mathbf{i}}_q,]\hfill \end{array}\mathrm{}\right)P\mathbf{V}\hfill \\ \hfill & =B^{-1}{I}_{n}P\mathbf{V}=B^{-1}P\mathbf{V},\hfill \end{array}$$ hence $\text{cov}({\mathbf{E}}_{{\mathbf{i}}_1,\dots ,{\mathbf{i}}_q})=B^{-1}{\mathit{\text{PAPB}}}^{-1}=B^{-1}P{B}^{-1}$. □

Remark 1.

The first part of Proposition 1 is purely deterministic, and hence applies to any function approximation methods that boils down to Gaussian conditioning equations (e.g., regularized quadratic risk minimization in Reproducing Kernel Hilbert Spaces).

Remark 2.

Note that for arbitary ${\mathbf{i}}_1,\dots ,{\mathbf{i}}_q\in \mathcal{S}$ that is without imposing ordering between them or that they form a partition, we would have a similar result yet without the above simplification, that is, $\text{cov}({\mathbf{E}}_{{\mathbf{i}}_1,\dots ,{\mathbf{i}}_q})=B^{-1}\Delta P{\Delta }^T{B}^{-1}$ with $\Delta ={[I_n{[{\mathbf{i}}_1,]}^{\top },\dots ,I_n{[{\mathbf{i}}_q,]}^{\top }]}^{\top }$. An extreme case would be to consider all possible non-empty subsets of $\{1,\dots ,n\}$, leading to $q=2^n-1$ and $n{2}^{n-1}$ lines for Δ.

Remark 3.

In the extreme case where q = n and the ${\mathbf{i}}_j$’s are set to (j) $(1\le j\le n)$, one recovers on the other hand fast leave-one-out cross-validation formulas, and we obtain as a by-product the covariance matrix of leave-one-out errors (28) $$\text{cov}({\mathbf{E}}_{(1),\dots ,(n)})=\text{diag}(P{\left[i\right]}^{-1})P\text{diag}(P{\left[i\right]}^{-1}).$$(28)

3.2 The Simple Kriging Case

We now come back to GP prediction settings and first focus on simple kriging, following the observation model $Z_i=\xi ({\mathbf{x}}_i)+{\epsilon }_i({\mathbf{x}}_i\in D,i\in \{1,\dots ,n\})$ from the previous section, where we recall that $\xi (\mathbf{x})=\mu (\mathbf{x})+\eta (\mathbf{x})$ ($\mathbf{x}\in D$) where $\eta$ is a centered GP with covariance kernel k and μ is a trend function, and $\mathit{\epsilon }\sim \mathcal{N}(\mathbf{0},{\Sigma }_{\epsilon })$ independently of $\eta$. μ is first assumed null, an assumption that we will relax later on. Like in the exposition of the previous section, we denote Z the observation vector and, for any $\mathbf{i}\in \mathcal{S},{\mathbf{E}}_{\mathbf{i}}=\mathbf{Z}[\mathbf{i}]-{\widehat{\mathbf{Z}}}^{(-\mathbf{i})}[\mathbf{i}]$ where ${\widehat{\mathbf{Z}}}^{(-\mathbf{i})}[\mathbf{i}]$ is the BLUP of $\mathbf{Z}[\mathbf{i}]$ based on $\mathbf{Z}[-\mathbf{i}]$. We assume throughout that $\Sigma =K+{\Sigma }_{\epsilon }$ is invertible, with $K={(k({\mathbf{x}}_i,{\mathbf{x}}_j))}_{i,j\in \{1,\dots ,n\}}$.

Corollary 1.

The ${\mathbf{E}}_{\mathbf{i}}$ $(\mathbf{i}\in \mathcal{S})$ are jointly Gaussian, centered, and (29) $$\text{cov}({\mathbf{E}}_{\mathbf{i}},{\mathbf{E}}_{\mathbf{j}})=Q{\left[\mathbf{i}\right]}^{-1}Q[\mathbf{i},\mathbf{j}]Q{\left[\mathbf{j}\right]}^{-1}(\mathbf{i},\mathbf{j}\in \mathcal{S}),$$(29) where $Q={\Sigma }^{-1}$. For the case of folds ${\mathbf{i}}_1,\dots ,{\mathbf{i}}_q$ for which concatenation gives $(1,\dots ,n)$, ${\mathbf{E}}_{{\mathbf{i}}_1,\dots ,{\mathbf{i}}_q}:={[{\mathbf{E}}_{{\mathbf{i}}_1}^{\top },\dots ,{\mathbf{E}}_{{\mathbf{i}}_q}^{\top }]}^{\top }$ has centered n-dimensional Gaussian distribution with covariance matrix (30) $$\text{cov}({\mathbf{E}}_{{\mathbf{i}}_1,\dots ,{\mathbf{i}}_q})=B^{-1}Q{B}^{-1}$$(30) where $B=\text{blockdiag}(Q[{\mathbf{i}}_1],\dots ,Q[{\mathbf{i}}_q])$.

Proof.

Directly follows from Lemma 1 using ${\mathbf{E}}_{\mathbf{i}}=e_{\mathbf{i}}(\mathbf{Z};\mathbf{0},\Sigma )$ with $\Sigma =K+{\Sigma }_{\epsilon }$. □

Remark 4.

Extending these results to the case of simple kriging with a known trend function μ directly follows from ${\mathbf{E}}_{\mathbf{i}}=e_{\mathbf{i}}(\mathbf{Z};\mathit{\mu },\Sigma )$ with $\mathit{\mu }={(\mu ({\mathbf{x}}_i))}_{1\le i\le n}$.

Remark 5.

In the case of Remark 3, we then obtain $\text{cov}(E_i,E_j)=Q{[i]}^{-1}Q[i,j]Q{[j]}^{-1}$ whereby $\text{corr}(E_i,E_j)=Q{[i]}^{-\frac{1}{2}}Q[i,j]Q{[j]}^{-\frac{1}{2}}$, leading in particular to a correlation matrix of LOO residuals (31) $$\text{corr}({\mathbf{E}}_{(1),\dots ,(n)})=\text{diag}\left(Q{\left[i\right]}^{-\frac{1}{2}}\right)Q\text{diag}\left(Q{\left[i\right]}^{-\frac{1}{2}}\right)=B^{\frac{1}{2}}Q{B}^{\frac{1}{2}}.$$(31)

As a by-product of Proposition 1 (in the limiting case of LOO, see Remark 3), we represent on Figure 2 the correlations between the LOO prediction residual at the leftmost location versus the LOO residuals at the other locations in the settings of the opening example from Figure 1.

Fig. 2 Correlations between LOO residuals at design points in the example displayed on Figure 1. This figure was realized with the R package corrplot Wei and Simko (2021).

A first interesting thing that we would like to stress here is that correlation between LOO prediction residuals at two successive locations is negative, with a value below–0.5 (Let us not that on Figure 2, the radii of the disks represent the magnitude of correlations, while their signs (resp. positive and negative) are represented via the black and white colors). An illustration of this effect can be seen on the upper panel of Figure 1, see in particular the fourth and fifth locations. Coming back to Figure 2, after the second location the correlation goes back to positive but with a smaller magnitude and subsequently appear to continue with a damped oscillation until stationing around zero. As we will see in Section 3, accounting for these correlations is instrumental in producing suitable QQ-plots from cross-validation residuals. We extend the analytical example in Section 3.3, with a focus on multiple-fold cross-validation.

3.3 Extension to Universal Kriging Case and More

Let us now formulate a corollary to Proposition 1 for the case of universal kriging such as presented throughout Section 2 in general and transductive settings.

Corollary 2.

For any $\mathbf{i}\in \mathcal{S}$ such that $F[-\mathbf{i},]$ is full column-ranked using the notation ${\widehat{\mathbf{Z}}}^{(-\mathbf{i})}[\mathbf{i}]$ to denote the BLUP of $\mathbf{Z}[\mathbf{i}]$ based on $\mathbf{Z}[-\mathbf{i}]$ in universal kriging settings, the residual ${\mathbf{E}}_{\mathbf{i}}=\mathbf{Z}[\mathbf{i}]-{\widehat{\mathbf{Z}}}^{(-\mathbf{i})}[\mathbf{i}]$ writes (32) $$\begin{array}{cc}{\mathbf{E}}_{\mathbf{i}}\hfill & ={(M^{-1}[\mathbf{i}])}^{-1}(M^{-1}[1:n]\mathbf{Z})\left[\mathbf{i}\right],\hfill \\ \hfill & ={(\tilde{Q}[\mathbf{i}])}^{-1}(\tilde{Q}\mathbf{Z})\left[\mathbf{i}\right],\hfill \end{array}$$(32) where M is defined in (21), and $\tilde{Q}=Q-\mathit{\text{QF}}{(F^{\top }\mathit{\text{QF}})}^{-1}{F}^{\top }Q$ with $Q={\Sigma }^{-1}$. Consequently, for any q > 1 and ${\mathbf{i}}_1,\dots ,{\mathbf{i}}_q\in \mathcal{S}$ such that the $F[-{\mathbf{i}}_j,]$’s are full column-ranked, the ${\mathbf{E}}_{{\mathbf{i}}_j}$ $(1\le j\le q)$ are jointly Gaussian, centered, and with covariance structure given by (33) $$\begin{array}{c}\text{cov}({\mathbf{E}}_{\mathbf{i}},{\mathbf{E}}_{\mathbf{j}})=\tilde{Q}{\left[\mathbf{i}\right]}^{-1}\tilde{Q}[\mathbf{i},\mathbf{j}]\tilde{Q}{\left[\mathbf{j}\right]}^{-1}\\ (\mathbf{i},\mathbf{j}\in \{{\mathbf{i}}_1,\dots ,{\mathbf{i}}_q\}).\end{array}$$(33)

In particular, for the case of index vectors ${\mathbf{i}}_1,\dots ,{\mathbf{i}}_q$ forming a partition of $\{1,\dots ,n\}$ and such that the concatenation of ${\mathbf{i}}_1,\dots ,{\mathbf{i}}_q$ gives $(1,\dots ,n)$, then the concatenated vector of cross-validation residuals ${\mathbf{E}}_{{\mathbf{i}}_1,\dots ,{\mathbf{i}}_q}:={[{\mathbf{E}}_{{\mathbf{i}}_1}^{\top },\dots ,{\mathbf{E}}_{{\mathbf{i}}_q}^{\top }]}^{\top }$ has a n-dimensional centered Gaussian distribution with covariance matrix (34) $$\text{cov}({\mathbf{E}}_{{\mathbf{i}}_1,\dots ,{\mathbf{i}}_q})={\tilde{B}}^{-1}\tilde{Q}{\tilde{B}}^{-1},$$(34) where $\tilde{B}=\text{blockdiag}(\tilde{Q}[{\mathbf{i}}_1],\dots ,\tilde{Q}[{\mathbf{i}}_q])$.

Proof.

A direct way to prove this property is to first generalize (23) to the case of generic folds $\mathbf{i}\in \mathcal{S}$ such that $F[-\mathbf{i},]$ is full column-ranked, which establishes (32), and then to follow the steps of Proposition 1 regarding the covariance structure and the joint Gaussian distribution of cross-validation residuals. Another proof taking a Bayesian route is sketched below. □

Remark 6.

As we have recalled in Section 2, it is known that universal kriging equations can be retrieved in a Bayesian setting by endowing $\mathit{\beta }$ with a Gaussian prior distribution $\mathcal{N}(\mathbf{0},{\Sigma }_{\beta })$ with invertible ${\Sigma }_{\beta }$ and letting ${\Sigma }_{\beta }^{-1}$ tend toward $\mathbf{0}\in {\mathbb{R}}^{p\times p}$. Under this model (at fixed invertible ${\Sigma }_{\beta }$), the random vector Z would be centered Gaussian with a covariance matrix $\Sigma +F{\Sigma }_{\beta }{F}^{\top }$. Then one may use Proposition 1 with $\mathbf{m}=0$ and $P=\Sigma +F{\Sigma }_{\beta }{F}^{\top }$ to find that the cross-validation errors $\mathbf{Z}[\mathbf{i}]-{\widehat{\mathbf{Z}}}_{{\Sigma }_{\beta }}^{(-\mathbf{i})}[\mathbf{i}]$ are jointly Gaussian with covariance structure $$\text{cov}(\mathbf{Z}[\mathbf{i}]-{\widehat{\mathbf{Z}}}_{{\Sigma }_{\beta }}^{(-\mathbf{i})}[\mathbf{i}],\mathbf{Z}[\mathbf{j}]-{\widehat{\mathbf{Z}}}_{{\Sigma }_{\beta }}^{(-\mathbf{j})}[\mathbf{i}])=P_{{\Sigma }_{\beta }}{\left[\mathbf{i}\right]}^{-1}{P}_{{\Sigma }_{\beta }}[\mathbf{i},\mathbf{j}]P_{{\Sigma }_{\beta }}{\left[\mathbf{j}\right]}^{-1}$$ where ${\widehat{\mathbf{Z}}}_{{\Sigma }_{\beta }}^{(-\mathbf{i})}[\mathbf{i}]$ denotes the BLUP (conditional expectation) of $\mathbf{Z}[\mathbf{i}]$ based on $\mathbf{Z}[-\mathbf{i}]$ under the prior distribution $\mathcal{N}(\mathbf{0},{\Sigma }_{\beta })$ on ${\Sigma }_{\beta }$, and $P_{{\Sigma }_{\beta }}={(\Sigma +F{\Sigma }_{\beta }{F}^{\top })}^{-1}$. Now by block inversion and Woodbury formula, we get $$\begin{array}{cc}{P}_{{\Sigma }_{\beta }}\hfill & ={\left(\begin{array}{cc}\Sigma \hfill & F\hfill \\ F^{\top }\hfill & -{\Sigma }_{\beta }^{-1}\hfill \end{array}\right)}^{-1}[1:n]\hfill \\ \hfill & =Q-\mathit{\text{QF}}{({\Sigma }_{\beta }^{-1}+F^{\top }\mathit{\text{QF}})}^{-1}{F}^{\top }Q.\hfill \end{array}$$

From there we can conclude since, for ${\Sigma }_{\beta }^{-1}\to \mathbf{0},P_{{\Sigma }_{\beta }}\to \tilde{Q}$.

Remark 7.

The latter corollary enables retrieving fast formulas for multiple-fold cross-validation residuals of linear regression models of Zhang (1993) as well as associated covariances. In fact, putting $K=\mathbf{0}$ and ${\Sigma }_{\epsilon }={\tau }^2{I}_n$ (where $\tau >0$) delivers ${\tau }^2\tilde{Q}=I_n-H$ with $H=F{(F^{\top }F)}^{-1}{F}^{\top }$, whereby ${\mathbf{E}}_{\mathbf{i}}=(I_n-H){[\mathbf{i}]}^{-1}(\mathbf{Z}-H\mathbf{Z})[\mathbf{i}]$ and $$\text{cov}({\mathbf{E}}_{\mathbf{i}},{\mathbf{E}}_{\mathbf{j}})=(I_n-H){\left[\mathbf{i}\right]}^{-1}{(I_n-H)}^{-1}[\mathbf{i},\mathbf{j}](I_n-H){\left[\mathbf{j}\right]}^{-1}(\mathbf{i},\mathbf{j}\in \mathcal{S}).$$

Similarly, in ridge regression, we obtain by putting ${\Sigma }_{\epsilon }={\tau }^2{I}_n$ and $K=\lambda I_n$, (35) $${\mathbf{E}}_{\mathbf{i}}=(I_n-H_{\lambda }){\left[\mathbf{i}\right]}^{-1}(\mathbf{Z}-H_{\lambda }\mathbf{Z})\left[\mathbf{i}\right],$$(35) where $H_{\lambda }=F{(F^{T}F+\lambda I_n)}^{-1}{F}^T$. Consequently, for any q > 1 and ${\mathbf{i}}_1,\dots ,{\mathbf{i}}_q\in \mathcal{S}$, the ${\mathbf{E}}_{{\mathbf{i}}_j,\lambda }$ $(1\le j\le q)$ are jointly Gaussian, centered, and with covariance structure given by (36) $$\begin{array}{c}\text{cov}({\mathbf{E}}_{\mathbf{i}},{\mathbf{E}}_{\mathbf{j}})=(I_n-H_{\lambda }){\left[\mathbf{i}\right]}^{-1}{(I_n-H_{\lambda })}^{-1}[\mathbf{i},\mathbf{j}]\\ (I_n-H_{\lambda }){\left[\mathbf{j}\right]}^{-1}(\mathbf{i},\mathbf{j}\in \mathcal{S}).\end{array}$$(36)

Many other models could enjoy similar results (encompassing kernel ridge regression and least squares support vector machines, see An, Liu, and Venkatesh (2007) where similar fast computations have been obtained for these model classes, yet without investigating the covariances and joint distribution of CV residuals) but in the following we will stick to GP as it is our main model of interest and the principles are readily applicable to other predictors falling under the settings of Proposition 1.

4 Some Consequences in GP Model Fitting

4.1 Graphical Diagnostics and Pivotal Statistics

We now wish to further explore the potential of the established formulas for diagnosing the quality of GP models. For simplicity, we will in what follows either remain generic or stick to simple kriging settings. A first immediate consequence of Proposition 1 and Corollary 2 for model fitting is that, under the hypothesis that observations are generated from the assumed model, cross-validation residuals can be transformed via appropriate matrix multiplications into a random vector with standard multivariate normal distribution. Assuming for simplicity that $\text{cov}({\mathbf{E}}_{{\mathbf{i}}_1,\dots ,{\mathbf{i}}_q})$ is full-ranked and denoting $t={\sum }_{j=1}^q\#{\mathbf{i}}_j$ it follows that $$\text{cov}{({\mathbf{E}}_{{\mathbf{i}}_1,\dots ,{\mathbf{i}}_q})}^{-1/2}{\mathbf{E}}_{{\mathbf{i}}_1,\dots ,{\mathbf{i}}_q}\sim \mathcal{N}\left(\mathbf{0},I_t\right).$$

If fact, any matrix $S\in {\mathbb{R}}^{t\times t}$ such that $S{B}^{-1}{\Sigma }^{-1}{B}^{-1}{S}^{\top }=I_n$ (remaining for simplicity under the final assumptions of Corollary 1 where $\Delta =I_n$) does the job. More specifically, one gets indeed with $S={\Sigma }^{1/2}B$ $$S{\mathbf{E}}_{{\mathbf{i}}_1,\dots ,{\mathbf{i}}_q}={\Sigma }^{-1/2}\mathbf{Z}\sim \mathcal{N}\left(\mathbf{0},I_n\right),$$ and hence the hypothesis that the model is correct can be questioned using any standard means relying on such a pivotal quantity with multivariate Gaussian distribution (e.g., by means of a ${\chi }^2$ test statistic or graphical diagnostics such as Q-Q plots). Similar considerations apply to the UK settings of Corollary 2.

Remark 8.

If Δ is not I_n but still an invertible matrix, then the above can be adapted with $S={\Sigma }^{1/2}{\Delta }^{-1}{B}^{-1}$. In case Δ is not invertible, analoguous approaches can be employed relying on the Moore-Penrose pseudo-inverse of $B^{-1}\Delta {\Sigma }^{-1}{\Delta }^{\top }{B}^{-1}$.

Coming back to our first one-dimensional example, Figure 3 represents Q-Q plots of the distribution of standardized LOO residuals (right panel) versus “transformed” LOO residuals (left panel). One may notice in particular, going from standardized to transformed LOO residuals, how the second point (corresponding to the second smallest theoretical quantile) switches from below the diagonal to above it.

Fig. 3 On the effect of accounting and correcting for correlation in Q-Q plots based on LOO residuals. Left panel: Q-Q plot against $\mathcal{N}(0,1)$ of LOO residuals merely divided by corresponding LOO standard deviations. Right panel: Q-Q plot against $\mathcal{N}(0,1)$ of duly transformed LOO residuals.

We now turn to the topic of covariance parameter estimation relying on cross-validation residuals. As we will see below, incorporating covariances between CV cross-validation residuals within such approaches leads to interesting developments.

4.2 Some Implications in Scale Parameter Estimation

Assuming a kernel of the form $k(\mathbf{x},\mathbf{x}\prime )={\sigma }^{2}r(\mathbf{x},\mathbf{x}\prime )$ with unknown ${\sigma }^2$ and $\Sigma =K$ (noiseless observations being for instance commonly encountered in the realm of computer experiments), the problem of estimating ${\sigma }^2$ by maximum likelihood versus cross-validation has been considered in earlier works, leading notably to results concerning the superiority of MLE in the well-specified case but also to situations when CV turned out to perform better than MLE in the misspecified case Bachoc (2013). Denoting by R the correlation matrix associated with the observation points ${\mathbf{x}}_1,\dots ,{\mathbf{x}}_n$, assumed full-rank, it is well-known (see for instance Santner, Williams, and Notz 2003) that the MLE of ${\sigma }^2$ can be written in closed form as (37) $${\widehat{\sigma }}_{\text{ML}}^2=\frac{1}{n}\mathbf{Z}{R}^{-1}\mathbf{Z}.$$(37)

In contrast, the leave-one-out-based estimator of ${\sigma }^2$ investigated in Bachoc (2013) reads (38) $${\widehat{\sigma }}_{\text{LOO}}^2=\frac{1}{n}\mathbf{Z}{R}^{-1}{(\text{diag}(R^{-1}))}^{-1}{R}^{-1}\mathbf{Z},$$(38) and originates from the idea, traced back by Bachoc (2013) to Cressie (1993), that the criterion (with notation $C_{\text{LOO}}$ inherited from Bachoc 2013) defined by (39) $$C_{\text{LOO}}^{(1)}({\sigma }^2)=\frac{1}{n}\sum _{i=1}^n\frac{{\mathbf{E}}_i^2}{{\sigma }^2{c}_{-i}^2},$$(39) should take a value close to one, where $c_{-i}^2=(s_{-i}^2)/{\sigma }^2$, leading to ${\widehat{\sigma }}_{\text{LOO}}^2=\frac{1}{n}{\sum }_{i=1}^n\frac{{\mathbf{E}}_i^2}{c_{-i}^2}$ and ultimately to (38). This estimator turns out to be unbiased as ${\widehat{\sigma }}_{\text{LOO}}^2$, yet with greater variance (see Bachoc 2013): (40) $$\begin{array}{c}\text{var}\left({\widehat{\sigma }}_{\text{LOO}}^2\right)=\frac{2{\sigma }^4\text{tr}({(R^{-1}{(\text{diag}(R^{-1}))}^{-1})}^2)}{n^2}\ge \frac{2{\sigma }^4}{n}\\ =\text{var}\left({\widehat{\sigma }}_{\text{ML}}^2\right).\end{array}$$(40)

In light of former considerations on the joint distribution of cross-validation residuals, it appears natural to revise $=\text{var}({\widehat{\sigma }}_{\text{ML}}^2).$ in order to correct for covariances between LOO residuals, resulting in (41) $$\begin{array}{c}{C}_{\widetilde{\text{LOO}}}^{(1)}({\sigma }^2)=\frac{1}{n}{\mathbf{E}}^{\top }{(B^{-1}Q{B}^{-1})}^{-1}\\ \mathbf{E}=\frac{1}{n{\sigma }^2}{\mathbf{E}}^{\top }\text{diag}(R^{-1})R\text{diag}(R^{-1})\\ \mathbf{E}=\frac{1}{n{\sigma }^2}{\mathbf{Z}}^{\top }{R}^{-1}\mathbf{Z},\end{array}$$(41) so that setting this modified criterion to 1 would plainly result in ${\widehat{\sigma }}_{\widetilde{\text{LOO}}}^2=\frac{1}{n}{\mathbf{Z}}^{\top }{R}^{-1}\mathbf{Z}={\widehat{\sigma }}_{\text{ML}}^2$ Hence, the advantages found for ${\widehat{\sigma }}_{\text{LOO}}^2$ in Bachoc (2013) in the case of model misspecification come at the price of neglecting covariances and hence possible redundancies between leave-one-out residuals, yet attempting to address this issue by naturally accounting for those covariances within the criterion leads to the MLE, known to enjoy optimality properties in the well-specified case of the considered settings but to be potentially suboptimal otherwise. We will see next that this is not the only instance where modifying LOO in order to account for residual covariances leads back to MLE.

4.3 On the Estimation of Further Covariance Parameters

Cross-validation has also been used to estimate covariance parameters beyond the scale parameter ${\sigma }^2$. We will now denote by θ those additional parameters (and by $\psi =({\sigma }^2,\theta )$ the concatenation of all covariance parameters) and review some approaches that have been used to estimate them via cross-validation. An important point compared to the last section on ${\sigma }^2$’s estimation is that Kriging predictors and CV residuals generally depend on θ, so that we will now stress it by noting $\mathbf{E}(\theta )$. Still in noiseless settings, the approach followed by Santner, Williams, and Notz (2003) and Bachoc (2013) to estimate θ based on leave-one-out cross-validation residuals is to minimize (42) $$C_{\text{LOO}}^{(2)}(\theta )=\sum _{i=1}^n{\mathbf{E}}_i{(\theta )}^2.$$(42)

Following Bachoc (2013), the latter criterion can be in turn expressed in closed form as (43) $$C_{\text{LOO}}^{(2)}(\theta )={\mathbf{Z}}^{\top }{R}_{\theta }^{-1}\text{diag}{(R_{\theta }^{-1})}^{-2}{R}_{\theta }^{-1}\mathbf{Z},$$(43) where $R_{\theta }$ stands for Z’s correlation matrix under correlation parameter θ. Based on our results, a natural extension of this criterion to multiple fold cross-validation can be similarly obtained. Considering q-fold settings such as in the final part of Corollary 1, $C_{\text{LOO}}^{(2)}$ becomes indeed (44) $$C_{\text{CV}}^{(2)}(\theta )=\sum _{j=1}^q{\left|\left|{\mathbf{E}}_{{\mathbf{i}}_j}(\theta )\right|\right|}^2={\left|\left|{\mathbf{E}}_{{\mathbf{i}}_1,\dots ,{\mathbf{i}}_q}(\theta )\right|\right|}^2.$$(44)

Then, building up upon our main theorem, we obtain that (45) $$\begin{array}{c}{C}_{\text{CV}}^{(2)}(\theta )={\mathbf{Z}}^{\top }{R}_{\theta }^{-1}\text{blockdiag}\\ \left({(R_{\theta }^{-1}[{\mathbf{i}}_1])}^{-1},\dots ,{(R_{\theta }^{-1}[{\mathbf{i}}_q])}^{-1}\right)R_{\theta }^{-1}\mathbf{Z},\end{array}$$(45) generalizing indeed (43). Note that applying a sphering transformation to the CV residuals previous to taking norms in (42) and (44) would lead to a criterion that boils down to $\theta \mapsto {\mathbf{Z}}^T{R}_{\theta }^{-1}\mathbf{Z}$ and hence appears as one of the two building blocks of the log-likelihood criterion.

On a different note and relaxing now the assumption of noiseless observations, (42) is only one among several possible ways to construct a criterion based on applying loss functions to leave-one-out residuals, as exposed in Rasmussen and Williams (2006). While (42) corresponds to the squared error loss, one arrives for instance by using instead the “negative log validation density loss” (after the terminology of Rasmussen and Williams (2006), sec. 5.4.2) at the “LOO log predictive probability,” also sometimes called pseudo-likelihood: (46) $$C_{\text{LOO}}^{(3)}(\psi )=\sum _{j=1}^n\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}\left(p_{Z\left[j\right]\left|\mathbf{Z}\right[-j]}(z_j|\mathbf{z}[-j];\psi )\right),$$(46) where $p_{Z[j]|\mathbf{Z}[-j]}(z[j]|\mathbf{z}[-j];\psi )$ denotes the conditional density of $Z[j]$ at the value $z[j]$ knowing that $\mathbf{Z}[-j]=\mathbf{z}[-j]$. Further scoring rules could be used as a base to define further criteria relying on cross-validation residuals, as illustrated for instance with the CRPS score in Petit et al. (2020) and Petit (2022).

Let us now focus on an extension of $C_{\text{LOO}}^{(3)}(\psi )$ to multiple fold cross-validation, and on the exploration of some links between the resulting class of criteria and the log-likelihood. First, $C_{\text{LOO}}^{(3)}(\psi )$ is straightforwardly adapted into (47) $$C_{\text{CV}}^{(3)}(\psi )=\sum _{j=1}^q\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}\left(p_{\mathbf{Z}\left[{\mathbf{i}}_j\right]\left|\mathbf{Z}\right[-{\mathbf{i}}_j]}(\mathbf{z}[{\mathbf{i}}_j\left]\right|\mathbf{z}[-{\mathbf{i}}_j];\psi )\right).$$(47)

Reformulating the sum of log terms as the log of a product, we see that (48) $$\begin{array}{c}{C}_{\text{CV}}^{(3)}(\psi )=\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}\left(\prod _{j=1}^q{p}_{\mathbf{Z}\left[{\mathbf{i}}_j\right]\left|\mathbf{Z}\right[-{\mathbf{i}}_j]}(\mathbf{z}[{\mathbf{i}}_j\left]\right|\mathbf{z}[-{\mathbf{i}}_j];\psi )\right)\\ =\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}\left(\prod _{j=1}^q{p}_{{\mathbf{E}}_{{\mathbf{i}}_j}}({\mathbf{e}}_{{\mathbf{i}}_j};\psi )\right)\end{array}$$(48) and so it appears that this approach amounts indeed in some sense to improperly assuming independence between the cross-validation residuals ${\mathbf{E}}_{{\mathbf{i}}_j}$ $(j=1\dots q)$.

We now examine in more detail (still in the settings of the final part of Corollary 1) when such cross-validation residuals are independent and show that in such case $C_{\text{CV}}^{(3)}(\psi )$ coincides with the log-likelihood criterion based on Z.

Corollary 3.

The following propositions are equivalent:

1. The cross-validation residuals ${\mathbf{E}}_{{\mathbf{i}}_j}$ $(j=1\dots q)$ are mutually independent,

2. $B=Q$,

3. the subvectors $\mathbf{Z}[{\mathbf{i}}_j]$ $(j=1\dots q)$ are mutually independent.

If they hold, then $C_{\text{CV}}^{(3)}(\psi )=\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}(p_{\mathbf{Z}}(\mathbf{z};\psi ))$ $(\psi \in \Psi )$.

Proof.

As we know from Corollary 1 that the cross-validation residuals ${\mathbf{E}}_{{\mathbf{i}}_j}$ $(j=1\dots q)$ form a Gaussian vector, they are mutually independent if and only if their cross-covariance matrices are null. Yet we also know by (34) from the same theorem that $\text{cov}({\mathbf{E}}_{\mathbf{i}},{\mathbf{E}}_{\mathbf{j}})={(Q[\mathbf{i}])}^{-1}Q[\mathbf{i},\mathbf{j}]{(Q[\mathbf{j}])}^{-1}(\mathbf{i},\mathbf{j}\in \mathcal{S})$. It follows that a) is equivalent to $Q[{\mathbf{i}}_j,{\mathbf{i}}_{j\prime }]=\mathbf{0}$ for $j\ne j\prime$ ($j,j\prime \in \{1,\dots ,q\}$), which coincides in turn with b). Looking finally at b) from the perspective of $\Sigma =Q^{-1}$ being block-diagonal, we obtain similarly that b) is equivalent to c). If these propositions hold, then ${\mathbf{E}}_{{\mathbf{i}}_j}=\mathbf{Z}[{\mathbf{i}}_j]$ $(j=1\dots q)$ from (25), and thus for $\psi \in \Psi$: $$\begin{array}{c}{C}_{\text{CV}}^{(3)}(\psi )=\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}\left(\prod _{j=1}^q{p}_{{\mathbf{E}}_{{\mathbf{i}}_j}}({\mathbf{e}}_{{\mathbf{i}}_j};\psi )\right)\\ =\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}\left(\prod _{j=1}^q{p}_{\mathbf{Z}\left[{\mathbf{i}}_j\right]}(\mathbf{z}[{\mathbf{i}}_j];\psi )\right)=\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}(p_{\mathbf{Z}}(\mathbf{z};\psi )).\end{array}$$

Also, one could have derived it from $p_{\mathbf{Z}[{\mathbf{i}}_j]|\mathbf{Z}[-{\mathbf{i}}_j]}(\mathbf{z}[{\mathbf{i}}_j]|\mathbf{z}[-{\mathbf{i}}_j];\psi )=p_{{\mathbf{Z}}_{{\mathbf{i}}_j}}(\mathbf{z}[{\mathbf{i}}_j];\psi )$. □

5 Computational Aspects

5.1 Accuracy and Speed-Up Checks for Fast Formulas

We have seen in previous sections that block matrix inversion enabled us to reformulate the equations of (multiple-fold) cross-validation residuals and associated covariances in compact form, from simple kriging to universal kriging settings. Here and in supplementary material, we briefly illustrate that our Schur-complement-based formuale do indeed reproduce results obtained by a straightforward (“naive”) computational approach, also keeping an eye on computational speed-ups and by-products. We focus below on how efficiency gains evolve depending on the number of folds. The next section focuses in more detail on the respective computational complexities of the two approaches. The settings are analogue to those presented in Section A.2. Accuracy and speed-up checks for LOO are presented in Section A.4 (also in supplementary material).

In the following experimental results and theoretical considerations, we focus on the particular case of homogeneous fold sizes and denote by $r\ge 1$ their common number of elements, so that $n=r\times q$. For convenience, we consider in experiments a design of size $n=2^{10}=1024$, and let the number of folds decrease from q = 1024 to q = 2 by successive halvings. This amounts to folds of sizes r ranging from 1 to 512, correspondingly. For each value of r, 50 random replications of the folds are considered; this is done by permuting the indices at random prior to arranging them in regularly constructed folds of successive r elements. The resulting speed-up distributions are represented in Figure 4, in function of r. With speed-up medians ranging from 1027.82 to 1.50 (means from 1429.63 to 1.28, all truncated after the second digit), we hence observed a speed-up whatever the value of r, yet with a substantial decrease when r increases compared to the LOO situation, as could be expected.

Fig. 4 Relative errors on CV means and covariances (between naive and fast methods, with the naive method as reference) measured in similar settings as on Figure 3.

Additionally, relative errors on cross-validation residuals and associated variances calculated using the fast versus naive approach are represented in Figure 5, also in function of r. It is remarkable that, while the relative error on the calculation of residuals increases moderately (with r) from a median of around $3.5\times {10}^{-14}$ to values nearing $4\times {10}^{-14}$, in the case of the covariances the variation is more pronounced, with an increase from approx. $2\times {10}^{-11}$ to $1.2\times {10}^{-10}$. Further analysis may shed light on the underlying phenomena, yet one can state that with a maximum relative error of magnitude ${10}^{-10}$, the observed differences remain of a negligible extent.

Fig. 5 Speed-up (ratio between times required to run the naive and fast methods) measured for q-fold CV, where q decreases from 1024 to 2 and 50 seeds are used that affect here both model fitting and the folds.

5.2 More on Complexity and Efficiency Gains

We now turn to a more specific examination on how our closed form formulas affect computational speed compared to a naive implementation of multiple fold cross-validation. As we will now develop, a first order reasoning on computational costs associated with naive versus closed multiple-fold cross-validation highlights the existence of two regimes depending on q: for large q it is faster to use close form formula yet for very small q one should better employ the naive approach. The rationale is as follows; in a naive implementation, the main computational cost is on inverting q covariance sub-matrices, all belonging to

${\mathbb{R}}^{(n-q)\times (n-q)}$. Denoting this (approximate) cost by ${\kappa }_{\text{naive}}$ and assuming a cubic cost for matrix inversion (with a multiplicative constant $\gamma >0$), we hence obtain that (49) $${\kappa }_{\text{naive}}=q\times \gamma {(n-r)}^3.$$(49)

On the other hand, the closed form approach requires one inversion of a n × n matrix and q inversions of r × r matrices, leading to an approximate cost (50) $${\kappa }_{\text{closed}}=\gamma (n^3+q{r}^3).$$(50)

We thus obtain an approximate speed up of $\frac{{\kappa }_{\text{naive}}}{{\kappa }_{\text{closed}}}=\frac{{(q-1)}^3}{(q^2+1)}$, implying that close form approach would be advantageous when ${(q-1)}^3\ge (q^2+1)$. As q is a positive integer, the latter condition be proved to be equivalent to $q\ge 4$, as $q\to {(q-1)}^3-(q^2+1)$ possesses a unique real root between 3 and 4, takes negative values for $q\in \{1,2,3\}$, and tends to $+\infty$ when $q\to +\infty$. In our implementation, however, the cost of the fast approach turns out to be smaller than $\gamma (n^3+q{r}^3)$ since the Cholesky factor is already pre-calculated. Modeling this with a damping factor of $\alpha \in (0,1]$ in front of the $\gamma n^3$ term, we end up with an approximate cost of ${\kappa }_{\text{close}}=\gamma (\alpha n^3+q{r}^3)$, leading via analogue calculations to an acceleration whenever ${(q-1)}^3\ge (\alpha q^2+1)$. Again, the polynomial involved possesses a single real-valued root, with a value shrinking toward 2 as α tends to 0.

In practice, speed-ups are already observed here from q = 2, which may be attributed to various reasons pertaining notably to implementation specifics and arbitrary settings used in the numerical experiments. Also, let us stress again that the reasoning done above is to be taken as an approximation as it does not account for matrix-vector multiplication costs, storage and communication overheads, and further auxiliary operations that could influence the actual running times.

On a different note, let us point out that the computational benefits highlighted above in the case of multiple-fold cross-validation with a unique partitioning do carry over and even further grow in the case of $N\ge 2$ replications with varying partitions. Sticking for simplicity to the previously considered settings (α = 1) but replicating the partitions, we find that a naive approach has a cost driven by (51) $${\kappa }_{\text{naive}}^{\text{rep}}=\mathit{\text{Nq}}\gamma {(n-r)}^3,$$(51) while a similar approach using the fast formulas would have a cost driven by (52) $${\kappa }_{\text{closed}}^{\text{rep}}=\gamma (n^3+{\mathit{\text{Nqr}}}^3).$$(52)

Then we obtain that ${\kappa }_{\text{closed}}^{\text{rep}}\le {\kappa }_{\text{naive}}^{\text{rep}}$ as soon as ${(q-1)}^3\ge \frac{q^2}{N}+1$, occurring already from q = 3.

6 Application Test Case

We now illustrate and investigate some applications of multiple-fold cross-validation on an application test case motivated by a contaminant source localization problem from hydrogeology Pirot et al. (2019). There the goal is to minimize an objective function that is quantifying the discrepancy between given contaminant concentrations at monitoring wells and analogue quantities obtained under varying scenarios regarding the potential source underlying this contamination. Solving this difficult optimization problem in order to localize the source is performed in Pirot et al. (2019) via Bayesian Optimization, and Gaussian process models are thus used to surrogate the objective function. An instance of such an objective function is presented in Figure 6. Our focus is on fitting a GP model to this function under a clustered experimental design, and to highlight notable differences arising here between leave-one-out and multiple-fold cross-validation using the cluster structure.

Fig. 6 Contaminant localization test function designed by summing dissimilarities between given concentrations at monitoring wells and corresponding simulation results when varying the candidate source localization (See Pirot et al. (2019) for more detail). The filled disks stand for the experimental design, consisting of 25 clover-shaped clusters of size 5. The small green circle stands for the actual (prescribed) contaminant source, and the red square indicates the grid point with the smallest misfit.

We wish to illustrate some differences between CV in the LOO versus in “Leave-One-Cluster-Out” settings (i.e., 25-fold CV with folds $\equiv$ clusters). Since we can evaluate the actual function on a fine grid, we can easily calculate the absolute prediction errors and compare them to absolute cross-validation residuals. The mean absolute error on the considered 2601-element grid is $2.27\times {10}^{-2}$. The mean absolute error of the 125 LOO residuals is $0.42\times {10}^{-2}$. In contrast, the mean absolute error of the 125 multiple-fold CV residuals is $3.49\times {10}^{-2}$. Thus, while absolute multiple-fold CV residuals appears to give a rather over-pessimistic picture of the prediction situation (as CV predictions are done only at points at points in left-out clusters), analogue LOO residuals give a very over-optimistic one (as CV predictions are done only at left-out points within clusters including four not left-out points).

Figure 7 represents LOO (upper panel) versus MFCV (medium panel) residual norms as a function of the range parameter and compares them to the true reconstruction error (lower panel). This example illustrates that the design of folds may crucially affect results. In the present case, it can be observed that, compared to the curve associated with LOO residuals, the curve associated with MFCV residuals is more similar to the actual (out-of-sample) reconstruction error. In contrast, the top curve associated with LOO features a flatter region near the optimum resulting in turn in an optimal range associated with a larger reconstruction error. Replications of the example with further realizations of the clustered designs delivered quite variable results (overall, neither the considered MFCV nor LOO stood up as evidently preferable to the other). Yet some recurring features emerged, in particular with the curves associated with LOO appearing more prone to decreasing to a plateau quite similarly as on Figure 7. This calls for more research on the impact of fold design on resulting indicators and estimators.

Fig. 7 Log square norm of LOO residuals (top), of MFCV residuals (center), and of reconstruction error (bottom) as a function of the range hyperparameter.

7 Discussion

Fast leave-one-out formulas for Gaussian process prediction that have been used for model diagnostics and hyperparameter fitting do carry over well to multiple-fold cross-validation, as our theoretical and numerical results in both simple and universal kriging settings suggest. The resulting formulas were found to allow substantial speed-ups, yet with a most favorable scenario in the leave-one-out case and decreasing advantages in cases with lesser folds of larger observation subsets. A first order analysis in terms of involved matrix inversion complexities corroborated the main trends observed in numerical experiments, although some room remains to analyze actual computation times more precisely depending on implementation specifics. In addition, established formulas enabled a closed-form calculation of the covariance structure of cross-validation residuals, and eventually correcting for an improper assumption of independence silently underlying QQ-plots with standardized LOO residuals.

As we established as well, the established formulas lend themselves well to generalizing existing covariance hyperparameter estimation approaches that are based on cross-validation residuals. Looking first at scale parameter estimation in the LOO case, we found that changing a criterion underlying the scale estimation approach used in Santner, Williams, and Notz (2003) and Bachoc (2013) by correcting for covariance between LOO residuals led back to the maximum likelihood estimator of scale. On a different note and with more general covariance hyperparameters in view, as established in Corollary 3 maximizing the pseudo-likelihood criterion mentioned in Rasmussen and Williams (2006) was found to coincide with MLE in the case of independent cross-validation residuals.

It is interesting to note as a potential starting point to future work on cross-validation-based parameter estimation that going beyond the independence assumption of Corollary 3 and accounting for the calculated residual covariances results in fact in a departure from basic Maximum Likelihood Estimation. Still assuming indeed that $\Delta =I_n$ and replacing the product in $C_{\text{CV}}^{(3)}(\psi )$ by the joint density of cross-validation residuals, we end up indeed with (53) $$\begin{array}{cc}{C}_{\widetilde{\text{CV}}}^{(3)}(\psi )\hfill & =\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}\left(p_{({\mathbf{E}}_{{\mathbf{i}}_1},\dots ,{\mathbf{E}}_{{\mathbf{i}}_q})}(({\mathbf{e}}_{{\mathbf{i}}_1},\dots ,{\mathbf{e}}_{{\mathbf{i}}_q});\psi )\right)\hfill \\ \hfill & =\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}\left(p_{B^{-1}Q\mathbf{Z}}(\mathbf{e};\psi )\right)\hfill \\ \hfill & =\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}\left(p_{B^{-1}Q\mathbf{Z}}(B^{-1}Q\mathbf{z};\psi )\right)\hfill \\ \hfill & =\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}\left(p_{\mathbf{Z}}(\mathbf{z};\psi )|\mathrm{\text{det}}(B^{-1}Q)|\right)\hfill \\ \hfill & =\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}\left(p_{\mathbf{Z}}(\mathbf{z};\psi )\right)-\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}\left(\mathrm{\text{det}}(B)\right)+\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}\left(\mathrm{\text{det}}(Q)\right),\hfill \end{array}$$(53) and we see that maximizing $C_{\widetilde{\text{CV}}}^{(3)}(\psi )$ generally departs indeed from MLE, as Q and B are functions of ψ, yet coincides with it when $\mathrm{\text{det}}(B)=\mathrm{\text{det}}(Q)$. While studying further this new criterion is out of scope of the present article, we believe that using the full distribution of cross-validation residuals could be helpful in designing novel criteria and procedures for covariance parameter estimation and more. Beyond this, efficient adjoint computation of gradients could be useful to speed-up covariance hyperparameter estimation procedures such as studied here, but also variations thereof relying on further scoring rules Petit et al. (2020), not only in LOO settings but also potentially in the case of multiple-fold cross-validation.

It is important to stress that while the presented results do not come with procedures for the design of folds, they might be of interest to create objective functions for their choice. Looking at the covariance matrix between cross-validation residuals might be a relevant ingredient to such procedures. Of course, designing the folds is not completely unrelated to choosing the ${\mathbf{x}}_i$’s and we expect interesting challenges to arise at this interface. Also, stochasic optimization under random folds is yet another approach of interest which could be eased thanks to fast computation of cross-validation residuals. Last but not least, as fast multiple-fold cross-validation carries over well to linear regression and GP modeling shares foundations with several other paradigms, it would be desirable to generalize presented results to broader model classes (see for instance Rabinowicz and Rosset (2020b), where mixed models are tackled).

Acknowledgments

Part of DG’s contributions have taken place within the Swiss National Science Foundation project number 178858, and he would like to thank Sam Allen, Athénaïs Gautier, Fabian Guignard, Philip Stange, Cédric Travelletti, Riccardo Turin as well as many attendees of events where this work was presented for stimulating scientific exchanges. Yves Deville should be warmly thanked too for insightful comments on an earlier version of this article. Warm thanks as well to the two referees and the members of the editorial team, whose comments have enabled improving the paper substantially. Calculations were performed on UBELIX, the High Performance Computing (HPC) cluster of the University of Bern. DG would also like to aknowledge support of Idiap Research Institute, his primary affiliation during earlier stages of this work.

Supplementary Materials

In the supplement, we provide extra detail about block inversion via Schur complements, as well as about the basics of universal kriging. Also, we give detail on the opening 1-dimensional example and present additional numerical experiments and results. The latter pertain on the one hand to the effect on cross-validation diagnostics of a change of design (when points are “doubled” by adding close neighbors and thus creating clusters) and, on the other hand, to LOO and MFCV speed-ups and relative errors depending on the size of considered designs.

Disclosure Statement

The authors report there are no competing interests to declare.