Skip to main content

Subsection 10.3.6 The Francis implicit QR Step

Starting with a tridiagonal Hermitian matrix, the Francis QR Step combines the steps

\begin{equation*} \begin{array}{l} A^{(k-1)} - \mu_k I \rightarrow Q^{(k)} R^{(k)} \\ A^{(k+1)} \becomes R^{(k)} Q^{(k)} + \mu_k I \end{array} \end{equation*}

into a single step.

Consider the \(4 \times 4 \) tridiagonal matrix

\begin{equation*} \left( \begin{array}{c c c c} \alpha_{0,0} \amp \alpha_{0,1} \amp 0 \amp 0 \\ \alpha_{1,0} \amp \alpha_{1,1} \amp \alpha_{1,2} \amp 0 \\ 0 \amp \alpha_{2,1} \amp \alpha_{2,2} \amp \alpha_{2,3} \\ 0 \amp 0 \amp \alpha_{3,2} \amp \alpha_{3,3} \\ \end{array} \right) - \mu I \end{equation*}

The first Givens' rotation is computed from \(\left( \begin{array}{c} \alpha_{0,0} - \mu\\ \alpha_{1,0} \end{array} \right) \text{,}\) yielding \(\gamma_{1,0} \) and \(\sigma_{1,0} \) so that

\begin{equation*} \left( \begin{array}{c c} \gamma_{1,0} \amp -\sigma_{1,0} \\ \sigma_{1,0} \amp \gamma_{1,0} \end{array} \right)^T \left( \begin{array}{c} \alpha_{0,0} - \mu I \\ \alpha_{1,0} \end{array} \right) \end{equation*}

has a zero second entry. Now, to preserve eigenvalues, any orthogonal matrix that is applied from the left must also have its transpose applied from the right. Let us compute

\begin{equation*} \begin{array}{l} \left( \begin{array}{| c c | c c} \hline \tilde{\alpha}_{0,0} \amp \widehat{\alpha}_{1,0} \amp \widehat{\alpha}_{2,0} \amp 0 \\ \widehat{\alpha}_{1,0} \amp \widehat{\alpha}_{1,1} \amp \widehat{\alpha}_{1,2} \amp 0 \\\hline \widehat{\alpha}_{2,0} \amp \widehat{\alpha}_{2,1} \amp \alpha_{2,2} \amp \alpha_{2,3} \\ 0 \amp 0 \amp \alpha_{3,2} \amp \alpha_{3,3} \end{array} \right)\\ ~~~ = \left( \begin{array}{| c c | c c} \hline \gamma_{1,0} \amp \sigma_{1,0} \amp 0 \amp 0 \\ - \sigma_{1,0} \amp \gamma_{1,0} \amp 0 \amp 0 \\ \hline 0 \amp 0 \amp 1 \amp 0 \\ 0 \amp 0 \amp 0 \amp 1 \end{array} \right) \left( \begin{array}{| c c | c c} \hline \alpha_{0,0} \amp \alpha_{0,1} \amp 0 \amp 0 \\ \alpha_{1,0} \amp \alpha_{1,1} \amp \alpha_{1,2} \amp 0 \\ \hline 0 \amp \alpha_{2,1} \amp \alpha_{2,2} \amp \alpha_{2,3} \\ 0 \amp 0 \amp \alpha_{3,2} \amp \alpha_{3,3} \end{array} \right) \left( \begin{array}{| c c | c c} \hline \gamma_{1,0} \amp - \sigma_{1,0} \amp 0 \amp 0 \\ \sigma_{1,0} \amp \gamma_{1,0} \amp 0 \amp 0 \\ \hline 0 \amp 0 \amp 1 \amp 0 \\ 0 \amp 0 \amp 0 \amp 1 \end{array} \right). \end{array} \end{equation*}

This is known as "introducing the bulge."

Next, from \(\left( \begin{array}{c} \widehat{\alpha}_{1,0} \\ \widehat{\alpha}_{2,0} \end{array} \right)\) one can compute \(\gamma_{2,0} \) and \(\sigma_{2,0} \) so that

