QR algorithm

1 Basic ideas and goals of the QR algorithm

The QR algorithm is an iterative method designed to compute the eigenvalues of a square matrix, and—through related procedures—its eigenvectors. It is especially valued in scientific computing because it is broadly reliable under floating-point arithmetic and scales well to large problems.

1.1 Eigenvalues and invariants under similarity transforms

Eigenvalues are preserved under similarity transformations. If a nonsingular matrix \(S\) is used to form \(B = S^{-1}AS\), then \(A\) and \(B\) share the same eigenvalues. The QR algorithm exploits this fact by constructing a sequence of similar matrices, so that every iterate has the same eigenvalues as the original matrix (up to numerical error).

1.2 QR decomposition (A = QR) and the “shifted” iteration concept

At each iteration, the current matrix \(A_k\) is factored as \(A_k = Q_kR_k\), where \(Q_k\) is orthogonal (or unitary in the complex case) and \(R_k\) is upper triangular. The next iterate is formed by reversing the factors: \[ A_{k+1} = R_kQ_k. \] A “shifted” variant instead factors \(A_k-\mu I\) (for a scalar \(\mu\)), applies the same idea, and then adds the shift back. Shifts are used to steer the iteration so that certain parts of the spectrum converge more quickly.

1.3 Convergence overview and what “converges” means

In the QR algorithm, convergence typically means that the iterates approach a form whose structure reveals the spectrum. For real symmetric problems, convergence often produces a nearly diagonal matrix. For general nonsymmetric matrices, convergence is commonly described via convergence to an upper triangular (or unitary similarity equivalent) structure, reflecting the eigenvalues along the diagonal (or in blocks in related decompositions). Practical implementations also use “deflation” to separate already-resolved eigencomponents.

2 QR decomposition fundamentals

The method’s foundation is the efficient and numerically stable computation of a QR decomposition, along with careful handling of orthogonality and triangular structure during iterative updates.

2.1 Householder reflections and orthogonal/unitary factors

A common way to compute \(A=QR\) is via Householder reflections. These are orthogonal transforms that systematically introduce zeros below the diagonal of \(R\). Accumulating the reflections yields \(Q\). Householder-based QR is widely used because it is stable and reduces sensitivity to rounding errors compared with naïve Gram–Schmidt approaches.

2.2 Givens rotations and sparse/structured updates

Givens rotations use plane rotations to eliminate specific entries. They can be more local than Householder transformations, which can be advantageous when matrices have sparsity or when one wants targeted updates. In QR iteration contexts, Givens-style updates may be used in specialized implementations or conceptual explanations to clarify how sparsity patterns evolve.

2.3 Numerical considerations (orthogonality, conditioning)

Numerical stability hinges on keeping \(Q\) close to orthogonal/unitary despite floating-point rounding. Well-conditioned problems usually behave better, but the QR algorithm is designed to be robust across a wide range of matrix structures. Conditioning affects how perturbations in the input translate into changes in computed eigenvalues and especially eigenvectors.

2.4 Complexity and implementation trade-offs

The cost of QR iteration depends heavily on matrix size and on preprocessing that reduces the matrix to a structured form (often Hessenberg for nonsymmetric and tridiagonal for symmetric cases). Trade-offs arise between update cost per iteration, number of iterations needed for convergence, and overhead for shifts and deflation logic. Implementations typically prefer algorithms that keep per-iteration work small while maintaining strong numerical behavior.

3 The unshifted QR algorithm

This section describes the core iteration without the acceleration benefits introduced by shifts, highlighting the basic similarity mechanism and typical convergence behavior.

3.1 Iteration definition and matrix similarity

Given \(A_k\), compute \(A_k = Q_kR_k\) and set \(A_{k+1} = R_kQ_k\). Because \(A_{k+1} = Q_k^{\!*}A_kQ_k\) (with appropriate notation for conjugate transpose), each step is a similarity transformation. Consequently, all iterates share the same eigenvalues as the original matrix, neglecting rounding errors.

3.2 Behavior on general matrices

For general nonsymmetric matrices, unshifted QR iteration may converge slowly or stall, particularly when eigenvalues are clustered or poorly separated. The iteration can still move the matrix toward an upper triangular-like structure, but the speed depends strongly on the spectral distribution and on how eigenvalues interact with the transformation process.

3.3 Convergence conditions and typical limitations

Unshifted convergence is not guaranteed to be fast in general. Certain matrices exhibit behavior where the diagonal entries do not meaningfully stabilize until many iterations have been performed. As a result, unshifted QR is often used mainly as a conceptual baseline or in special cases, while practical eigenvalue solvers almost always incorporate shifts and structured reductions.

