A Second Course in Linear Algebra

Subsection 2.2.1 QR Decomposition via Gram-Schmidt

The Gram-Schmidt procedure is based on Theorem GSP. We begin with a set of linearly independent vectors and progressively convert them into a new set of nonzero vectors that form an orthogonal set, which can be easily converted to an orthonormal set. An orthonormal set is a wondrous thing, so the other portion of the conclusion is often overlooked. The new set of vectors spans the same subspace as the space spanned by the original set. This is half the reason that the Gram-Schmidt procedure is so useful.

As a preview of our main theorem, let's convert this idea into the language of matrices for a special case. Let \(A\) be a square matrix of size \(n\text{.}\) Take the columns of \(A\) as a set of vectors, which will form a basis of \(\complex{n}\text{.}\) Apply the Gram-Schmidt procedure to this set of \(n\) vectors. The manufacture of the \(i\)-th vector of this new set is the \(i\)-th vector of the original set, added to a linear combination of vectors \(1\) through \(i-1\) of the new set. If we recursively unpack these linear combinations, we can express each new vector as a linear combination of vectors \(1\) through \(i\) of the original set, where the \(i\)-th vector has a coefficient of \(1\text{.}\) Record the scalars of this linear combination in a column vector, whose last nonzero entry is \(1\text{.}\) Make these column vectors the columns of a square matrix, \(R^\prime\text{,}\) of size \(n\text{.}\) Define \(Q=AR^\prime\text{.}\)

By the Gram-Schmidt procedure, the columns of \(Q\) are an orthogonal set of nonzero vetors, and so \(\adjoint{Q}Q\) will be a diagonal matrix with nonzero entries. The matrix \(R^\prime\) is square, upper-triangular, and each diagonal entry is \(1\text{.}\) Hence \(R^\prime\) is invertible, so let \(R\) denote the inverse, which is again upper-triangular with diagonal entries equal to \(1\text{.}\) We then obtain \(A=QR\text{.}\) It is a simple matter to scale the columns of \(Q\) to form an orthonormal set, and the requisite scaling of the columns of \(R^\prime\) will not impede the existence of \(R\text{,}\) though we can only claim diagonal entries are nonzero. In this way, we can claim that \(Q\) is a unitary matrix.

A QR decomposition can be created for any matrix — it need not be square and it need not have full rank. The matrix \(Q\) is unitary, and \(R\) is upper triangular. Thus, each column of \(A\) can be expressed as a linear combination of the columns of \(Q\text{,}\) which form an orthonormal basis. So the column space of \(A\) is spanned by an orthonormal subset of the columns of \(Q\text{,}\) giving us the essence of the Gram-Schmidt procedure without the hypothesis that our original set is linearly independent. For the statement of Theorem GSP, it was a convenience to hypothesize that \(S\) is linearly independent. Can you examine the proof and see what changes are required if we lift this hypothesis?

We now state, and prove, a sequence of theorems which solidify the discussion above and generalize to rectangular matrices that may not have full rank.

Notice that the rank \(n\) condition of Theorem Theorem 2.2.1 necessarily implies \(m\geq n\text{.}\) We can expand \(Q\) and \(R\) simultaneously to get a decomposition where \(Q\) is a unitary matrix.

The column space of a matrix is an important property of a matrix. For example, the column space of the coefficient matrix of a system of equations is the set of vectors that can be used to form a consistent system when paired with the coefficient matrix (Theorem CSCS). Not only does \(Q\) have a column space equal to that of \(A\text{,}\) the first \(i\) columns of \(Q\) are a basis for the space spanned by the first \(i\) columns of \(A\text{,}\) for \(1\leq i\leq n\text{.}\)

The decomposition of Theorem Theorem 2.2.1 is referred to as a thin decomposition, while the decomposition of Theorem Theorem 2.2.2 is referred to as a full decomposition.

What happens if \(A\) does not have full rank? Then there will be some relation of linear dependence on the columns of the matrix. In the course of working the Gram-Schmidt procedure, as in the construction given for the proof of Theorem Theorem 2.2.1, this will be discovered when the newly created vector is the zero vector. No matter, we simply add any vector as the next column of \(Q\) which will preserve the columns as an orthonormal set. This vector can be determined in a fashion entirely similar to the device used in the proof of Theorem Theorem 2.2.2. However, this nearly arbitrary choice outside the span of the current set of columns, requires that we add a row of zeros in \(R\) and lose our positive diagonal entry. The ultimate price for this is that certain uniqueness results [provisional cross-reference: uniqueness of QR] no longer hold.