\begin{equation*} \left( \begin{array}{c c} \gamma_{2,0} \amp -\sigma_{2,0} \\ \sigma_{2,0} \amp \gamma_{2,0} \end{array} \right)^T \left( \begin{array}{c} \widehat{\alpha}_{1,0} \\ \widehat{\alpha}_{2,0} \end{array} \right) = \left( \begin{array}{c} \tilde{\alpha}_{1,0} \\ 0 \end{array} \right). \end{equation*}

Then

\begin{equation*} \begin{array}{l} \left( \begin{array}{| c c c | c } \hline \tilde{\alpha}_{0,0} \amp \tilde{\alpha}_{1,0} \amp 0 \amp 0 \\ \tilde{\alpha}_{1,0} \amp \tilde{\alpha}_{1,1} \amp \widehat{\widehat{\alpha}}_{2,1} \amp \widehat{\alpha}_{3,1} \\ 0 \amp \widehat{\widehat{\alpha}}_{2,1} \amp \widehat{\alpha}_{2,2} \amp\widehat{\alpha}_{2,3} \\ \hline 0 \amp \widehat{\alpha}_{3,1} \amp \widehat{\alpha}_{3,2} \amp {\alpha}_{3,3} \\ \end{array} \right) \\ ~~~= \left( \begin{array}{| c c c | c} \hline 1 \amp 0 \amp 0 \amp 0 \\ 0 \amp \gamma_{2,0} \amp \sigma_{2,0} \amp 0 \\ 0 \amp - \sigma_{2,0} \amp \gamma_{2,0} \amp 0 \\ \hline 0 \amp 0 \amp 0 \amp 1 \end{array} \right) \left( \begin{array}{| c c c | c } \hline \tilde{\alpha}_{0,0} \amp \widehat{\alpha}_{1,0} \amp \widehat{\alpha}_{2,0} \amp 0 \\ \widehat{\alpha}_{1,0} \amp \widehat{\alpha}_{1,1} \amp \widehat{\alpha}_{1,2} \amp 0 \\ \widehat{\alpha}_{2,0} \amp \widehat{\alpha}_{2,1} \amp \alpha_{2,2} \amp \alpha_{2,3} \\ \hline 0 \amp 0 \amp \alpha_{3,2} \amp \alpha_{3,3} \\ \end{array} \right) \left( \begin{array}{| c c c | c} \hline 1 \amp 0 \amp 0 \amp 0 \\ 0 \amp \gamma_{2,0} \amp - \sigma_{2,0} \amp 0 \\ 0 \amp \sigma_{2,0} \amp \gamma_{2,0} \amp 0 \\ \hline 0 \amp 0 \amp 0 \amp 1 \end{array} \right) \end{array} \end{equation*}

again preserves eigenvalues. Finally, from

\begin{equation*} \left( \begin{array}{c} \widehat{\alpha}_{2,1} \\ \widehat{\alpha}_{3,1} \end{array} \right) \end{equation*}

one can compute \(\gamma_{3,1} \) and \(\sigma_{3,1} \) so that

\begin{equation*} \begin{array}{l} \left( \begin{array}{c c} \gamma_{3,1} \amp -\sigma_{3,1} \\ \sigma_{3,1} \amp \gamma_{3,1} \end{array} \right)^T \left( \begin{array}{c} \widehat{\alpha}_{2,1} \\ \widehat{\alpha}_{3,1} \end{array} \right) \\ ~~~ = \left( \begin{array}{c} \tilde{\alpha}_{2,1} \\ 0 \end{array} \right) . \end{array} \end{equation*}

Then

