QR algorithm

The QR algorithm is a numerical method for calculating all eigenvalues and possibly the eigenvectors of a square matrix . The QR method or QR iteration, also known as the QR method, is based on the QR decomposition and was introduced in 1961–1962 independently of one another by John GF Francis and Wera Nikolajewna Kublanowskaja . A forerunner was the LR algorithm by Heinz Rutishauser (1958), but it is less stable and is based on the LR decomposition . Often the iterates from the QR algorithm converge on the Schur formthe matrix. The original method is quite complex and therefore impractical - even on today's computers - for matrices with hundreds of thousands of rows and columns.

Derived variants such as the Multishift method by Z. Bai and James Demmel 1989 and the numerically more stable variant by K. Braman, R. Byers and R. Mathias 2002 have practical running times that are cubic in the size of the matrix. The latter procedure is implemented in the numerical software library LAPACK , which in turn is used in many computer algebra systems (CAS) for the numerical matrix algorithms.

Description of the QR algorithm

Goal of the calculation process

Starting from a real or complex matrix or a sequence of orthogonal or unitary matrices is determined, so that the recursive matrix sequence largely converges to an upper triangular matrix. Since all transformations in the recursion are similarity transformations, all matrices of the matrix sequence have the same eigenvalues ​​with the same multiplicities. ${\displaystyle A\in \mathbb {R} ^{n\times n}}$ ${\displaystyle A\in \mathbb {C} ^{n\times n}}$ ${\displaystyle Q_{1},Q_{2},\dots ,}$ ${\displaystyle A_{1}=A,\,A_{2}=Q_{1}^{-1}A_{1}Q_{1},\dots ,A_{k+1}=Q_{k}^{-1}A_{k}Q_{k},\dots }$ ${\displaystyle A_{1},A_{2},\dots ,}$

If an upper triangular matrix is ​​reached in the limit value, a Schur form of , the eigenvalues ​​can be read from the diagonal elements. In the case of a real matrix with orthogonal transformations, however, there can also be complex eigenvalues. Then the result of the method is an upper block triangle matrix, the diagonal blocks of which have the size for the real eigenvalues ​​or the size for complex eigenvalues. The diagonal block of the corresponding rotational extension corresponds to an eigenvalue and its complex conjugate eigenvalue . ${\displaystyle A}$ ${\displaystyle 1}$ ${\displaystyle 2}$ ${\displaystyle \lambda =a+ib}$ ${\displaystyle \left({\begin{smallmatrix}a&-b\\b&a\end{smallmatrix}}\right)}$

General scheme of the procedure

Starting from a matrix (or ), a sequence of matrices is determined according to the following rule: ${\displaystyle A_{1}=A\in \mathbb {R} ^{n\times n}}$ ${\displaystyle A\in \mathbb {C} ^{n\times n}}$ ${\displaystyle A_{i}}$

1. for do ${\displaystyle i\in \mathbb {N} }$
2. Find a polynomial and evaluate it in the matrix . ${\displaystyle p_{i}}$ ${\displaystyle A_{i}}$
3. Compute the QR decomposition of . ${\displaystyle p_{i}(A_{i})=Q_{i}R_{i}}$
4. Calculate the new iterate: . ${\displaystyle A_{i+1}=Q_{i}^{-1}A_{i}Q_{i}}$
5. end for

The variants with (implicit) shifts and multiple shifts can also be displayed and analyzed in this general form.

The matrix function is often a polynomial (here a linear combination of matrix powers) with real (or complex) coefficients. In the simplest case, a linear polynomial is chosen, which then has the shape . In this case, the algorithm is simplified to the classic version (with single shift): ${\displaystyle p_{i}}$ ${\displaystyle p_{i}(A_{i})=A_{i}-\kappa _{i}I}$

1. for do ${\displaystyle i\in \mathbb {N} }$
2. Find an appropriate number ${\displaystyle \kappa _{i}}$
3. Compute the QR decomposition of ${\displaystyle A_{i}-\kappa _{i}I=Q_{i}R_{i}}$
4. Compute the new iterate: ${\displaystyle A_{i+1}=Q_{i}^{-1}A_{i}Q_{i}=Q_{i}^{-1}(Q_{i}R_{i}+\kappa _{i}I)Q_{i}=R_{i}Q_{i}+\kappa _{i}I}$
5. end for