4 Shifted QR algorithms

Shifts are central to the practical success of QR iteration. They modify the iteration so that the algorithm behaves like a refined process targeting specific eigenvalues.

4.1 Single-shift strategy and the role of spectral information

In single-shift QR, a scalar \(\mu\) is chosen and the factorization is applied to \(A_k-\mu I\). The shift effectively changes the dynamics of the iteration so that eigenvalues near \(\mu\) typically converge more rapidly. Conceptually, the method uses partial spectral information to focus progress.

4.2 Rayleigh quotient shifts

A Rayleigh quotient shift uses an estimate derived from the current iterate, often tied to a vector or a local substructure. When the shift approximates an eigenvalue well, the method can exhibit faster convergence, particularly in favorable scenarios where the iterate aligns with the corresponding invariant subspace.

4.3 Wilkinson shift for improved practical convergence

The Wilkinson shift selects \(\mu\) using information from the trailing \(2\times 2\) principal submatrix. This choice is designed to approximate an eigenvalue likely to govern the convergence of the bottom-right portion of the matrix. In practice, the Wilkinson shift is known for producing consistently strong convergence behavior in many real symmetric and related settings.

4.4 Multi-shift approaches (high-level coverage)

Multi-shift strategies apply several shifts per iteration or use polynomial filtering ideas to accelerate convergence across multiple eigenvalues. High-level schemes can reduce the number of matrix-vector-like update stages by packing several acceleration steps into a structured transformation. While the details vary across algorithms, the overarching goal is the same: accelerate convergence without sacrificing stability.

5 Reduction to Hessenberg form

Before applying QR iteration efficiently, matrices are often transformed into forms that preserve eigenvalues while enabling cheap updates.

5.1 Why Hessenberg reduction is used

A matrix in Hessenberg form has zeros below the first subdiagonal, meaning only a narrow band of entries can be nonzero. This structure reduces the work needed in each QR iteration step, because orthogonal transformations only need to affect a limited region.

5.2 Algorithms for transforming to (upper) Hessenberg matrices

Common procedures use similarity transformations built from Householder reflections. The goal is to eliminate entries below the first subdiagonal while maintaining similarity to the original matrix. The resulting Hessenberg matrix has the same eigenvalues as the original matrix, and the orthogonal similarity factors can be optionally stored for later eigenvector recovery.

5.3 Effect on computational cost and convergence behavior

Reduction to Hessenberg form typically lowers the cost from what would be needed for dense QR factorizations at every step. While the convergence characteristics depend on spectral properties and the use of shifts, the Hessenberg structure makes each iteration computationally feasible for large matrices.

5.4 Back-transformation and eigenvector recovery

If eigenvectors are desired, the transformations used during reduction must be applied (directly or implicitly) to the eigenvectors of the reduced matrix. Because QR iteration may produce eigenvectors for the Hessenberg form (or for an even more structured final form), back-transformation ensures the eigenvectors correspond to the original matrix.

6 The QR algorithm for symmetric matrices

Symmetric matrices offer a particularly favorable environment: real eigenvalues, strong stability properties, and well-behaved convergence to a nearly diagonal form.

6.1 Real symmetric case and guaranteed stability

For a real symmetric matrix, QR iteration can be carried out using orthogonal matrices, producing iterates that remain symmetric. This structure supports numerical stability and makes it easier to interpret convergence. The eigenvalues are real, and the algorithm’s progression toward diagonalization aligns with the symmetry-driven invariant subspaces.

6.2 Tridiagonal reduction and accelerated iterations

Symmetric matrices are often reduced further to tridiagonal form, which has nonzeros only on the main diagonal and the first sub- and super-diagonals. Because tridiagonal QR iteration can be updated using only short-range operations, each iteration becomes much cheaper. Additionally, convergence is typically accelerated because the algorithm’s transformations preserve the tridiagonal structure.

6.3 Interlacing properties and practical interpretation

In many symmetric eigenproblems, eigenvalues associated with leading principal submatrices interlace with those of the full matrix. This provides interpretive guidance about how convergence proceeds as the iteration deflates the matrix into smaller blocks.

6.4 Eigenvalue extraction from tridiagonal/triangular forms

As QR iteration converges, the off-diagonal entries decay toward zero, leaving an approximately diagonal matrix (or a block diagonal structure in certain extensions). The diagonal entries of the converged form then serve as the computed eigenvalues. In implementations, eigenvalues are often extracted as deflation identifies portions of the matrix that have effectively decoupled.

7 The QR algorithm for general (nonsymmetric) matrices

