QR Factorization

One of the main virtues of orthogonal matrices is that they can be easily inverted—the transpose is the inverse. This fact, combined with the factorization theorem in this section, provides a useful way to simplify many matrix calculations (for example, in least squares approximation).

The importance of the factorization lies in the fact that there are computer algorithms that accomplish it with good control over round-off error, making it particularly useful in matrix calculations. The factorization is a matrix version of the Gram-Schmidt process.

Suppose \(A = \left [ \begin{array}{cccc} |&|& &| \\ \vec {c}_{1} & \vec {c}_{2} & \cdots & \vec {c}_{n}\\ |&|& &| \end{array}\right ]\) is an \(m \times n\) matrix with linearly independent columns \(\vec {c}_{1}, \vec {c}_{2}, \dots , \vec {c}_{n}\). The Gram-Schmidt algorithm can be applied to these columns to provide orthogonal columns \(\vec {f}_{1}, \vec {f}_{2}, \dots , \vec {f}_{n}\) where \(\vec {f}_{1} = \vec {c}_{1}\) and

\begin{equation*} \vec {f}_{k} = \vec {c}_{k} - \frac {\vec {c}_{k} \dotp \vec {f}_{1}}{\norm { \vec {f}_{1} }^2}\vec {f}_{1} + \frac {\vec {c}_{k} \dotp \vec {f}_{2}}{\norm { \vec {f}_{2} }^2}\vec {f}_{2} - \dots - \frac {\vec {c}_{k} \dotp \vec {f}_{k-1}}{\norm { \vec {f}_{k-1} }^2}\vec {f}_{k-1} \end{equation*}

for each \(k = 2, 3, \dots , n\). Now write \(\vec {q}_{k} = \frac {1}{\norm { \vec {f}_{k} }}\vec {f}_{k}\) for each \(k\). Then \(\vec {q}_{1}, \vec {q}_{2}, \dots , \vec {q}_{n}\) are orthonormal columns, and the above equation becomes

\begin{equation*} \norm { \vec {f}_{k} } \vec {q}_{k} = \vec {c}_{k} - (\vec {c}_{k} \dotp \vec {q}_{1})\vec {q}_{1} - (\vec {c}_{k} \dotp \vec {q}_{2})\vec {q}_{2} - \dots - (\vec {c}_{k} \dotp \vec {q}_{k-1})\vec {q}_{k-1} \end{equation*}

Using these equations, express each \(\vec {c}_{k}\) as a linear combination of the \(\vec {q}_{i}\):

\begin{equation*} \begin{array}{ccl} \vec {c}_{1} &=& \norm { \vec {f}_{1} } \vec {q}_{1} \\ \vec {c}_{2} &=& (\vec {c}_{2} \dotp \vec {q}_{1})\vec {q}_{1} + \norm { \vec {f}_{2} } \vec {q}_{2} \\ \vec {c}_{3} &=& (\vec {c}_{3} \dotp \vec {q}_{1})\vec {q}_{1} + (\vec {c}_{3} \dotp \vec {q}_{2})\vec {q}_{2} + \norm { \vec {f}_{3} } \vec {q}_{3} \\ \vdots && \vdots \\ \vec {c}_{n} &=& (\vec {c}_{n} \dotp \vec {q}_{1})\vec {q}_{1} + (\vec {c}_{n} \dotp \vec {q}_{2})\vec {q}_{2} + (\vec {c}_{n} \dotp \vec {q}_{3})\vec {q}_{3} + \dots + \norm { \vec {f}_{n} } \vec {q}_{n} \end{array} \end{equation*}

These equations have a matrix form that gives the required factorization:

\begin{align} A & = \left [ \begin{array}{ccccc} |&|&|& &| \\ \vec {c}_{1} & \vec {c}_{2} & \vec {c}_{3} &\cdots & \vec {c}_{n}\\ |&|&|& &| \end{array}\right ] \nonumber \\ &= \left [ \begin{array}{ccccc} |&|&|& &| \\ \vec {q}_{1} & \vec {q}_{2} & \vec {q}_{3} & \cdots & \vec {q}_{n}\\ |&|&|& &| \end{array}\right ] \left [ \begin{array}{ccccc} \norm { \vec {f}_{1} } & \vec {c}_{2} \dotp \vec {q}_{1} & \vec {c}_{3} \dotp \vec {q}_{1} & \cdots & \vec {c}_{n} \dotp \vec {q}_{1} \\ 0 & \norm { \vec {f}_{2} } & \vec {c}_{3} \dotp \vec {q}_{2} & \cdots & \vec {c}_{n} \dotp \vec {q}_{2} \\ 0 & 0 & \norm { \vec {f}_{3} } & \cdots & \vec {c}_{n} \dotp \vec {q}_{3} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & \norm { \vec {f}_{n} } \end{array} \right ] \label {matrixFactEq} \end{align}

