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:

[latex display=”true”] H = I – 2vv^T [/latex]

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:

[latex display=”true”] u = \begin{bmatrix} x_1 + sign(x_1) ||x_1|| \\\ x_2 \\\ \vdots \\\ x_n \end{bmatrix} [/latex]

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

[latex display=”true”] H(x) = I – 2vv^T = I – 2 \frac{uu^T}{u^Tu} [/latex]

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.

[latex display=”true”] Q = H_1 H_2 \cdots H_{m-2}H_{m-1}[/latex]

Consider the matrix $A$.

[latex display=”true”]A = \begin{bmatrix} 2 & – 2 & 18 \\\ 2 & 1 & 0 \\\ 1 & 2 & 0 \end{bmatrix}[/latex]

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$.

[latex display=”true”] 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} [/latex]

With a corresponding Householder matrix:

[latex display=”true”] 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}} = [/latex]
[latex display=”true”] 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} [/latex]

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

[latex display=”true”] 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} [/latex]
[latex display=”true”] H_1 A = \begin{bmatrix} -3 & 0 & -12 \\\ 0 & 1.8 & 12 \\\ 0 & 2.4 & -6 \end{bmatrix} [/latex]

$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$:

[latex display=”true”] A^{(1)} = \begin{bmatrix} 1.8 & 12 \\\ 2.4 & -6 \end{bmatrix} [/latex]

Thus, $v_2$ is equal to:

[latex display=”true”] \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}[/latex]

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

[latex display=”true”] 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}} [/latex]
[latex display=”true”] H_2 = \begin{bmatrix} 1 & 0 & 0 \\\ 0 & -0.6 & -0.8 \\\ 0 & -0.8 & 0.6 \end{bmatrix} [/latex]

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$:

[latex display=”true”] \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} [/latex]
[latex display=”true”] H_2 A = \begin{bmatrix} -3 & 0 & -12 \\\ 0 & -3 & 12 \\\ 0 & 0 & 6 \end{bmatrix}[/latex]

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

[latex display=”true”] \begin{bmatrix} 6 \end{bmatrix} – \sqrt{\sum^m_{j=1} a_3^2} \begin{bmatrix} 1 \end{bmatrix} = 12 [/latex]

The corresponding Householder matrix $H_3$ is then:

[latex display=”true”] 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}[/latex]
[latex display=”true”]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} [/latex]
[latex display=”true”] H_3 A = R = \begin{bmatrix} -3 & 0 & -12 \\\ 0 & -3 & 12 \\\ 0 & 0 & -6 \end{bmatrix} [/latex]

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.

[latex display=”true”] H_1 H_2 H_3 = Q [/latex]
[latex display=”true”] 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} [/latex]
[latex display=”true”] 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} [/latex]

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

[latex display=”true”] 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} [/latex]

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.