In each iteration is set, the method is consistent with the on subspaces (here the full vector space) extended power method match. ${\displaystyle \kappa _{i}=0}$

The polynomial controlling the QR decomposition can be fixed from the start or dynamically adapted to the current matrix in each iteration step . There are also attempts to use for rational matrix functions, i. H. Functions of the form with polynomials and . ${\displaystyle A_{i}}$ ${\displaystyle p_{i}}$ ${\displaystyle p(A)=h(A)^{-1}g(A)}$ ${\displaystyle g}$ ${\displaystyle h}$

deflation

If it turns out in an iteration step that a lower left block of the columns and their rows in the amounts of all of its components falls below a predetermined threshold, the process is continued on the two diagonal square sub-blocks of the rows / columns and separately. If the matrices resulting from the division are small enough, the determination of the eigenvalues ​​can e.g. B. be terminated by calculating the characteristic polynomial and its zeros. ${\displaystyle 1,\dots ,i}$ ${\displaystyle i+1,\dots ,n}$ ${\displaystyle 1,\dots ,i}$ ${\displaystyle i+1,\dots ,n}$

Transformation to Hessenberg form

The goal of the QR algorithm is to create an upper triangle shape, or block triangle shape, in which blocks of the size correspond to complex eigenvalues. This can be done almost in a simple manner by a similarity transformation to Hessenberg form , i. H. on a matrix with only one (not identically disappearing) lower secondary diagonal. ${\displaystyle 2}$

• Set ${\displaystyle {\tilde {A}}=A}$
• for k = 1 to n-1 do

1. Form the column segment ${\displaystyle u={\tilde {A}}_{k+1:n,\,k}}$
2. Find the householder reflection that maps to the first canonical basis vector, ${\displaystyle {\tilde {P}}_{k}=I-2vv^{T}}$ ${\displaystyle u}$
3. Use the block matrix to replace with the similar matrix . ${\displaystyle P_{k}=\left({\begin{smallmatrix}I_{k}&0\\0&{\tilde {P}}_{k}\end{smallmatrix}}\right)}$ ${\displaystyle {\tilde {A}}}$ ${\displaystyle P_{k}^{-1}{\tilde {A}}P_{k}=P_{k}{\tilde {A}}P_{k}}$

• Note the overall transformation matrix . ${\displaystyle Q_{0}=P_{1}P_{2}\cdots P_{n-1}}$
• ${\displaystyle A_{1}={\tilde {A}}=Q_{0}^{-1}AQ_{0}}$ is now in Hessenberg form.

The Hessenberg form makes it easier to determine the characteristic polynomials of partial matrices . The Hessenberg form of a symmetrical matrix is a tridiagonal form , which considerably speeds up further calculations.

Furthermore, the Hessenberg form is obtained in every step of the QR algorithm. The basis for this is the arithmetic of generalized triangular matrices in which all entries disappear below a lower secondary diagonal. A Hessenberg matrix is ​​therefore a generalized triangular matrix with a secondary diagonal. With multiplication, the numbers of non-vanishing secondary diagonals add up, with addition, the larger number is usually retained.

Therefore, like the orthogonal matrix, have the number of lower secondary diagonals. Well, because of , too ${\displaystyle p_{i}(A_{i})}$ ${\displaystyle Q_{i}=p_{i}(A_{i})R_{i}^{-1}}$ ${\displaystyle m=\deg p_{i}}$ ${\displaystyle A_{i+1}=Q_{i}^{-1}A_{i}Q_{i}}$

${\displaystyle p_{i}(A_{i+1})=Q_{i}^{-1}p_{i}(A_{i})Q_{i}=R_{i}Q_{i}}$,

and the latter product also has secondary diagonals. This is generally only possible if it has exactly one secondary diagonal, i.e. is again in the Hessenberg shape. From the decomposition of into linear factors it follows (see below) that this is always the case. ${\displaystyle m}$ ${\displaystyle A_{i+1}}$ ${\displaystyle P-i(x)}$

