Subsection10.2.3A simple QR algorithm

We now morph the subspace iteration discussed in the last unit into a simple algorithm known as the QR algorithm. It starts by performing subspace iteration with an $m \times m$ (square) matrix so that the method finds all eigenvectors simultaneously (under mild conditions). Rather than starting with a random matrix $V \text{,}$ we now start with the identity matrix. This yields the algorithm on the left:

\begin{equation*} \begin{array}{l} \mbox{Subspace iteration with }\widehat V = I \\ \mbox{(a square matrix)} \\ \widehat V := I \\ {\bf for~} k:=0, \ldots \\ ~~~( \widehat V, R ) := {\rm QR}( A \widehat V ) \\ {\bf endfor} \\ \widehat \Lambda = \widehat V^H A \widehat V \\ \end{array} \end{equation*}

\begin{equation*} \begin{array}{l} \mbox{A simple QR algorithm} \\ \\ V = I \\ {\bf for~} k:=0, \ldots \\ ~~~( Q, R ) := {\rm QR}( A ) \\ ~~~ A = R Q \\ ~~~ V = V Q \\ {\bf endfor} \\ \end{array} \end{equation*}

We contrast this subspace iteration algorithm with the algorithm on the right, which is a simple incarnation of an algorithm known as the QR algorithm.

The magic lies in the fact that the matrices computed by the QR algorithm are identical to those computed by the subspace iteration: Upon completion $\widehat \Lambda = A$ and $\widehat V = V \text{.}$ To be able to prove this, we annotate the algorithm so we can reason about the contents of the matrices for each iteration:

\begin{equation*} \begin{array}{l} \mbox{Subspace iteration with }\widehat V = I \\ \mbox{(a square matrix)} \\ \widehat A^{(0)} = A \\ \widehat V^{(0)} := I \\ \widehat R^{(0)} := I \\ {\bf for~} k:=0, \ldots \\ ~~~( \widehat V^{(k+1)}, \widehat R^{(k+1)} ) := {\rm QR}( A \widehat V^{(k)} ) \\ ~~~ \widehat A^{(k+1)} = { \widehat V^{(k+1)}~}^H A \widehat V^{(k+1)} \\ {\bf endfor} \end{array} \end{equation*}

\begin{equation*} \begin{array}{l} \mbox{A simple QR algorithm} \\ \\ A^{(0)} = A \\ V^{(0)} = I \\ R^{(0)} = I \\ {\bf for~} k:=0, \ldots \\ ~~~( Q^{(k+1)}, R^{(k+1)} ) := {\rm QR}( A^{(k)} ) \\ ~~~ A^{(k+1)} = R^{(k+1)} Q^{(k+1)} \\ ~~~ V^{(k+1)} = V^{(k)} Q^{(k+1)} \\ {\bf endfor} \end{array} \end{equation*}

Let's start by showing how the QR algorithm applies unitary equivalence transformations to the matrices $A^{(k)} \text{:}$

Homework10.2.3.1.
Show that for the algorithm on the right $A^{(k+1)} = {Q^{(k+1)}~}^H A^{(k)} Q^{(k+1)} \text{.}$
Solution

The algorithm computes the QR factorization of $A^{(k)}$

\begin{equation*} A^{(k)} = Q^{(k+1)} R^{(k+1)} \end{equation*}

after which

\begin{equation*} A^{(k+1)} := R^{(k+1)} Q^{(k+1)} \end{equation*}

Hence

\begin{equation*} A^{(k+1)} = R^{(k+1)} Q^{(k+1)} = {Q^{(k+1)}~}^H A^{(k)} Q^{(k+1)}. \end{equation*}

We now prove these algorithms mathematically equivalent.

We will employ a proof by induction.

• Base case: $k = 0$

This is trivially true.

• Inductive step: Assume that $\widehat A^{(k)} = A^{(k)}\text{,}$ $\widehat R^{(k)} = R^{(k)}\text{,}$ $\widehat V^{(k)} = V^{(k)} \text{,}$ $V^{(k)} = Q^{(0)} \cdots Q^{(k)}\text{,}$ and $A^{k-1} = V^{(k)} R^{(k)} \cdots R^{(0)} \text{.}$ Show that $\widehat A^{(k+1)} = A^{(k+1)} \text{,}$ $\widehat R^{(k+1)} = R^{(k+1)}\text{,}$ $\widehat V^{(k+1)} = V^{(k+1)}\text{,}$ $V^{(k+1)} = Q^{(0)} \cdots Q^{(k+1)}\text{,}$ and $A^{k} = V^{(k+1)} R^{(k+1)} \cdots R^{(0)} \text{.}$

