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

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.

• can you add your code. Something like that should work. see hastie.su.domains/Papers/ESLII.pdf section 3.2.3 Multiple regression from simple univariate regression page 52 Sep 11, 2023 at 16:07
• Yes, there is a sequential algorithm for OLS, e.g. Algorithm 3.1 in Hastie et al, but it's complexity is the same. In step 2, there are up to (p - 1) predictors. Since I need only RSS, I hoped I could get it by running p simple linear regressions as I described above, but it doesn't work. Sep 11, 2023 at 17:09
• What is the motivation? Regression is pretty easy for computers so why not do it all at once? I had a case a while ago with many variables which had so many rows that it could not be loaded into the memory. In that case I found that it was sufficient to compute all pairwise covariances, then do SVD on the design matrix and RSS could be computed from that, without any actual regression step. Is this what you are after?
– Cryo
Sep 11, 2023 at 19:19
• The reason is that p is too large in the original model. The OLS complexity is O(p^3), so it takes too much time. Sep 11, 2023 at 19:26
• Libraries for neural networks routinely operate with millions of features, and have tools to schedule and distribute those computations. Can your use case be solved in full, but with these libraries? OLS is simply few matrix operations + solving a linear system afterall
– Cryo
Sep 12, 2023 at 4:47

## 3 Answers

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

• Thank you, that's exactly what I was after. That also explains why there are simple, formula-based solutions for n-way ANOVA, but they only applied to balanced designs, i.e. when each column or at least the groups of columns corresponding to different terms are orthogonal to one another. Sep 13, 2023 at 2:31

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(n*p), nrow = n)
BETA  <- rnorm(p)
ERROR <- sigma*rnorm(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.

• This should be the accepted answer. (Minor, you write you create $n=10$ mock data points.) Sep 12, 2023 at 11:27
• Thanks a lot for replying. Your solution is similar to one in Hastie, and its complexity is the same as that of original model. I hoped there was a way around it to get just the RSS, but it seems impossible without getting the regression coefficients into the bargain, hence the complexity stays the same. Sep 12, 2023 at 21:45
• @ChristophHanck: Thanks for that --- edited. Please feel free to edit to correct simple typos.
– Ben
Sep 12, 2023 at 23:15
• Re complexity, the simplest way to compute $RSS$ is to bypass the model-fit altogether and first compute the $R^2$ using the linked formula. This is a single matrix computation that involves solving a system of $p$ linear equations. It is then simple to obtain $RSS = (1-R^2) ||\mathbb{y}||^2$.
– Ben
Sep 14, 2023 at 22:18

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