Chapter 94. QR Algorithm

The QR algorithm is an iterative method for computing eigenvalues of a square matrix. It is one of the central algorithms of numerical linear algebra.

The method is based on a simple idea. Factor the current matrix into an orthogonal factor and an upper triangular factor,

A_k = Q_kR_k,

then reverse the factors:

A_{k+1}=R_kQ_k.

Since

A_{k+1} = R_kQ_k = Q_k^TA_kQ_k,

the new matrix is orthogonally similar to the old one. Therefore it has the same eigenvalues. Repeating this process transforms the matrix toward a form where the eigenvalues can be read from the diagonal or from small diagonal blocks. The practical QR algorithm is usually implemented after reduction to Hessenberg form and with shifts and deflation.

94.1 The Eigenvalue Problem

Let

A\in \mathbb{R}^{n\times n}.

The eigenvalue problem asks for scalars $\lambda$ and nonzero vectors $v$ such that

Av=\lambda v.

Equivalently, the eigenvalues are the roots of

\det(A-\lambda I)=0.

For small matrices, this equation may be solved directly. For large matrices, computing and solving the characteristic polynomial is numerically poor. The QR algorithm avoids explicit polynomial root-finding.

Its goal is to transform $A$ by similarity transformations into a nearly triangular matrix. Since triangular matrices have eigenvalues on the diagonal, the eigenvalues then become visible.

94.2 Similarity Transformations

Two matrices $A$ and $B$ are similar if

B=S^{-1}AS

for some invertible matrix $S$ .

Similar matrices have the same eigenvalues.

The QR algorithm uses orthogonal similarity transformations:

A_{k+1}=Q_k^TA_kQ_k.

This is numerically important because orthogonal transformations preserve vector norms and are stable in floating point arithmetic.

The algorithm therefore changes the representation of the matrix without changing its eigenvalues.

94.3 QR Factorization Step

At iteration $k$ , compute a QR factorization

A_k=Q_kR_k,

where:

Symbol	Meaning
$Q_k$	orthogonal matrix
$R_k$	upper triangular matrix

Then define

A_{k+1}=R_kQ_k.

Because

R_k=Q_k^TA_k,

we get

A_{k+1}=Q_k^TA_kQ_k.

Thus every QR step is an orthogonal similarity transformation.

The sequence

A_0,A_1,A_2,\ldots

preserves all eigenvalues of the original matrix.

94.4 Basic QR Algorithm

The unshifted QR algorithm is:

qr_algorithm(A, max_iter):
    A_k = A

    for k = 0 to max_iter - 1:
        Q, R = qr_factorization(A_k)
        A_k = R * Q

    return A_k

If the method converges, $A_k$ approaches an upper triangular matrix:

A_k \to T.

The eigenvalues are then approximately

t_{11},t_{22},\ldots,t_{nn}.

For symmetric matrices, the limit is diagonal rather than merely triangular.

94.5 Why Reversing the Factors Works

The QR step starts from

A_k=Q_kR_k.

Multiplying on the left by $Q_k^T$ gives

R_k=Q_k^TA_k.

Then

A_{k+1}=R_kQ_k = Q_k^TA_kQ_k.

Thus the operation $A_k\mapsto A_{k+1}$ is not an arbitrary multiplication. It is a change of orthonormal basis.

The eigenvalues remain fixed, while the matrix representation changes.

The iteration attempts to choose these basis changes so that the matrix becomes closer to triangular form.

94.6 Schur Form

The Schur decomposition states that a square complex matrix can be written as

A=QTQ^*,

where $Q$ is unitary and $T$ is upper triangular.

For real matrices, the real Schur form is quasi-upper triangular. It has $1\times 1$ blocks for real eigenvalues and $2\times 2$ blocks for complex conjugate pairs.

The QR algorithm may be understood as a method for computing Schur form. Once Schur form is obtained, eigenvalues are read from the diagonal entries or from the $2\times 2$ diagonal blocks.

94.7 Symmetric Case

A=A^T,

then all eigenvalues are real, and the Schur form is diagonal.

The QR algorithm attempts to compute

A=Q\Lambda Q^T,

where

