https://rstudio-pubs-static.s3.amazonaws.com/260581_5014d1d133354f5084d2ea67b0aee649.html

QR Decomposition with the Gram-Schmidt Procedure

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 A
, into two components, Q, and R.

A=QR

Where Q
is an orthogonal matrix, and R is an upper triangular matrix. Recall an orthogonal matrix is a square matrix with orthonormal row and column vectors such that QTQ=I, where I is the identity matrix. The term orthonormal implies the vectors are of unit length and are perpendicular (orthogonal) to each other.
QR decomposition is often used in linear least squares estimation and is, in fact, the method used by R in its lm() function. Signal processing and MIMO systems 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.

The Gram-Schmidt Process

The Gram-Schmidt process is used to find an orthogonal basis from a non-orthogonal basis. An orthogonal basis 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:

||an||=an⋅an−−−−−√=aTnan−−−−√=1

Where an
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 qn for each column vector an and then subtracting its projections onto the previous projections (qj). The resulting vector is then divided by the length of that vector to produce a unit vector.
Consider a matrix A
with n column vectors such that:

A=[a1|a2|⋯|an]

The Gram-Schmidt process proceeds by finding the orthogonal projection of the first column vector a1
.

v1=a1,e1=v1||v1||

Because a1
is the first column vector, there is no preceeding projections to subtract. The second column a2 is subtracted by the previous projection on the column vector:

v2=a2−projv1(a2)=a2−(a2⋅e1)e1,e2=v2||v2||

This process continues up to the n
column vectors, where each incremental step k+1 is computed as:

vk+1=ak+1−(ak+1⋅e1)e1−⋯−(ak+1⋅ek)ek,ek+1=uk+1||uk+1||

The ||⋅||
is the L2 norm which is defined as:

∑j=1mv2k−−−−−⎷

The projection can also be defined by:
Thus the matrix A
can be factorized into the QR matrix as the following:

A=[a1|a2|⋯|an]=[e1|e2|⋯|en]⎡⎣⎢⎢⎢⎢a1⋅e10⋮0a2⋅e1a2⋅e2⋮0⋯⋯⋯an⋅e1an⋅e2⋮an⋅en⎤⎦⎥⎥⎥⎥=QR

Gram-Schmidt Process Example

Consider the matrix A

:

⎡⎣⎢221−2121800⎤⎦⎥

We would like to orthogonalize this matrix using the Gram-Schmidt process. The resulting orthogonalized vector is also equivalent to Q
in the QR decomposition.
The Gram-Schmidt process on the matrix A
proceeds as follows:

v1=a1=⎡⎣⎢221⎤⎦⎥e1=v1||v1||=⎡⎣⎢221⎤⎦⎥∑⎡⎣⎢221⎤⎦⎥2−−−−−−−⎷

e1=⎡⎣⎢⎢232313⎤⎦⎥⎥

v2=a2−(a2⋅e1)e1=⎡⎣⎢−212⎤⎦⎥−⎛⎝⎜⎜⎡⎣⎢−212⎤⎦⎥,⎡⎣⎢⎢232313⎤⎦⎥⎥⎞⎠⎟⎟⎡⎣⎢⎢232313⎤⎦⎥⎥

v2=⎡⎣⎢−212⎤⎦⎥e2=v2||v2||=⎡⎣⎢−212⎤⎦⎥∑⎡⎣⎢−212⎤⎦⎥2−−−−−−−−−⎷

e2=⎡⎣⎢⎢−231323⎤⎦⎥⎥

v3=a3−(a3⋅e1)e1−(a3⋅e2)e2

v3=⎡⎣⎢1800⎤⎦⎥−⎛⎝⎜⎜⎡⎣⎢1800⎤⎦⎥,⎡⎣⎢⎢232313⎤⎦⎥⎥⎞⎠⎟⎟⎡⎣⎢⎢232313⎤⎦⎥⎥−⎛⎝⎜⎜⎡⎣⎢1800⎤⎦⎥,⎡⎣⎢⎢−231323⎤⎦⎥⎥⎞⎠⎟⎟⎡⎣⎢⎢−231323⎤⎦⎥⎥