For nonsymmetric matrices, eigenvalues can be complex and the algorithm’s convergence is typically described via triangular (or Schur-like) forms rather than diagonalization.

7.1 Complex eigenvalues and real Schur form

In the complex setting, QR iteration can converge to an upper triangular matrix whose diagonal contains the eigenvalues. In real arithmetic with real inputs, convergence may produce a real Schur form, where complex conjugate eigenpairs appear as \(2\times 2\) blocks on the diagonal. This representation preserves real computations while still capturing complex spectral information.

7.2 Bulge chasing (conceptual description)

When working with Hessenberg matrices, QR iteration with shifts affects only a local region near the trailing part of the matrix, but that local action creates nonzeros outside the intended band. “Bulge chasing” is a conceptual way to describe how subsequent orthogonal transformations restore the Hessenberg structure by sweeping the introduced nonzeros along the matrix. The name reflects the transient “bulge” that moves toward the bottom.

7.3 Robustness considerations and deflation

Nonsymmetric problems can be more sensitive than symmetric ones, especially regarding convergence rate and eigenvector conditioning. Deflation remains important: once subdiagonal entries become sufficiently small, the matrix effectively splits, and the iteration can focus on smaller subproblems to improve efficiency and reduce numerical work.

7.4 Computing a Schur decomposition using QR iteration

A standard use of QR iteration for nonsymmetric matrices is to compute a Schur decomposition, where \(A = UTU^{\!*}\), with \(U\) unitary/orthogonal and \(T\) upper triangular (or block upper triangular in real arithmetic). The eigenvalues appear in \(T\), and \(U\) provides an orthonormal basis associated with the Schur vectors. This is often more stable and reliable than directly attempting full eigen-decomposition.

8 Deflation and stopping criteria

Deflation is the practical mechanism that turns an iterative process into a sequence of smaller problems, accelerating both speed and accuracy.

8.1 Deflation concept and identifying nearly decoupled subproblems

Deflation occurs when an entry that couples two parts of the matrix becomes negligibly small. For Hessenberg or tridiagonal forms, this often corresponds to subdiagonal elements approaching zero. When such an element is tiny relative to the scale of the matrix, the matrix behaves like a block diagonal system, and the QR iteration can proceed independently on each block.

8.2 Thresholding and practical convergence checks

Implementations use thresholds to decide when an entry is “small enough” to treat as zero. These checks account for both absolute and relative magnitudes, and they also consider machine precision. The objective is to avoid declaring convergence too early (which could impair accuracy) while preventing unnecessary iterations.

8.3 Balancing of matrices (when applicable)

For nonsymmetric problems, scaling or balancing can reduce the spread of magnitudes among rows and columns. Although “balancing” can vary across libraries, the general idea is to improve numerical conditioning and make subdiagonal decay more interpretable. This can enhance convergence behavior and reduce roundoff sensitivity.

8.4 Storing and updating subdiagonal elements

Efficient stopping and deflation depend on monitoring specific entries, such as subdiagonal elements in Hessenberg or tridiagonal forms. These values must be updated accurately as transformations are applied. Data handling and update ordering influence both runtime and the reliability of convergence tests.

9 Eigenvectors and accumulation of transformations

Eigenvalues alone can be extracted from triangular-like limits, but eigenvectors require additional structure and bookkeeping.

9.1 Accumulating orthogonal/unitary factors

To compute eigenvectors, one must track the orthogonal (or unitary) transformations generated during reduction and QR iteration. Some algorithms accumulate these factors explicitly, while others apply them implicitly to reduce memory and computational overhead. The final eigenvectors are obtained by applying the accumulated transforms to vectors from the converged form.

9.2 Backward error considerations (high-level)

Even with a stable eigenvalue method, eigenvectors may exhibit larger relative errors because eigenvectors can change substantially under small perturbations if eigenvalues are close together or multiple. Backward error perspectives evaluate how much perturbation of the original matrix would be needed to make the computed vectors exact, offering a framework for understanding observed accuracy.

9.3 Refinement and normalization of computed eigenvectors

After eigenvectors are computed, normalization ensures consistent scaling. Refinement procedures may be used to improve accuracy, such as applying additional iteration or solving small linear systems related to residual minimization. These steps aim to reduce the norm of \(Av-\lambda v\), improving practical correctness.

10 Variants and related algorithms

QR iteration has many closely related forms and conceptual relatives that differ mainly in how shifts are applied or how triangular structure is managed.

10.1 Implicit QR iteration (high-level explanation)

