<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/qr-decomposition-with-householder-reflections/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.aaronschlegel.com/qr-decomposition-householder-reflections/"> R – Aaron Schlegel</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>The more common approach to QR decomposition is employing Householder reflections rather than utilizing Gram-Schmidt. In practice, the <a target="_blank" rel="nofollow" href="http://wp.me/p4aZEo-5R0">Gram-Schmidt procedure</a> is not recommended as it can lead to cancellation that causes inaccuracy of the computation of <span data-display="false" class="wp-katex-eq">q_j</span>, which may result in a non-orthogonal <span data-display="false" class="wp-katex-eq">Q</span> matrix. Householder reflections are another method of orthogonal transformation that transforms a vector <span data-display="false" class="wp-katex-eq">x</span> into a unit vector <span data-display="false" class="wp-katex-eq">y</span> parallel with <span data-display="false" class="wp-katex-eq">x</span>. The Householder reflection matrix with normal vector <span data-display="false" class="wp-katex-eq">v</span> takes the form:</p>
<p><span data-display="true" class="wp-katex-eq"> H = I – 2vv^T </span></p>
<p>Thus we need to build <span data-display="false" class="wp-katex-eq">H</span> so that <span data-display="false" class="wp-katex-eq">Hx = \alpha e_1</span> for some constant <span data-display="false" class="wp-katex-eq">\alpha</span> and <span data-display="false" class="wp-katex-eq">e_1 = \begin{bmatrix} 1 &amp; 0 &amp; 0 \end{bmatrix}^T</span>.</p>
<p>Since <span data-display="false" class="wp-katex-eq">H</span> is orthogonal, <span data-display="false" class="wp-katex-eq">||Hx|| = ||x||</span> and <span data-display="false" class="wp-katex-eq">||\alpha e_1|| = |\alpha| ||e_1|| = |\alpha|</span>. Therefore, <span data-display="false" class="wp-katex-eq">\alpha = \pm ||x||</span>. The sign is selected so it has the opposite sign of <span data-display="false" class="wp-katex-eq">x_1</span> (we’ll use <span data-display="false" class="wp-katex-eq">+</span> for the remaining definitions). The vector <span data-display="false" class="wp-katex-eq">u</span> we seek is thus:</p>
<p><span data-display="true" class="wp-katex-eq"> u = \begin{bmatrix} x_1 + sign(x_1) ||x_1|| \\\ x_2 \\\ \vdots \\\ x_n \end{bmatrix} </span></p>
<p>With the unit vector <span data-display="false" class="wp-katex-eq">v</span> defined as <span data-display="false" class="wp-katex-eq">u = \frac{v}{||v||}</span>. The corresponding Householder reflection is then:</p>
<p><span data-display="true" class="wp-katex-eq"> H(x) = I – 2vv^T = I – 2 \frac{uu^T}{u^Tu} </span></p>
<h6>QR Decomposition with Householder Reflections</h6>
<p>The Householder reflection method of QR decomposition works by finding appropriate <span data-display="false" class="wp-katex-eq">H</span> matrices and multiplying them from the left by the original matrix <span data-display="false" class="wp-katex-eq">A</span> to construct the upper triangular matrix <span data-display="false" class="wp-katex-eq">R</span>. As we saw earlier, unlike the Gram-Schmidt procedure, the Householder reflection approach does not explicitly form the <span data-display="false" class="wp-katex-eq">Q</span> matrix. However, the <span data-display="false" class="wp-katex-eq">Q</span> matrix can be found by taking the dot product of each successively formed Householder matrix.</p>
<p><span data-display="true" class="wp-katex-eq"> Q = H_1 H_2 \cdots H_{m-2}H_{m-1}</span></p>
<p>Consider the matrix <span data-display="false" class="wp-katex-eq">A</span>.</p>
<p><span data-display="true" class="wp-katex-eq">A = \begin{bmatrix} 2 &amp; – 2 &amp; 18 \\\ 2 &amp; 1 &amp; 0 \\\ 1 &amp; 2 &amp; 0 \end{bmatrix}</span></p>
<p>We find the reflection of the first column vector <span data-display="false" class="wp-katex-eq">a_1 = \begin{bmatrix}2 &amp; 2 &amp; 1 \end{bmatrix}^T</span>, <span data-display="false" class="wp-katex-eq">v_1 = a_1 + sign(a_{11})||a_1||e_1</span>.</p>
<p><span data-display="true" class="wp-katex-eq"> 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} </span></p>
<p>With a corresponding Householder matrix:</p>
<p><span data-display="true" class="wp-katex-eq"> H_1 = \begin{bmatrix} 1 &amp; 0 &amp; 0 \\\ 0 &amp; 1 &amp; 0 \\\ 0 &amp; 0 &amp; 1 \end{bmatrix} – 2 \frac{\begin{bmatrix} 2 &amp; 2 &amp; 1 \end{bmatrix}\begin{bmatrix} 2 \\\ 2 \\\ 1 \end{bmatrix}}{\begin{bmatrix} 2 \\\ 2 \\\ 1 \end{bmatrix}\begin{bmatrix} 2 &amp; 2 &amp; 1 \end{bmatrix}} = </span><br>
<span data-display="true" class="wp-katex-eq"> H_1 = \begin{bmatrix} -\frac{2}{3} &amp; -\frac{2}{3} &amp; -\frac{1}{3} \\\ -\frac{2}{3} &amp; 0.7333 &amp; -0.1333 \\\ -\frac{1}{3} &amp; -0.1333 &amp; 0.9333 \end{bmatrix} </span></p>
<p>With <span data-display="false" class="wp-katex-eq">H_1</span>, we then find <span data-display="false" class="wp-katex-eq">R</span> by <span data-display="false" class="wp-katex-eq">H_1 A</span>:</p>
<p><span data-display="true" class="wp-katex-eq"> H_1 A = \begin{bmatrix} -\frac{2}{3} &amp; -\frac{2}{3} &amp; -\frac{1}{3} \\\ -\frac{2}{3} &amp; 0.7333 &amp; -0.1333 \\\ -\frac{1}{3} &amp; -0.1333 &amp; 0.9333 \end{bmatrix} \begin{bmatrix} 2 &amp; – 2 &amp; 18 \\\ 2 &amp; 1 &amp; 0 \\\ 1 &amp; 2 &amp; 0 \end{bmatrix} </span><br>
<span data-display="true" class="wp-katex-eq"> H_1 A = \begin{bmatrix} -3 &amp; 0 &amp; -12 \\\ 0 &amp; 1.8 &amp; 12 \\\ 0 &amp; 2.4 &amp; -6 \end{bmatrix} </span></p>
<p><span data-display="false" class="wp-katex-eq">H_1A</span> replaces <span data-display="false" class="wp-katex-eq">A</span> in the next iteration. Now we move to <span data-display="false" class="wp-katex-eq">a_2</span> and <span data-display="false" class="wp-katex-eq">H_2</span>. To this end we only consider the submatrix of <span data-display="false" class="wp-katex-eq">A</span>:</p>
<p><span data-display="true" class="wp-katex-eq"> A^{(1)} = \begin{bmatrix} 1.8 &amp; 12 \\\ 2.4 &amp; -6 \end{bmatrix} </span></p>
<p>Thus, <span data-display="false" class="wp-katex-eq">v_2</span> is equal to:</p>
<p><span data-display="true" class="wp-katex-eq"> \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}</span></p>
<p>Therefore, the corresponding Householder matrix <span data-display="false" class="wp-katex-eq">H_2</span> is equal to:</p>
<p><span data-display="true" class="wp-katex-eq"> H_2 = \begin{bmatrix} 1 &amp; 0 \\\ 0 &amp; 1 \end{bmatrix} – 2 \frac{\begin{bmatrix} 4.8 &amp; 2.4 \end{bmatrix} \begin{bmatrix} 4.8 \\\ 2.4 \end{bmatrix}}{\begin{bmatrix} 4.8 \\\ 2.4 \end{bmatrix} \begin{bmatrix} 4.8 &amp; 2.4 \end{bmatrix}} </span><br>
<span data-display="true" class="wp-katex-eq"> H_2 = \begin{bmatrix} 1 &amp; 0 &amp; 0 \\\ 0 &amp; -0.6 &amp; -0.8 \\\ 0 &amp; -0.8 &amp; 0.6 \end{bmatrix} </span></p>
<p>The first column <span data-display="false" class="wp-katex-eq">\begin{bmatrix}1 \\\ 0 \\\ 0 \end{bmatrix}</span> and first row <span data-display="false" class="wp-katex-eq">\begin{bmatrix}1 &amp; 0 &amp; 0 \end{bmatrix}</span> are added to the resulting <span data-display="false" class="wp-katex-eq">H_2</span> matrix to keep it <span data-display="false" class="wp-katex-eq">n \times n</span>.</p>
<p>Then we find <span data-display="false" class="wp-katex-eq">H_2 A</span>:</p>
<p><span data-display="true" class="wp-katex-eq"> \begin{bmatrix} 1 &amp; 0 &amp; 0 \\\ 0 &amp; -0.6 &amp; -0.8 \\\ 0 &amp; -0.8 &amp; 0.6 \end{bmatrix} \begin{bmatrix} -3 &amp; 0 &amp; -12 \\\ 0 &amp; 1.8 &amp; 12 \\\ 0 &amp; 2.4 &amp; -6 \end{bmatrix} </span><br>
<span data-display="true" class="wp-katex-eq"> H_2 A = \begin{bmatrix} -3 &amp; 0 &amp; -12 \\\ 0 &amp; -3 &amp; 12 \\\ 0 &amp; 0 &amp; 6 \end{bmatrix}</span></p>
<p>Moving to <span data-display="false" class="wp-katex-eq">a_3</span> and <span data-display="false" class="wp-katex-eq">H_3</span>, the submatrix of <span data-display="false" class="wp-katex-eq">H_2 A</span> is thus <span data-display="false" class="wp-katex-eq">[6]</span>. Therefore, <span data-display="false" class="wp-katex-eq">v_3</span> is equal to:</p>
<p><span data-display="true" class="wp-katex-eq"> \begin{bmatrix} 6 \end{bmatrix} – \sqrt{\sum^m_{j=1} a_3^2} \begin{bmatrix} 1 \end{bmatrix} = 12 </span></p>
<p>The corresponding Householder matrix <span data-display="false" class="wp-katex-eq">H_3</span> is then:</p>
<p><span data-display="true" class="wp-katex-eq"> 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 &amp; 0 &amp; 0 \\\ 0 &amp; 1 &amp; 0 \\\ 0 &amp; 0 &amp; -1 \end{bmatrix}</span><br>
<span data-display="true" class="wp-katex-eq">H_3 A = \begin{bmatrix}1 &amp; 0 &amp; 0 \\\ 0 &amp; 1 &amp; 0 \\\ 0 &amp; 0 &amp; -1 \end{bmatrix} \begin{bmatrix} -3 &amp; 0 &amp; -12 \\\ 0 &amp; -3 &amp; 12 \\\ 0 &amp; 0 &amp; 6 \end{bmatrix} </span><br>
<span data-display="true" class="wp-katex-eq"> H_3 A = R = \begin{bmatrix} -3 &amp; 0 &amp; -12 \\\ 0 &amp; -3 &amp; 12 \\\ 0 &amp; 0 &amp; -6 \end{bmatrix} </span></p>
<p>Which is the <span data-display="false" class="wp-katex-eq">R</span> factorization in the QR decomposition method. The <span data-display="false" class="wp-katex-eq">Q</span> factorization of QR decomposition is found by multiplying all the <span data-display="false" class="wp-katex-eq">H</span> matrices together as mentioned earlier.</p>
<p><span data-display="true" class="wp-katex-eq"> H_1 H_2 H_3 = Q </span><br>
<span data-display="true" class="wp-katex-eq"> Q = \begin{bmatrix} -\frac{2}{3} &amp; -\frac{2}{3} &amp; -\frac{1}{3} \\\ -\frac{2}{3} &amp; 0.7333 &amp; -0.1333 \\\ -\frac{1}{3} &amp; -0.1333 &amp; 0.9333 \end{bmatrix} \begin{bmatrix} 1 &amp; 0 &amp; 0 \\\ 0 &amp; -0.6 &amp; -0.8 \\\ 0 &amp; -0.8 &amp; 0.6 \end{bmatrix} \begin{bmatrix}1 &amp; 0 &amp; 0 \\\ 0 &amp; 1 &amp; 0 \\\ 0 &amp; 0 &amp; -1 \end{bmatrix} </span><br>
<span data-display="true" class="wp-katex-eq"> Q = \begin{bmatrix} -\frac{2}{3} &amp; \frac{2}{3} &amp; -\frac{1}{3} \\\ -\frac{2}{3} &amp; -\frac{1}{3} &amp; \frac{2}{3} \\\ -\frac{1}{3} &amp; -\frac{2}{3} &amp; -\frac{2}{3} \end{bmatrix} </span></p>
<p>Thus we obtain the same result as we did utilizing the Gram-Schmidt procedure.</p>
<p><span data-display="true" class="wp-katex-eq"> QR = \begin{bmatrix} -\frac{2}{3} &amp; \frac{2}{3} &amp; -\frac{1}{3} \\\ -\frac{2}{3} &amp; -\frac{1}{3} &amp; \frac{2}{3} \\\ -\frac{1}{3} &amp; -\frac{2}{3} &amp; -\frac{2}{3} \end{bmatrix} \begin{bmatrix} -3 &amp; 0 &amp; -12 \\\ 0 &amp; -3 &amp; 12 \\\ 0 &amp; 0 &amp; -6 \end{bmatrix} </span></p>
<h6>Householder Reflection QR Decomposition in R</h6>
<p>The following function implements the Householder reflections approach to QR decomposition. The <code>bdiag()</code> function in the <a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/Matrix/index.html">Matrix package</a> is used in constructing the <span data-display="false" class="wp-katex-eq">H</span> matrices as seen above in the calculation of <span data-display="false" class="wp-katex-eq">H_2</span>.</p>
<pre>
qr.householder &lt;- function(A) {
  require(Matrix)
  
  R &lt;- as.matrix(A) # Set R to the input matrix A
  
  n &lt;- ncol(A)
  m &lt;- nrow(A)
  H &lt;- list() # Initialize a list to store the computed H matrices to calculate Q later
  
  if (m > n) {
    c &lt;- n
  } else {
    c &lt;- m
  }
  
  for (k in 1:c) {
    x &lt;- R[k:m,k] # Equivalent to a_1
    e &lt;- as.matrix(c(1, rep(0, length(x)-1)))
    vk &lt;- sign(x[1]) * sqrt(sum(x^2)) * e + x
    
    # Compute the H matrix
    hk &lt;- diag(length(x)) - 2 * as.vector(vk %*% t(vk)) / (t(vk) %*% vk)
    if (k > 1) {
      hk &lt;- bdiag(diag(k-1), hk)
    }
    
    # Store the H matrix to find Q at the end of iteration
    H[[k]] &lt;- hk
    
    R &lt;- hk %*% R
  }

  Q &lt;- Reduce(&quot;%*%&quot;, H) # Calculate Q matrix by multiplying all H matrices
  res &lt;- list('Q'=Q,'R'=R)
  return(res)
}
</pre>
<p>Use the function to compute the QR decomposition of the following matrix <span data-display="false" class="wp-katex-eq">A</span> with Householder reflections.</p>
<pre>
A &lt;- 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
</pre>
<pre>
qr.householder(A)
## $Q
## 3 x 3 Matrix of class &quot;dgeMatrix&quot;
##            [,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 &quot;dgeMatrix&quot;
##               [,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
</pre>
<p>The only package that I found that directly implements the Householder reflection approach to QR is the <a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/pracma/index.html">pracma package</a>.</p>
<pre>
library(pracma)
</pre>
<pre>
house &lt;- 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
</pre>
<h6>References</h6>
<p>En.wikipedia.org. (2017). QR decomposition. [online] Available at: <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/QR_decomposition#Using_Householder_reflections">https://en.wikipedia.org/wiki/QR_decomposition#Using_Householder_reflections</a> [Accessed 10 Apr. 2017].</p>
<p><a target="_blank" rel="nofollow" href="http://www.math.iit.edu/~fass/477577_Chapter_4.pdf">http://www.math.iit.edu/~fass/477577_Chapter_4.pdf</a></p>
<p>Trefethen, L. and Bau, D. (1997). Numerical linear algebra. 1st ed. Philadelphia: SIAM.</p>
<p>The post <a target="_blank" rel="nofollow" href="http://www.aaronschlegel.com/qr-decomposition-householder-reflections/">QR Decomposition with Householder Reflections</a> appeared first on <a target="_blank" rel="nofollow" href="http://www.aaronschlegel.com/">Aaron Schlegel</a>.</p>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://www.aaronschlegel.com/qr-decomposition-householder-reflections/"> R – Aaron Schlegel</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/fNB2g7EtKQc"></a></div>