Building on this, one can show (Francis 1962 after Bai / Demmel) that the first column of , which is also proportional to the first column of , completely determines the following matrix . More precisely, an orthogonal matrix whose first column is proportional to is created by bringing the transformed matrix back into Hessenberg form. Since only the components of the rows differ from zero in, a modification of the identity matrix can be in the upper left -block with a . ${\displaystyle x}$ ${\displaystyle p_{i}(A_{i})}$ ${\displaystyle Q_{i}}$ ${\displaystyle A_{i+1}=Q_{i}^{-1}A_{i}Q_{i}}$ ${\displaystyle U}$ ${\displaystyle x}$ ${\displaystyle A_{i+1}}$ ${\displaystyle U^{-1}A_{i}U}$ ${\displaystyle x}$ ${\displaystyle 1,\ldots ,m+1}$ ${\displaystyle U}$ ${\displaystyle s\times s}$ ${\displaystyle s>m}$

Variants of the QR algorithm

Simple QR iteration

It is elected. The method can then be given briefly as QR decomposition followed by multiplying together in the reverse order. This procedure is the direct generalization of the simultaneous power method for determining the first eigenvalues ​​of a matrix. This connection is derived from the subspace iteration . A simultaneous inverse power method is also performed indirectly. ${\displaystyle p_{i}(A)=A}$ ${\displaystyle A_{i}=Q_{i}R_{i}}$ ${\displaystyle A_{i+1}=R_{i}Q_{i}}$ ${\displaystyle k}$

QR algorithm with simple shifts

It is set. This means that the method can alternatively also be displayed as a QR decomposition followed by multiplying together corrected for the shift . Usually an attempt is made to approximate the smallest eigenvalue with the shift . The last diagonal element can be selected for this. The simple QR iteration results from setting all shifts to zero. ${\displaystyle p_{i}(A)=A-\kappa _{i}I}$ ${\displaystyle A_{i}-\kappa _{i}I=Q_{i}R_{i}}$ ${\displaystyle A_{i+1}=R_{i}Q_{i}+\kappa _{i}I}$ ${\displaystyle \kappa _{i}}$ ${\displaystyle \kappa _{i}=(A_{i})_{n,n}}$

Has so must Hessenberg form as a product of a matrix and a matrix without secondary diagonals also have Hessenberg form. The same therefore also applies to . So if in preparation of the QR algorithm in Hessenberg form, this is retained throughout the algorithm. ${\displaystyle A_{i}}$ ${\displaystyle Q_{i}=(A_{i}-\kappa _{i}I)R_{i}^{-1}}$ ${\displaystyle A_{i+1}=R_{i}Q_{i}+\kappa _{i}I}$ ${\displaystyle A_{1}=Q_{0}^{-1}AQ_{0}}$

Simple shifts for symmetric matrices

A symmetric real matrix has only real eigenvalues. The symmetry is preserved in all of them during the QR algorithm . For symmetric matrices, Wilkinson (1965) suggested that the eigenvalue of the submatrix on the right ${\displaystyle A_{1}=A}$ ${\displaystyle A_{i}}$ ${\displaystyle 2\times 2}$

${\displaystyle {\begin{pmatrix}a_{n-1,n-1}&a_{n-1,n}\\a_{n,n-1}&a_{n,n}\end{pmatrix}}}$

to choose as a shift that is closer to . Wilkinson showed that the matrix sequence determined in this way converges to a diagonal matrix whose diagonal elements are the eigenvalues ​​of . The speed of convergence is quadratic. ${\displaystyle a_{n,n}}$ ${\displaystyle (A_{i})_{i\in \mathbb {N} }}$ ${\displaystyle A=A_{1}}$

QR algorithm with double shifts

A pair of simple shifts can be combined in one iteration step. As a consequence, this means that complex shifts can be dispensed with for real matrices. In the above notation is

${\displaystyle (A_{1}-\kappa _{2}I)(A_{1}-\kappa _{1}I)=(A_{1}-\kappa _{2}I)Q_{1}R_{1}=Q_{1}(A_{2}-\kappa _{2}I)R_{1}=Q_{1}Q_{2}R_{2}R_{1}}$

a QR decomposition for the quadratic polynomial , evaluated in . The coefficients of this polynomial are real even for a conjugate pair of complex shifts. Thus, complex eigenvalues ​​of real matrices can also be approximated without complex numbers appearing in the calculation. ${\displaystyle p_{1}(x)=x^{2}-(\kappa _{1}+\kappa _{2})x+\kappa _{1}\kappa _{2}}$ ${\displaystyle A_{1}}$