Implicit QR iteration avoids explicitly forming the full \(Q_k\) and \(R_k\) matrices. Instead, it applies the necessary effect of \(R_kQ_k\) to the current structured matrix using sequences of orthogonal transformations. This reduces both cost and memory, and it is particularly effective on Hessenberg or tridiagonal forms.

10.2 Francis QR algorithm connection (overview)

The Francis QR algorithm is a well-known implementation strategy for shifted QR iteration on Hessenberg or tridiagonal matrices. It provides an efficient mechanism for applying shifts implicitly while maintaining structure and enabling bulge chasing. The method is often associated with the practical use of shifts such as the Wilkinson shift.

10.3 Comparison to power iteration and divide-and-conquer

Power iteration targets the dominant eigenvalue and eigenvector and converges depending on eigenvalue separation, but it is not generally suited to computing the full spectrum. Divide-and-conquer approaches, especially for symmetric tridiagonal matrices, can be faster for some cases. QR iteration remains attractive because of its general applicability, stability, and effectiveness for producing all eigenvalues with controlled numerical behavior.

10.4 When QR is preferred in practice

QR-based methods are commonly preferred when a robust general-purpose eigenvalue computation is needed, including for nonsymmetric problems. In symmetric cases, QR may be chosen for its simplicity and reliability, while specialized alternatives may outperform it in certain structured regimes. The choice depends on matrix type, desired outputs, and performance requirements.

11 Applications and practical use

QR iteration appears across numerical linear algebra, serving as a core component in many eigenvalue solver workflows.

11.1 Standard workflows in eigenvalue solvers

A typical workflow is: preprocess the matrix into a structured form (Hessenberg or tridiagonal), apply shifted QR iteration with deflation until convergence, extract eigenvalues from the final form, and optionally recover eigenvectors using accumulated transformations. For large problems, the method is often paired with sparsity-aware preprocessing or efficient dense-to-structured transformations.

11.2 Software/library perspectives (conceptual, not vendor-specific)

Most mature numerical libraries implement QR-based eigenvalue routines with careful attention to stability, performance, and edge cases such as clustered eigenvalues and nearly defective matrices. They also manage memory usage for accumulating transforms and provide options controlling shifts, tolerances, and whether eigenvectors are computed.

11.3 Common benchmarks and performance metrics

Performance is commonly assessed by runtime and by the accuracy of returned eigenpairs, often using residual norms or backward error estimates. Benchmarks may include random dense matrices, structured test suites, and matrices with challenging eigenvalue distributions to evaluate convergence behavior under realistic conditions.

12 Numerical stability and error analysis (applied focus)

QR’s reputation for stability comes from the orthogonal transformations it uses and from its compatibility with floating-point arithmetic.

12.1 Backward stability intuition

Backward stability suggests that the computed results are close to the exact results for a slightly perturbed input matrix. Since QR factorization using Householder reflections is backward stable in standard settings, the iterative process tends to inherit good error control, especially for structured matrices like symmetric ones.

12.2 Sources of floating-point error in QR iteration

Errors can accumulate from repeated orthogonal transformations and from finite precision when forming and updating \(A_k-\mu I\), computing shifts, and applying bulge-chasing transformations. Cancellation can occur when subtracting nearly equal quantities, and rounding affects the smallest subdiagonal entries most directly connected to deflation logic.

12.3 Mitigation strategies in implementations

Implementations mitigate error using orthogonality-preserving transformations, careful scaling and balancing, robust shift selection, and adaptive stopping thresholds. When eigenvectors are required, residual-based checks and optional refinement can improve the final quality.

13 Worked examples and intuition-building

Small-scale demonstrations help clarify how QR iteration moves a matrix toward a form where eigenvalues become visible.

13.1 Small 2×2 and 3×3 demonstrations (conceptual)

For a \(2\times 2\) matrix, one iteration of QR (especially with a shift) can quickly reveal how diagonal entries relate to eigenvalues, even if the intermediate matrices are less interpretable. A \(3\times 3\) example highlights that convergence often proceeds from the bottom-right corner under shifted strategies, reflecting how trailing principal submatrices influence the shift.

13.2 Tridiagonal toy problems

Tridiagonal toy matrices make the banded structure explicit. As iterations progress, the subdiagonal entries shrink, and deflation can separate the matrix into smaller blocks. This provides an intuitive picture of how QR iteration reduces the problem size over time.

13.3 Reading eigenvalues from the converged form

Once convergence has occurred, eigenvalues can be read from the diagonal of the final triangular form (or from block diagonals in real Schur representations). The residual magnitude provides a way to interpret numerical reliability: if the matrix is nearly in its converged structure and the residuals are small, the extracted eigenvalues are typically accurate.