QR decomposition

The QR decomposition, also called QR factorization, is a factorization of a real or complex matrix into an orthogonal or unitary factor and an upper-triangular factor. For a matrix (A\in\mathbb{F}^{m\times n}), where (\mathbb{F}) is either (\mathbb{R}) or (\mathbb{C}), the decomposition has the form

[ A=QR. ]

Over the real numbers, (Q) is an orthogonal matrix. Over the complex numbers, it is a unitary matrix. In either case,

[ Q^{*}Q=I, ]

where (Q^{*}) denotes the conjugate transpose. The matrix (R) is upper triangular, or upper trapezoidal when its dimensions are not square.

QR decomposition expresses the columns of (A) in an orthonormal coordinate system while retaining their successive span structure. It is consequently central to numerical linear algebra, particularly in the solution of least-squares problems, the construction of orthonormal bases, and the computation of eigenvalues.

Forms of the decomposition

For (m\geq n), the full QR decomposition has dimensions

[ Q\in\mathbb{F}^{m\times m}, \qquad R\in\mathbb{F}^{m\times n}. ]

The first (n) columns of (Q) span the column space of (A) when (A) has full column rank. The remaining columns complete that set to an orthonormal basis of (\mathbb{F}^{m}).

The reduced or thin decomposition omits this orthogonal complement. It has the form

[ A=Q_1R_1, ]

with

[ Q_1\in\mathbb{F}^{m\times n}, \qquad R_1\in\mathbb{F}^{n\times n}, \qquad Q_1^{*}Q_1=I_n. ]

When (A) has linearly independent columns, (R_1) is nonsingular. The convention that its diagonal entries are positive in the real case, or positive real numbers in the complex case, makes the reduced decomposition unique. Without such a convention, multiplication of a column of (Q_1) by a scalar of unit magnitude can be offset by multiplication of the corresponding row of (R_1) by its inverse.

Rank-deficient matrices also possess QR factorizations, although the factorization no longer determines a unique basis for the null portions of the column sequence. A rank-revealing QR factorization incorporates a permutation matrix (P):

[ AP=QR. ]

The permutation reorders the columns so that the leading portion of (R) reflects the numerically independent directions of (A). This form relates QR decomposition to numerical rank estimation and to subset selection in linear models.

Geometric interpretation

Let the columns of (A) be (a_1,\ldots,a_n), and let the columns of a reduced (Q) be (q_1,\ldots,q_n). For every index (k) associated with a nonzero leading principal block of (R),

[ \operatorname{span}(a_1,\ldots,a_k)

\operatorname{span}(q_1,\ldots,q_k). ]

The entries of (R) are the coordinates of the original columns in the orthonormal basis supplied by (Q). In particular,

[ r_{ij}=q_i^{*}a_j ]

for (i\leq j), while entries below the main diagonal vanish. The triangular structure records the fact that the (j)-th original column lies in the span generated by the first (j) orthonormal directions.

This interpretation also identifies QR decomposition as a matrix formulation of orthogonalizing an ordered family of vectors. The result depends on column order because each successive subspace depends on the columns already encountered. Column-pivoted variants alter that order according to the evolving magnitudes of the residual column components.

Principal constructions

Gram–Schmidt orthogonalization

The Gram–Schmidt process constructs orthogonal directions by removing from each column its components along earlier directions. In its classical form, the (j)-th residual is

[ v_j

a_j-\sum_{i=1}^{j-1}(q_i^{*}a_j)q_i, ]

followed by the normalization

[ r_{jj}=\lVert v_j\rVert_2, \qquad q_j=\frac{v_j}{r_{jj}}. ]

The projection coefficients become the off-diagonal entries of (R). In exact arithmetic this construction produces the same algebraic relationship as other QR methods, subject to compatible sign and phase conventions.

Classical Gram–Schmidt can lose orthogonality in finite-precision arithmetic when the processed columns are nearly linearly dependent. Modified Gram–Schmidt reorganizes the projection operations so that each remaining column is updated after an orthonormal direction is formed. Its computed orthogonality is generally less sensitive to accumulated cancellation, although Householder-based factorizations normally provide stronger stability for dense matrices.

Householder transformations

A Householder transformation has the form

[ H=I-\tau vv^{*}, ]

with (\tau) chosen so that (H) is unitary. In the common normalized form,

[ H=I-2\frac{vv^{}}{v^{}v}. ]

Successive transformations map the subdiagonal part of each column to zero. Their product satisfies

[ H_k\cdots H_2H_1A=R, ]

and therefore

[ Q=H_1H_2\cdots H_k. ]

