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.

The post QR Decomposition with Householder Reflections appeared first on Aaron Schlegel.