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 [latex]q_j[/latex], which may result in a non-orthogonal [latex]Q[/latex] matrix. Householder reflections are another method of orthogonal transformation that transforms a vector [latex]x[/latex] into a unit vector [latex]y[/latex] parallel with [latex]x[/latex]. The Householder reflection matrix with normal vector [latex]v[/latex] takes the form:

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

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

Since [latex]H[/latex] is orthogonal, [latex]||Hx|| = ||x||[/latex] and [latex]||\alpha e_1|| = |\alpha| ||e_1|| = |\alpha|[/latex]. Therefore, [latex]\alpha = \pm ||x||[/latex]. The sign is selected so it has the opposite sign of [latex]x_1[/latex] (we’ll use [latex]+[/latex] for the remaining definitions). The vector [latex]u[/latex] 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 [latex]v[/latex] defined as [latex]u = \frac{v}{||v||}[/latex]. 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 [latex]H[/latex] matrices and multiplying them from the left by the original matrix [latex]A[/latex] to construct the upper triangular matrix [latex]R[/latex]. As we saw earlier, unlike the Gram-Schmidt procedure, the Householder reflection approach does not explicitly form the [latex]Q[/latex] matrix. However, the [latex]Q[/latex] 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 [latex]A[/latex].

[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 [latex]a_1 = \begin{bmatrix}2 & 2 & 1 \end{bmatrix}^T[/latex], [latex]v_1 = a_1 + sign(a_{11})||a_1||e_1[/latex].

[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 [latex]H_1[/latex], we then find [latex]R[/latex] by [latex]H_1 A[/latex]:

[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]

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

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

Thus, [latex]v_2[/latex] 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 [latex]H_2[/latex] 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 [latex]\begin{bmatrix}1 \\\ 0 \\\ 0 \end{bmatrix}[/latex] and first row [latex]\begin{bmatrix}1 & 0 & 0 \end{bmatrix}[/latex] are added to the resulting [latex]H_2[/latex] matrix to keep it [latex]n \times n[/latex].

Then we find [latex]H_2 A[/latex]:

[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 [latex]a_3[/latex] and [latex]H_3[/latex], the submatrix of [latex]H_2 A[/latex] is thus [latex][6][/latex]. Therefore, [latex]v_3[/latex] 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 [latex]H_3[/latex] 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 [latex]R[/latex] factorization in the QR decomposition method. The [latex]Q[/latex] factorization of QR decomposition is found by multiplying all the [latex]H[/latex] 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 [latex]H[/latex] matrices as seen above in the calculation of [latex]H_2[/latex].

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 [latex]A[/latex] 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.

2 Comments
  • Drew Schmidt

    April 13, 2017 at 6:38 pm Reply

    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

  • […] article was first published on R – Aaron Schlegel, and kindly contributed to […]

Post a Comment