[This article was first published on R – Aaron Schlegel, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

The more common approach to QR decomposition is employing Householder reflections rather than utilizing Gram-Schmidt. In practice, the Gram-Schmidt procedure is not recommended as it can lead to cancellation that causes inaccuracy of the computation of $$q_j$$, which may result in a non-orthogonal $$Q$$ matrix. Householder reflections are another method of orthogonal transformation that transforms a vector $$x$$ into a unit vector $$y$$ parallel with $$x$$. The Householder reflection matrix with normal vector $$v$$ takes the form:

$$H = I – 2vv^T$$

Thus we need to build $$H$$ so that $$Hx = \alpha e_1$$ for some constant $$\alpha$$ and $$e_1 = \begin{bmatrix} 1 & 0 & 0 \end{bmatrix}^T$$.

Since $$H$$ is orthogonal, $$||Hx|| = ||x||$$ and $$||\alpha e_1|| = |\alpha| ||e_1|| = |\alpha|$$. Therefore, $$\alpha = \pm ||x||$$. The sign is selected so it has the opposite sign of $$x_1$$ (we’ll use $$+$$ for the remaining definitions). The vector $$u$$ we seek is thus:

$$u = \begin{bmatrix} x_1 + sign(x_1) ||x_1|| \\\ x_2 \\\ \vdots \\\ x_n \end{bmatrix}$$

With the unit vector $$v$$ defined as $$u = \frac{v}{||v||}$$. The corresponding Householder reflection is then:

$$H(x) = I – 2vv^T = I – 2 \frac{uu^T}{u^Tu}$$
QR Decomposition with Householder Reflections

The Householder reflection method of QR decomposition works by finding appropriate $$H$$ matrices and multiplying them from the left by the original matrix $$A$$ to construct the upper triangular matrix $$R$$. As we saw earlier, unlike the Gram-Schmidt procedure, the Householder reflection approach does not explicitly form the $$Q$$ matrix. However, the $$Q$$ matrix can be found by taking the dot product of each successively formed Householder matrix.

$$Q = H_1 H_2 \cdots H_{m-2}H_{m-1}$$

Consider the matrix $$A$$.

$$A = \begin{bmatrix} 2 & – 2 & 18 \\\ 2 & 1 & 0 \\\ 1 & 2 & 0 \end{bmatrix}$$

We find the reflection of the first column vector $$a_1 = \begin{bmatrix}2 & 2 & 1 \end{bmatrix}^T$$, $$v_1 = a_1 + sign(a_{11})||a_1||e_1$$.

$$v_1 = \begin{bmatrix}2 \\\ 2 \\\ 1 \end{bmatrix} + \sqrt{\sum^m_{k=1}{a_1}^2} \begin{bmatrix}1 \\\ 0 \\\ 0 \end{bmatrix} = \begin{bmatrix} 5 \\\ 2 \\\ 1 \end{bmatrix}$$

With a corresponding Householder matrix:

$$H_1 = \begin{bmatrix} 1 & 0 & 0 \\\ 0 & 1 & 0 \\\ 0 & 0 & 1 \end{bmatrix} – 2 \frac{\begin{bmatrix} 2 & 2 & 1 \end{bmatrix}\begin{bmatrix} 2 \\\ 2 \\\ 1 \end{bmatrix}}{\begin{bmatrix} 2 \\\ 2 \\\ 1 \end{bmatrix}\begin{bmatrix} 2 & 2 & 1 \end{bmatrix}} =$$
$$H_1 = \begin{bmatrix} -\frac{2}{3} & -\frac{2}{3} & -\frac{1}{3} \\\ -\frac{2}{3} & 0.7333 & -0.1333 \\\ -\frac{1}{3} & -0.1333 & 0.9333 \end{bmatrix}$$

With $$H_1$$, we then find $$R$$ by $$H_1 A$$:

$$H_1 A = \begin{bmatrix} -\frac{2}{3} & -\frac{2}{3} & -\frac{1}{3} \\\ -\frac{2}{3} & 0.7333 & -0.1333 \\\ -\frac{1}{3} & -0.1333 & 0.9333 \end{bmatrix} \begin{bmatrix} 2 & – 2 & 18 \\\ 2 & 1 & 0 \\\ 1 & 2 & 0 \end{bmatrix}$$
$$H_1 A = \begin{bmatrix} -3 & 0 & -12 \\\ 0 & 1.8 & 12 \\\ 0 & 2.4 & -6 \end{bmatrix}$$

$$H_1A$$ replaces $$A$$ in the next iteration. Now we move to $$a_2$$ and $$H_2$$. To this end we only consider the submatrix of $$A$$:

$$A^{(1)} = \begin{bmatrix} 1.8 & 12 \\\ 2.4 & -6 \end{bmatrix}$$

Thus, $$v_2$$ is equal to:

$$\begin{bmatrix} 1.8 \\\ 2.4 \end{bmatrix} + \sqrt{\sum^m_{j=1} a_2^2} \begin{bmatrix} 1 \\\ 0 \end{bmatrix} = \begin{bmatrix} 4.8 \\\ 2.4 \end{bmatrix}$$

Therefore, the corresponding Householder matrix $$H_2$$ is equal to:

$$H_2 = \begin{bmatrix} 1 & 0 \\\ 0 & 1 \end{bmatrix} – 2 \frac{\begin{bmatrix} 4.8 & 2.4 \end{bmatrix} \begin{bmatrix} 4.8 \\\ 2.4 \end{bmatrix}}{\begin{bmatrix} 4.8 \\\ 2.4 \end{bmatrix} \begin{bmatrix} 4.8 & 2.4 \end{bmatrix}}$$
$$H_2 = \begin{bmatrix} 1 & 0 & 0 \\\ 0 & -0.6 & -0.8 \\\ 0 & -0.8 & 0.6 \end{bmatrix}$$

The first column $$\begin{bmatrix}1 \\\ 0 \\\ 0 \end{bmatrix}$$ and first row $$\begin{bmatrix}1 & 0 & 0 \end{bmatrix}$$ are added to the resulting $$H_2$$ matrix to keep it $$n \times n$$.

Then we find $$H_2 A$$:

$$\begin{bmatrix} 1 & 0 & 0 \\\ 0 & -0.6 & -0.8 \\\ 0 & -0.8 & 0.6 \end{bmatrix} \begin{bmatrix} -3 & 0 & -12 \\\ 0 & 1.8 & 12 \\\ 0 & 2.4 & -6 \end{bmatrix}$$
$$H_2 A = \begin{bmatrix} -3 & 0 & -12 \\\ 0 & -3 & 12 \\\ 0 & 0 & 6 \end{bmatrix}$$

Moving to $$a_3$$ and $$H_3$$, the submatrix of $$H_2 A$$ is thus $$[6]$$. Therefore, $$v_3$$ is equal to:

$$\begin{bmatrix} 6 \end{bmatrix} – \sqrt{\sum^m_{j=1} a_3^2} \begin{bmatrix} 1 \end{bmatrix} = 12$$

The corresponding Householder matrix $$H_3$$ is then:

$$H_3 = \begin{bmatrix} 1 \end{bmatrix} – 2 \frac{\begin{bmatrix} 12 \end{bmatrix}\begin{bmatrix} 12 \end{bmatrix}}{\begin{bmatrix} 12 \end{bmatrix}\begin{bmatrix} 12 \end{bmatrix}} = \begin{bmatrix}1 & 0 & 0 \\\ 0 & 1 & 0 \\\ 0 & 0 & -1 \end{bmatrix}$$
$$H_3 A = \begin{bmatrix}1 & 0 & 0 \\\ 0 & 1 & 0 \\\ 0 & 0 & -1 \end{bmatrix} \begin{bmatrix} -3 & 0 & -12 \\\ 0 & -3 & 12 \\\ 0 & 0 & 6 \end{bmatrix}$$
$$H_3 A = R = \begin{bmatrix} -3 & 0 & -12 \\\ 0 & -3 & 12 \\\ 0 & 0 & -6 \end{bmatrix}$$

Which is the $$R$$ factorization in the QR decomposition method. The $$Q$$ factorization of QR decomposition is found by multiplying all the $$H$$ matrices together as mentioned earlier.

$$H_1 H_2 H_3 = Q$$
$$Q = \begin{bmatrix} -\frac{2}{3} & -\frac{2}{3} & -\frac{1}{3} \\\ -\frac{2}{3} & 0.7333 & -0.1333 \\\ -\frac{1}{3} & -0.1333 & 0.9333 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\\ 0 & -0.6 & -0.8 \\\ 0 & -0.8 & 0.6 \end{bmatrix} \begin{bmatrix}1 & 0 & 0 \\\ 0 & 1 & 0 \\\ 0 & 0 & -1 \end{bmatrix}$$
$$Q = \begin{bmatrix} -\frac{2}{3} & \frac{2}{3} & -\frac{1}{3} \\\ -\frac{2}{3} & -\frac{1}{3} & \frac{2}{3} \\\ -\frac{1}{3} & -\frac{2}{3} & -\frac{2}{3} \end{bmatrix}$$

Thus we obtain the same result as we did utilizing the Gram-Schmidt procedure.

$$QR = \begin{bmatrix} -\frac{2}{3} & \frac{2}{3} & -\frac{1}{3} \\\ -\frac{2}{3} & -\frac{1}{3} & \frac{2}{3} \\\ -\frac{1}{3} & -\frac{2}{3} & -\frac{2}{3} \end{bmatrix} \begin{bmatrix} -3 & 0 & -12 \\\ 0 & -3 & 12 \\\ 0 & 0 & -6 \end{bmatrix}$$
Householder Reflection QR Decomposition in R

The following function implements the Householder reflections approach to QR decomposition. The bdiag() function in the Matrix package is used in constructing the $$H$$ matrices as seen above in the calculation of $$H_2$$.

qr.householder <- function(A) {
require(Matrix)

R <- as.matrix(A) # Set R to the input matrix A

n <- ncol(A)
m <- nrow(A)
H <- list() # Initialize a list to store the computed H matrices to calculate Q later

if (m > n) {
c <- n
} else {
c <- m
}

for (k in 1:c) {
x <- R[k:m,k] # Equivalent to a_1
e <- as.matrix(c(1, rep(0, length(x)-1)))
vk <- sign(x[1]) * sqrt(sum(x^2)) * e + x

# Compute the H matrix
hk <- diag(length(x)) - 2 * as.vector(vk %*% t(vk)) / (t(vk) %*% vk)
if (k > 1) {
hk <- bdiag(diag(k-1), hk)
}

# Store the H matrix to find Q at the end of iteration
H[[k]] <- hk

R <- hk %*% R
}

Q <- Reduce("%*%", H) # Calculate Q matrix by multiplying all H matrices
res <- list('Q'=Q,'R'=R)
return(res)
}


Use the function to compute the QR decomposition of the following matrix $$A$$ with Householder reflections.

A <- rbind(c(2,-2,18),c(2,1,0),c(1,2,0))
A
##      [,1] [,2] [,3]
## [1,]    2   -2   18
## [2,]    2    1    0
## [3,]    1    2    0

qr.householder(A)
## $Q ## 3 x 3 Matrix of class "dgeMatrix" ## [,1] [,2] [,3] ## [1,] -0.6666667 0.6666667 -0.3333333 ## [2,] -0.6666667 -0.3333333 0.6666667 ## [3,] -0.3333333 -0.6666667 -0.6666667 ## ##$R
## 3 x 3 Matrix of class "dgeMatrix"
##               [,1]          [,2] [,3]
## [1,] -3.000000e+00  2.220446e-16  -12
## [2,] -1.165734e-16 -3.000000e+00   12
## [3,]  1.554312e-16  0.000000e+00   -6


The only package that I found that directly implements the Householder reflection approach to QR is the pracma package.

library(pracma)

house <- householder(A)
house
## $Q ## [,1] [,2] [,3] ## [1,] -0.6666667 0.6666667 0.3333333 ## [2,] -0.6666667 -0.3333333 -0.6666667 ## [3,] -0.3333333 -0.6666667 0.6666667 ## ##$R
##               [,1]          [,2] [,3]
## [1,] -3.000000e+00  2.220446e-16  -12
## [2,] -1.165734e-16 -3.000000e+00   12
## [3,] -1.554312e-16  0.000000e+00    6

References

En.wikipedia.org. (2017). QR decomposition. [online] Available at: https://en.wikipedia.org/wiki/QR_decomposition#Using_Householder_reflections [Accessed 10 Apr. 2017].

http://www.math.iit.edu/~fass/477577_Chapter_4.pdf

Trefethen, L. and Bau, D. (1997). Numerical linear algebra. 1st ed. Philadelphia: SIAM.

The post QR Decomposition with Householder Reflections appeared first on Aaron Schlegel.

To leave a comment for the author, please follow the link and comment on their blog: R – Aaron Schlegel.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update! Subscribe to R-bloggers to receive e-mails with the latest R posts.(You will not see this message again.)

Click here to close (This popup will not appear again)