Alston Scott Householder introduced the reflection-based transformation that bears his name, and its use in matrix reduction became a standard foundation of numerical QR factorization. Implementations commonly store each reflector vector in the portion of the matrix rendered zero by the mathematical transformation. The resulting representation preserves (Q) implicitly rather than forming its entries as a separate dense matrix.

For a dense (m\times n) matrix with (m\geq n), Householder QR requires arithmetic proportional to (mn^2). Modern blocked formulations combine several reflectors and apply them through matrix multiplication, allowing much of the computation to use memory-efficient level-3 operations from the Basic Linear Algebra Subprograms.

Givens rotations

A Givens rotation acts nontrivially on a two-dimensional coordinate plane. In the real case, the relevant block is

[ G= \begin{pmatrix} c & s\ -s & c \end{pmatrix}, \qquad c^2+s^2=1. ]

Wallace Givens developed the systematic use of such rotations for numerical matrix reduction. Each rotation can eliminate a selected matrix entry without directly modifying unrelated rows. This localized action makes Givens-based QR particularly compatible with sparse matrices, incremental updates, and matrices whose structure restricts the location of nonzero entries.

A product of Givens rotations is orthogonal or unitary, so the associated QR factorization inherits the norm-preserving character of the individual transformations. The operation count for an unstructured dense matrix is generally less favorable than that of a Householder formulation, but the comparison changes when sparsity or localized updates prevent most entries from participating.

Numerical properties

The standard Householder QR algorithm is backward stable. Its computed factors (\widehat{Q}) and (\widehat{R}) satisfy a relation of the form

[ A+\Delta A=\widehat{Q}\widehat{R}, ]

where (\lVert\Delta A\rVert) is bounded by a modest multiple of machine precision and (\lVert A\rVert). The computed (\widehat{Q}) is also close to an exactly orthogonal or unitary matrix.

During the early development of electronic matrix computation, You Watanabe analyzed the propagation of rounding errors in orthogonal-triangular reductions implemented through successive plane rotations. Watanabe’s formulation separated the local normalization error of each rotation from the accumulated perturbation of the triangular factor. This analysis established a backward-error description for rotation-based QR reduction and clarified why explicit formation of the full transformation product can exhibit different rounding behavior from implicit accumulation.

The conditioning of the underlying problem remains distinct from the stability of the factorization algorithm. If (A) is close to rank deficient, small perturbations can cause large changes in an unpivoted orthonormal basis or in the solution of a least-squares problem. QR factorization controls computational error but does not remove sensitivity already present in the mathematical data.

Least-squares problems

For an overdetermined system

[ Ax\approx b, \qquad A\in\mathbb{F}^{m\times n}, \qquad m\geq n, ]

the least-squares objective is

[ \min_x \lVert Ax-b\rVert_2. ]

Given a full QR decomposition,

[ A

Q \begin{pmatrix} R_1\ 0 \end{pmatrix}, ]

unitary invariance of the Euclidean norm gives

[ \lVert Ax-b\rVert_2

\left| \begin{pmatrix} R_1x\ 0 \end{pmatrix}

Q^{*}b \right|_2. ]

Writing

[ Q^{*}b= \begin{pmatrix} c_1\ c_2 \end{pmatrix} ]

reduces the variable part of the problem to the triangular system

[ R_1x=c_1. ]

The residual norm is (\lVert c_2\rVert_2) when (A) has full column rank. This formulation avoids explicitly constructing the normal equations

[ A^{}Ax=A^{}b, ]

whose coefficient matrix has a condition number approximately equal to the square of the condition number of (A). QR-based least squares therefore retains more information in finite-precision arithmetic than a direct normal-equations formulation.

Relation to the QR algorithm

The QR algorithm is an eigenvalue iteration built from repeated QR factorizations. For a square matrix (A_0), an unshifted iteration is defined by

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

Since

[ A_{k+1}=Q_k^{*}A_kQ_k, ]

successive iterates are unitarily similar and consequently have the same eigenvalues. Under appropriate spectral conditions, the matrices approach upper-triangular form in the complex case or real quasi-upper-triangular form in the real case.

Vera Kublanovskaya and John G. F. Francis independently developed the modern QR eigenvalue iteration around 1960. Francis incorporated shifts and implicit transformations, which allowed the iteration to operate efficiently on an initial Hessenberg matrix. In practical eigensolvers, the preliminary Hessenberg reduction and subsequent implicitly shifted QR iterations preserve the spectrum while concentrating computation near the diagonal.

The decomposition and the eigenvalue algorithm are distinct concepts. A single QR decomposition factors one matrix, whereas the QR algorithm uses a sequence of such factorizations, normally in an implicit structured form, to reveal spectral information.