QR algorithm

The QR algorithm is a numerical method in linear algebra designed to solve the complete eigenvalue problem, that is, to find all eigenvalues and eigenvectors of the matrix . It was developed in the late 1950s independently by V.N. Kublanovskaya and J. Francis.

Let A be a real matrix for which we want to find eigenvalues and vectors. Put A ₀ = A. At the kth step (starting from k = 0), we calculate the QR decomposition A _k = Q _k R _k , where Q _k is the orthogonal matrix (i.e., Q _k ^T = Q _k ⁻¹ ), and R _k is the upper triangular matrix . Then we define A _{k +1} = R _k Q _k .

notice, that

A_{k+1}=R_{k}Q_{k}=Q_{k}^{-1}Q_{k}R_{k}Q_{k}=Q_{k}^{-1}A_{k}Q_{k}=Q_{k}^{T}A_{k}Q_{k},

{\ displaystyle A_ {k + 1} = R_ {k} Q_ {k} = Q_ {k} ^ {- 1} Q_ {k} R_ {k} Q_ {k} = Q_ {k} ^ {- 1} A_ {k} Q_ {k} = Q_ {k} ^ {T} A_ {k} Q_ {k},}

{\ displaystyle A_ {k + 1} = R_ {k} Q_ {k} = Q_ {k} ^ {- 1} Q_ {k} R_ {k} Q_ {k} = Q_ {k} ^ {- 1} A_ {k} Q_ {k} = Q_ {k} ^ {T} A_ {k} Q_ {k},}

that is, all matrices A _k are similar , that is, their eigenvalues are equal.

Suppose that all diagonal minors of the matrix A are not degenerate . Then the sequence of matrices A _k as k → ∞ converges in shape to the cellular right triangular form corresponding to cells with identical eigenvalues modulo same. ^[one]

In order to obtain the eigenvectors of the matrix, we need to multiply all the matrices Q _k .

Aglorhythm is considered computationally stable , because produced by orthogonal similarity transformations.

Proof for a symmetric positive definite matrix

We assume that the eigenvalues of a positive definite matrix A are ordered in descending order:

\lambda _{1}>\lambda _{2}>...>\lambda _{n}>0.

{\ displaystyle \ lambda _ {1}> \ lambda _ {2}> ...> \ lambda _ {n}> 0.}

\lambda _{1}>\lambda _{2}>...>\lambda _{n}>0.

Let be

\Lambda =\mathrm {diag} \left(\lambda _{1},...,\lambda _{n}\right),

{\ displaystyle \ Lambda = \ mathrm {diag} \ left (\ lambda _ {1}, ..., \ lambda _ {n} \ right),}

\Lambda =\mathrm {diag} \left(\lambda _{1},...,\lambda _{n}\right),

and S is a matrix composed of eigenvectors of the matrix A. Then, matrix A can be written as spectral decomposition

A=S\Lambda S^{T}.

{\ displaystyle A = S \ Lambda S ^ {T}.}

A=S\Lambda S^{T}.

We find an expression for the powers of the original matrix in terms of the matrices Q _k and R _k . On the one hand, by the definition of the QR algorithm:

A^{k}=A_{1}^{k}=\left(Q_{1}R_{1}\right)^{k}=Q_{1}\left(R_{1}Q_{1}\right)^{k-1}R_{1}=Q_{1}A_{2}^{k-1}R_{1}.

{\ displaystyle A ^ {k} = A_ {1} ^ {k} = \ left (Q_ {1} R_ {1} \ right) ^ {k} = Q_ {1} \ left (R_ {1} Q_ { 1} \ right) ^ {k-1} R_ {1} = Q_ {1} A_ {2} ^ {k-1} R_ {1}.}

Applying this relation recursively, we obtain:

A^{k}=Q_{1}\cdot ...\cdot Q_{k}\cdot R_{k}\cdot ...\cdot R_{1}

{\ displaystyle A ^ {k} = Q_ {1} \ cdot ... \ cdot Q_ {k} \ cdot R_ {k} \ cdot ... \ cdot R_ {1}}

By entering the following notation:

S_{k}=Q_{1}\cdot ...\cdot Q_{k},

{\ displaystyle S_ {k} = Q_ {1} \ cdot ... \ cdot Q_ {k},}

T_{k}=R_{k}\cdot ...\cdot R_{1},

{\ displaystyle T_ {k} = R_ {k} \ cdot ... \ cdot R_ {1},}

we get

A^{k}=S_{k}T_{k}.

{\ displaystyle A ^ {k} = S_ {k} T_ {k}.}

On the other hand:

A^{k}=S\Lambda ^{k}S^{T}.

