A.1 Leave-one-out cross validation for linear regression in one step

Now, in each step of the leave-one-out procedure, one column of $\mathbf{X}$ and the corresponding element of $\bm{y}$ are omitted. If we call the missing column $\bm{x}_{n}$, this omission evidently turns $\mathbf{X}\bm{y}$, the second factor in the normal equations (Eq. A.1), into $(\mathbf{X}\bm{y}-\bm{x}_{n}y_{n})$. The first term, $\mathbf{X}\mathbf{X}^{\text{T}}$, can be expressed as a sum of outer products, $\sum_{n}\bm{x}_{n}\bm{x}_{n}^{\text{T}}$, so losing $\bm{x}_{n}$ simply omits one outer product: $\mathbf{X}\mathbf{X}^{\text{T}}-\bm{x}_{n}\bm{x}_{n}^{\text{T}}$. Inverting this quantity with the Sherman-Morrison formula, Eq. B.15, and bringing together with the second term, gives:

for the regression coefficients at the $n^{\text{th}}$ step.

Let us write the prediction at the $n^{\text{th}}$ step for the held out datum, namely ${\hat{y}}_{\text{LOO},n}\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.} }}=\bm{x}_{n}^{\text{T}}\bm{w}_{\text{LOO},n}$. It is easily seen that $\bm{x}_{n}^{\text{T}}{\mathopen{}\mathclose{{}\left(\mathbf{X}\mathbf{X}^{\text{T}}}\right)}^{-1}\bm{x}_{n}=\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.} }}h_{nn}$ is the $n^{\text{th}}$ diagonal element of the hat matrix, Eq. A.3, which allows us to simplify the expression for the prediction:

where on the final line we use the fact, from Eq. A.4, that the diagonal entries of the residual matrix are $m_{nn}=1-h_{nn}$. Collecting together the predictions for the held-out data at all the steps, we have

We can clean up this equation by giving names to the diagonal matrices constructed from the diagonals of $\mathbf{M}$ and $\mathbf{H}$; that is,

in which case

From the leave-out-out predictions it is simple to calculate the leave-one-out residuals:

Eqs. A.6 and A.7 are pleasingly elegant, and yield intuitive interpretations. To begin with, note that in a standard linear regression, the element $h_{nn}$ of the hat matrix $\mathbf{H}$ tells us how much weight the observation $y_{n}$ gets in predicting itself under the normal equations.¹¹ 1 Put this way, the claim may sound strange—shouldn’t $y_{n}$ get a weight of 1, and all other samples be weighted with zeros? Remember: the model requires that the vector of samples be “compressed” into the space of the inputs, $\bm{x}$, before expanding back out into predictions. As long as there are fewer dimensions than samples (as is typically the case in regression), the hat matrix cannot be the identity matrix: $\mathbf{X}^{\text{T}}$ will be “tall,” in which case (see Eq. A.3) there does not exist a matrix $\mathbf{P}$ such that $\mathbf{X}^{\text{T}}\mathbf{P}=\mathbf{I}$. Intuitively, we expect this weight to be between 0 and 1, and indeed we prove this in the next section. Consequently, although perhaps less intuitively, the diagonal elements of the residual matrix, $m_{nn}$, are also between 0 and 1 (since $m_{nn}=1-h_{nn}$). They tell us what fraction of $y_{n}$ turns into error in the prediction of itself.

Eq. A.6 tells us how to transform the standard predictions ($\bm{\hat{y}}$) into leave-one-out predictions ($\bm{\hat{y}}_{\text{LOO}}$). First, we remove from each prediction $\bm{\hat{y}}_{n}$ the contribution from $y_{n}$ (to which by assumption we have no access!), namely $h_{nn}y_{n}$. The vector of all contributions to be removed is $\mathbf{H}_{\text{d}}\bm{y}$, so the complete adjustment is $\bm{\hat{y}}-\mathbf{H}_{\text{d}}\bm{y}$. Since $0\leq h_{nn}\leq 1$ for all $n$, this shrinks $\bm{\hat{y}}$, and therefore we will have to scale it back up to the right size. This is accomplished by the pre-multiplication with $\mathbf{M}_{\text{d}}^{-1}$, which scales the $n^{\text{th}}$ entry of $\bm{\hat{y}}-\mathbf{H}_{\text{d}}\bm{y}$ by $1/m_{nn}$. Since $0\leq m_{nn}\leq 1$, this scaling will indeed increase the magnitude of the prediction. And the amount of scaling applied is sensible: If, for example, the original (non-cross-validated) prediction of $y_{n}$ relied heavily on $y_{n}$—that is, if $h_{nn}$ is close to 1—then the corresponding $m_{nn}$ would be close to zero, which makes sense: the prediction would be scaled up considerably to make up for the missing $y_{n}$ from the inputs. At the other extreme, if the original prediction of $y_{n}$ relied mostly on the other samples $y_{m\neq n}$—that is, if $h_{nn}$ is close to zero—then the corresponding $m_{nn}$ would be close to 1: very little scaling would be applied to ${\hat{y}}_{n}-h_{nn}y_{n}$, which is already close to ${\hat{y}}_{n}$.