\Lambda= \operatorname{diag}(\lambda_1,\lambda_2,\ldots,\lambda_n).

In this case, the iterates tend toward a diagonal matrix under suitable conditions.

The diagonal entries approximate the eigenvalues. The accumulated orthogonal transformations approximate the eigenvectors.

94.8 Relation to Power Iteration

Power iteration repeatedly applies $A$ to one vector and normalizes:

x,\quad Ax,\quad A^2x,\ldots.

The QR algorithm may be viewed as a simultaneous version of power iteration.

Instead of tracking one vector, it tracks a full orthonormal basis.

Repeated multiplication by $A$ tends to align basis vectors with invariant directions. QR factorization reorthogonalizes the basis at each step.

This explains the connection between eigenvalue separation and convergence. When eigenvalues are well separated, directions separate faster.

94.9 Accumulating Eigenvectors

If eigenvectors are desired, one accumulates the orthogonal transformations.

Let

Z_0=I.

At each step,

A_k=Q_kR_k, \qquad A_{k+1}=R_kQ_k,

and update

Z_{k+1}=Z_kQ_k.

Then

A_k=Z_k^TAZ_k.

If $A_k$ converges to a diagonal matrix $\Lambda$ , then

A\approx Z_k\Lambda Z_k^T.

The columns of $Z_k$ approximate eigenvectors in the symmetric case.

94.10 Cost of the Naive Algorithm

A direct QR factorization of a dense $n\times n$ matrix costs

O(n^3)

operations per iteration.

Many iterations may be required.

Therefore the naive algorithm is too expensive for practical dense eigenvalue computation.

The practical QR algorithm first reduces the matrix to a simpler form:

Matrix type	Reduced form
General square matrix	upper Hessenberg form
Symmetric matrix	tridiagonal form

The QR iteration then preserves this reduced structure and becomes much cheaper.

94.11 Hessenberg Form

An upper Hessenberg matrix has zeros below the first subdiagonal:

H= \begin{bmatrix} h_{11} & h_{12} & h_{13} & \cdots & h_{1n}\\ h_{21} & h_{22} & h_{23} & \cdots & h_{2n}\\ 0 & h_{32} & h_{33} & \cdots & h_{3n}\\ 0 & 0 & h_{43} & \cdots & h_{4n}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & h_{n,n-1} & h_{nn} \end{bmatrix}.

Every square matrix can be reduced to Hessenberg form by orthogonal similarity transformations.

For a general matrix, practical QR iteration is applied to this Hessenberg matrix rather than to the original dense matrix. This reduction is a standard first stage of the QR algorithm.

94.12 Tridiagonal Form for Symmetric Matrices

If $A$ is symmetric, Hessenberg reduction preserves symmetry.

A symmetric Hessenberg matrix is tridiagonal:

T= \begin{bmatrix} d_1 & e_1 & 0 & \cdots & 0\\ e_1 & d_2 & e_2 & \cdots & 0\\ 0 & e_2 & d_3 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & e_{n-1}\\ 0 & 0 & 0 & e_{n-1} & d_n \end{bmatrix}.

Thus the symmetric QR algorithm usually has two stages:

Stage	Operation
Reduction	$A$ is reduced to tridiagonal form
Iteration	QR steps are applied to the tridiagonal matrix

This structure makes the symmetric eigenvalue problem much cheaper than the fully general case.

94.13 Preservation of Hessenberg Structure

If $H_k$ is upper Hessenberg, then a QR step produces another upper Hessenberg matrix:

H_{k+1}=R_kQ_k.

This is crucial.

It means the expensive reduction is performed once. Later iterations stay within the Hessenberg class.

For a Hessenberg matrix, one QR step can be computed in

O(n^2)

operations.

For a symmetric tridiagonal matrix, one step can be computed in

O(n)

operations.

94.14 Shifted QR Algorithm

Unshifted QR iteration may converge slowly.

A shifted QR step chooses a scalar $\mu_k$ , factors

A_k-\mu_k I=Q_kR_k,

then defines

A_{k+1}=R_kQ_k+\mu_k I.

This is still a similarity transformation:

A_{k+1} = Q_k^TA_kQ_k.

The shift changes the factorization step but preserves eigenvalues.

The purpose of $\mu_k$ is to accelerate convergence by choosing a value close to an eigenvalue. Practical QR algorithms use shifts extensively.

94.15 Shifted QR Pseudocode

shifted_qr_algorithm(A, max_iter):
    A_k = hessenberg_reduction(A)

    for k = 0 to max_iter - 1:
        mu = choose_shift(A_k)

        Q, R = qr_factorization(A_k - mu * I)
        A_k = R * Q + mu * I

        deflate if a subdiagonal entry is small

    return A_k

This pseudocode hides important implementation details. Modern algorithms usually use implicit shifts rather than explicitly forming a full QR factorization at every step.

94.16 Rayleigh Quotient Shift

A simple shift choice is

\mu_k=(A_k)_{nn}.

This uses the bottom-right diagonal entry as an eigenvalue estimate.

For some matrices, this accelerates convergence substantially.

However, more refined shifts are usually preferred.

94.17 Wilkinson Shift

For symmetric tridiagonal matrices, the Wilkinson shift is a standard choice.

It uses the eigenvalue of the trailing $2\times 2$ principal submatrix that is closer to the bottom-right diagonal entry.

If the trailing block is

\begin{bmatrix} a & b\\ b & d \end{bmatrix},

then the shift is chosen from its two eigenvalues.

This often gives very rapid convergence of the bottom subdiagonal entry to zero.

94.18 Deflation

Deflation occurs when a subdiagonal entry becomes sufficiently small.

Suppose

|(A_k)_{i+1,i}|

is tiny.

Then $A_k$ is nearly block upper triangular:

A_k \approx \begin{bmatrix} A_{11} & A_{12}\\ 0 & A_{22} \end{bmatrix}.

The eigenvalues of $A_k$ are then approximately the union of the eigenvalues of $A_{11}$ and $A_{22}$ .

The problem can be split into smaller eigenvalue problems.

In practice, deflation is essential for efficiency. It allows one eigenvalue or one conjugate pair to be isolated and removed from the active iteration.

94.19 Practical Deflation Test

A common deflation test checks whether

|a_{i+1,i}| \le \tau \left(|a_{ii}|+|a_{i+1,i+1}|\right),

where $\tau$ is a tolerance related to machine precision.

If this inequality holds, the subdiagonal entry is set to zero.

This produces a block split.

The QR algorithm is then applied separately to each active block.

94.20 Real Schur Form and Complex Eigenvalues

A real matrix may have complex eigenvalues.

Since the QR algorithm can be performed using real arithmetic, complex conjugate eigenvalues appear as $2\times 2$ diagonal blocks in real Schur form.

A block

\begin{bmatrix} a & b\\ c & d \end{bmatrix}

represents the eigenvalues of that $2\times 2$ matrix.

If its characteristic polynomial has negative discriminant, the corresponding eigenvalues are complex conjugates.

Thus real QR iteration does not need to store complex numbers to represent complex eigenvalue pairs.

94.21 Implicit QR Algorithm

Modern implementations usually use the implicit QR algorithm.

Instead of explicitly computing

Q_kR_k=A_k-\mu I

and multiplying

R_kQ_k+\mu I,

the algorithm applies a sequence of orthogonal transformations directly to the Hessenberg matrix.

These transformations introduce a small nonzero pattern called a bulge and then chase it down the matrix until Hessenberg form is restored.

This process is called bulge chasing.

Implicit QR is more efficient, more stable, and better suited to multiple shifts.

94.22 Francis QR Step

The implicit shifted QR step is often called the Francis QR step.

It was developed independently by John Francis and Vera Kublanovskaya in the late 1950s and early 1960s. The QR algorithm is historically one of the major developments in numerical eigenvalue computation.

In modern dense eigenvalue solvers, the Francis step is the practical core of the QR algorithm.

94.23 Double Shifts

For real matrices, complex eigenvalues occur in conjugate pairs.

A double-shift QR step uses a quadratic polynomial

p(t)=(t-\mu)(t-\overline{\mu})

so that the computation can remain real even when the shifts are complex.

