Leave-one-out cross-validation and error correction

Getting the inverse of a submatrix

Suppose we have a matrix K and know its inverse K⁻¹. Suppose that K has block structure $$ K = \begin{bmatrix} A~B \\ C~D \end{bmatrix}$$ Now we want to find out how to find A⁻¹ using K⁻¹ instead of doing the full inverse. We can write K⁻¹ in block structure $$K^{-1} = \begin{bmatrix} E~F \\ G~H \end{bmatrix}$$

Now we use the fact that KK⁻¹ = I $$ \begin{bmatrix} I~0 \\ 0~I \end{bmatrix} = \begin{bmatrix} A~B \\ C~D \end{bmatrix}\begin{bmatrix} E~F \\ G~H \end{bmatrix} $$

This gives the equations AE + BG = I AF + BH = 0 CE + DG = 0 CF + DH = I

Solving the first equation gives that A = (I − BG)E⁻¹ or A⁻¹ = E(I − BG)⁻¹

Leave-one-out covariance matrix inverse for Gaussian processes

For Gaussian processes we can consider the block matrix for the covariance (or correlation) matrix where a single row and its corresponding column is being removed. Let the first n − 1 rows and columns be the covariance of the points in design matrix X, while the last row and column are the covariance for the vector x with X and x. Then we can have

$$ K = \begin{bmatrix} C(X,X)~ C(X,\mathbf{x}) \\ C(\mathbf{x},X)~C(\mathbf{x},\mathbf{x}) \end{bmatrix}$$

Using the notation from the previous subsection we have A = C(X, X) and B = C(X, x), and E and G will be submatrices of the full K⁻¹. B is a column vector, so I’ll write it as a vector b, and G is a row vector, so I’ll write it as a vector g^T. So we have C(X, X)⁻¹ = E(I − C(X, x)G)⁻¹ So if we want to calculate A⁻¹ = E(I − bg^T)⁻¹ we still have to invert I − BG, which is a large matrix. However this can be done efficiently since it is a rank one matrix using the Sherman-Morrison formula. $$ (I - \mathbf{b}\mathbf{g}^T)^{-1} = I^{-1} - \frac{I^{-1}\mathbf{b}\mathbf{g}^TI^{-1}}{1+\mathbf{g}^TI^{-1}\mathbf{b}} = I - \frac{\mathbf{b}\mathbf{g}^T}{1+\mathbf{g}^T\mathbf{b}} $$ Thus we have the shortcut for A⁻¹ that is only multiplication $$ A^{-1} = E (I - \frac{\mathbf{b}\mathbf{g}^T}{1+\mathbf{g}^T\mathbf{b}})$$

To speed this up it should be calculated as

$$ A^{-1} = E - \frac{(E\mathbf{b})\mathbf{g}^T}{1+\mathbf{g}^T\mathbf{b}}$$ Below demonstrates that we get a speedup of almost twenty by using this shortcut.

set.seed(0)
corr <- function(x,y) {exp(sum(-30*(x-y)^2))}
n <- 200
d <- 2
X <- matrix(runif(n*d),ncol=2)
R <- outer(1:n,1:n, Vectorize(function(i,j) {corr(X[i,], X[j,])}))
Rinv <- solve(R)
A <- R[-n,-n]
Ainv <- solve(A)
E <- Rinv[-n, -n]
b <- R[n,-n]
g <- Rinv[n,-n]
Ainv_shortcut <- E + E %*% b %*% g / (1-sum(g*b))
summary(c(Ainv - Ainv_shortcut))

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -30245.9     -0.1      0.0      0.0      0.1  43135.9

if (requireNamespace("microbenchmark", quietly = TRUE)) {
  microbenchmark::microbenchmark(solve(A), E + E %*% b %*% g / (1-sum(g*b)))
}

## Unit: microseconds
##                                expr     min       lq      mean   median
##                            solve(A) 914.578 936.3385 1188.4202 985.4555
##  E + E %*% b %*% g/(1 - sum(g * b)) 142.787 144.3145  189.9717 147.8365
##        uq      max neval
##  1395.751 3343.372   100
##   285.213  332.751   100

In terms of the covariance matrices already calculated, this is the following, where M_−i is the matrix M with the ith row and column removed, and M_i, −i is the ith row of the matrix M with the value from the ith column removed.

$$ R(X_{-i})^{-1} = R(X)_{-i} - \frac{(R(X)_{-i}R(X)_{-i,i}) (R(X)^{-1})_{i,-i}^T }{1 + (R(X)^{-1})_{i,-i}^T R(X)_{-i,i}}$$

Leave-one-out prediction

Recall that the predicted mean at a new point is ŷ = μ̂ + R(x, X₂)R(X₂)⁻¹(y₂ − μ1_n))

ŷ = μ̂ + R(x_i, X₋₁)R(X_−i)⁻¹(y_−i − μ1_n))