Residual Sum of Squares from Generalized Least Squares (GLS) always 0?

Question

Standard set-up: $$ Y=X\beta+e $$ where $e\sim N(0,\Sigma)$

We know that the GLS estimate of $\beta$ is $\hat{\beta}=(X'\Sigma^{-1}X)^{-1}X'\Sigma^{-1}Y$

The (generalized) residual sum of square is then: \begin{align} SS(\hat{\beta})&=(Y-X\hat{\beta})'\Sigma^{-1}(Y-X\hat{\beta})\\ &=(Y-X(X'\Sigma^{-1}X)^{-1}X'\Sigma^{-1}Y)'\Sigma^{-1}(Y-X(X'\Sigma^{-1}X)^{-1}X'\Sigma^{-1}Y)\\ &=Y'\Sigma^{-1}Y-Y'\Sigma^{-1}X(X'\Sigma^{-1}X)^{-1}X'\Sigma^{-1}Y \end{align}

The second equality is by substituting $\hat{\beta}$ into the LHS, the third equality is by brute force expanding the brackets but I think it is just standard algebra stuff one do with least square problems. In fact, if $\Sigma$ is the Identity matrix it becomes $Y'(I-H)Y$ where $H$ is the 'hat- matrix'

The thing I got quite confused now is that, if I am correct, $$ \Sigma^{-1}X(X'\Sigma^{-1}X)^{-1}X'\Sigma^{-1}=\Sigma^{-1}\quad(*) $$ which would make $SS(\hat{\beta})=0$.

Which does not make sense, as $SS(\hat{\beta})$ is the objective function we try to minimize. Can someone tell me what I did wrong?

The way I got $(*)$ is by a bunch of factorization (inspired by this post):

We first use cholesky to get $\Sigma=LL'$ And QR to get $L^{-1}X=QR$ Therefore \begin{align} &\Sigma^{-1}X(X'\Sigma^{-1}X)^{-1}X'\Sigma^{-1}\\ =&(LL')^{-1}X(X'(LL')^{-1}X)^{-1}X'(LL')^{-1}\\ =&(LL')^{-1}X(R'Q'QR)^{-1}X'(LL')^{-1}\\ =&L'^{-1}QR(R'R)^{-1}R'Q'L^{-1}\\ =&L'^{-1}L^{-1}\\ =&\Sigma^{-1} \end{align}

Answer 1

I found by problem by generating random matrices and checking every equality. The problem is in the second last equality. Turns out, $$ L'^{-1}QR(R'R)^{-1}R'Q'L^{-1}\neq L'^{-1}L^{-1} $$ because $$ R(R'R)^{-1}R' =\begin{bmatrix} I_p&0\\ 0&0 \end{bmatrix} \neq I $$ where $I_p$ is the identity matrix of dimension $p<N$, where $N$ is the dimension of the full matrix.

Answer 2

There is often some ambiguity when people talk about the QR factorization of a rectangular matrix. Do we mean the full QR factorization or the "economy" size? The full QR factorization can be written as: $$ X = QR = Q \begin{bmatrix} R_1 \\ 0 \end{bmatrix} = \begin{bmatrix} Q_1 & Q_2 \end{bmatrix} \begin{bmatrix} R_1 \\ 0 \end{bmatrix} = Q_1 R_1, $$ where $X$ is $m \times n$, $R_1$ is an $n \times n$ upper triangular matrix, $0$ is an $(m − n) \times n$ zero matrix, $Q_1$ is $m \times n$, $Q_2$ is $m \times (m − n)$, and $Q_1$ and $Q_2$ both have orthogonal columns. For purposes of computation, we typically only care about $Q_1$ and $R_1$ (the "economy" size). It's important to note that, while $Q$ is orthogonal, $Q_1$ is not (it's not even square; though it is true that $Q_1^T Q_1 = I$).

With that in mind, let's work through the computation of the residual. I'll use the subscripts to specify that I'm using the "economy" QR factorization. And I'll use some results from the previous post, namely $\Sigma = L L^T$, $L^{-1} X = Q_1 R_1$ and $\beta = R_1^{-1} Q_1^T L^{-1} Y$ (and their transposes). $$ \begin{array}{} (X \beta - Y)^T \Sigma^{-1} (X \beta - Y) & = & \beta^T X^T \Sigma^{-1} X \beta - 2 Y^T \Sigma^{-1} X \beta + Y^T \Sigma^{-1} Y \\ & = & Y^T L^{-T} Q_1 R_1^{-T} X^T L^{-T} L^{-1} X R_1^{-1} Q_1^T L^{-1} Y - 2 Y^T L^{-T} L^{-1} X R_1^{-1} Q_1^T L^{-1} Y + Y^T L^{-T} L^{-1} Y \\ & = & Y^T L^{-T} (Q_1 R_1^{-T} X^T L^{-T} L^{-1} X R_1^{-1} Q_1^T - 2 L^{-1} X R_1^{-1} Q_1^T + I) L^{-1} Y \\ & = & Y^T L^{-T} (Q_1 R_1^{-T} R_1^T Q_1^T Q_1 R_1 R_1^{-1} Q_1^T - 2 Q_1 R_1 R_1^{-1} Q_1^T + I) L^{-1} Y \\ & = & Y^T L^{-T} (Q_1 Q_1^T - 2 Q_1 Q_1^T + I) L^{-1} Y \\ & = & Y^T L^{-T} (-Q_1 Q_1^T + I) L^{-1} Y \\ \end{array} $$ At this point, you might think we're stuck. But there is one more trick up our sleeve. Recall that $Q = \begin{bmatrix} Q_1 & Q_2 \end{bmatrix}$ is orthogonal. So we have: $$ \begin{array}{} Q Q^T & = & \begin{bmatrix} Q_1 & Q_2 \end{bmatrix} \begin{bmatrix} Q_1 & Q_2 \end{bmatrix}^T \\ & = & \begin{bmatrix} Q_1 & Q_2 \end{bmatrix} \begin{bmatrix} Q_1^T \\ Q_2^T \end{bmatrix} \\ & = & Q_1 Q_1^T + Q_2 Q_2^T \\ & = & I \end{array} $$ Thus $I - Q_1 Q_1^T = Q_2 Q_2^T$ and we can write: $$ \begin{array}{} (X \beta - Y)^T \Sigma^{-1} (X \beta - Y) & = & Y^T L^{-T} Q_2 Q_2^T L^{-1} Y \\ & = & ||Q_2^T L^{-1} Y||_2^2 \end{array} $$ Computationally, I'm not sure this is very useful. If you want to compute the residual, you're probably better off first computing $\beta$ and then just plugging it in. However, this form can offer you some insight as to where your residual is coming from. In the full QR factorization $X = \begin{bmatrix} Q_1 & Q_2 \end{bmatrix} \begin{bmatrix} R_1 \\ 0 \end{bmatrix}$, the columns of $Q_2$ form a basis for the null space of $X$. In the linear least squares application, the null space of $X$ is all the stuff your model does not explain. As you add more terms to your model, the null space of $X$ gets smaller and smaller (assuming the terms you're adding are linearly independent).