A common choice for this double shift consists of the eigenvalues ​​of the sub-matrix at the bottom right , i.e. H. the quadratic polynomial is the characteristic of this block, ${\displaystyle 2\times 2}$

${\displaystyle p_{i}(x)=x^{2}-(a_{n-1,n-1}+a_{n,n})x+a_{n-1,n-1}a_{n,n}-a_{n-1,n}a_{n,n-1}}$.

QR algorithm with multiple shifts

A number is specified that is larger but significantly smaller than the size of the matrix . The polynomial can be chosen as the characteristic polynomial of the bottom right quadratic sub-matrix of the current matrix . Another strategy is the eigenvalues of the very lowest to determine -Teilmatrix and the amount smallest eigenvalues select among them. With these, a QR decomposition of ${\displaystyle m}$ ${\displaystyle 2}$ ${\displaystyle n}$ ${\displaystyle A}$ ${\displaystyle p_{i}(x)}$ ${\displaystyle m\times m}$ ${\displaystyle A_{i}}$ ${\displaystyle s>m}$ ${\displaystyle s\times s}$ ${\displaystyle m}$ ${\displaystyle s_{1},\dots ,s_{m}}$

${\displaystyle p_{i}(A_{i})=(A_{i}-s_{1}I)(A_{i}-s_{2}I)\cdots (A_{i}-s_{m}I)=Q_{i}R_{i}}$ and ${\displaystyle A_{i+1}=Q_{i}^{-1}A_{i}Q_{i}}$

certainly.

With a multishift method, it is often achieved that the components of the lower left -block become small particularly quickly in the sequence of the iterated matrices, thus splitting the eigenvalue problem. ${\displaystyle (n-m+1,\dots ,n)\times (1,\dots ,n-m)}$

Implicit multishift iteration

Combining multiple shifts in the general form is very time-consuming. As mentioned above, the effort can be reduced by producing the Hessenberg shape in a preparatory step . Since each multishift step can be composed of individual (even complex) shifts, the Hessenberg form is retained throughout the entire algorithm. ${\displaystyle A_{1}=Q_{0}^{-1}AQ_{0}}$

This allows the QR algorithm to be converted into a “bulge-chasing” algorithm that creates a dent in the Hessenberg shape at the upper end of the diagonal and then “chases” it down the diagonal and out of the matrix at the lower end.

1. for do ${\displaystyle j\in \mathbb {N} _{0}}$
2. Calculate the polynomial according to one of the given variants, ${\displaystyle p(x)}$
3. Find the vector . ${\displaystyle x=p(A_{j})e_{1}}$
4. Find a reflection from onto the first unit vector. Since only the first components do not disappear, this reflection has a simple block shape. ${\displaystyle x}$ ${\displaystyle x}$ ${\displaystyle m+1}$
5. Form the matrix and transform it so that it is in Hessenberg form again. ${\displaystyle {\tilde {A}}=U^{-1}A_{j}U}$ ${\displaystyle A_{j+1}=Q^{-1}U^{-1}A_{j}UQ}$
6. end for

If householder reflections are put together , they have a block diagonal shape . ${\displaystyle Q}$ ${\displaystyle Q=P_{1}P_{2}\cdots P_{n-1}}$ ${\displaystyle P_{k}=\mathrm {diag} (I_{k},{\tilde {P}}_{k},I_{\max(0,n-k-m)})}$

Notes on functionality

Similarity transformations

The matrices calculated in the QR algorithm are unitarily similar to one another, since due to ${\displaystyle A_{i}-\kappa _{i}I=Q_{i}R_{i}\quad \Leftrightarrow \quad R_{i}=Q_{i}^{H}(A_{i}-\kappa _{i}I)}$

${\displaystyle A_{i+1}=R_{i}Q_{i}+\kappa _{i}I=Q_{i}^{H}(A_{i}-\kappa _{i}I)Q_{i}+\kappa _{i}I=Q_{i}^{H}A_{i}Q_{i}-\kappa _{i}Q_{i}^{H}Q_{i}+\kappa _{i}I=Q_{i}^{H}A_{i}Q_{i}}$

applies. Thus all matrices have the same eigenvalues ​​(counted with the algebraic and geometric multiplicity). ${\displaystyle A_{i}}$