\begin{equation*} \begin{array}{l} \left( \begin{array}{c | c c c |} \tilde{\alpha}_{0,0} \amp \tilde{\alpha}_{1,0} \amp 0 \amp 0 \\ \hline \tilde{\alpha}_{1,0} \amp \tilde{\alpha}_{1,1} \amp \tilde{\alpha}_{2,1} \amp 0 \\ 0 \amp \tilde{\alpha}_{2,1} \amp \tilde{\alpha}_{2,2} \amp \tilde{\alpha}_{2,3} \\ 0 \amp 0 \amp \tilde{\alpha}_{3,2} \amp \tilde{\alpha}_{3,3} \\ \hline \end{array} \right) \\ ~~~ = \left( \begin{array}{c | c c c |} 1 \amp 0 \amp 0 \amp 0 \\\hline 0 \amp 1 \amp 0 \amp 0 \\ 0 \amp 0 \amp \gamma_{3,2} \amp \sigma_{3,2} \\ 0 \amp 0 \amp - \sigma_{3,2} \amp \gamma_{3,2} \\ \hline \end{array} \right) \left( \begin{array}{c | c c c |} \tilde{\alpha}_{0,0} \amp \tilde{\alpha}_{1,0} \amp 0 \amp 0 \\ \hline \tilde{\alpha}_{1,0} \amp \tilde{\alpha}_{1,1} \amp \widehat{\widehat{\alpha}}_{2,1} \amp \widehat{\alpha}_{3,1} \\ 0 \amp \widehat{\widehat{\alpha}}_{2,1} \amp \widehat{\alpha}_{2,2} \amp\widehat{\alpha}_{2,3} \\ 0 \amp \widehat{\alpha}_{3,1} \amp \widehat{\alpha}_{3,2} \amp {\alpha}_{3,3} \\ \hline \end{array} \right) \left( \begin{array}{c | c c c |} 1 \amp 0 \amp 0 \amp 0 \\ \hline 0 \amp 1 \amp 0 \amp 0 \\ 0 \amp 0 \amp \gamma_{3,1} \amp - \sigma_{3,1} \\ 0 \amp 0 \amp \sigma_{3,1} \amp \gamma_{3,1} \\ \hline \end{array} \right) \end{array} \end{equation*}

The process of transforming the matrix that results from introducing the bulge (the nonzero element \(\widehat \alpha_{2,0}\)) back into a tridiagonal matrix is commonly referred to as "chasing the bulge."

The matrix \(Q \) is the orthogonal matrix that results from multiplying the different Givens' rotations together:

\begin{equation*} Q = \left( \begin{array}{| c c | c c} \hline \gamma_{1,0} \amp - \sigma_{1,0} \amp 0 \amp 0 \\ \sigma_{1,0} \amp \gamma_{1,0} \amp 0 \amp 0 \\ \hline 0 \amp 0 \amp 1 \amp 0 \\ 0 \amp 0 \amp 0 \amp 1 \end{array} \right) \left( \begin{array}{c | c c | c} 1 \amp 0 \amp 0 \amp 0 \\ \hline 0 \amp\gamma_{2,0} \amp - \sigma_{2,0} \amp 0 \\ 0 \amp \sigma_{2,0} \amp \gamma_{2,0} \amp 0 \\ \hline 0 \amp 0 \amp 0 \amp 1 \end{array} \right) \left( \begin{array}{c c | c c |} 1 \amp 0 \amp 0 \amp 0 \\ 0 \amp 1 \amp 0 \amp 0 \\ \hline 0 \amp 0 \amp \gamma_{3,1} \amp - \sigma_{3,1} \\ 0 \amp 0 \amp \sigma_{3,1} \amp \gamma_{3,1} \\ \hline \end{array} \right). \end{equation*}

Importantly, the first column of \(Q \) is given by

\begin{equation*} \left( \begin{array}{c} \gamma_{1,0} \\ \sigma_{1,0} \\ 0 \\ 0 \end{array} \right), \end{equation*}

is exactly the same first column had \(Q \) been computed as in Subsection 10.3.4 (10.3.1). Thus, by the Implicit Q Theorem, the tridiagonal matrix that results from this approach is equal to the tridiagonal matrix that would be computed by applying the QR factorization from Section~\ref{sec:QRTriDiag} with \(A - \mu I \text{,}\) \(A - \mu I \rightarrow QR \) followed by the formation of \(R Q + \mu I \) using the algorithm for computing \(R Q \) in Subsection 10.3.4.