Here the first factor \(Q = \left [ \begin{array}{ccccc} |&|&|& &| \\ \vec {q}_{1} & \vec {q}_{2} & \vec {q}_{3} & \cdots & \vec {q}_{n}\\ |&|&|& &| \end{array}\right ]\) has orthonormal columns, and the second factor is an \(n \times n\) upper triangular matrix \(R\) with positive diagonal entries (and so is invertible). We record this in the following theorem.

The matrices \(Q\) and \(R\) in Theorem th:QR-025133 are uniquely determined by \(A\); we return to this below.

If a matrix \(A\) has independent rows and we apply QR-factorization to \(A^{T}\), the result is:

Since a square matrix with orthonormal columns is orthogonal, we have

We now take the time to prove the uniqueness of the QR-factorization.

Write \(Q = \left [ \begin{array}{cccc} |&|& &| \\ \vec {c}_{1} & \vec {c}_{2} & \cdots & \vec {c}_{n}\\ |&|& &| \end{array}\right ]\) and \(Q_{1} = \left [ \begin{array}{cccc} |&|& &| \\ \vec {d}_{1} & \vec {d}_{2} & \cdots & \vec {d}_{n}\\ |&|& &| \end{array}\right ]\) in terms of their columns, and observe first that \(Q^TQ = I_{n} = Q_{1}^TQ_{1}\) because \(Q\) and \(Q_{1}\) have orthonormal columns. Hence it suffices to show that \(Q_{1} = Q\) (then \(R_{1} = Q_{1}^TA = Q^TA = R\)). Since \(Q_{1}^TQ_{1} = I_{n}\), the equation \(QR = Q_{1}R_{1}\) gives \(Q_{1}^TQ = R_{1}R^{-1}\); for convenience we write this matrix as

\begin{equation*} Q_{1}^TQ = R_{1}R^{-1} = \left [ \begin{array}{c} t_{ij} \end{array}\right ] \end{equation*}
This matrix is upper triangular with positive diagonal elements (since this is true for \(R\) and \(R_{1}\)), so \(t_{ii} > 0\) for each \(i\) and \(t_{ij} = 0\) if \(i > j\). On the other hand, the \((i, j)\)-entry of \(Q_{1}^TQ\) is \(\vec {d}_{i}^T\vec {c}_{j} = \vec {d}_{i} \dotp \vec {c}_{j}\), so we have \(\vec {d}_{i} \dotp \vec {c}_{j} = t_{ij}\) for all \(i\) and \(j\). But each \(\vec {c}_{j}\) is in \(\mbox {span}\{\vec {d}_{1}, \vec {d}_{2}, \dots , \vec {d}_{n}\}\) because \(Q = Q_{1}(R_{1}R^{-1})\). We know how to write a vector as a linear combination of an orthonormal basis (using Corollary cor:orthonormal from Orthogonality and Projections):
\begin{equation*} \vec {c}_{j} = (\vec {d}_{1} \dotp \vec {c}_{j})\vec {d}_{1} + (\vec {d}_{2} \dotp \vec {c}_{j})\vec {d}_{2} + \dots + (\vec {d}_{n} \dotp \vec {c}_{j})\vec {d}_{n} = t_{1j}\vec {d}_{1} + t_{2j}\vec {d}_{2} + \dots + t_{jj}\vec {d}_{i}, \end{equation*}
because \(\vec {d}_{i} \dotp \vec {c}_{j} = t_{ij} = 0\) if \(i > j\). The first few equations here are
\begin{equation*} \begin{array}{ccl} \vec {c}_{1} &=& t_{11}\vec {d}_{1} \\ \vec {c}_{2} &=& t_{12}\vec {d}_{1} + t_{22}\vec {d}_{2} \\ \vec {c}_{3} &=& t_{13}\vec {d}_{1} + t_{23}\vec {d}_{2} + t_{33}\vec {d}_{3} \\ \vec {c}_{4} &=& t_{14}\vec {d}_{1} + t_{24}\vec {d}_{2} + t_{34}\vec {d}_{3} + t_{44}\vec {d}_{4} \\ \vdots && \vdots \end{array} \end{equation*}
The first of these equations gives \(1 = \norm { \vec {c}_{1} } = \norm { t_{11}\vec {d}_{1} } = | t_{11} | \norm { \vec {d}_{1} } = t_{11}\), whence \(\vec {c}_{1} = \vec {d}_{1}\). But then we have \(t_{12} = \vec {d}_{1} \dotp \vec {c}_{2} = \vec {c}_{1} \dotp \vec {c}_{2} = 0\), so the second equation becomes \(\vec {c}_{2} = t_{22}\vec {d}_{2}\). Now a similar argument gives \(\vec {c}_{2} = \vec {d}_{2}\), and then \(t_{13} = 0\) and \(t_{23} = 0\) follows in the same way. Hence \(\vec {c}_{3} = t_{33}\vec {d}_{3}\) and \(\vec {c}_{3} = \vec {d}_{3}\). Continue in this way to get \(\vec {c}_{i} = \vec {d}_{i}\) for all \(i\). This proves that \(Q_{1} = Q\).

