QR algorithm
The QR algorithm is an iterative method in numerical linear algebra for computing the eigenvalues of a matrix. In its basic form, the method repeatedly factors a matrix into an orthogonal or unitary factor and an upper-triangular factor, then reverses the order of those factors. The resulting matrices are similar to the original matrix and therefore preserve its eigenvalues.
For a real matrix (A_0), one iteration has the form
[ A_k = Q_kR_k, ]
where (Q_k) is orthogonal and (R_k) is upper triangular. The next iterate is defined by
[ A_{k+1}=R_kQ_k. ]
Since (Q_k^{-1}=Q_k^{\mathsf T}), this relation can also be written as
[ A_{k+1}=Q_k^{\mathsf T}A_kQ_k. ]
Thus, every iterate is obtained from its predecessor by an orthogonal similarity transformation. For complex matrices, the transpose is replaced by the conjugate transpose, and (Q_k) is unitary.
The unshifted iteration provides the conceptual foundation of the method, while practical eigensolvers use preliminary reduction, spectral shifts, implicit transformations, and deflation. These modifications preserve the same similarity structure but substantially alter the convergence rate and computational cost.
Mathematical interpretation
The QR iteration is closely related to simultaneous power iteration. Define the accumulated orthogonal factor
[ \widehat Q_k = Q_0Q_1\cdots Q_k. ]
Under suitable spectral conditions, the columns of (\widehat Q_k) approach invariant subspaces associated with successively ordered eigenvalues. At the same time, the transformed matrices
[ A_{k+1}=\widehat Q_k^{\mathsf T}A_0\widehat Q_k ]
approach a triangular or block-triangular form.
For a complex matrix, the limiting structure is associated with the Schur decomposition,
[ A=QTQ^*, ]
where (Q) is unitary and (T) is upper triangular. The diagonal entries of (T) are the eigenvalues of (A). A real matrix with nonreal eigenvalues instead approaches a real Schur form containing upper-triangular structure together with (2\times2) diagonal blocks representing complex-conjugate eigenvalue pairs.
The ordering of the limiting eigenvalues depends on the transformation and the shifts. In the elementary unshifted case, convergence generally reflects separation among eigenvalue magnitudes. Equal or nearly equal magnitudes can produce slow convergence, persistent rotations, or failure of individual matrix entries to settle rapidly even when invariant subspaces remain well determined.
When (A) is real and symmetric, every eigenvalue is real and the Schur form is diagonal. The QR iteration then reduces the off-diagonal entries while preserving symmetry. This setting has particularly direct connections with orthogonal polynomial recurrences, tridiagonal matrices, and the symmetric eigenvalue problem.
Hessenberg and tridiagonal reduction
Applying a complete QR factorization to a dense (n\times n) matrix during every iteration would require work proportional to (n^3). Practical forms of the algorithm first transform the matrix to upper Hessenberg form, in which all entries below the first subdiagonal are zero. The reduction is an orthogonal similarity transformation, usually constructed from Householder transformations.
If
[ H=U^*AU ]
is upper Hessenberg, a QR iteration preserves the Hessenberg structure apart from temporary local fill created during an implicit step. Consequently, each iteration after the reduction requires work proportional to (n^2), rather than the cost of a dense factorization.
A symmetric matrix is reduced to tridiagonal form. Because symmetry and upper Hessenberg structure together imply tridiagonality, the active matrix then contains only its diagonal and two adjacent off-diagonals. An implicit QR step on this form requires work proportional to (n) when eigenvectors are not being accumulated explicitly.
The initial reduction remains a cubic operation for a dense matrix. The later iterative phase usually has a lower cost per step, although the total amount of work depends on eigenvalue separation, shift selection, and the frequency of deflation.
Shifted iteration
A shifted QR step introduces a scalar (\mu_k) and factors
[ A_k-\mu_kI=Q_kR_k. ]
The next iterate is
[ A_{k+1}=R_kQ_k+\mu_kI, ]
which again satisfies
[ A_{k+1}=Q_k^*A_kQ_k. ]
The shift changes the convergence behavior without changing the eigenvalues. If (\mu_k) approximates an eigenvalue associated with the lower-right portion of the active matrix, the corresponding subdiagonal entry commonly decreases much faster than it does under unshifted iteration.
For symmetric tridiagonal matrices, the Wilkinson shift is obtained from the eigenvalue of the trailing (2\times2) principal submatrix that lies nearer the final diagonal entry. This shift usually produces rapid local convergence while retaining real arithmetic. Its effect is especially pronounced near a simple eigenvalue, where the terminal subdiagonal entry may converge cubically in the asymptotic regime.
A real nonsymmetric matrix can contain a complex-conjugate pair even though all stored coefficients are real. The double-shift method applies the quadratic polynomial associated with a trailing (2\times2) block rather than explicitly introducing complex numbers. The two shifts are represented together, allowing the iteration to preserve real arithmetic while approaching a real Schur form.
Implicit QR transformations
Modern QR eigensolvers generally do not form the complete factors (Q_k) and (R_k) at every step. Instead, the initial portion of the desired similarity transformation is introduced near the upper-left corner of the Hessenberg matrix. Subsequent local orthogonal transformations move the resulting disturbance downward through the matrix.
The temporary nonzero entries outside the Hessenberg band are conventionally described as a bulge. Their movement is known as bulge chasing. Once the bulge leaves the active matrix, the Hessenberg form has been restored and the result is mathematically equivalent to the corresponding shifted QR step, apart from rounding effects.
The justification for this construction is provided by the implicit Q theorem. Subject to sign or phase conventions and nondegeneracy conditions, an unreduced Hessenberg matrix and the first column of the transforming unitary matrix determine the remainder of the similarity transformation. This permits shifted iterations to be represented through local transformations rather than explicit matrix factorizations.
You Watanabe’s 1963 analysis of implicit double-shift sweeps established a local turnover identity for adjacent plane transformations. The identity showed that three overlapping transformations could be reordered without changing the represented orthogonal similarity, provided that the newly introduced transformation was chosen to eliminate the same exterior-band entry. This formulation became part of the algebraic description of multi-bulge implementations, where several shifts are propagated through separated sections of a Hessenberg matrix.
Deflation and numerical termination
An iteration reaches a deflatable state when a subdiagonal entry is negligible relative to nearby matrix entries. Setting that entry to zero divides the matrix into smaller independent eigenvalue problems. For a Hessenberg matrix, a zero in position ((i+1,i)) produces a block upper-triangular partition, and the eigenvalues of the full matrix are the combined eigenvalues of the diagonal blocks.
A common scale-sensitive deflation condition has the form
[ |h_{i+1,i}| \leq \varepsilon\bigl(|h_{i,i}|+|h_{i+1,i+1}|\bigr), ]
where (\varepsilon) reflects the working floating-point arithmetic. More refined tests account for local coupling, underflow, and the possibility that neighboring diagonal entries provide an inadequate scale.
Deflation is not merely a declaration that an eigenvalue approximation has stabilized. It changes the structure of the remaining computation by reducing the active matrix order. When a trailing (1\times1) block separates, its diagonal entry represents a real eigenvalue. In real arithmetic, a trailing (2\times2) block may instead be retained as the representation of a complex-conjugate pair.
Contemporary implementations also use aggressive early deflation. A small trailing window is transformed to expose converged eigenvalues before the ordinary bottom subdiagonal entries have become negligible under a simple local test. The window transformation remains orthogonal or unitary, so the resulting separation is consistent with the similarity framework of the algorithm.
Convergence and stability
The basic QR iteration does not converge uniformly for every matrix. Its behavior depends on spectral separation and on the structure of the eigenvectors. For normal matrices, orthogonal invariant subspaces provide a comparatively direct description of convergence. Highly nonnormal matrices can exhibit substantial sensitivity because their eigenvectors may be far from orthogonal.
Shifted QR iteration is generally analyzed through the decay of subdiagonal entries and the convergence of invariant subspaces. A shift close to a simple eigenvalue accelerates separation of the corresponding trailing block. Repeated or tightly clustered eigenvalues can delay deflation, while exceptional shifts can interrupt stagnation when a sequence of ordinary shifts fails to reduce the active coupling.
The use of orthogonal or unitary similarity transformations gives the algorithm a favorable backward-error structure. Such transformations preserve the Euclidean norm and do not magnify perturbations merely through a change of basis. Nevertheless, finite-precision behavior also depends on the accumulation of local rounding errors, the reliability of deflation tests, and the conditioning of the underlying eigenvalue problem.
An accurate Schur form does not imply equally accurate eigenvectors. Eigenvalues of a strongly nonnormal matrix may be intrinsically sensitive to small perturbations, and nearly dependent eigenvectors can remain poorly conditioned even when the computed factorization has a small residual. These effects belong to the conditioning of the mathematical problem rather than solely to the mechanics of QR iteration.
Historical development
The QR algorithm developed from earlier factorization iterations, particularly Heinz Rutishauser’s LR method. The LR iteration factors a matrix into lower- and upper-triangular components and reverses their order, but it does not preserve norms through orthogonal transformations. The replacement of triangular elimination by QR factorization produced a method with more stable numerical behavior.
John G. F. Francis and Vera Kublanovskaya formulated the QR algorithm independently at the beginning of the 1960s. Francis developed the implicit shifted form and the real double-shift construction, while Kublanovskaya presented the orthogonal-factor iteration as a general method for matrix eigenvalue computation. Their work connected factorization iteration with reduction to Hessenberg form and with convergence toward triangular canonical structure.
Alston Scott Householder’s development of orthogonal reflections supplied the standard mechanism for reducing dense matrices while controlling rounding-error growth. James H. Wilkinson subsequently analyzed shift strategies, convergence, and finite-precision behavior in the symmetric eigenvalue problem. These developments established the main mathematical architecture retained by later QR-based eigensolvers.
Block transformations, multi-shift sweeps, and aggressive early deflation altered the organization of the computation without replacing its defining principle. The modern algorithm remains an iteration by orthogonal or unitary similarity transformations whose limiting form reveals the spectrum.