Choice of shifts, convergence

A simple but not very effective choice is to choose shifts equal to zero. The iterates of the resulting algorithm, the QR algorithm in its basic form, partially converge, if all eigenvalues ​​differ in their absolute value, to an upper triangular matrix with the eigenvalues ​​on the diagonal. Partial convergence here means that the elements of the lower triangle go from towards zero and the diagonal elements towards the eigenvalues. So nothing is said about the elements in the upper triangle . ${\displaystyle A_{i}}$ ${\displaystyle A_{i}}$

If the shifts are chosen as the matrix element at the bottom right of the current iterate, that is , the algorithm converges squarely or even cubically under suitable circumstances. This shift is known as the Rayleigh quotient shift because it is related to a Rayleigh quotient via the inverse iteration . ${\displaystyle \kappa _{i}=a_{nn}^{(i)}}$

When calculating in the real ( ) one would like to calculate the real Schur form and still be able to work with conjugate complex eigenvalues. There are various shift strategies for this. ${\displaystyle A\in \mathbb {R} ^{n\times n}}$

An extension of simple shifts is the Wilkinson shift named after James Hardy Wilkinson . For the Wilkinson Shift, the eigenvalue closer to the last matrix element becomes the last main subsection matrix ${\displaystyle 2\times 2}$

${\displaystyle {\begin{pmatrix}a_{n-1,n-1}^{(i)}&a_{n-1,n}^{(i)}\\a_{n,n-1}^{(i)}&a_{n,n}^{(i)}\end{pmatrix}}}$

used.

The QR algorithm as an extension of the power method

For the analysis of the algorithm, the additional matrix effects of the products are accumulated and , defined. The products of orthogonal or unitary matrices are again orthogonal or unitary matrices, and the products of right-upper triangular matrices are again right-upper triangular matrices. The matrices of the QR iteration result from similarity transformations , because ${\displaystyle {\tilde {Q}}_{k}=Q_{1}Q_{2}\cdots Q_{k}}$ ${\displaystyle {\tilde {R}}_{k}=R_{k}R_{k-1}\cdots R_{2}R_{1}}$ ${\displaystyle k\in \mathbb {N} }$ ${\displaystyle {\tilde {Q}}_{k}}$ ${\displaystyle {\tilde {R}}_{k}}$ ${\displaystyle A}$

${\displaystyle A_{1}=A,\;A_{2}=R_{1}Q_{1}=Q_{1}^{-1}A_{1}Q_{1},\;A_{3}=Q_{2}^{-1}A_{2}Q_{2}={\tilde {Q}}_{2}^{-1}A{\tilde {Q}}_{2},\dots ,\;A_{k}={\tilde {Q}}_{k}^{-1}A{\tilde {Q}}_{k},\dots }$.

From this one deduces QR-decomposition of the powers of : ${\displaystyle A}$

${\displaystyle {\begin{aligned}A=&&&Q_{1}R_{1}&=&{\tilde {Q}}_{1}{\tilde {R}}_{1}\\A^{2}=&AQ_{1}R_{1}=Q_{1}A_{2}R_{1}&=&Q_{1}Q_{2}R_{2}R_{1}&=&{\tilde {Q}}_{2}{\tilde {R}}_{2}\\\vdots &\\A^{k}=&A{\tilde {Q}}_{k-1}{\tilde {R}}_{k-1}={\tilde {Q}}_{k-1}A_{k}{\tilde {R}}_{k-1}&=&{\tilde {Q}}_{k-1}Q_{k}R_{k}{\tilde {R}}_{k-1}&=&{\tilde {Q}}_{k}{\tilde {R}}_{k}\end{aligned}}}$

In the course of the algorithm, QR-decompositions of the powers of the matrix are determined implicitly . Due to the form of these decompositions, it applies that for each the first columns of the matrix span the same subspace as the first columns of the matrix ; in addition, the columns are orthonormal to each other. However, this is exactly the situation after the -th step of a simultaneous power iteration. The diagonal elements of are the approximations of the eigenvalues ​​of . Therefore the convergence properties of the power iteration can be carried over: ${\displaystyle A}$ ${\displaystyle j=1,2,\ldots ,n}$ ${\displaystyle j}$ ${\displaystyle A^{k}}$ ${\displaystyle j}$ ${\displaystyle {\tilde {Q}}_{k}}$ ${\displaystyle {\tilde {Q}}_{k}}$ ${\displaystyle k}$ ${\displaystyle R_{k}}$ ${\displaystyle A{\tilde {Q}}_{k-1}={\tilde {Q}}_{k}R_{k}}$ ${\displaystyle A}$