QR-Algorithm for approximating eigenvalues

In The Power Method and the Dominant Eigenvalue, we learned about an iterative method for computing eigenvalues. We also mentioned that a better method for approximating the eigenvalues of an invertible matrix \(A\) depends on the QR-factorization of \(A\). While it is beyond the scope of this book to pursue a detailed discussion of this method, we give an example and conclude with some remarks on the QR-algorithm. The interested reader is referred to J. M. Wilkinson, The Algebraic Eigenvalue Problem (Oxford, England: Oxford University Press, 1965) or G. W. Stewart, Introduction to Matrix Computations (New York: Academic Press, 1973).

The QR-algorithm uses QR-factorization repeatedly to create a sequence of matrices \(A_{1} =A, A_{2}, A_{3}, \dots ,\) as follows:

Shifting.

We first learned about the concept of shifting in The Power Method and the Dominant Eigenvalue. Convergence is accelerated if, at stage \(k\) of the algorithm, a number \(\tau _{k}\) is chosen and \(A_{k} - \tau _{k}I\) is factored in the form \(Q_{k}R_{k}\) rather than \(A_{k}\) itself. Then

\begin{equation*} Q_{k}^{-1}A_{k}Q_{k} = Q_{k}^{-1}(Q_{k}R_{k} + \tau _{k}I)Q_{k} = R_{k}Q_{k} + \tau _{k}I \end{equation*}

so we take \(A_{k+1} = R_{k}Q_{k} + \tau _{k}I\). If the shifts \(\tau _{k}\) are carefully chosen, convergence can be greatly improved.

Preliminary Preparation.

A matrix such as

\begin{equation*} \left [ \begin{array}{rrrrr} \ast & \ast & \ast & \ast & \ast \\ \ast & \ast & \ast & \ast & \ast \\ 0 & \ast & \ast & \ast & \ast \\ 0 & 0 & \ast & \ast & \ast \\ 0 & 0 & 0 & \ast & \ast \\ \end{array}\right ] \end{equation*}