v3=⎡⎣⎢2−44⎤⎦⎥e3=v3||v3||=⎡⎣⎢2−44⎤⎦⎥∑⎡⎣⎢2−44⎤⎦⎥2−−−−−−−−−⎷

e3=⎡⎣⎢⎢13−2323⎤⎦⎥⎥

Thus, the orthogonalized matrix resulting from the Gram-Schmidt process is:

⎡⎣⎢⎢232313−23131313−2323⎤⎦⎥⎥

The component R
of the QR decomposition can also be found from the calculations made in the Gram-Schmidt process as defined above.

R=⎡⎣⎢⎢⎢⎢a1⋅e10⋮0a2⋅e1a2⋅e2⋮0⋯⋯⋯an⋅e1an⋅e2⋮an⋅en⎤⎦⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡⎣⎢221⎤⎦⎥⋅⎡⎣⎢⎢232313⎤⎦⎥⎥00⎡⎣⎢−212⎤⎦⎥⋅⎡⎣⎢⎢232313⎤⎦⎥⎥⎡⎣⎢−212⎤⎦⎥⋅⎡⎣⎢⎢−231323⎤⎦⎥⎥0⎡⎣⎢1800⎤⎦⎥⋅⎡⎣⎢⎢232313⎤⎦⎥⎥⎡⎣⎢1800⎤⎦⎥⋅⎡⎣⎢⎢−231323⎤⎦⎥⎥⎡⎣⎢1800⎤⎦⎥⋅⎡⎣⎢⎢13−2323⎤⎦⎥⎥⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

R=⎡⎣⎢30003012−126⎤⎦⎥

The Gram-Schmidt Algorithm in R

We use the same matrix A

to verify our results above.

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

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 can be found here. The Modified Gram-Schmidt algorithm was used above due to its improved numerical stability, which results in more orthogonal columns over the Classical algorithm.

gramschmidt <- function(x) {
  x <- as.matrix(x)
  # Get the number of rows and columns of the matrix
  n <- ncol(x)
  m <- nrow(x)
  
  # Initialize the Q and R matrices
  q <- matrix(0, m, n)
  r <- 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] <- 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 <- v - r[i,j] * q[,i] 
      }      
    }
    # Find the L2 norm of the jth diagonal of R
    r[j,j] <- sqrt(sum(v^2))
    # The orthogonalized result is found and stored in the ith column of Q.
    q[,j] <- v / r[j,j]
  }
  
  # Collect the Q and R matrices into a list and return
  qrcomp <- list('Q'=q, 'R'=r)
  return(qrcomp)
}

Perform the Gram-Schmidt orthogonalization process on the matrix A
using our function.

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

The results of our function match those of our manual calculations!
The qr() function performs QR decomposition using Householder reflections, but we can use it to verify our results in this case. Householder reflections is the more common approach to performing QR decomposition as it is more numerically stable. The Gram-Schmidt process can result in a non-orthogonal Q
matrix due to rounding errors. The qr() function does not output the Q and R

matrices, which must be found by calling qr.Q() and qr.R(), respectively, on the qr object.

A.qr <- qr(A)
A.qr.out <- 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

Thus the qr() function in R matches our function and manual calculations as well.

References

http://www.calpoly.edu/~jborzell/Courses/Year%2005-06/Spring%202006/304Gram_Schmidt_Exercises.pdf
http://cavern.uark.edu/~arnold/4353/CGSMGS.pdf
https://www.math.ucdavis.edu/~linear/old/notes21.pdf
http://www.math.ucla.edu/~yanovsky/Teaching/Math151B/handouts/GramSchmidt.pdf