QR Algorithm

The QR algorithm is an iterative method in numerical linear algebra for computing the eigenvalues of a square matrix and, optionally, a Schur decomposition. Its basic step factors a matrix into an orthogonal or unitary factor $Q$ and an upper triangular factor $R$, then multiplies these factors in the reverse order. Practical implementations combine this principle with preliminary matrix reduction, shifts, implicit transformations, and deflation. These refinements make QR a standard method for the complete eigenvalue problem of general dense matrices. (cs.cornell.edu)

Basic iteration

Let $A_0=A$ be an $n\times n$ matrix. The unshifted QR iteration is

where the first equation is a QR decomposition. For real matrices, $Q_k$ is an orthogonal matrix; for complex matrices, it is a unitary matrix. Since $Q_k^*Q_k=I$,

where $^{*}$ denotes the conjugate transpose. Thus each step is a similarity transformation, preserving the eigenvalues in exact arithmetic. (cs.cornell.edu)

If the transformations are accumulated as

then $A_k=Z_k^*AZ_k$. Under suitable spectral and initial-subspace conditions, this process approaches Schur form. The unshifted iteration is closely related to orthogonal iteration: it implicitly performs simultaneous subspace iteration rather than finding one eigenvector at a time. Its convergence depends on eigenvalue-modulus ratios and may be slow or fail to produce a triangular limit when moduli coincide. (cs.cornell.edu)

Schur form and eigenvectors

For a complex matrix, the desired decomposition is

with $Z$ unitary and $T$ upper triangular. The eigenvalues are the diagonal entries of $T$. For a real matrix, real arithmetic produces a real Schur form: $Z$ is orthogonal and $T$ is upper quasi-triangular, with $1\times1$ blocks representing real eigenvalues and $2\times2$ blocks representing complex-conjugate eigenvalue pairs. (netlib.org)

The columns of $Z$, called Schur vectors, are not generally individual eigenvectors. They provide orthonormal bases for suitable invariant subspaces. To recover a right eigenvector, one solves

using the triangular or block-triangular structure and then forms $x=Zy$. LAPACK provides separate routines for this recovery. Consequently, the QR algorithm is more precisely a Schur-form method than a direct diagonalization method. (netlib.org)

Hessenberg reduction

Factoring a general dense matrix at every iteration costs $O(n^3)$ operations per step. Practical QR algorithms therefore first reduce $A$ to an upper Hessenberg matrix:

The preliminary reduction takes $O(n^3)$ operations, usually using Householder transformations. Hessenberg structure is preserved by appropriately implemented QR steps and reduces the cost of a single-shift or double-shift sweep to $O(n^2)$. (cs.cornell.edu)

For a real symmetric matrix, a symmetric Hessenberg matrix is tridiagonal: entries can be nonzero only on the diagonal and its two neighboring diagonals. This stronger structure supports much cheaper iterations and specialized symmetric eigensolvers. (cs.cornell.edu)

Shifts

A shifted QR step replaces the basic iteration by

where $I$ is the identity matrix and $\mu_k$ is a selected scalar. The similarity identity remains

A shift near an eigenvalue can accelerate its separation from the remaining active matrix. Shift selection is therefore central to practical convergence. (cs.cornell.edu)

For symmetric tridiagonal matrices, a widely used choice is the Wilkinson shift: the eigenvalue of the trailing $2\times2$ principal submatrix closest to the bottom-right diagonal entry. With appropriate tie handling, this strategy typically gives rapid convergence of the trailing subdiagonal toward zero. Its local convergence is often cubic, but cubic convergence is not universal; particular symmetric matrices exhibit strictly quadratic convergence. (arxiv.org)

Francis double shifts

A real nonsymmetric matrix can have nonreal eigenvalues, making a real single-shift strategy insufficient. The Francis double-shift step combines two shifts, commonly obtained from the trailing $2\times2$ block. If they are a conjugate pair $\mu,\bar{\mu}$, the associated quadratic is

Its coefficients are real, so the combined transformation can be performed entirely in real arithmetic. Conceptually, it uses the orthogonal factor of $p(H)$, but an efficient implementation does not explicitly form the full matrix polynomial. (cs.cornell.edu)

Implicit QR and bulge chasing

The implicit QR algorithm constructs the required similarity transformation without explicitly factoring the entire shifted matrix. For an unreduced Hessenberg matrix—one whose subdiagonal entries are nonzero—the implicit Q theorem establishes that the first column of the transforming matrix, together with the requirement to restore Hessenberg form, determines the transformation up to sign or phase conventions. (cs.cornell.edu)

In a double-shift step, a small Householder transformation is constructed from the first column of $p(H)$. Applying it on both sides creates a small cluster of entries below the Hessenberg band, called a bulge. Subsequent local transformations move the bulge downward until it exits the matrix. This process, bulge chasing, realizes the intended shifted similarity transformation while retaining the low cost afforded by Hessenberg structure. (cs.cornell.edu)

Deflation and termination

Deflation separates a converged part of the spectrum from the remaining problem. If a Hessenberg subdiagonal entry becomes exactly zero, the matrix splits into block upper-triangular form:

The eigenvalues of the full matrix are then the eigenvalues of the two diagonal blocks. In floating-point arithmetic, sufficiently small subdiagonal entries can be set to zero using scale-aware tests. Production implementations include safeguards beyond a simple absolute threshold. (cs.cornell.edu)

For real nonsymmetric problems, an isolated $2\times2$ block representing a complex-conjugate pair is a valid converged block and need not become diagonal. Successful termination produces the appropriate triangular or quasi-triangular Schur form. Implementations also use exceptional shifts when ordinary shift choices lead to prolonged stagnation. (netlib.org)

Symmetric problems and computational cost

For a real symmetric matrix, the QR process begins with orthogonal reduction to symmetric tridiagonal form. An implicit QR sweep on the tridiagonal representation costs $O(n)$ when only eigenvalues are required. Assuming roughly a bounded number of sweeps per deflated eigenvalue, finding all eigenvalues of the tridiagonal matrix costs $O(n^2)$; the initial reduction of a general dense symmetric matrix still costs $O(n^3)$. Accumulating all eigenvectors adds substantial work. (cs.cornell.edu)

For a general dense nonsymmetric matrix, preliminary Hessenberg reduction and subsequent QR iteration are conventionally assigned an overall $O(n^3)$ operation count. This estimate assumes a modest average number of sweeps per eigenvalue or pair, not a uniform bound for every input. Actual iteration counts depend on the matrix and shift strategy. Full dense matrix storage requires $O(n^2)$ space. (netlib.org)

Specialized symmetric alternatives include divide-and-conquer eigensolvers, Jacobi methods, and bisection with inverse iteration. Depending on the requested eigenvalues or eigenvectors, these can be preferable to QR. (cs.cornell.edu)

Numerical stability and limitations

The principal numerical advantage of QR is its use of orthogonal or unitary transformations. Properly implemented implicit QR iteration, with suitably controlled deflation, is backward stable: its computed Schur decomposition corresponds closely to an exact decomposition of a slightly perturbed input matrix. This is a statement about numerical stability, not a guarantee that every computed eigenvalue has a small relative error. (cs.cornell.edu)

Eigenvalue and eigenvector accuracy also depends on the sensitivity of the original problem, quantified by suitable condition numbers. LAPACK therefore provides separate facilities for estimating eigenpair and invariant-subspace conditioning. A stable algorithm cannot eliminate sensitivity intrinsic to the input matrix. (netlib.org)

Optional matrix balancing can rescale a nonsymmetric problem before Hessenberg reduction. It may improve accuracy, but indiscriminate balancing can also worsen eigenvector errors; the scaling and its reversal require care. (arxiv.org)

Implementations and related methods

LAPACK separates the nonsymmetric computation into stages: xGEHRD reduces a matrix to Hessenberg form, xHSEQR computes its eigenvalues and optional Schur form, and routines such as xTREVC recover eigenvectors. The prefix x denotes a family of routines for different numerical types and precisions. (netlib.org)

Advanced implementations use small-bulge multishift QR, which processes multiple shifts while organizing updates as efficient block matrix operations. Aggressive early deflation examines a trailing window to identify converged eigenvalues before conventional subdiagonal tests would detect them. These techniques improve practical performance without changing the underlying similarity-transformation principle. (netlib.org)

Related QR iterations also occur in singular value decomposition, after reduction to bidiagonal form. They apply transformations from the left and right to compute singular values, rather than similarity transformations to compute eigenvalues. The related QZ algorithm handles generalized eigenvalue problems for a matrix pair through generalized Schur form. (cs.cornell.edu)

History

John G. F. Francis and Vera Kublanovskaya independently developed QR methods around 1959–1961. Francis submitted his first theoretical paper on October 29, 1959; Kublanovskaya submitted an initial summary on July 5, 1960. Their early publications established QR as a method for the complete matrix eigenvalue problem, with Francis's work particularly associated with the implicit double-shift construction. (atm.org.uk)

Subsequent development concentrated on reliable shift strategies, more stringent deflation tests, multishift sweeps, and efficient blocked implementations. LAPACK 3.1 incorporated small-bulge multishift QR with aggressive early deflation, extending the original method into a substantially more elaborate production eigensolver. (netlib.org)

References

1. CS 4220: Numerical Analysis — March 6, 2026cs.cornell.edu
2. A Parallel QR Algorithm for the Nonsymmetric Eigenvalue Problemnetlib.org
3. LAPACK: CHSEQRnetlib.org
4. LAPACK Users' Guidenetlib.org
5. CS 6210: Matrix Computations — October 29, 2025cs.cornell.edu
6. CS 6210: Matrix Computations — November 5, 2025cs.cornell.edu
7. The Asymptotics of Wilkinson's Iteration: Loss of Cubic Convergencearxiv.org
8. LAPACK 3.1 xHSEQR: Tuning and Implementation Notes on the Small Bulge Multi-shift QR Algorithm with Aggressive Early Deflationnetlib.org
9. A Parallel Implementation of the Nonsymmetric QR Algorithmnetlib.org
10. LAPACK Working Note 13netlib.org
11. On Matrix Balancing and Eigenvector Computationarxiv.org
12. The QR Algorithm: 50 Years Later—Its Genesis by John Francis and Vera Kublanovskaya and Subsequent Developmentsatm.org.uk