is said to be in upper Hessenberg form, and the QR-factorizations of such matrices are greatly simplified. Given an \(n \times n\) matrix \(A\), a series of orthogonal matrices \(H_{1}, H_{2}, \dots , H_{m}\) (called Householder matrices can be easily constructed such that

\begin{equation*} B = H_{m}^T \cdots H_{1}^TAH_{1} \cdots H_{m} \end{equation*}

is in upper Hessenberg form. Then the QR-algorithm can be efficiently applied to \(B\) and, because \(B\) is similar to \(A\), it produces the eigenvalues of \(A\).

Complex Eigenvalues.

If some of the eigenvalues of a real matrix \(A\) are not real, the QR-algorithm converges to a block upper triangular matrix where the diagonal blocks are either \(1 \times 1\) (the real eigenvalues) or \(2 \times 2\) (each providing a pair of conjugate complex eigenvalues of \(A\)).

Practice Problems

In each case find the QR-factorization of \(A\).
(a)
\(A = \left [ \begin{array}{rr} 1 & -1 \\ -1 & 0 \end{array}\right ]\)
(b)
\(A = \left [ \begin{array}{rr} 2 & 1 \\ 1 &1 \end{array}\right ]\)
(c)
\(A = \left [ \begin{array}{rrr} 1 & 1 & 1 \\ 1 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{array}\right ]\)
(d)
\(A = \left [ \begin{array}{rrr} 1 & 1 & 0 \\ -1 & 0 & 1 \\ 0 & 1 & 1 \\ 1 & -1 & 0 \end{array}\right ]\)

Click the arrow to see answer.

(a)
\(Q = \frac {1}{\sqrt {5}}\left [ \begin{array}{rr} 2 & -1 \\ 1 & 2 \end{array}\right ]\), \(R = \frac {1}{\sqrt {5}}\left [ \begin{array}{rr} 5 & 3 \\ 0 & 1 \end{array}\right ]\)
(b)
\(Q = \frac {1}{\sqrt {3}}\left [ \begin{array}{rrr} 1 & 1 & 0 \\ -1 & 0 & 1 \\ 0 & 1 & 1 \\ 1 & -1 & 1 \end{array}\right ]\),
\(R = \frac {1}{\sqrt {3}}\left [ \begin{array}{rrr} 3 & 0 & -1 \\ 0 & 3 & 1 \\ 0 & 0 & 2 \end{array}\right ]\)
If \(R\) is upper triangular and invertible, show that there exists a diagonal matrix \(D\) with diagonal entries \(\pm 1\) such that \(R_{1} = DR\) is invertible, upper triangular, and has positive diagonal entries.
If \(A\) has independent columns, let
\(A = QR\) where \(Q\) has orthonormal columns and \(R\) is invertible and upper triangular. (Some authors do not require a \(QR\)-factorization to have positive diagonal entries.) Show that there is a diagonal matrix \(D\) with diagonal entries \(\pm 1\) such that \(A = (QD)(DR)\) is the QR-factorization of \(A\).
See Practice Problem
In each case, find the exact eigenvalues and then approximate them using the QR-algorithm.
(a)
\(A = \left [ \begin{array}{rr} 1 & 1 \\ 1 & 0 \end{array}\right ]\)
(b)
\(A = \left [ \begin{array}{rr} 3 & 1 \\ 1 & 0 \end{array}\right ]\)

Click the arrow to see answer.

(a)
Eigenvalues \(\lambda _{1} = \frac {1}{2}(3 + \sqrt {13}) = 3.302776\), \(\lambda _{2} = \frac {1}{2}(3 - \sqrt {13}) = -0.302776\)

\(A_{1} = \left [ \begin{array}{rr} 3 & 1 \\ 1 & 0 \end{array}\right ]\), \(Q_{1} = \frac {1}{\sqrt {10}}\left [ \begin{array}{rr} 3 & -1 \\ 1 & 3 \end{array}\right ]\), \(R_{1} = \frac {1}{\sqrt {10}}\left [ \begin{array}{rr} 10 & 3 \\ 0 & -1 \end{array}\right ]\)
\(A_{2} = \frac {1}{10}\left [ \begin{array}{rr} 33 & -1 \\ -1 & -3 \end{array}\right ]\),
\(Q_{2} = \frac {1}{\sqrt {1090}}\left [ \begin{array}{rr} 33 & 1 \\ -1 & 33 \end{array}\right ]\),
\(R_{2} = \frac {1}{\sqrt {1090}}\left [ \begin{array}{rr} 109 & -3 \\ 0 & -10 \end{array}\right ]\)
\(A_{3} = \frac {1}{109}\left [ \begin{array}{rr} 360 & 1 \\ 1 & -33 \end{array}\right ]\)
\({} = \left [ \begin{array}{rr} 3.302775 & 0.009174 \\ 0.009174 & -0.302775 \end{array}\right ]\)

If \(A\) is symmetric, show that each matrix \(A_{k}\) in the QR-algorithm is also symmetric. Deduce that they converge to a diagonal matrix.

Click the arrow to see answer.

Use induction on \(k\). If \(k = 1\), \(A_{1} = A\). In general \(A_{k+1} = Q_{k}^{-1}A_{k}Q_{k} = Q_{k}^{T}A_{k}Q_{k}\), so the fact that \(A_{k}^{T} = A_{k}\) implies \(A_{k+1}^{T} = A_{k+1}\). The eigenvalues of \(A\) are all real (Theorem ??), so the \(A_{k}\) converge to an upper triangular matrix \(T\). But \(T\) must also be symmetric (it is the limit of symmetric matrices), so it is diagonal.
Apply the QR-algorithm to
\(A = \left [ \begin{array}{rr} 2 & -3 \\ 1 & -2 \end{array}\right ]\). Explain.
Given a matrix \(A\), let \(A_{k}\), \(Q_{k}\), and \(R_{k}\), \(k \geq 1\), be the matrices constructed in the QR-algorithm. Show that \(A_{k} = (Q_{1}Q_{2} \cdots Q_{k})(R_{k} \cdots R_{2}R_{1})\) for each \(k \geq 1\) and hence that this is a QR-factorization of \(A_{k}\).
Show that \(Q_{k}R_{k} = R_{k-1}Q_{k-1}\) for each \(k \geq 2\), and use this equality to compute \((Q_{1}Q_{2} \cdots Q_{k})(R_{k} \cdots R_{2}R_{1})\) “from the centre out.” Use the fact that \((AB)^{n+1} = A(BA)^{n}B\) for any square matrices \(A\) and \(B\).

Text Source

This section was adapted from Sections 8.4 and 8.5 of Keith Nicholson’s Linear Algebra with Applications. (CC-BY-NC-SA)

W. Keith Nicholson, Linear Algebra with Applications, Lyryx 2018, Open Edition, pp. 437–445.