Method of least squares, first order condition and QR decomposition

Question

When you use the method of least squares you estimate the parameters in the following way:

$$\min_{\mathbf{b}} (\mathbf{y} - \mathbf{X}\mathbf{b})^T(\mathbf{y} - \mathbf{X}\mathbf{b})$$

Where $\mathbf{y}_{n \times 1}$, $\mathbf{X}_{n \times (p + 1)}$ and $\mathbf{b}_{(p+1) \times 1}$

If you solve the problem you obtain the following first order condition:

$$\mathbf{X}^T\mathbf{X}\mathbf{b} = \mathbf{X}^T\mathbf{y}$$

According to R version 4.3.1 to find $\mathbf{b}$ QR decomposition is used (see ?lm in relation to method = "qr"). Thefore we have that if $\mathbf{X}$ is of full column rank it can be expressed as $\mathbf{X} = \mathbf{QR}$ where $\mathbf{Q}$ is a orthogonal matrix with dimensions $n \times (p+1)$ and $\mathbf{R}$ is a upper triangular matrix with dimensions $(p + 1) \times (p + 1)$. According to wikipedia (here, here and here) we have the following:

• $\mathbf{Q}\mathbf{Q}^T = \mathbf{Q}^T\mathbf{Q} = \mathbf{I}$
• For $\mathbf{R}$ we have that $r_{ij} = 0$ for $i > j$

Then applying this to the first order condition we have that:

$$\mathbf{X}^T\mathbf{X}\mathbf{b} = \mathbf{X}^T\mathbf{y}$$ $$(\mathbf{QR})^T\mathbf{QR}\mathbf{b} = (\mathbf{QR})^T\mathbf{y}$$ $$\mathbf{R}^T\mathbf{Q}^T\mathbf{QR}\mathbf{b} = \mathbf{R}^T\mathbf{Q}^T\mathbf{y}$$ $$\mathbf{Q}^T\mathbf{QR}\mathbf{b} = \mathbf{Q}^T\mathbf{y}$$ $$\mathbf{Rb} = \mathbf{Q}^T\mathbf{y}$$ $$\mathbf{b} = \mathbf{R}^{-1}\mathbf{Q}^T\mathbf{y}$$

Using the following reproducible example you can see that the above result applies:

model <- lm(mpg ~ wt + hp, data = mtcars)
X <- model.matrix(object = model)
QR <- qr(x = X)
R <- qr.R(qr = QR)
Q <- qr.Q(qr = QR)
y <- as.matrix(mtcars['mpg'])
b <- solve(R) %*% t(Q) %*% y
b
#>                     mpg
#> (Intercept) 37.22727012
#> wt          -3.87783074
#> hp          -0.03177295
coefficients(model)
#> (Intercept)          wt          hp 
#> 37.22727012 -3.87783074 -0.03177295

^{Created on 2023-11-04 with reprex v2.0.2}

Then we have the following result:

$$\mathbf{\hat{y}} = \mathbf{Xb} = \mathbf{X}\mathbf{R}^{-1}\mathbf{Q}^T\mathbf{y} = \mathbf{QR}\mathbf{R}^{-1}\mathbf{Q}^T\mathbf{y} = \mathbf{Q}\mathbf{Q}^T\mathbf{y} = \mathbf{y}$$

Which is totally false because you dont have always that $\mathbf{\hat{y}} = \mathbf{y}$.

Taking into account this I have the following questions:

• ¿Is really the definition of a orthogonal matrix this $\mathbf{Q}\mathbf{Q}^T = \mathbf{Q}^T\mathbf{Q} = \mathbf{I}$ or it is only this $\mathbf{Q}^T\mathbf{Q} = \mathbf{I}$?
• If the definition of orthogonal matrix is $\mathbf{Q}\mathbf{Q}^T = \mathbf{Q}^T\mathbf{Q} = \mathbf{I}$ why using the following reproducible example you get this in R version 4.3.1:

model <- lm(mpg ~ wt + hp, data = mtcars)
X <- model.matrix(object = model)
QR <- qr(x = X)
Q <- qr.Q(qr = QR)
dim(Q)
#> [1] 32  3
# t(Q) %*% Q is an identity matrix 
identity_matrix <- t(Q) %*% Q
# Q %*% t(Q) is not the identity matrix
not_an_identity_matrix <- Q %*% t(Q)

^{Created on 2023-11-04 with reprex v2.0.2}

Where Q %*% t(Q) is not an identity matrix but t(Q) %*% Q is an identity matrix

Answer 1

Taking into account the comments by @PBulls, @Zhanxiong and @whuber I was able to understand the problem I have and I was able to find the following reference The QR decomposition of a matrix that corresponds to the Hyper-Textbook: Optimization Models and Applications by Laurent El Ghaoui.

So there are 2 types of QR decomposition:

• Case when $X$ is full column rank
• Full QR decomposition

In the first case for $X$ you have only that $\mathbf{Q}^T\mathbf{Q} = \mathbf{I}$ so you get only the following result: $\mathbf{\hat{y}} = \mathbf{Q}\mathbf{Q}^T\mathbf{y}$

Using the following reproducible example you can see that the above result applies:

model <- lm(mpg ~ wt + hp, data = mtcars)
X <- model.matrix(object = model)
QR <- qr(x = X)
Q <- qr.Q(qr = QR)
y <- as.matrix(mtcars['mpg'])
y_hat <- Q %*% t(Q) %*% y 
fitted_values <- as.matrix(model$fitted.values)

^{Created on 2023-11-04 with reprex v2.0.2}

Where y_hat and fitted_values are the same

The case of the Full QR decomposition is not used in this particular context but you can check this topic for other purposes in The QR decomposition of a matrix > Full QR decomposition.