{\ displaystyle A ^ {k} = S \ Lambda ^ {k} S ^ {T}.}

Equating the right sides of the last two formulas, we obtain:

S\Lambda ^{k}S^{T}=S_{k}T_{k}.

{\ displaystyle S \ Lambda ^ {k} S ^ {T} = S_ {k} T_ {k}.}

Suppose that there is an LU decomposition of the matrix S ^T :

S^{T}=LU,

{\ displaystyle S ^ {T} = LU,}

then

S\Lambda ^{k}LU=S_{k}T_{k}.

{\ displaystyle S \ Lambda ^ {k} LU = S_ {k} T_ {k}.}

We multiply on the right by the matrix inverse to U , and then by the inverse to Λ ^k :

S\Lambda ^{k}L=S_{k}T_{k}U^{-1},

{\ displaystyle S \ Lambda ^ {k} L = S_ {k} T_ {k} U ^ {- 1},}

S\Lambda ^{k}L\Lambda ^{-k}=S_{k}T_{k}U^{-1}\Lambda ^{-k}.

{\ displaystyle S \ Lambda ^ {k} L \ Lambda ^ {- k} = S_ {k} T_ {k} U ^ {- 1} \ Lambda ^ {- k}.}

It can be shown that

\Lambda ^{k}L\Lambda ^{-k}\to \mathrm {diag} \left(l_{11},...,l_{nn}\right)=L^{\prime }

{\ displaystyle \ Lambda ^ {k} L \ Lambda ^ {- k} \ to \ mathrm {diag} \ left (l_ {11}, ..., l_ {nn} \ right) = L ^ {\ prime} }

at $k\to \infty ,$ ${\ displaystyle k \ to \ infty,}$ without loss of generality, we can assume that there are units on the diagonal of the matrix L ; therefore,

S_{k}T_{k}U^{-1}\Lambda ^{-k}\to S.

{\ displaystyle S_ {k} T_ {k} U ^ {- 1} \ Lambda ^ {- k} \ to S.}

We denote

P_{k}=T_{k}U^{-1}\Lambda ^{-k},

{\ displaystyle P_ {k} = T_ {k} U ^ {- 1} \ Lambda ^ {- k},}

moreover, the matrix P _k is upper triangular, as the product of upper triangular and diagonal matrices.

Thus, we have proved that

S_{k}P_{k}\to S

{\ displaystyle S_ {k} P_ {k} \ to S}

.

It follows from the uniqueness of the QR decomposition that if the product of the orthogonal and triangular matrix converges to the orthogonal matrix, then the triangular matrix converges to the identity matrix . From the foregoing it follows that

S_{k}=Q_{1}\cdot ...\cdot Q_{k}\to S.

{\ displaystyle S_ {k} = Q_ {1} \ cdot ... \ cdot Q_ {k} \ to S.}

That is, the matrices S _k converge to the eigenvector matrix of the matrix A.

Because

A_{k+1}=Q_{k}^{T}A_{k}Q_{k}=...=\left(Q_{k}^{T}\cdot ...\cdot Q_{1}^{T}\right)A_{1}\left(Q_{1}\cdot ...\cdot Q_{k}\right)=\left(Q_{1}\cdot ...\cdot Q_{k}\right)^{T}A\left(Q_{1}\cdot ...\cdot Q_{k}\right),

{\ displaystyle A_ {k + 1} = Q_ {k} ^ {T} A_ {k} Q_ {k} = ... = \ left (Q_ {k} ^ {T} \ cdot ... \ cdot Q_ {1} ^ {T} \ right) A_ {1} \ left (Q_ {1} \ cdot ... \ cdot Q_ {k} \ right) = \ left (Q_ {1} \ cdot ... \ cdot Q_ {k} \ right) ^ {T} A \ left (Q_ {1} \ cdot ... \ cdot Q_ {k} \ right),}

then

A_{k+1}=S_{k}^{T}AS_{k}.

{\ displaystyle A_ {k + 1} = S_ {k} ^ {T} AS_ {k}.}

Passing to the limit, we get:

\lim _{k\to \infty }A_{k}=\lim _{k\to \infty }A_{k+1}=S^{T}AS=S^{T}S\Lambda S^{T}S=\Lambda .

{\ displaystyle \ lim _ {k \ to \ infty} A_ {k} = \ lim _ {k \ to \ infty} A_ {k + 1} = S ^ {T} AS = S ^ {T} S \ Lambda S ^ {T} S = \ Lambda.}

So, we have proved that the QR algorithm allows us to solve the complete eigenvalue problem for a symmetric positive definite matrix.

QR algorithm implementation

