QR Decomposition with Householder Reflections

[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 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=I2vvT

Thus 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)=I2vvT=I2uuTuTu
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=H1H2Hm2Hm1

Consider the matrix A.

A=[2218 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=[232313 230.73330.1333 130.13330.9333]

With H1, we then find R by H1A:

H1A=[232313 230.73330.1333 130.13330.9333][2218 210 120]
H1A=[3012 01.812 02.46]

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.46]

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 00.60.8 00.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 00.60.8 00.80.6][3012 01.812 02.46]
H2A=[3012 0312 006]

Moving to a3 and H3, the submatrix of H2A is thus [6]. Therefore, v3 is equal to:

[6]mj=1a23[1]=12

The corresponding Householder matrix H3 is then:

H3=[1]2[12][12][12][12]=[100 010 001]
H3A=[100 010 001][3012 0312 006]
H3A=R=[3012 0312 006]

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=[232313 230.73330.1333 130.13330.9333][100 00.60.8 00.80.6][100 010 001]
Q=[232313 231323 132323]

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

QR=[232313 231323 132323][3012 0312 006]
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.

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)