The shifts are often chosen from the trailing $2\times 2$ block.

Double shifts are standard in real nonsymmetric QR algorithms.

94.24 Numerical Stability

The QR algorithm uses orthogonal transformations.

Orthogonal transformations preserve norms:

\|Qx\|_2=\|x\|_2.

This gives the method strong numerical stability compared with methods based on nonorthogonal transformations.

The LR algorithm, an older method based on LU factorization, is generally less stable. QR replaced it in practical eigenvalue computation largely because orthogonal transformations behave better numerically.

94.25 Eigenvectors

The QR algorithm naturally computes eigenvalues.

Eigenvectors require additional work.

For symmetric matrices, if the accumulated orthogonal matrix $Z$ is stored, then

A\approx Z\Lambda Z^T,

and the columns of $Z$ approximate eigenvectors.

For nonsymmetric matrices, eigenvectors can be computed from Schur form by solving triangular systems.

Often numerical software first computes Schur form, then optionally computes eigenvectors.

94.26 Relation to Hessenberg and Schur Decompositions

The practical dense eigenvalue pipeline is:

Step	Purpose
Balance, optional	improve scaling
Hessenberg reduction	reduce cost of iteration
Implicit shifted QR	converge toward Schur form
Deflation	split completed blocks
Back transformation	recover Schur vectors or eigenvectors

For symmetric matrices, the corresponding pipeline is:

Step	Purpose
Tridiagonal reduction	reduce dense symmetric problem
Shifted QR or related method	compute eigenvalues
Accumulate transformations	compute eigenvectors

The QR algorithm is therefore part of a larger eigenvalue solver.

94.27 Comparison with Power Iteration

Feature	Power iteration	QR algorithm
Goal	dominant eigenpair	full spectrum
Main operation	matrix-vector product	orthogonal similarity steps
Memory	low	higher
Best use	large sparse leading eigenvalue	dense full eigenvalue problem
Convergence	depends on dominant gap	improved by shifts and deflation
Eigenvectors	one or few	many or all

Power iteration is simple and sparse-friendly. QR is more complete and is the standard dense eigenvalue method.

94.28 Comparison with Krylov Methods

Krylov methods such as Arnoldi and Lanczos are preferred for large sparse matrices when only a few eigenvalues are needed.

The QR algorithm is preferred for dense matrices when most or all eigenvalues are needed.

Setting	Typical choice
Dense matrix, all eigenvalues	QR algorithm
Symmetric dense matrix	tridiagonal reduction plus QR or related methods
Sparse matrix, few eigenvalues	Lanczos or Arnoldi
Matrix-free operator	Krylov method
Full Schur form	QR algorithm

Thus QR and Krylov methods serve different computational regimes.

94.29 Failure Modes and Limitations

The QR algorithm is robust, but practical implementations must handle several issues.

Issue	Handling
Slow convergence	use shifts
Complex eigenvalues	use real Schur $2\times 2$ blocks
High per-step cost	reduce to Hessenberg or tridiagonal form
Completed eigenvalues	deflate
Roundoff	use orthogonal transformations
Poor scaling	balance matrix before iteration

The naive unshifted QR algorithm is mainly pedagogical. Practical QR algorithms use Hessenberg reduction, implicit shifts, and deflation.

94.30 Summary

The QR algorithm computes eigenvalues by repeatedly applying QR factorization and reversing the factors:

A_k=Q_kR_k, \qquad A_{k+1}=R_kQ_k.

Each step is an orthogonal similarity transformation:

A_{k+1}=Q_k^TA_kQ_k.

Therefore the eigenvalues are preserved.

Under suitable conditions, the iterates approach Schur form. For symmetric matrices, they approach diagonal form. In practical implementations, the algorithm is made efficient by:

Technique	Role
Hessenberg reduction	lowers cost for general matrices
Tridiagonal reduction	lowers cost for symmetric matrices
Shifts	accelerate convergence
Deflation	splits completed blocks
Implicit QR	avoids explicit expensive steps
Orthogonal transformations	provide numerical stability

The QR algorithm is the classical dense-matrix method for computing the full eigenvalue spectrum.