The simple QR algorithm only converges if all eigenvalues ​​differ from one another in their amounts. Are the eigenvalues ​​after

${\displaystyle |\lambda _{1}|<|\lambda _{2}|<\dots <|\lambda _{n}|}$

sorted, the speed of convergence is linear with a contraction factor that corresponds to the minimum of the quotients . ${\displaystyle {\tfrac {|\lambda _{k}|}{|\lambda _{k+1}|}}}$

For real matrices in particular, this algorithm can only converge if all eigenvalues ​​are real (since otherwise complex conjugate pairs would exist with the same amount). This situation is given for all symmetric matrices.

The QR algorithm as a simultaneous inverse power iteration

If it is invertible, it holds for the transpose (for complex matrices the Hermitian adjoint ) the inverse of and all their powers ${\displaystyle A}$ ${\displaystyle A}$

${\displaystyle (A^{-T})^{k}=(A^{k})^{-T}=({\tilde {R}}_{k}^{-1}{\tilde {Q}}_{k}^{-1})^{T}={\tilde {Q}}_{k}{\tilde {R}}_{k}^{-T}}$

The inverse of an upper right triangular matrix is ​​again an upper right triangular matrix, the transpose of which is a lower left triangular matrix. The QR algorithm thus also determines a QL decomposition of the powers of . From the form of this decomposition it is evident that for each the last columns of span the same subspace as the last columns of . In the last column of , an inverse power iteration for is carried out, so this column converges towards the dual eigenvector to the smallest eigenvalue of . In the product , the lower left component is the so-called Rayleigh quotient of the inverse power iteration. The convergence properties are analogous to those given above. ${\displaystyle A^{-T}}$ ${\displaystyle j=1,2,\ldots ,n}$ ${\displaystyle j}$ ${\displaystyle {\tilde {Q}}_{k}}$ ${\displaystyle j}$ ${\displaystyle (A^{-T})^{k}}$ ${\displaystyle {\tilde {Q}}_{k}}$ ${\displaystyle A^{T}}$ ${\displaystyle A}$ ${\displaystyle {\tilde {Q}}_{k}^{T}A{\tilde {Q}}_{k}=A_{k}}$ ${\displaystyle (A_{k})_{nn}}$

Web links

• LP: QR procedure for general EWP , Georg-August-Universität Göttingen

literature

• JGF Francis (1961): The QR Transformation: A Unitary Analogue to the LR Transformation — Part 1 . The Computer Journal Vol. 4 (3), pp. 265-271. doi : 10.1093 / comjnl / 4.3.265
• JGF Francis (1962): The QR Transformation — Part 2 . The Computer Journal 1962 4 (4): 332-345; ( online )
• David S. Watkins (1982): Understanding the QR algorithm , SIAM Review, Vol. 24, pp. 427-440 ( JSTOR )
• Z. Bai; J. Demmel (1989): On a block implementation of Hessenberg multishift QR iteration , International Journal of High Speed ​​Computing, Vol. 1 (1), pp. 97-112. (see LAPACK Working Notes )
• AA Dubrulle; GH Golub (1994): A multishift QR iteration without computation of the shifts . Numerical Algorithms, Vol 7, pp 173-181
• K. Braman; R. Byers; R. Mathias (2002): The Multishift QR Algorithm. Part I: Maintaining Well-Focused Shifts and Level 3 Performance (PDF; 224 kB). SIAM Journal on Matrix Analysis and Applications, Vol. 23, No. 4, pp. 929-947
• K. Braman; R. Byers; R. Mathias (2002): The Multishift QR Algorithm. Part II: Aggressive Early Deflation (PDF; 265 kB). SIAM Journal on Matrix Analysis and Applications, Vol. 23, No. 4, pp. 948-989
• David S. Watkins: The QR algorithm revisited (PDF; 417 kB), SIAM Review, Vol. 50, No. 1, pp. 133-145