Residual Sum of squares in Weighted regression

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

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.

• Could you please provide a brief mathematical interpretation on what this code is doing – Wis Sep 9 '16 at 4:58
• The simple way to think about it is that for the unweighted case you have an initial problem like $X\beta\approx y$. For the weighted version, if you "divide" both sides of the above by the "matrix square root" of the covariance, then everything else is exactly the same. Specifically you have $\hat{X}\beta\approx\hat{y}$, where $\hat{X}=C^{-1/2}X$ and $\hat{y}=C^{-1/2}y$. The Cholesky factor is the simplest way to get a "matrix square root" ($C=LL^T$, where $L$ is a lower triangular matrix, the "square root".) – GeoMatt22 Sep 9 '16 at 5:31
• @GeoMatt22 What is the computational complexity of a backsolve operator ? – Wis Sep 16 '16 at 1:46
• @raw5 not sure what the context is? For the above answer, the Cholesky factor $L=R^T$ is triangular, so $Lx=b$ can be solved in a single pass, i.e. "complexity" of O[# non-zeros in $L$]. Wait, did you mean to ask the other Matt this, as a comment to his answer? (This comment thread is on Gordon's answer!) – GeoMatt22 Sep 16 '16 at 1:56
• @GeoMatt22 I am sorry I meant to ask Gordon , sorry for that and thanks for your input – Wis Sep 16 '16 at 4:20

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.

• I think there was a query for you in a comment to Gordon's answer? (It was directed to me by accident I think, but I answered.) – GeoMatt22 Sep 16 '16 at 2:01