Obtaining Residual Sum of Squares from a large OLS problem using a naive sequential approach - why doesn't it work?

Question

Suppose we have an OLS problem with a large number of predictors: $Y = X_1 + X_2 + \cdots + X_p$. I want to obtain its RSS. I don't need to know the regression coefficients or individual residuals, just RSS.

A naive way to go about it is first to regress $Y$ on $X_1$, take the residuals, regress them on $X_2$, take the residuals, regress them on $X_3$, and so on. The final regression of residuals on $X_p$ is expected to produce the same RSS as the original model. I experimented in R and that way does not result in the correct RSS.

I understand why this approach doesn't work for getting the regression coefficients: if it were to work, then the time complexity of OLS would be much smaller because it's much faster to fit 100 regressions with one predictor each than one regression with 100 predictors. However, I still don't understand the theoretical reason why we can't get the RSS that way.

Answer 1

Note: It is possible to compute the residual sum-of-squares directly, without fitting the regression model. To do this, you first compute the coefficient of determination directly from the correlations between the relevant vectors (see e.g., here) and compute the total sum-of-squares from the response vector, then you compute the residual sum-of-squares from these quantities. In the present answer I will ignore this method and deal with the problem you are having with your iterative method.

This approach is an example of iterative regression modelling, and the resulting phenomena you are observing is related to the Frisch-Waugh-Lovell theorem in regression. At the moment, your problem is that you are adding each new explanatory variable in as a raw vector, without first modifying it to remove the parts that are already attributable to the previously fit explanatory variables. If you want to do this iterative procedure correctly, you would need to regress each new explanatory variable on the previous explanatory variable in addition to updating the residuals for the response variable in this way. Here is how a corrected version of the procedure would work:

Start the procedure by regressing $\mathbf{y}$ on $\mathbf{x}^{(1)} = \mathbf{x}_1$ to get residuals $\mathbf{r}_1$.

For $t=2,...,p$ do the following:

• Regress $\mathbf{x}_t$ on $\mathbf{x}^{(t-1)}$ to get residuals $\mathbf{s}_t$,
• Set the new design matrix $\mathbf{x}^{(t)} = [\mathbf{x}^{(t-1)}, \mathbf{s}_{t}]$,
• Regress $\mathbf{r}_{t-1}$ on $\mathbf{s}_t$ to get residuals $\mathbf{r}_{t}$.

Finally, compute $\text{RSS} = ||\mathbf{r}_p||^2$.

This algorithm gives an iterative update of the residuals adding one explanatory variable to the model at a time. As you can see, in addition to updating the residuals from the response variable, the algorithm also keeps track of an iterated design matrix that consists of the initial explanatory vector and then the residuals from regressing the later explanatory vectors against each previous design matrix.

---

Computational Implementation: We can implement this iterative algorithm in R to onfirm that it gives the correct RSS. Below we create mock regression data with $n=100$ data points and $p=10$ explanatory variables and we then compute the RSS using the full regression and the iterative method. (In the present case we use a model that includes an intercept term, but the code can be altered to remove this term if preferred.)

#Set parameters
n <- 100
p <- 10
sigma <- 3
#Generate some regression data
set.seed(1)
X     <- matrix(rnorm(np), nrow = n)
BETA  <- rnorm(p)
ERROR <- sigmarnorm(n)
Y     <- X %*% BETA + ERROR
#Fit full regression and compute RSS
MODEL.FULL  <- lm(Y ~ X)
RSS.FULL    <- anova(MODEL.FULL)[2,2]
#Fit iterative regression and compute RSS
XX <- X[,1]-mean(X[,1])
RESID <- vector(mode = "list", length = p)
RESID[[1]] <- lm(Y ~ XX)$residuals
for (t in 2:p) {
  SS <- lm(X[,t] ~ XX)$residuals
  XX <- cbind(XX, SS)
  RESID[[t]] <- lm(RESID[[t-1]] ~ SS)$residuals }
RSS.ITER    <- sum(RESID[[p]]^2)

Here is the outcome of the computation:

RSS.FULL
[1] 789.7609
RSS.ITER
[1] 789.7609
RSS.FULL - RSS.ITER
[1] 0

This confirms that the iterative procedure yields the correct RSS value.

Answer 2

So the conceptual question you are asking is: I understand that to get the linear coefficients, I need to orthogonalise the X variables. But I am only interested in the RSS. So if the residual has the component along $X_i$ removed, then the residual is orthogonal to each component up to column i.

The issue is that removing the component along 1 nonorthogonal direction introduces a residual with a component in the previous components

[Please just try drawing this example]

consider Y = (1,0) and $X_1 = (1,1)$ and $X_2=(0,1)$

$resid_1 = (.5, -0.5)$

but $resid_2 = (0.5,0)$ whereas Y is in the span of X so the residual should be zero

Answer 3

Here is how you can try to get RSS without finding the regression coefficients

Lets say the problem is:

$$ \begin{align} \mathbf{Y}&=\mathbf{X}\cdot\boldsymbol{\beta}+\boldsymbol{\rho}\\ \left(\begin{array}\\Y_1\\Y_2\\\vdots\\Y_n\end{array}\right)&= \left(\begin{array}\\ X_{11} & X_{12} & \dots & X_{1p} \\ X_{21} & X_{22} & \dots & X_{2p} \\ \vdots \\ X_{n1} & X_{n2} & \dots & X_{np} \\ \end{array}\right)\cdot\left(\begin{array}\\\beta_1\\\vdots\\\beta_p\end{array}\right)+\left(\begin{array}\\\rho_1\\\rho_2\\\vdots\\\rho_n\end{array}\right) \end{align} $$

Where $\boldsymbol{\rho}$ are your residuals. You can decompose $\mathbf{X}$ using singular value decomposition into an $n\times n$ orthogonal matrix $\mathbf{U}$, a $p\times p$ orthogonal matrix $\mathbf{V}$ and a diagonal matrix $\mathbf{D}$. It is a good idea to impose regularization condition where you define certain diagonal elements small enough to be zero:

$$ \mathbf{X}=\mathbf{U}\cdot \left(\begin{array}\\ d_1 & 0 & 0& 0 & 0 &\dots & 0 \\ 0 & d_2 & 0& 0 & 0 &\dots & 0 \\ 0 & 0 & \ddots& 0 & 0 &\dots & 0 \\ 0 & 0 & 0& d_s & 0 &\dots & 0 \\ 0 & 0 & 0& 0 & 0 &\dots & 0 \\ \dots\\ 0 & 0 & 0& 0 & 0 &\dots & 0 \\ \end{array}\right) \cdot \mathbf{V} $$

One can then apply an orthogonal transformation:

$$ \begin{align} \mathbf{Z}&=\mathbf{U}^T \cdot \mathbf{Y} \\ \mathbf{r}&=\mathbf{U}^T \cdot \boldsymbol{\rho} \\ \boldsymbol{\alpha}&=\mathbf{V}^T\cdot\boldsymbol{\beta} \end{align} $$

Then we have:

$$ \begin{align} \left(\begin{array}\\Z_1\\Z_2\\\vdots\\Z_s\end{array}\right)&= \left(\begin{array}\\ d_1 & 0 & \dots & 0\\ 0 & d_2 & \dots & 0\\ 0 & 0 & \ddots & 0 \\ 0 & 0 & \dots& d_s\\ \end{array}\right)\cdot\left(\begin{array}\\\alpha_1\\\vdots\\\alpha_s\end{array}\right)+\left(\begin{array}\\r_1\\\vdots\\r_s\end{array}\right)\\ \left(\begin{array}\\Z_{s+1}\\\vdots\\Z_{n}\end{array}\right)&= \left(\begin{array}\\r_{s+1}\\\vdots\\r_{n}\end{array}\right) \end{align} $$

You have a freedom to choose any $\alpha_{\dots}$ and, by construction $d_{1\dots s}\neq 0$, so residuals $r_1\dots r_s$ can all be set to zero. Due to orthogonality of $\mathbf{V}$ then:

$$ \boldsymbol{\rho}^T\cdot\boldsymbol{\rho}=\rho_1^2+\dots+\rho^2_n=\mathbf{r}^T\cdot\mathbf{r}=Z_{s+1}^2+\dots+Z_{n}^2=\sum_{i=1\dots n} Y_i^2-\left(Z_1^2+\dots + Z_s^2\right) $$

Let $\mathbf{v}_j$ be the j-th column vector of matrix $\mathbf{V}$. It then follows that:

$$ \mathbf{X}^T\cdot\mathbf{Y}=\left(\mathbf{v}_1\,\mathbf{v}_2\,\dots\,\mathbf{v}_s\right)\cdot \left(\begin{array}\\ d_1 & 0 & \dots & 0\\ 0 & d_2 & \dots & 0\\ 0 & 0 & \ddots & 0 \\ 0 & 0 & \dots& d_s\\ \end{array}\right)\cdot \left(\begin{array}\\Z_1\\Z_2\\\vdots\\Z_s\end{array}\right) $$

So that:

$$ Z_{i=1\dots s}=\frac{\mathbf{v}_i^T\cdot \mathbf{X}^T\cdot \mathbf{Y}}{d_i} $$

Thus:

$$ \rho_1^2+\dots+\rho^2_n=\sum_{i=1\dots n} Y_i^2-\sum_{j=1\dots s}\frac{\left(\mathbf{v}_j^T\cdot \mathbf{X}^T\cdot \mathbf{Y}\right)^2}{d_j^2} $$

This is what you are after, I believe. You still need a way of extracting $\mathbf{v}$ and $d$

Finally, the non-normalized covariance matrix is:

$$ \left(n-1\right)\mathbf{C}=\mathbf{X}^T\cdot\mathbf{X}=\mathbf{V}^T\cdot \left(\begin{array}\\ d_1^2 & 0 & 0& 0 & 0 &\dots & 0 \\ 0 & d_2^2 & 0& 0 & 0 &\dots & 0 \\ 0 & 0 & \ddots& 0 & 0 &\dots & 0 \\ 0 & 0 & 0& d_s^2 & 0 &\dots & 0 \\ 0 & 0 & 0& 0 & 0 &\dots & 0 \\ \dots\\ 0 & 0 & 0& 0 & 0 &\dots & 0 \\ \end{array}\right)\cdot\mathbf{V} $$

So all one needs to do is diagonalize the covariance matrix, and work with its non-zero eigenvalues and eigenvectors