<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-the-gram-schmidt-algorithm/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.aaronschlegel.com/qr-decomposition-gram-schmidt-algorithm/"> R – Aaron Schlegel</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>QR decomposition is another technique for decomposing a matrix into a form that is easier to work with in further applications. The QR decomposition technique decomposes a square or rectangular matrix, which we will denote as <span data-display="false" class="wp-katex-eq">A</span>, into two components, <span data-display="false" class="wp-katex-eq">Q</span>, and <span data-display="false" class="wp-katex-eq">R</span>.<br>
<span data-display="true" class="wp-katex-eq"> A = QR </span></p>
<p>Where <span data-display="false" class="wp-katex-eq">Q</span> is an orthogonal matrix, and <span data-display="false" class="wp-katex-eq">R</span> is an upper triangular matrix. Recall an orthogonal matrix is a square matrix with orthonormal row and column vectors such that <span data-display="false" class="wp-katex-eq">Q^T Q = I</span>, where <span data-display="false" class="wp-katex-eq">I</span> is the identity matrix. The term orthonormal implies the vectors are of unit length and are perpendicular (orthogonal) to each other.</p>
<p>QR decomposition is often used in linear least squares estimation and is, in fact, the method used by R in its <a target="_blank" rel="nofollow" href="https://stat.ethz.ch/R-manual/R-devel/library/stats/html/lm.html"><code>lm()</code> function</a>. <a target="_blank" rel="nofollow" href="http://www.auburn.edu/~tamtiny/Xiamen%20U%202010%20Dec.pdf">Signal processing and MIMO systems</a> also employ QR decomposition. There are several methods for performing QR decomposition, including the Gram-Schmidt process, Householder reflections, and Givens rotations. This post is concerned with the Gram-Schmidt process.</p>
<h6>The Gram-Schmidt Process</h6>
<p>The Gram-Schmidt process is used to find an orthogonal basis from a non-orthogonal basis. An <a target="_blank" rel="nofollow" href="https://www.math.ucdavis.edu/~linear/old/notes21.pdf">orthogonal basis</a> has many properties that are desirable for further computations and expansions. As noted previously, an orthogonal matrix has row and column vectors of unit length:</p>
<p><span data-display="true" class="wp-katex-eq"> ||a_n|| = \sqrt{a_n \cdot a_n} = \sqrt{a_n^T a_n} = 1 </span></p>
<p>Where <span data-display="false" class="wp-katex-eq">a_n</span> is a linearly independent column vector of a matrix. The vectors are also perpendicular in an orthogonal basis. The Gram-Schmidt process works by finding an orthogonal projection <span data-display="false" class="wp-katex-eq">q_n</span> for each column vector <span data-display="false" class="wp-katex-eq">a_n</span> and then subtracting its projections onto the previous projections <span data-display="false" class="wp-katex-eq">(q_j)</span>. The resulting vector is then divided by the length of that vector to produce a unit vector.</p>
<p>Consider a matrix <span data-display="false" class="wp-katex-eq">A</span> with <span data-display="false" class="wp-katex-eq">n</span> column vectors such that:</p>
<p><span data-display="true" class="wp-katex-eq"> A = \left[ a_1 | a_2 | \cdots | a_n \right] </span></p>
<p>The Gram-Schmidt process proceeds by finding the orthogonal projection of the first column vector <span data-display="false" class="wp-katex-eq">a_1</span>.</p>
<p><span data-display="true" class="wp-katex-eq"> v_1 = a_1, \qquad e_1 = \frac{v_1}{||v_1||} </span></p>
<p>Because <span data-display="false" class="wp-katex-eq">a_1</span> is the first column vector, there is no preceeding projections to subtract. The second column <span data-display="false" class="wp-katex-eq">a_2</span> is subtracted by the previous projection on the column vector:</p>
<p><span data-display="true" class="wp-katex-eq"> v_2 = a_2 – proj_{v_1} (a_2) = a_2 – (a_2 \cdot e_1) e_1, \qquad e_2 = \frac{v_2}{||v_2||} </span></p>
<p>This process continues up to the <span data-display="false" class="wp-katex-eq">n</span> column vectors, where each incremental step <span data-display="false" class="wp-katex-eq">k + 1</span> is computed as:</p>
<p><span data-display="true" class="wp-katex-eq"> v_{k+1} = a_{k+1} – (a_{k+1} \cdot e_{1}) e_1 – \cdots – (a_{k+1} \cdot e_k) e_k, \qquad e_{k+1} = \frac{u_{k+1}}{||u_{k+1}||} </span></p>
<p>The <span data-display="false" class="wp-katex-eq">|| \cdot ||</span> is the <span data-display="false" class="wp-katex-eq">L_2</span> norm which is defined as:</p>
<p><span data-display="true" class="wp-katex-eq"> \sqrt{\sum^m_{j=1} v_k^2} </span> The projection can also be defined by:</p>
<p>Thus the matrix <span data-display="false" class="wp-katex-eq">A</span> can be factorized into the <span data-display="false" class="wp-katex-eq">QR</span> matrix as the following:</p>
<p><span data-display="true" class="wp-katex-eq"> A = \left[a_1 | a_2 | \cdots | a_n \right] = \left[e_1 | e_2 | \cdots | e_n \right] \begin{bmatrix}a_1 \cdot e_1 &amp; a_2 \cdot e_1 &amp; \cdots &amp; a_n \cdot e_1 \\\ 0 &amp; a_2 \cdot e_2 &amp; \cdots &amp; a_n \cdot e_2 \\\ \vdots &amp; \vdots &amp; &amp; \vdots \\\ 0 &amp; 0 &amp; \cdots &amp; a_n \cdot e_n\end{bmatrix} = QR</span></p>
<h6>Gram-Schmidt Process Example</h6>
<p>Consider the matrix <span data-display="false" class="wp-katex-eq">A</span>:</p>
<p><span data-display="true" class="wp-katex-eq">\begin{bmatrix} 2 &amp; – 2 &amp; 18 \\\ 2 &amp; 1 &amp; 0 \\\ 1 &amp; 2 &amp; 0 \end{bmatrix}</span></p>
<p>We would like to orthogonalize this matrix using the Gram-Schmidt process. The resulting orthogonalized vector is also equivalent to <span data-display="false" class="wp-katex-eq">Q</span> in the <span data-display="false" class="wp-katex-eq">QR</span> decomposition.</p>
<p>The Gram-Schmidt process on the matrix <span data-display="false" class="wp-katex-eq">A</span> proceeds as follows:</p>
<p><span data-display="true" class="wp-katex-eq"> v_1 = a_1 = \begin{bmatrix}2 \\\ 2 \\\ 1\end{bmatrix} \qquad e_1 = \frac{v_1}{||v_1||} = \frac{\begin{bmatrix}2 \\\ 2 \\\ 1\end{bmatrix}}{\sqrt{\sum{\begin{bmatrix}2 \\\ 2 \\\ 1\end{bmatrix}^2}}}</span> <span data-display="true" class="wp-katex-eq"> e_1 = \begin{bmatrix} \frac{2}{3} \\\ \frac{2}{3} \\\ \frac{1}{3} \end{bmatrix}</span><br>
<span data-display="true" class="wp-katex-eq"> v_2 = a_2 – (a_2 \cdot e_1) e_1 = \begin{bmatrix}-2 \\\ 1 \\\ 2\end{bmatrix} – \left(\begin{bmatrix}-2 \\\ 1 \\\ 2\end{bmatrix}, \begin{bmatrix} \frac{2}{3} \\\ \frac{2}{3} \\\ \frac{1}{3} \end{bmatrix}\right)\begin{bmatrix} \frac{2}{3} \\\ \frac{2}{3} \\\ \frac{1}{3} \end{bmatrix} </span><br>
<span data-display="true" class="wp-katex-eq"> v_2 = \begin{bmatrix}-2 \\\ 1 \\\ 2\end{bmatrix} \qquad e_2 = \frac{v_2}{||v_2||} = \frac{\begin{bmatrix}-2 \\\ 1 \\\ 2\end{bmatrix}}{\sqrt{\sum{\begin{bmatrix}-2 \\\ 1 \\\ 2\end{bmatrix}^2}}}</span><br>
<span data-display="true" class="wp-katex-eq"> e_2 = \begin{bmatrix} -\frac{2}{3} \\\ \frac{1}{3} \\\ \frac{2}{3} \end{bmatrix}</span><br>
<span data-display="true" class="wp-katex-eq"> v_3 = a_3 – (a_3 \cdot e_1) e_1 – (a_3 \cdot e_2) e_2 </span><br>
<span data-display="true" class="wp-katex-eq"> v_3 = \begin{bmatrix}18 \\\ 0 \\\ 0\end{bmatrix} – \left(\begin{bmatrix}18 \\\ 0 \\\ 0\end{bmatrix}, \begin{bmatrix} \frac{2}{3} \\\ \frac{2}{3} \\\ \frac{1}{3} \end{bmatrix}\right)\begin{bmatrix} \frac{2}{3} \\\ \frac{2}{3} \\\ \frac{1}{3} \end{bmatrix} – \left(\begin{bmatrix}18 \\\ 0 \\\ 0\end{bmatrix}, \begin{bmatrix} -\frac{2}{3} \\\ \frac{1}{3} \\\ \frac{2}{3} \end{bmatrix} \right)\begin{bmatrix} -\frac{2}{3} \\\ \frac{1}{3} \\\ \frac{2}{3} \end{bmatrix}</span><br>
<span data-display="true" class="wp-katex-eq"> v_3 = \begin{bmatrix}2 \\\ – 4 \\\ 4 \end{bmatrix} \qquad e_3 = \frac{v_3}{||v_3||} = \frac{\begin{bmatrix}2 \\\ -4 \\\ 4\end{bmatrix}}{\sqrt{\sum{\begin{bmatrix}2 \\\ -4 \\\ 4\end{bmatrix}^2}}}</span><br>
<span data-display="true" class="wp-katex-eq"> e_3 = \begin{bmatrix} \frac{1}{3} \\\ -\frac{2}{3} \\\ \frac{2}{3} \end{bmatrix} </span></p>
<p>Thus, the orthogonalized matrix resulting from the Gram-Schmidt process is:</p>
<p><span data-display="true" class="wp-katex-eq"> \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{1}{3} &amp; \frac{2}{3} \end{bmatrix} </span></p>
<p>The component <span data-display="false" class="wp-katex-eq">R</span> of the QR decomposition can also be found from the calculations made in the Gram-Schmidt process as defined above.</p>
<p><span data-display="true" class="wp-katex-eq">R = \begin{bmatrix}a_1 \cdot e_1 &amp; a_2 \cdot e_1 &amp; \cdots &amp; a_n \cdot e_1 \\\ 0 &amp; a_2 \cdot e_2 &amp; \cdots &amp; a_n \cdot e_2 \\\ \vdots &amp; \vdots &amp; &amp; \vdots \\\ 0 &amp; 0 &amp; \cdots &amp; a_n \cdot e_n \end{bmatrix} = \begin{bmatrix} \begin{bmatrix} 2 \\\ 2 \\\ 1 \end{bmatrix} \cdot \begin{bmatrix} \frac{2}{3} \\\ \frac{2}{3} \\\ \frac{1}{3} \end{bmatrix} &amp; \begin{bmatrix} -2 \\\ 1 \\\ 2 \end{bmatrix} \cdot \begin{bmatrix} \frac{2}{3} \\\ \frac{2}{3} \\\ \frac{1}{3} \end{bmatrix} &amp; \begin{bmatrix} 18 \\\ 0 \\\ 0 \end{bmatrix} \cdot \begin{bmatrix} \frac{2}{3} \\\ \frac{2}{3} \\\ \frac{1}{3} \end{bmatrix} \\\ 0 &amp; \begin{bmatrix} -2 \\\ 1 \\\ 2 \end{bmatrix} \cdot \begin{bmatrix} -\frac{2}{3} \\\ \frac{1}{3} \\\ \frac{2}{3} \end{bmatrix} &amp; \begin{bmatrix} 18 \\\ 0 \\\ 0 \end{bmatrix} \cdot \begin{bmatrix} -\frac{2}{3} \\\ \frac{1}{3} \\\ \frac{2}{3} \end{bmatrix} \\\ 0 &amp; 0 &amp; \begin{bmatrix} 18 \\\ 0 \\\ 0 \end{bmatrix} \cdot \begin{bmatrix} \frac{1}{3} \\\ -\frac{2}{3} \\\ \frac{2}{3} \end{bmatrix}\end{bmatrix}</span><br>
<span data-display="true" class="wp-katex-eq"> R = \begin{bmatrix} 3 &amp; 0 &amp; 12 \\\ 0 &amp; 3 &amp; -12 \\\ 0 &amp; 0 &amp; 6 \end{bmatrix} </span></p>
<h6>The Gram-Schmidt Algorithm in R</h6>
<p>We use the same matrix <span data-display="false" class="wp-katex-eq">A</span> to verify our results above.</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>
<p>The following function is an implementation of the Gram-Schmidt algorithm using the modified version of the algorithm. A good comparison of the classical and modified versions of the algorithm <a target="_blank" rel="nofollow" href="http://cavern.uark.edu/~arnold/4353/CGSMGS.pdf">can be found here</a>. The Modified Gram-Schmidt algorithm was used above due to its improved numerical stability, which results in more orthogonal columns over the Classical algorithm.</p>
<pre>
gramschmidt &lt;- function(x) {
  x &lt;- as.matrix(x)
  # Get the number of rows and columns of the matrix
  n &lt;- ncol(x)
  m &lt;- nrow(x)
  
  # Initialize the Q and R matrices
  q &lt;- matrix(0, m, n)
  r &lt;- matrix(0, n, n)
  
  for (j in 1:n) {
    v = x[,j] # Step 1 of the Gram-Schmidt process v1 = a1
    # Skip the first column
    if (j > 1) {
      for (i in 1:(j-1)) {
        r[i,j] &lt;- t(q[,i]) %*% x[,j] # Find the inner product (noted to be q^T a earlier)
        # Subtract the projection from v which causes v to become perpendicular to all columns of Q
        v &lt;- v - r[i,j] * q[,i] 
      }      
    }
    # Find the L2 norm of the jth diagonal of R
    r[j,j] &lt;- sqrt(sum(v^2))
    # The orthogonalized result is found and stored in the ith column of Q.
    q[,j] &lt;- v / r[j,j]
  }
  
  # Collect the Q and R matrices into a list and return
  qrcomp &lt;- list('Q'=q, 'R'=r)
  return(qrcomp)
}
</pre>
<p>Perform the Gram-Schmidt orthogonalization process on the matrix <span data-display="false" class="wp-katex-eq">A</span> using our function.</p>
<pre>
gramschmidt(A)
## $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    0   12
## [2,]    0    3  -12
## [3,]    0    0    6
</pre>
<p>The results of our function match those of our manual calculations!</p>
<p>The <code>qr()</code> function in R also performs the Gram-Schmidt process. The <a target="_blank" rel="nofollow" href="https://stat.ethz.ch/R-manual/R-devel/library/base/html/qr.html">qr() function</a> does not output the <span data-display="false" class="wp-katex-eq">Q</span> and <span data-display="false" class="wp-katex-eq">R</span> matrices, which must be found by calling <code>qr.Q()</code> and <code>qr.R()</code>, respectively, on the <code>qr</code> object.</p>
<pre>A.qr &lt;- qr(A)
A.qr.out &lt;- list('Q'=qr.Q(A.qr), 'R'=qr.R(A.qr))
A.qr.out
## $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    0  -12
## [2,]    0   -3   12
## [3,]    0    0    6
</pre>
<p>Thus the <code>qr()</code> function in R matches our function and manual calculations as well.</p>
<h6>References</h6>
<p><a target="_blank" rel="nofollow" href="http://www.calpoly.edu/~jborzell/Courses/Year%2005-06/Spring%202006/304Gram_Schmidt_Exercises.pdf">http://www.calpoly.edu/~jborzell/Courses/Year%2005-06/Spring%202006/304Gram_Schmidt_Exercises.pdf</a></p>
<p><a target="_blank" rel="nofollow" href="http://cavern.uark.edu/~arnold/4353/CGSMGS.pdf">http://cavern.uark.edu/~arnold/4353/CGSMGS.pdf</a></p>
<p><a target="_blank" rel="nofollow" href="https://www.math.ucdavis.edu/~linear/old/notes21.pdf">https://www.math.ucdavis.edu/~linear/old/notes21.pdf</a></p>
<p><a target="_blank" rel="nofollow" href="http://www.math.ucla.edu/~yanovsky/Teaching/Math151B/handouts/GramSchmidt.pdf">http://www.math.ucla.edu/~yanovsky/Teaching/Math151B/handouts/GramSchmidt.pdf</a></p>
<p>The post <a target="_blank" rel="nofollow" href="http://www.aaronschlegel.com/qr-decomposition-gram-schmidt-algorithm/">QR Decomposition with the Gram-Schmidt Algorithm</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-gram-schmidt-algorithm/"> 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/OYzBi4i-yPc"></a></div>