In
numerical linear algebra
Numerical linear algebra, sometimes called applied linear algebra, is the study of how matrix operations can be used to create computer algorithms which efficiently and accurately provide approximate answers to questions in continuous mathematics. ...
, the QR algorithm or QR iteration is an
eigenvalue algorithm
In numerical analysis, one of the most important problems is designing efficient and stable algorithms for finding the eigenvalues of a matrix. These eigenvalue algorithms may also find eigenvectors.
Eigenvalues and eigenvectors
Given an square ...
: that is, a procedure to calculate the
eigenvalues and eigenvectors
In linear algebra, an eigenvector () or characteristic vector of a linear transformation is a nonzero vector that changes at most by a scalar factor when that linear transformation is applied to it. The corresponding eigenvalue, often denoted b ...
of a
matrix
Matrix most commonly refers to:
* ''The Matrix'' (franchise), an American media franchise
** ''The Matrix'', a 1999 science-fiction action film
** "The Matrix", a fictional setting, a virtual reality environment, within ''The Matrix'' (franchis ...
. The QR algorithm was developed in the late 1950s by
John G. F. Francis and by
Vera N. Kublanovskaya, working independently. The basic idea is to perform a
QR decomposition
In linear algebra, a QR decomposition, also known as a QR factorization or QU factorization, is a decomposition of a matrix ''A'' into a product ''A'' = ''QR'' of an orthogonal matrix ''Q'' and an upper triangular matrix ''R''. QR decompo ...
, writing the matrix as a product of an
orthogonal matrix
In linear algebra, an orthogonal matrix, or orthonormal matrix, is a real square matrix whose columns and rows are orthonormal vectors.
One way to express this is
Q^\mathrm Q = Q Q^\mathrm = I,
where is the transpose of and is the identity ma ...
and an upper
triangular matrix
In mathematics, a triangular matrix is a special kind of square matrix. A square matrix is called if all the entries ''above'' the main diagonal are zero. Similarly, a square matrix is called if all the entries ''below'' the main diagonal are ...
, multiply the factors in the reverse order, and iterate.
The practical QR algorithm
Formally, let ''A'' be a real matrix of which we want to compute the eigenvalues, and let ''A''
0:=''A''. At the ''k''-th step (starting with ''k'' = 0), we compute the
QR decomposition
In linear algebra, a QR decomposition, also known as a QR factorization or QU factorization, is a decomposition of a matrix ''A'' into a product ''A'' = ''QR'' of an orthogonal matrix ''Q'' and an upper triangular matrix ''R''. QR decompo ...
''A''
''k''=''Q''
''k''''R''
''k'' where ''Q''
''k'' is an
orthogonal matrix
In linear algebra, an orthogonal matrix, or orthonormal matrix, is a real square matrix whose columns and rows are orthonormal vectors.
One way to express this is
Q^\mathrm Q = Q Q^\mathrm = I,
where is the transpose of and is the identity ma ...
(i.e., ''Q''
''T'' = ''Q''
−1) and ''R''
''k'' is an upper triangular matrix. We then form ''A''
''k''+1 = ''R''
''k''''Q''
''k''. Note that
:
so all the ''A''
''k'' are
similar and hence they have the same eigenvalues. The algorithm is
numerically stable
In the mathematical subfield of numerical analysis, numerical stability is a generally desirable property of numerical algorithms. The precise definition of stability depends on the context. One is numerical linear algebra and the other is algorit ...
because it proceeds by ''orthogonal'' similarity transforms.
Under certain conditions,
the matrices ''A''
''k'' converge to a triangular matrix, the
Schur form of ''A''. The eigenvalues of a triangular matrix are listed on the diagonal, and the eigenvalue problem is solved. In testing for convergence it is impractical to require exact zeros, but the
Gershgorin circle theorem
In mathematics, the Gershgorin circle theorem may be used to bound the spectrum of a square matrix. It was first published by the Soviet mathematician Semyon Aronovich Gershgorin in 1931. Gershgorin's name has been transliterated in several diffe ...
provides a bound on the error.
In this crude form the iterations are relatively expensive. This can be mitigated by first bringing the matrix ''A'' to upper
Hessenberg form (which costs
arithmetic operations using a technique based on
Householder reduction), with a finite sequence of orthogonal similarity transforms, somewhat like a two-sided QR decomposition.
(For QR decomposition, the Householder reflectors are multiplied only on the left, but for the Hessenberg case they are multiplied on both left and right.) Determining the QR decomposition of an upper Hessenberg matrix costs
arithmetic operations. Moreover, because the Hessenberg form is already nearly upper-triangular (it has just one nonzero entry below each diagonal), using it as a starting point reduces the number of steps required for convergence of the QR algorithm.
If the original matrix is
symmetric
Symmetry (from grc, συμμετρία "agreement in dimensions, due proportion, arrangement") in everyday language refers to a sense of harmonious and beautiful proportion and balance. In mathematics, "symmetry" has a more precise definiti ...
, then the upper Hessenberg matrix is also symmetric and thus
tridiagonal
In linear algebra, a tridiagonal matrix is a band matrix that has nonzero elements only on the main diagonal, the subdiagonal/lower diagonal (the first diagonal below this), and the supradiagonal/upper diagonal (the first diagonal above the main di ...
, and so are all the ''A''
''k''. This procedure costs
arithmetic operations using a technique based on Householder reduction.
[ Determining the QR decomposition of a symmetric tridiagonal matrix costs operations.
The rate of convergence depends on the separation between eigenvalues, so a practical algorithm will use shifts, either explicit or implicit, to increase separation and accelerate convergence. A typical symmetric QR algorithm isolates each eigenvalue (then reduces the size of the matrix) with only one or two iterations, making it efficient as well as robust.
]
Visualization
The basic QR algorithm can be visualized in the case where ''A'' is a positive-definite symmetric matrix. In that case, ''A'' can be depicted as an ellipse
In mathematics, an ellipse is a plane curve surrounding two focus (geometry), focal points, such that for all points on the curve, the sum of the two distances to the focal points is a constant. It generalizes a circle, which is the special ty ...
in 2 dimensions or an ellipsoid
An ellipsoid is a surface that may be obtained from a sphere by deforming it by means of directional scalings, or more generally, of an affine transformation.
An ellipsoid is a quadric surface; that is, a surface that may be defined as the ...
in higher dimensions. The relationship between the input to the algorithm and a single iteration can then be depicted as in Figure 1 (click to see an animation). Note that the LR algorithm is depicted alongside the QR algorithm.
A single iteration causes the ellipse to tilt or "fall" towards the x-axis. In the event where the large semi-axis of the ellipse is parallel to the x-axis, one iteration of QR does nothing. Another situation where the algorithm "does nothing" is when the large semi-axis is parallel to the y-axis instead of the x-axis. In that event, the ellipse can be thought of as balancing precariously without being able to fall in either direction. In both situations, the matrix is diagonal. A situation where an iteration of the algorithm "does nothing" is called a fixed point. The strategy employed by the algorithm is iteration towards a fixed-point. Observe that one fixed point is stable while the other is unstable. If the ellipse were tilted away from the unstable fixed point by a very small amount, one iteration of QR would cause the ellipse to tilt away from the fixed point instead of towards. Eventually though, the algorithm would converge to a different fixed point, but it would take a long time.
Finding eigenvalues versus finding eigenvectors
It's worth pointing out that finding even a single eigenvector of a symmetric matrix is not computable (in exact real arithmetic according to the definitions in computable analysis
In mathematics and computer science, computable analysis is the study of mathematical analysis from the perspective of computability theory. It is concerned with the parts of real analysis and functional analysis that can be carried out in a compu ...
). This difficulty exists whenever the multiplicities of a matrix's eigenvalues are not knowable. On the other hand, the same problem does not exist for finding eigenvalues. The eigenvalues of a matrix are always computable.
We will now discuss how these difficulties manifest in the basic QR algorithm. This is illustrated in Figure 2. (Remember to click the thumbnail). Recall that the ellipses represent positive-definite symmetric matrices. As the two eigenvalues of the input matrix approach each other, the input ellipse changes into a circle. A circle corresponds to a multiple of the identity matrix. A near-circle corresponds to a near-multiple of the identity matrix whose eigenvalues are nearly equal to the diagonal entries of the matrix. Therefore the problem of approximately finding the eigenvalues is shown to be easy in that case. But notice what happens to the semi-axes of the ellipses. An iteration of QR (or LR) tilts the semi-axes less and less as the input ellipse gets closer to being a circle. The eigenvectors can only be known when the semi-axes are parallel to the x-axis and y-axis. The number of iterations needed to achieve near-parallelism increases without bound as the input ellipse becomes more circular.
While it may be impossible to compute the eigendecomposition
In linear algebra, eigendecomposition is the factorization of a matrix into a canonical form, whereby the matrix is represented in terms of its eigenvalues and eigenvectors. Only diagonalizable matrices can be factorized in this way. When the matri ...
of an arbitrary symmetric matrix, it is always possible to perturb the matrix by an arbitrarily small amount and compute the eigendecomposition of the resulting matrix. In the case when the matrix is depicted as a near-circle, the matrix can be replaced with one whose depiction is a perfect circle. In that case, the matrix is a multiple of the identity matrix, and its eigendecomposition is immediate. Be aware though that the resulting eigenbasis
In linear algebra, an eigenvector () or characteristic vector of a linear transformation is a nonzero vector that changes at most by a scalar factor when that linear transformation is applied to it. The corresponding eigenvalue, often denoted ...
can be quite far from the original eigenbasis.
The implicit QR algorithm
In modern computational practice, the QR algorithm is performed in an implicit version which makes the use of multiple shifts easier to introduce. The matrix is first brought to upper Hessenberg form as in the explicit version; then, at each step, the first column of is transformed via a small-size Householder similarity transformation to the first column of (or ), where , of degree , is the polynomial that defines the shifting strategy (often , where and are the two eigenvalues of the trailing principal submatrix of , the so-called ''implicit double-shift''). Then successive Householder transformations of size are performed in order to return the working matrix to upper Hessenberg form. This operation is known as ''bulge chasing'', due to the peculiar shape of the non-zero entries of the matrix along the steps of the algorithm. As in the first version, deflation is performed as soon as one of the sub-diagonal entries of is sufficiently small.
Renaming proposal
Since in the modern implicit version of the procedure no QR decomposition
In linear algebra, a QR decomposition, also known as a QR factorization or QU factorization, is a decomposition of a matrix ''A'' into a product ''A'' = ''QR'' of an orthogonal matrix ''Q'' and an upper triangular matrix ''R''. QR decompo ...
s are explicitly performed, some authors, for instance Watkins, suggested changing its name to ''Francis algorithm''. Golub and Van Loan use the term ''Francis QR step''.
Interpretation and convergence
The QR algorithm can be seen as a more sophisticated variation of the basic "power" eigenvalue algorithm. Recall that the power algorithm repeatedly multiplies ''A'' times a single vector, normalizing after each iteration. The vector converges to an eigenvector of the largest eigenvalue. Instead, the QR algorithm works with a complete basis of vectors, using QR decomposition to renormalize (and orthogonalize). For a symmetric matrix ''A'', upon convergence, ''AQ'' = ''QΛ'', where ''Λ'' is the diagonal matrix of eigenvalues to which ''A'' converged, and where ''Q'' is a composite of all the orthogonal similarity transforms required to get there. Thus the columns of ''Q'' are the eigenvectors.
History
The QR algorithm was preceded by the ''LR algorithm'', which uses the LU decomposition
In numerical analysis and linear algebra, lower–upper (LU) decomposition or factorization factors a matrix as the product of a lower triangular matrix and an upper triangular matrix (see matrix decomposition). The product sometimes includes a pe ...
instead of the QR decomposition. The QR algorithm is more stable, so the LR algorithm is rarely used nowadays. However, it represents an important step in the development of the QR algorithm.
The LR algorithm was developed in the early 1950s by Heinz Rutishauser
Heinz Rutishauser (30 January 1918 – 10 November 1970) was a Swiss mathematician and a pioneer of modern numerical mathematics and computer science.
Life
Rutishauser's father died when he was 13 years old and his mother died three years lat ...
, who worked at that time as a research assistant of Eduard Stiefel
Eduard L. Stiefel (21 April 1909 – 25 November 1978) was a Swiss mathematician. Together with Cornelius Lanczos and Magnus Hestenes, he invented the conjugate gradient method, and gave what is now understood to be a partial construction of the ...
at ETH Zurich
(colloquially)
, former_name = eidgenössische polytechnische Schule
, image = ETHZ.JPG
, image_size =
, established =
, type = Public
, budget = CHF 1.896 billion (2021)
, rector = Günther Dissertori
, president = Joël Mesot
, ac ...
. Stiefel suggested that Rutishauser use the sequence of moments ''y''0T ''A''''k'' ''x''0, ''k'' = 0, 1, … (where ''x''0 and ''y''0 are arbitrary vectors) to find the eigenvalues of ''A''. Rutishauser took an algorithm of Alexander Aitken
Alexander Craig "Alec" Aitken (1 April 1895 – 3 November 1967) was one of New Zealand's most eminent mathematicians. In a 1935 paper he introduced the concept of generalized least squares, along with now standard vector/matrix notation fo ...
for this task and developed it into the ''quotient–difference algorithm'' or ''qd algorithm''. After arranging the computation in a suitable shape, he discovered that the qd algorithm is in fact the iteration ''A''''k'' = ''L''''k''''U''''k'' (LU decomposition), ''A''''k''+1 = ''U''''k''''L''''k'', applied on a tridiagonal matrix, from which the LR algorithm follows.
Other variants
One variant of the ''QR algorithm'', ''the Golub-Kahan-Reinsch'' algorithm starts with reducing a general matrix into a bidiagonal one. This variant of the ''QR algorithm'' for the computation of singular values In mathematics, in particular functional analysis, the singular values, or ''s''-numbers of a compact operator T: X \rightarrow Y acting between Hilbert spaces X and Y, are the square roots of the (necessarily non-negative) eigenvalues of the self- ...
was first described by . The LAPACK
LAPACK ("Linear Algebra Package") is a standard software library for numerical linear algebra. It provides routines for solving systems of linear equations and linear least squares, eigenvalue problems, and singular value decomposition. It also ...
subroutin
DBDSQR
implements this iterative method, with some modifications to cover the case where the singular values are very small . Together with a first step using Householder reflections and, if appropriate, QR decomposition
In linear algebra, a QR decomposition, also known as a QR factorization or QU factorization, is a decomposition of a matrix ''A'' into a product ''A'' = ''QR'' of an orthogonal matrix ''Q'' and an upper triangular matrix ''R''. QR decompo ...
, this forms th
DGESVD
routine for the computation of the singular value decomposition
In linear algebra, the singular value decomposition (SVD) is a factorization of a real or complex matrix. It generalizes the eigendecomposition of a square normal matrix with an orthonormal eigenbasis to any \ m \times n\ matrix. It is related ...
. The QR algorithm can also be implemented in infinite dimensions with corresponding convergence results.
References
External links
*
Notes on orthogonal bases and the workings of the QR algorithm
by Peter J. Olver
Module for the QR Method
C++ Library
{{DEFAULTSORT:Qr Algorithm
Numerical linear algebra