QR Decomposition with Householder Reflections
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 qj, 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–2vvTThus we need to build H so that Hx=αe1 for some constant α and e1=[100]T.
Since H is orthogonal, ||Hx||=||x|| and ||αe1||=|α|||e1||=|α|. Therefore, α=±||x||. The sign is selected so it has the opposite sign of x1 (we’ll use + for the remaining definitions). The vector u we seek is thus:
u=[x1+sign(x1)||x1|| x2 ⋮ xn]With the unit vector v defined as u=v||v||. The corresponding Householder reflection is then:
H(x)=I–2vvT=I–2uuTuTuQR 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=H1H2⋯Hm−2Hm−1Consider the matrix A.
A=[2–218 210 120]We find the reflection of the first column vector a1=[221]T, v1=a1+sign(a11)||a1||e1.
v1=[2 2 1]+√∑mk=1a12[1 0 0]=[5 2 1]With a corresponding Householder matrix:
H1=[100 010 001]–2[221][2 2 1][2 2 1][221]=
H1=[−23−23−13 −230.7333−0.1333 −13−0.13330.9333]
With H1, we then find R by H1A:
H1A=[−23−23−13 −230.7333−0.1333 −13−0.13330.9333][2–218 210 120]
H1A=[−30−12 01.812 02.4−6]
H1A replaces A in the next iteration. Now we move to a2 and H2. To this end we only consider the submatrix of A:
A(1)=[1.812 2.4−6]Thus, v2 is equal to:
[1.8 2.4]+√∑mj=1a22[1 0]=[4.8 2.4]Therefore, the corresponding Householder matrix H2 is equal to:
H2=[10 01]–2[4.82.4][4.8 2.4][4.8 2.4][4.82.4]
H2=[100 0−0.6−0.8 0−0.80.6]
The first column [1 0 0] and first row [100] are added to the resulting H2 matrix to keep it n×n.
Then we find H2A:
[100 0−0.6−0.8 0−0.80.6][−30−12 01.812 02.4−6]
H2A=[−30−12 0−312 006]
Moving to a3 and H3, the submatrix of H2A is thus [6]. Therefore, v3 is equal to:
[6]–√∑mj=1a23[1]=12The corresponding Householder matrix H3 is then:
H3=[1]–2[12][12][12][12]=[100 010 00−1]
H3A=[100 010 00−1][−30−12 0−312 006]
H3A=R=[−30−12 0−312 00−6]
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.
H1H2H3=Q
Q=[−23−23−13 −230.7333−0.1333 −13−0.13330.9333][100 0−0.6−0.8 0−0.80.6][100 010 00−1]
Q=[−2323−13 −23−1323 −13−23−23]
Thus we obtain the same result as we did utilizing the Gram-Schmidt procedure.
QR=[−2323−13 −23−1323 −13−23−23][−30−12 0−312 00−6]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 H2.
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.
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.