Residual Sum of squares in Weighted regression

Question

The Residual Sum of squares (RSS) in Weighted regression is written as $$(\mathbf{y-X\hat{\boldsymbol\beta}})^{'}\mathbf{C}^{-1}(\mathbf{y-X\hat{\boldsymbol\beta}})$$ Where $$\hat{\boldsymbol\beta}=(\mathbf{X^{'}C^{-1}X})^{-1}\mathbf{X^{'}C^{-1}y}$$

I am trying to write the RSS in an efficient manner which reduces computational complexity, for example I am able to write the RSS as follows $$(\mathbf{y-X\hat{\boldsymbol\beta}})^{'}\mathbf{C}^{-1}(\mathbf{y-X\hat{\boldsymbol\beta}})=tr(\mathbf{e^{'}C^{-1}e})=tr(\mathbf{e^{'}eC^{-1}})=tr(\mathbf{EC^{-1}})$$ where $tr=trace, \mathbf{e}=(\mathbf{y-X\hat{\boldsymbol\beta}})$ and $E=\mathbf{e^{'}e}$ Although this expression seems mathematically simple however the computational complexity is the same

I hope anyone can help me find an efficient way to write and code the RSS in weighted regression. Also, I would appreciate a reference to an $R$ function that can find this RSS so that I can take a look how the expression is written.

PS : $C$ need not to be a diagonal matrix, however it is symmetric positive semi-definite e.g. a covariance matrix

Answer 1

There are several good ways to do this using R. One classical method is to compute the Choleski factor of the covariance matrix:

R <- chol(C)
yc <- backsolve(R, y, transpose=TRUE)
Xc <- backsolve(R, X, transpose=TRUE)
fit <- lm.fit(Xc, yc)
RSS <- sum(fit$effects[-(1:fit$rank)]^2)

This code requires C to be a strictly positive definite matrix. C has to be positive definite anyway in order to guarantee that the RSS is finite.

If you want to make the computation even more explicit, you could replace the last two lines with this:

QR <- qr(Xc)
e <- qr.qty(QR, yc)
e <- e[-(1:QR$rank)]
RSS <- sum(e^2)

The above code is performing the following mathematical steps. First, we factorize $$C = R^TR$$ where $R$ is an upper triangular matrix. Then we solve the linear systems $$R^Ty_c=y$$ and $$R^TX_c=X$$ for $y_c$ and $X_c$ using an efficient forward substitution algorithm. Note that the above two steps are far more efficient than inverting $C$.

From this point we can view this as an unweighted regression problem with $y_c$ and $X_c$. Amongst other things, the lm.fit function uses the QR decomposition of $X_c$ to find an $n\times(n-p)$ matrix $Q$ such that $Q^TQ=I$ and $Q^TX=0$. Here, $p$ is the column rank of $X$. The orthogonal residuals (or effects) can then be computed as $$e=Q^Ty$$ and finally the RSS is $e^Te$. Actually the function computed $Q^Ty$, where $Q$ was $n\times n$, and stored this vector in fit$effects. We then threw away the first $p$ values to get $e$.

You might have been hoping for a simpler mathematical formula, but efficient computation requires that one avoids evaluating mathematical entities such as inverse matrices or ordinary residuals.

Answer 2

Big idea: solve linear systems rather than compute matrix inverse and multiply

Computing the inverse of a large matrix is slow and imprecise. If you have a huge matrix, don't compute the inverse. Instead:

• Observe $A^{-1} \mathbf{b}$ is the solution to $A\mathbf{x} = \mathbf{b}$. More generally, $A^{-1}B$ is the solution to $AX=B$. Anytime you see, $A^{-1}B$ you should think, "Hey, maybe I should solve linear systems rather than compute the matrix inverse?" (Note: for small problems, it probably doesn't matter.)

• Another idea: solving triangular systems is super fast! $$\left[\begin{array}{ccc} r_{11}&r_{12} & r_{13}\\0&r_{22}&r_{23}\\0&0&r_{33} \end{array} \right] \left[ \begin{array}{c}x_1\\x_2\\x_3\end{array} \right] = \left[ \begin{array}{c}b_1\\b_2\\b_3\end{array} \right] $$ Notice you immediately have $x_{3} = b_3 / r_{33}$. Then substitute $x_3$ and second row $x_2r_{22}+x_3r_{23}=b_2$ immediately gives you $x_2$. Computers can solve triangular systems at lightning speed applying this backsolve algorithm. (Don't reinvent the wheel, call Linear Algebra libraries to do this.) This is what Gordon is getting at when advocating taking a Cholesky decomposition on $A$. Once you have the Cholesky Decomposition $A=RR'$ where $R$ is lower triangular, you can solve $A\mathbf{x}=\mathbf{b}$ with two back substitutions, which is nearly free.

$$\hat{\boldsymbol\beta}=(\mathbf{X^{'}C^{-1}X})^{-1}\mathbf{X^{'}C^{-1}y}$$

If you Cholesky decomposition for $C$ such that $ C = RR' $ and $C^{-1} = R'^{-1}R^{-1}$ then:

$$\hat{\boldsymbol\beta}=(X^{'}R'^{-1}R^{-1} X)^{-1}X^{'}R'^{-1}R^{-1}y$$

So you want to solve:

• $R^{-1}X$ with a backsolve
• $R^{-1}y$ with a backsolve
• and then one final linear system

Practical coding way forward

Since there's already an answer in R, I'll give some code in Matlab. In Matlab, there's the operator \ which basically examines your matrix and does the right thing to solve a linear system fast.

C_inv_X = C \ X;                     % Solves C^{-1}X
bhat = (X'*C_inv_X) \ (C_inv_X'*y);  % Final linear system solve

You could also:

R       = chol(C);  %Careful... does your library give you an UT or LT?
                    %chol gives you an upper triangular in Matlab
Rp_inv_X = R'\X; %mldivide auto-detects triangular system
                 % in other environments, you may have to declare it to
                 % get full speed
Rp_inv_y = R'\y;
bhat    = (Rp_inv_X'*Rp_inv_X) \ (Rp_inv_X'*Rp_inv_y);  % Final linear system solve

And then you have the triangular matrix $R$ around to use later.