Figure 10.3.6.1. Illustration of the "chasing of the bulge". (Note to self: grab a higher resolution version of this.)
\begin{equation*} % Define the operation to be computed \renewcommand{\routinename}{ T \becomes \mbox{ChaseBulge}( T )} % Step 3: Loop-guard \renewcommand{\guard}{ m( T_{BR} ) > 0 } % Step 8: Insert the updates required to change the % state from that given in Step 6 to that given in Step 7 % Note: The below needs editing!!! \renewcommand{\update}{ \begin{array}{l} \mbox{Compute } ( \gamma,\sigma ) \mbox{ s.t. } G_{\gamma, \sigma}^T t_{21} = \left( \begin{array}{c} \tau_{21} \\ 0 \end{array} \right), \mbox{ and assign } t_{21} = \left( \begin{array}{c} \tau_{21} \\ 0 \end{array} \right) \\ T_{22} = G_{\gamma, \sigma}^T T_{22} G_{\gamma, \sigma} \\[2mm] t_{32}^T = t_{32}^T G_{\gamma, \sigma} {\rm ~~~~~~(not~performed~during~final~step)}\\ \end{array} } % Step 4: Redefine Initialize \renewcommand{\partitionings}{ T \rightarrow \left( \begin{array}{ c | c | c } T_{TL} \amp \star \amp \star \\ \hline T_{ML} \amp T_{MM} \amp \star \\ \hline 0 \amp T_{BM} \amp T_{BR} \end{array} \right) } \renewcommand{\partitionsizes}{ T_{TL} \mbox{ is } 0 \times 0 \mbox{ and } T_{MM} \mbox{ is } 3 \times 3 } % Step 5a: Repartition the operands \renewcommand{\repartitionings}{ \left( \begin{array}{ c | c | c } T_{TL} \amp \star \amp 0 \\ \hline T_{ML} \amp T_{MM} \amp \star \\ \hline 0 \amp T_{BM} \amp T_{BR} \end{array} \right) \rightarrow \left( \begin{array}{ c | c c | c c } T_{00} \amp \star \amp 0 \amp 0 \amp 0 \\ \hline t_{10}^T \amp \tau_{11} \amp \star \amp 0 \amp 0 \\ 0 \amp t_{21} \amp T_{22} \amp \star \amp 0 \\ \hline 0 \amp 0 \amp t_{32}^T \amp \tau_{33} \amp \star \\ 0 \amp 0 \amp 0 \amp t_{43} \amp T_{44} \end{array} \right) } \renewcommand{\repartitionsizes}{ \tau_{11} \mbox{ and } \tau_{33} \mbox{ are scalars} \\ \mbox{(during final step,} \tau_{33} \mbox{ is } 0 \times 0 \mbox{)} } % Step 5b: Move the double lines \renewcommand{\moveboundaries}{ \left( \begin{array}{ c | c | c } T_{TL} \amp \star \amp 0 \\ \hline T_{ML} \amp T_{MM} \amp \star \\ \hline 0 \amp T_{BM} \amp T_{BR} \end{array} \right) \leftarrow \left( \begin{array}{ c c | c c | c } T_{00} \amp \star \amp 0 \amp 0 \amp 0 \\ t_{10}^T \amp \tau_{11} \amp \star \amp 0 \amp 0 \\ \hline 0 \amp t_{21} \amp T_{22} \amp \star \amp 0 \\ 0 \amp 0 \amp t_{32}^T \amp \tau_{33} \amp \star \\ \hline 0 \amp 0 \amp 0 \amp t_{43} \amp T_{44} \end{array} \right) } \FlaAlgorithm \end{equation*}
Figure 10.3.6.2. Algorithm for "chasing the bulge" that, given a tridiagonal matrix with an additional nonzero \(\alpha_{2,0}\) element, reduces the given matrix back to a tridiagonal matrix.

The successive elimination of elements \(\widehat{\alpha}_{i+1,i} \) is often referred to as {\em chasing the bulge} while the entire process that introduces the bulge and then chases it is known as a Francis Implicit QR Step. Obviously, the method generalizes to matrices of arbitrary size, as illustrated in Figure 10.3.6.1. An algorithm for the chasing of the bulge is given in Figure 10.3.6.2. (Note that in those figures \(T \) is used for \(A \text{,}\) something that needs to be made consistent in these notes, eventually.) In practice, the tridiagonal matrix is not stored as a matrix. Instead, its diagonal and subdiagonal are stored as vectors.