Since $\widehat A^{(k)}$ is computed as $\widehat A^{(k)} = { \widehat V^{(k)}~}^H A \widehat V^{(k)}$ by the algorithm on the left, the inductive hypothesis implies that

\begin{equation*} A^{(k)} = \widehat A^{(k)} = {\widehat V^{(k)}~}^H A \widehat V^{(k)} = {V^{(k)}}^H A V^{(k)}. \end{equation*}

Now, in the algorithm on the left,

$$A \widehat V^{(k)} = \widehat V^{(k+1)} \widehat R^{(k+1)}\label{chapter10-QR-eq-1}\tag{10.2.1}$$

and in the algorithm on the right

\begin{equation*} A^{(k)}= Q^{(k+1)} R^{(k+1)}. \end{equation*}

Substituting ${V^{(k)}}^H A V^{(k)}$ for $A^{(k)}$ yields

\begin{equation*} {V^{(k)}}^H A V^{(k)} = Q^{(k+1)} R^{(k+1)}. \end{equation*}

and hence

$$A V^{(k)} = V^{(k)} Q^{(k+1)} R^{(k+1)}.\label{chapter10-QR-eq-2}\tag{10.2.2}$$

Comparing (10.2.1) and (10.2.2), and remembering that a QR factorization is unique if its diagonal elements are positive, we conclude that

\begin{equation*} \widehat V^{(k+1)} = \begin{array}[t]{c} \underbrace{ V^{(k)} Q^{(k+1)}} \\ V^{(k+1)} \end{array} \quad \mbox{and} \quad \widehat R^{(k+1)} = R^{(k+1)}. \end{equation*}
• By the Principle of Mathematical Induction, the result holds for all $k \text{.}$

Homework10.2.3.2.

In the above algorithms, show that for all $k$

• $V^{(k)} = Q^{(0)} Q^{(1)} \cdots Q^{(k)}\text{.}$

• $A^{k-1} = V^{(k)} R^{(k)} \cdots R^{(1)} R^{(0)} \text{.}$

Solution

We will employ a proof by induction.

• Base case: $k = 0$

This is trivially true.

• Inductive step: Assume that $V^{(k)} = Q^{(0)} \cdots Q^{(k)}$ and $A^{k-1} = V^{(k)} R^{(k)} \cdots R^{(0)} \text{.}$ Show that $V^{(k+1)} = Q^{(0)} \cdots Q^{(k+1)}$ and $A^{k} = V^{(k+1)} R^{(k+1)} \cdots R^{(0)} \text{.}$

\begin{equation*} V^{(k+1)} = V^{(k)} Q^{(k+1)} = Q^{(0)} \cdots Q^{(k)} Q^{(k+1)}. \end{equation*}

by the inductive hypothesis.

Also,

\begin{equation*} \begin{array}{l} A^k \\ ~~~=~~~~ \lt \mbox{ definition } \gt \\ A A^{k-1} \\ ~~~=~~~~ \lt \mbox{ inductive hypothesis } \gt \\ A V^{(k)} R^{(k)} \cdots R^{(0)} \\ ~~~=~~~~ \lt \mbox{ inductive hypothesis } \gt \\ A \widehat V^{(k)} R^{(k)} \cdots R^{(0)} \\ ~~~=~~~~ \lt \mbox{ left algorithm } \gt \\ \widehat V^{(k+1)} \widehat R^{(k+1)} R^{(k)}\cdots R^{(0)} \\ ~~~=~~~~ \lt V^{(k+1)} = \widehat V^{(k+1)}; R^{(k+1)} = \widehat R^{(k+1)} \gt \\ V^{(k+1)} R^{(k+1)} R^{(k)}\cdots R^{(0)} . \end{array} \end{equation*}
• By the Principle of Mathematical Induction, the result holds for all $k \text{.}$

This last exercise also shows that

\begin{equation*} A^{k-1} = \begin{array}[t]{c} \underbrace{ Q^{(0)} Q^{(1)} \cdots Q^{(k)} } \\ \mbox{unitary} \end{array} \begin{array}[t]{c} \underbrace{ R^{(k)} \cdots R^{(1)} R^{(0)} } \\ \mbox{upper triangular} \end{array} \end{equation*}

which exposes a QR factorization of $A^{k-1} \text{.}$