It is easily verified that $\mathbf{H}^{\text{T}}\mathbf{H}=\mathbf{H}$, and therefore

Letting $\bm{v}$ be a “one-hot” vector, so that the first term just picks out a diagonal element, we see that $h_{nn}\geq 0$. Furthermore, for all vectors $\bm{v}$ of length 1, the quadratic form $\bm{v}^{\text{T}}\mathbf{H}\bm{v}$ is maximal when $\bm{v}$ is the leading eigenvector, $\lambda_{\text{max}}$, in which case $\bm{v}^{\text{T}}\mathbf{H}\bm{v}=\lambda_{\text{max}}$. Since one-hot vectors have length 1, $\bm{v}^{\text{T}}\mathbf{H}\bm{v}\leq\lambda_{\text{max}}$ for any one-hot vector (with equality only when the leading eigenvector is precisely this one-hot vector). It only remains to show that $\lambda_{\text{max}}$ is no more than 1, which follows from the fact that $\mathbf{H}$ is idempotent, i.e. $\mathbf{H}\mathbf{H}=\mathbf{H}$. For any eigenvector $\bm{u}$, we have

which can only be satisfied when $\lambda$ is -1, 0, or 1 (assuming non-trivial eigenvectors). Therefore, $0\leq h_{nn}\leq 1$.

Since the diagonal elements of the residual matrix, $m_{nn}$, are just 1 minus the diagonal elements of hat matrix, we see that $0\leq m_{nn}\leq 1$ as well.

With a little effort, the analysis can be extended to more general folded cross validation, in which the folds can consist of more than one sample, although must still be disjoint. The left out samples from fold $k$ then form a matrix which we will call $\mathbf{X}_{k}$, and which for concreteness we imagine (without loss of generality) to comprise consecutive columns of the data matrix $\mathbf{X}$. Note that the derivation will not assume that $\mathbf{X}_{k}$ contains the same number of columns as block $\mathbf{X}_{j}$ for any other $j$, although usually one lets all folds except the last contain an equal number of samples. Now, to invert a matrix involving this matrix $\mathbf{X}_{k}$ rather than vector $\bm{x}_{n}$, we will need the full-blown Woodbury inversion lemma (Eq. B.14):

where we have given a name to the $k^{\text{th}}$ block on the diagonal of the hat matrix:

Now, it is easily verified that the bracketed quantity on the final line is merely $(\mathbf{I}-\mathbf{H}_{k})^{-1}$ (since multiplying it by $\mathbf{I}-\mathbf{H}_{k}$ yields the identity), so it follows that the $k^{\text{th}}$ leave-out-fold-out predictions are

Correspondingly, the $k^{\text{th}}$ leave-out-fold-out residuals are

where we have named the $k^{\text{th}}$ block $\mathbf{I}-\mathbf{H}_{k}=\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.} }}\mathbf{M}_{k}$, analogously to the leave-one-out residual matrix. Now if we define $\mathbf{M}_{\text{bd}}$ to be the block diagonal matrix with the blocks $\mathbf{M}_{k}$ on the diagonal, we can write all the residuals at once as

since a block diagonal matrix can be inverted by inverting its blocks.

Eq. A.9 looks reassuringly like Eq. A.7, but the resemblance is somewhat misleading: In practice we would not (typically) invert $\mathbf{M}_{\text{bd}}$, but would rather invert each block individually. Thus although it is possible to compute these residuals in “one step,” one typically takes $K$ steps, like the naïve calculation. But Eq. A.9 will still be cheaper than the naïve calculation, which inverts (much) larger matrices at each step.