Definition 6.3.1.
The QR factorization of a matrix \(A\) is given by
\begin{equation*}
A = QR
\end{equation*}
where \(Q\) is an orthogonal matrix and \(R\) is an upper triangular matrix.
Section 6.3 QR Factorization
The QR factorization of a matrix is a decomposition of a matrix into the product of an orthogonal matrix and an upper triangular matrix. This factorization is useful in solving systems of linear equations, computing eigenvalues, and in numerical linear algebra.
Subsection 6.3.1 QR Factorization Algorithm
The QR-factorization can be performed using the Gram-Schmidt orthogonalization seen in Subsection 4.4.4. We will first show this with an example and then explain why this works.
Example 6.3.2.
Perform Gram-Schmidt Orthogonalization to find an orthonormal set of vectors of the columns of
\begin{equation*}
A = \begin{bmatrix} -9 \amp 4 \amp 7 \\ -18 \amp 15 \amp -28 \\ 6 \amp 2 \amp 49 \\ \end{bmatrix}
\end{equation*}
and then use the resulting matrix \(Q\) to find \(R\) using
\begin{equation*}
R = Q^{\intercal} A
\end{equation*}
Solution.
Let \(\boldsymbol{a}_j\) be the columns of \(A\) and we will call the resulting orthogonal vectors as \(\boldsymbol{q}_j\text{.}\)
\begin{equation*}
\boldsymbol{q}_1 = \frac{\boldsymbol{a}_1}{||\boldsymbol{a}_1||} = \frac{1}{\sqrt{9^2+18^2+6^2}} \begin{bmatrix} -9 \\ -18 \\ 6 \end{bmatrix} = \begin{bmatrix} -3/7 \\ -6/7 \\ 2/7 \end{bmatrix}
\end{equation*}
For the second vector,
\begin{equation*}
\begin{aligned} \boldsymbol{u}_2 \amp = \boldsymbol{a}_2 - \langle \boldsymbol{q}_1, \boldsymbol{a}_2 \rangle \boldsymbol{q}_1\\ \amp = \begin{bmatrix} 4 \\ 15 \\ 2 \end{bmatrix} - (-14) \begin{bmatrix} -3/7 \\ -6/7 \\ 2/7 \end{bmatrix} = \begin{bmatrix} 4 \\ 15 \\ 2 \end{bmatrix} + \begin{bmatrix} -6 \\ -12 \\ 4 \end{bmatrix} = \begin{bmatrix} -2 \\ 3 \\ 6 \end{bmatrix} \end{aligned}
\end{equation*}
and then this vectors is normalized with
\begin{equation*}
\boldsymbol{q}_2 = \frac{1}{||\boldsymbol{u}_2||} \boldsymbol{u_2} = \frac{1}{\sqrt{2^2+3^2+6^2}} \begin{bmatrix} -2 \\ 3 \\ 6 \end{bmatrix} = \frac{1}{7} \begin{bmatrix} -2 \\ 3 \\ 6 \end{bmatrix} = \begin{bmatrix} -2/7 \\ 3/7 \\ 6/7 \end{bmatrix}
\end{equation*}
Lastly, the third column of \(Q\) can be found with
\begin{equation*}
\begin{aligned} \boldsymbol{u}_3 \amp = \boldsymbol{a}_3 - \langle \boldsymbol{q}_2, \boldsymbol{a}_3 \rangle \boldsymbol{q}_2 - \langle \boldsymbol{q}_1, \boldsymbol{a}_3 \rangle \boldsymbol{q}_1 \\ \amp = \begin{bmatrix} 7 \\ -28 \\ 49 \end{bmatrix} - 28 \begin{bmatrix} -2/7 \\ 3/7 \\ 6/7 \end{bmatrix} - 35 \begin{bmatrix} -3/7 \\ -6/7 \\ 2/7 \end{bmatrix} \\ \amp = \begin{bmatrix} 7 \\ -28 \\ 49 \end{bmatrix} + \begin{bmatrix} 8 \\ -12 \\ -24 \end{bmatrix} + \begin{bmatrix} 15 \\ 30 \\ -10 \end{bmatrix} = \begin{bmatrix} 30 \\ -10 \\ 15 \end{bmatrix} \end{aligned}
\end{equation*}
and lastly, we normalize this vector as
\begin{equation*}
\boldsymbol{q3} = \frac{1}{||\boldsymbol{u}_2||} \boldsymbol{u_2} = \frac{1}{\sqrt{30^2+10^2+15^2}} \begin{bmatrix} 30 \\ -10 \\ 15 \end{bmatrix} = \frac{1}{35} \begin{bmatrix} 30 \\ -10 \\ 15 \end{bmatrix} = \begin{bmatrix} 6/7 \\ -2/7 \\ 3/7 \end{bmatrix}
\end{equation*}
and thus the matrix \(Q\) is the matrix with these three vectors or
\begin{equation*}
Q = \begin{bmatrix} -3/7 \amp -2/7 \amp 6/7 \\ -6/7 \amp 3/7 \amp -2/7 \\ 2/7 \amp 6/7 \amp 3/7 \end{bmatrix}
\end{equation*}
Lastly, we find \(R\) with
\begin{equation*}
\begin{aligned} R \amp = Q^{\intercal} A \\ \amp = \begin{bmatrix} -3/7 \amp -6/7 \amp 2/7 \\ -2/7 \amp 3/7 \amp 6/7 \\ 6/7 \amp -2/7 \amp 3/7 \end{bmatrix} \begin{bmatrix} -9 \amp 4 \amp 7 \\ -18 \amp 15 \amp -28 \\ 6 \amp 2 \amp 49 \end{bmatrix} \\ \amp = \begin{bmatrix} 21 \amp -14 \amp 35 \\ 0 \amp 7 \amp 28 \\ 0 \amp 0 \amp 35 \\ \end{bmatrix} \end{aligned}
\end{equation*}
A natural question to ask is why does this work? Why does the Gram-Schmidt process give us the QR factorization? The answer is that the Gram-Schmidt process gives us an orthonormal basis for the column space of \(A\) and thus the matrix \(Q\) is an orthogonal matrix. Then, the matrix \(R\) is the change of basis matrix from the standard basis to the orthonormal basis given by the columns of \(Q\text{.}\)
Examining the result of \(R = Q^{\intercal}A \text{,}\) the matrix \(Q^{\intercal}\) projects the columns of \(A\) onto \(Q\text{.}\) Because \(Q\) is constructed as the
Subsection 6.3.2 Solving Linear Systems with QR factoring.
We can solving linear systems with the QR factorization and is often done if the original matrix is difficult to solve (Section where this is explained??). Let’s look the example in Subsection 6.1.4.
-
Write \(A = QR\)
-
Solve \(R\boldsymbol{x} = Q^{\intercal}\boldsymbol{b}\text{.}\) Note that \(R\) is upper triangular and thus this can be solved with back substitution.
Example 6.3.3.
Solve the linear system \(A\boldsymbol{x}= \boldsymbol{b}\) using the \(QR\)-decomposition where
\begin{equation*}
A = \begin{bmatrix} 3 \amp -2 \amp 1 \\ 6 \amp 0 \amp 4 \\ -6 \amp -8 \amp -7 \end{bmatrix} \qquad \boldsymbol{b} = \begin{bmatrix} 14 \\ 22 \\ -9 \end{bmatrix}
\end{equation*}
Solution.
First, we need to find the factorization. We can follow the same steps using Gram-Schmidt Orthogonalization as in Example 6.3.2, however, generally the resulting elements of \(Q\) are not rational, so we will resort to software to do this. Using Matlab or Julia, the result in
\begin{equation*}
\begin{aligned} Q \amp = \begin{bmatrix} -0.333333 \amp -0.522976 \amp 0.784465 \\ -0.666667 \amp -0.457604 \amp -0.588348 \\ 0.666667 \amp -0.719092 \amp -0.196116 \\ \end{bmatrix} \\ R \amp = \begin{bmatrix} -9.0 \amp -4.66667 \amp -7.66667 \\ 0.0 \amp 6.79869 \amp 2.68025 \\ 0.0 \amp 0.0 \amp -0.196116 \\ \end{bmatrix} \end{aligned}
\end{equation*}
and this just shows the first 5 or so digits of the two matrices.
The next step is to find \(Q^{\intercal} \boldsymbol{b}\) or
\begin{equation*}
Q^{\intercal} \boldsymbol{b} = \begin{bmatrix} -25.33333333333333 \\ -10.917131522692245 \\ -0.19611613513818504 \\ \end{bmatrix}
\end{equation*}
And lastly, we solve \(R\boldsymbol{x} = Q^{\intercal} \boldsymbol{b}\text{.}\) This is done with back-substitution. The result of this is:
\begin{equation*}
\boldsymbol{x} = \begin{bmatrix} 2.999999999999997 \\ -2.000000000000001 \\ 1.0000000000000033 \\ \end{bmatrix}
\end{equation*}
which is close to the answer found in Subsection 6.1.4, the difference being that here we used floating-point numbers which have round off errors.
Subsection 6.3.3 Using QR factorization to find Eigenvalues
Another use of \(QR\)-factorization is that of numerically estimating the eigenvalues of a matrix.
Lemma 6.3.4.
The matrix \(A = QR\) for orthogonal matrix \(Q\) and upper triangular matrix \(R\) has the same eigenvalues as \(RQ\text{.}\)
Proof.
This amounts to showing that the matrices \(QR\) and \(RQ\) are similar.
\begin{equation*}
QR = S^{-1} RQ S
\end{equation*}
Let \(S = Q^{\intercal}\) and since \(Q\) is orthogonal, then \(S^{-1} = Q\)
\begin{equation*}
QR = Q RQ Q^{\intercal} = Q R
\end{equation*}
We use this fact to generate the following algorithm
Algorithm 6.3.5.
The following steps generates a sequence of matrices \(A_k\text{.}\) In the limit \(A_k\) goes to the diagonal matrix with the eigenvalues of \(A\) on the diagonal.
-
Let \(A_{k+1} = R_k Q_k\text{.}\) That is, reverse the order of the matrix multiplies.
-
If the maximum of the absolute value of the lower left elements of \(A_{k+1}\) is within \(\epsilon\text{,}\) then stop. Otherwise increment \(k\) by 1 and return to step 2.
We will show an example of this in action and then prove later why this converges to the correct answer.
Example 6.3.6.
Estimate the eigenvalues of the matrix
\begin{equation*}
A = \begin{bmatrix} 1 \amp 0 \amp 0 \\ 3 \amp -2 \amp 6 \\ 2 \amp 1 \amp 3 \end{bmatrix}
\end{equation*}
as found in Example 5.1.10
Solution.
First let \(A_0=A\text{.}\) Then factor \(A_0\) into \(Q_0\) and \(R_0\text{,}\) where \(A_0=Q_0R_0\text{.}\) This is done with software. If using Julia,
QR₀ = qr(A₀)
and the matrices are then located in or
Matrix(QR₀.Q)1
Note: if you leave off the
Matrix, then the matrix QR₀.Q is in a compact form.3×3 Matrix{Float64}:
-0.267261 0.145479 -0.952579
-0.801784 -0.581914 0.136083
-0.534522 0.800132 0.272166
and
QR₀.R or3×3 Matrix{Float64}:
-3.74166 1.06904 -6.41427
0.0 1.96396 -1.09109
0.0 0.0 1.63299
Then the order of the matrices is reversed to
A₁ = QR₀.R * QR₀.Q, resulting in3×3 Matrix{Float64}:
3.57143 -6.29869 1.96396
-0.99146 -2.01587 -0.0296957
-0.872872 1.30661 0.444444
Recall that the algorithm is supposed to result in an upper triangular matrix \(A_k\text{.}\) We’ll do one more step by hand. Let
QR₁ = qr(A₁). This gives us the two matrices (which we won’t show here), but the next \(A\) matrix is found by reversing the \(R\) and \(Q\) with A₂ = QR₁.R * QR₁.Q resulting in3×3 Matrix{Float64}:
4.65025 -4.45991 -2.88562
0.830799 -3.49591 -0.700881
0.199152 0.000981054 0.845661
Notice that the numbers in the lower right are getting smaller, but not very quickly. We’ll do a number of steps with
Aₖ = A
for i=1:25
QRₖ = qr(Aₖ)
Aₖ = QRₖ.R * QRₖ.Q
end
Aₖ
and the result is
3×3 Matrix{Float64}:
3.99934 -5.00092 3.53563
-0.000921692 -2.99934 0.706641
-2.34445e-15 1.27349e-12 1.0
The terms in the lower left are quite close to zero, so we can estimate that the eigenvalues are the diagonal elements or
3.99934, -2.9934 and 1.0.
We now discuss why this works.