Under certain conditions, a sequence of matrices $A_{k}$ ${\ displaystyle A_ {k}}$ converges to a triangular matrix, Schur decomposition of the matrix $A$ ${\ displaystyle A}$ . In this case, the eigenvalues of the triangular matrix are located on its diagonal, and the problem of finding the eigenvalues is considered solved. In convergence tests, it is not practical to require exact zeros in the zero part of the matrix, but you can use sets the error limits.

In the initial state of the matrix (without additional transformations), the cost of iterations is relatively high. The cost of the algorithm can be reduced by first bringing the matrix $A$ ${\ displaystyle A}$ to the form of the upper Hessenberg matrix (the cost of which is estimated as ${10 \over 3}n^{3}+O(n^{2})$ ${\ displaystyle {10 \ over 3} n ^ {3} + O (n ^ {2})}$ arithmetic operations using a method based on the Householder transform ), and using a finite sequence of orthogonal similarity transformations. This algorithm is somewhat similar to a two-way QR decomposition. (In a conventional QR decomposition, the Householder reflection matrix is multiplied only on the left, when using the Hessenberg form the reflection matrix is multiplied on both the left and right.) Finding the QR decomposition of the upper Hessenberg matrix is estimated as $6n^{2}+O(n)$ ${\ displaystyle 6n ^ {2} + O (n)}$ arithmetic operations. Due to the fact that the shape of the Hessenberg matrix is almost upper triangular (it has only one sub-diagonal element that is not equal to zero), it is possible to immediately reduce the number of iterations required for converging the QR algorithm.

If the original matrix is symmetric, the upper Hessenberg matrix is also symmetric and therefore is three-diagonal. The whole sequence of matrices has the same property. $A_{k}$ ${\ displaystyle A_ {k}}$ . In this case, the cost of the procedure is estimated as ${4 \over 3}n^{3}+O(n^{2})$ ${\ displaystyle {4 \ over 3} n ^ {3} + O (n ^ {2})}$ arithmetic operations using the Householder Transformation Method. Finding a QR decomposition of a symmetric tridiagonal matrix is estimated as $O(n)$ ${\ displaystyle O (n)}$ operations.

The rate of convergence depends on the degree of separation of the eigenvalues, and in practical implementation will be used "shifts", explicitly or implicitly, to enhance the separation of eigenvalues and to accelerate convergence. In a typical form for symmetric matrices, the QR algorithm accurately finds one eigenvalue (decreasing the dimension of the matrix) in one or two iterations, making this approach both effective and reliable.

Implicit QR code implementation

In modern computational practice, the QR algorithm is performed using its implicit version, which makes it easier to add multiple “shifts”. Initially, the matrix is reduced to the form of the upper Hessenberg matrix $A_{0}=QAQ^{T}$ ${\ displaystyle A_ {0} = QAQ ^ {T}}$ as well as in the explicit version. Then, at every step, the first column $A_{k}$ ${\ displaystyle A_ {k}}$ converted through a small-sized transformation of Householder similarity to the first column $p(A_{k})$ ${\ displaystyle p (A_ {k})}$ (or $p(A_{k})e_{1}$ ${\ displaystyle p (A_ {k}) e_ {1}}$ ), where $p(A_{k})$ ${\ displaystyle p (A_ {k})}$ is a polynomial of degree $r$ ${\ displaystyle r}$ which defines a strategy for “shifts” (usually $p(x)=(x-\lambda )(x-{\bar {\lambda }})$ ${\ displaystyle p (x) = (x- \ lambda) (x - {\ bar {\ lambda}})}$ where $\lambda$ ${\ displaystyle \ lambda}$ and ${\bar {\lambda }}$ ${\ displaystyle {\ bar {\ lambda}}}$ these are two eigenvalues of the residual submatrix $A_{k}$ ${\ displaystyle A_ {k}}$ size 2 x 2, this is the so-called implicit double shift). Then consecutive Householder transformations of dimension $r+1$ ${\ displaystyle r + 1}$ are made in order to return the working matrix $A_{k}$ ${\ displaystyle A_ {k}}$ to the form of the upper Hessenberg matrix.

Notes

↑ Numerical methods / N. S. Bakhvalov, N. P. Zhidkov, G. M. Kobelkov. - 3rd ed. - M: BINOM, Laboratory of knowledge, 2004. - S. 321. - 636 p. - ISBN 5-94774-175-X .

Links

Notes by Peter Olver on Orthogonal Bases and the QR Algorithm

[1] Numerical methods / N. S. Bakhvalov, N. P. Zhidkov, G. M. Kobelkov. - 3rd ed. - M: BINOM, Laboratory of knowledge, 2004. - S. 321. - 636 p. - ISBN 5-94774-175-X .