Gaussian process regression: leave-one-out prediction According to Dubrule's Cross validation of kriging in a unique neighborhood, it is possible to compute leave-one-out the gaussian process prediction $\hat{Y}_{-i}(x_i)$ at a point $x_i$ from the design of experiments $\mathcal{X}$, where the kriging prediction is done by removing $x_i$ from $\mathcal{X}$. The formula for the prediction mean is 
$$
\hat{Y}_{-i}(x_i) = y_i - [K^{-1}y]_i/[K^{-1}]_{ii}
$$
assuming $y$ is the vector of measurements and $K$ the covariance matrix $R(X,X)$, $R$ being the covariance function. I was thinking of computing this quantity in the case of noisy measurements i.e. $y_i = f(x_i) + \epsilon_i$, $\epsilon_i \sim_{i.i.d.} \mathcal{N}(0,\sigma_\epsilon^2) $ independent of $x$. 
Should there be any significant modification or does it suffice it to replace $K$ by $K + \sigma_\epsilon^2I_n$. 
I ask this because all gaussian process toolboxes seem to ignore the case of noisy measurements. Perhaps because the above generalization does not hold any longer ?
 A: In the general noisy or ``signal $+$ noise'' framework $y_i = f(\mathbf{x}_i)
+ \epsilon_i$, several observations can be done at the same location
$\mathbf{x}_i$, so the notations $Y(\mathbf{x}_i)$ and $f_{-i}$ 
can be misleading then. 
Suppose first that the $n$ locations $\mathbf{x}_i$ are distinct, so
that deleting an observation and deleting a location are the same
thing.  Using gaussian conditional expectation, the predicted (or
smoothed) value for the signal vector $\mathbf{f}$ is
$\widehat{\mathbf{f}} = \mathbf{H} \mathbf{y}$ where $\mathbf{H}:=
\mathbf{K}(\mathbf{K} + \sigma^2_\epsilon \mathbf{I}_n)^{-1}$.  Then we can use
Cook and Weisberg's formula for a linear smoother
$$
    y_i - \widehat{f}_{-i} = \frac{y_i - \widehat{f}_i }{1 - H_{ii}}
$$
which requires $H_{ii} \neq 1$.  This rewrites in a form similar to
yours $\widehat{f}_{-i} = y_i - [\mathbf{B}\mathbf{y}]_i /
[\mathbf{B}]_{ii}$ with
$
\mathbf{B} := \mathbf{I}_n - \mathbf{H}  =
 \sigma^2_\epsilon (\mathbf{K} + \sigma^2_\epsilon \mathbf{I}_n)^{-1}
$ 
and we get your formula after simplifying $\sigma^2_\epsilon$ in the fraction.
So in this case it is enough to replace the covariance $\mathbf{K}$ by
$\mathbf{K} + \sigma^2_\epsilon \mathbf{I}_n$.
In the general case, the $n$ observations $y_i$ are related to $p
\leqslant n$ signal values $f_k:=f(\mathbf{x}_k)$ at $p$ distinct
points.  Most likely, a cross-validation will then leave out all the
$y_i$ at the same location $\mathbf{x}_k$ (and not only one $y_i$).  Then we are back to the first case by averaging the
observations at the same location and by modifying consequently the
variance. For the location $\mathbf{x}_k$, the variance
$\sigma^2_\epsilon$ must be replaced by $\sigma^2_\epsilon/n_k$ where
$n_k$ is the number of observations at site $k$, hence $\sigma^2_\epsilon
\mathbf{I}_n$ above will be replaced by a diagonal matrix.  This general
context can be regarded as a Bayes linear regression where the
parameter $\mathbf{f}$ has prior
$\text{Norm}(\mathbf{0},\,\mathbf{K})$ and $\sigma^2_\epsilon$ is
known, so conjugacy holds and we have update formulas that can be used
for the deletion of one observation if needed.
The Cook and Weisberg's formula above (p. 33 of their famous 1982 book
Residuals and Influence in Regression) holds for smoothing splines
as well as for linear regression. It can be proved by using the 
Woodbury matrix identiy 
or by purely statistical arguments involving optimality (e.g. minimal
variance). In
$$
  \widehat{f}_i = H_{ii} \,y_i + (1- H_{ii}) \sum_{j \neq i} \frac{H_{ij}}{1- H_{ii}} \,y_j,
$$
the sum $\sum$ at right hand side can be identified as the optimal
prediction of $f_i$ based on the $y_j$ with $ j\neq i$.  Indeed it is
unbiased, and it is also optimal because $\widehat{f}_i$ would
otherwise no longer be such. The formula does not hold as such for
non-noisy (interpolating) kriging because then $H_{ii} = 1$.
