Calculation of leave-one-out cross-validation statistic in linear regression when design matrix is square

Question

Suppose we are performing a least-squares multiple linear regression of the form $$ \mathbf{Y}=\mathbf{X}\boldsymbol{\beta}+\boldsymbol{\varepsilon}\,, $$ where $\mathbf{Y}=(y_1,y_2,...,y_N)$ are the response variables, $\mathbf{X}$ is the $N \times P$ design matrix, $\boldsymbol{\beta}=(\beta_1,\beta_2,...,\beta_P)$ are the parameters and $\boldsymbol{\varepsilon}$ is the error.

The leave-one-out cross-validation statistic $CV$ is defined as $$ CV=\frac{1}{N}\sum_{i=1}^N y_i-\hat{y}_{[i]}\,,\tag{1}\label{1} $$ where $\hat{y}_{[i]}$ is the predicted value of $y_i$ obtained by fitting the model to all observations except the $i$th one.

According to various sources—I'm following Rob Hyndman's notation here, but it's also described in e.g. An Introduction to Statistical Learning—a well-known "shortcut" to calculating $CV$ in the context of ordinary least squares without having to re-fit the model $N$ times is

$$ CV=\frac{1}{N}\sum_{i=1}^{N}\left(\frac{y_i-\hat{y}_i}{1-h_i}\right)^2\,,\tag{2}\label{2} $$

where $\hat{y}_i$ is the predicted value of $y_i$ from the "full" fitted model using all observations. Here, $h_i$ is the leverage statistic for observation $i$, which Hyndman writes is equal to the $i$th entry on the diagonal of the "hat-matrix" $\mathbf{H}=\mathbf{X} (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'$.

In the case where $N=P$ and therefore $\mathbf{X}$ is a square matrix, all $y_i=\hat{y}_i$ and all $h_i=1$, and the "shortcut" equation above is undefined. However, $CV$ can still be obtained using equation \eqref{1} above.

Is there a "shortcut" to obtain $CV$ without refitting the model $N$ times like in equation \eqref{2}, but which is defined when $\mathbf{X}$ is a square matrix?

Answer 1

You're quite right that when $N = P$, all $y_i = \hat y_i$. However, in this case, $\hat y_{[i]}$ is actually not defined, so you will need to (slightly!) modify your problem for it to make sense. To see why, recall that linear regression estimates the coefficients by minimizing the loss $$L(\beta) = \|\mathbf{Y} - \mathbf{X} \beta\|_2^2 = \sum_{i=1}^N (y_i - X_i^T \beta)^2$$ where here $X_i$ are the rows of $\mathbf{X}$. When $N < P$, there are many possible choices of $\beta$ which ensure that $X_i^T \beta = y_i$ for each observed datapoint and thus achieve a loss of zero. Thus, there is no unique solution. Of course, as you've pointed out, in this case $\hat y_i = y_i$ even when $N < P$ for the observed data-points, but for new datapoints, different choices of $\beta$ yield different predictions. Thus, if we leave the $i$th data-point out, there is no unique choice of the estimated coefficients based on the first $N-1$ datapoints, and therefore no unique choice of $\hat y_{[i]}$.

In general, the standard multiple linear regression technique will perform poorly when $N$ is close to $P$. A better approach in high dimensions is to regularize linear regression, for example by using a ridge estimator: $$\hat \beta_{\mathrm{ridge}} = \arg\min_{\beta} \sum_{i=1}^N (y_i - X_i^T \beta)^2 + \lambda \|\beta\|_2^2$$ which penalizes choices of $\beta$ that have large norms. Ridge regression is very nice, and there are some ways to get efficient expressions of the leave-one-out cross validation error; however, depending how large $N$ and $P$ are, this is an active area of research, so perhaps the best thing for me to do is refer you to two references. First, An Introduction to Statistical Learning by James, Witten, Hastie, and Tibshirani is a classical textbook on this subject and covers ridge regression in detail in Chapter 6 (it is freely available online). More specifically, Lecture notes on ridge regression by WN van Wieringen addresses the leave-one-out cross validation question; it is
linked here.