## QR Decomposition with Householder Reflections

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.

Base R’s qr() function uses Householder reflectors. By default it uses a modified version of the old LINPACK subroutine DQRDC (which they call DQRDC2), but you can optionally use the modern LAPACK function DGEQP3. Neither of these uses Gram-Schmidt. Both are rank-revealing however, as they use column pivoting. For an approach that avoids column pivoting, see https://github.com/wrathematics/lqr