For the most general case of a QR decomposition for a rectangular matrix of arbitrary rank, we could fashion a proof based on the discussion of the previous paragraph. However, while Gram-Schmidt provides a good theoretical grounding for the QR decomposition, in numerical contexts its performance is weak ([provisional cross-reference: TB, Lecture 7 to 10]). Instead we give a constructive proof based on Householder reflections.

Subsection 2.2.2 QR Decomposition via Householder Reflections

Householder reflections provide a usefull tool for creating the zero entries required in each column of \(R\text{,}\) and will provide an algorithm with better performance (speed and accuracy), compared to the more intuitive approach of the Gram-Schmidt procedure.

The inductive proof of Theorem Theorem 2.2.3 will automatically provide a recipe for a recursive procedure to compute a QR decomposition. But recursion is rarely, if ever, a good idea for efficiency. Fortunately, the proof suggests a better procedure. Work on \(A\) column-by-column, progressively using Householder matrices to “zero out” each column below the diagonal entry. The construction of such a Householder matrix will require a nonzero entry in the diagonal entry. A row swap, accomplished by a (unitary) permutation matrix, can move a nonzero entry from elsewhere in the column. Of course, if the whole remainder of the column is all zeros, then we can just move on to the next column.

Notice how the unitary matrices change over the course of these iterations. In later steps each unitary matrix has a larger identity matrix as the block in the upper-left corner of the matrix. So if a product of these matrices is required, it can be more efficient to begin building up the product starting with the last unitary matrix. Of course, order matters in matrix multiplication, and maybe you need the product in the reverse order. No matter, compute the transpose of the desired matrix, which necessarily will reverse the order of the product. Once the product is concluded, then transpose the result.

Subsection 2.2.3 Solving Systems with a QR Decomposition

Consider the problem of solving a linear system \(A\vect{x}=\vect{b}\text{.}\) Replace \(A\) by a QR decomposition, to obtain \(QR\vect{x}=\vect{b}\text{.}\) Now, the inverse of the unitary matrix \(Q\) is its adjoint, so we have the new system \(R\vect{x}=\adjoint{Q}\vect{b}\text{.}\) Since \(R\) is upper-triangular, the new system can be solved via back-solving, much as we did with an LU decomposition (Subsection Subsection 2.1.2).

A primary advantage of this approach is that a unitary matrix is generally well-behaved in numerical computations. Since the columns have unit norm, the entries can never be too large. Many of the problems of floating-point arithmetic come from combining numbers of grossly different magnitudes, and this is avoided in computations with a unitary matrix, such as the matrix-vector product in the rearrangement of a linear system in the previous paragraph.

What about the number of operations required? For a QR decomposition of an \(m\times n\) matrix via the Gram-Schmidt procedure the operation count is \(\orderof{2mn^2}\text{.}\) When \(m\geq n\) the determination of \(R\) can be accomplished in \(\orderof{2mn^2-\frac{2}{3}n^3}\) operations when using a sequence of Householder reflections ([provisional cross-reference: TB Lecture 10]). This includes storing the Householder vector for each iteration, but not the computation of the Householder matrix itself, or the accumulated product of all the Householder matrices into one grand unitary matrix \(Q\text{.}\) Notice that for a square matrix, when \(m=n\text{,}\) the count is \(\orderof{\frac{4}{3}n^3}\text{,}\) which is twice the cost of an LU factorization

If we use Householder reflections, then we must decide what we want to do with Householder matrices. If we are solving a linear system, we can successively multiply each new Householder matrix times the column vector on the right hand side of the equation. Thus, we modify both sides of the equation in the same way, producing equivalent systems as we go. In practice, the product of a column vector by a Householder matrix can be accomplished very efficiently, and not by explicitly forming the matrix. (See Exercise Checkpoint 1.5.8.) It should never be necessary to explicitly form a Householder matrix, but instead the Householder vector should be enough information to perform whatever computation is at hand.

Subsection 2.2.4 Uniqueness of QR Decompositions

Generally, when \(A\) has full rank and we require that diagonal entries of \(R\) be positive, we get a unique QR decomposition. Proving this will be easier once we learn about Cholesky decompositions (Section [provisional cross-reference: a section about Cholesky decompositions]), so we will defer a proof until then.

Subsection 2.2.5 Final Thoughts