QR factorization

The original problem: imperfect observations share unknowns

In his 1805 appendix on least squares, Legendre considers observation equations whose known coefficients vary while the unknown quantities remain shared. Usually the equations cannot all hold exactly. He proposes minimizing the sum of squared errors. See the translated original, Smith, pp. 576–577. This supplies the historical fitting problem used here; it is not a claim that Legendre introduced QR.

Original symbolic problem, modern matrix re-derivation. Collect $m$ observation equations for $n$ unknown parameters into $Ax\approx b$. The real matrix $A\in\mathbb R^{m\times n}$ stores the known coefficients, $x\in\mathbb R^n$ the parameters, and $b\in\mathbb R^m$ the observations. The residual $r=b-Ax$ contains the discrepancies. For equal observation weights, solve

Legendre writes an error as $E=a+bx+cy+fz+\cdots$. Our matrix notation moves his known constant term to the other side; our residual has the opposite sign, which leaves its square unchanged. The symbolic model specifies no numerical observation table here. The numerical example below is a labeled teaching reconstruction, not Legendre’s measured data.

The LU page develops the normal equations $A^TAx=A^Tb$. They characterize the best fit, but forming $A^TA$ can worsen numerical sensitivity. A second route first changes to perpendicular model directions. Gram–Schmidt constructs these directions; QR records them as a matrix product.

Build the factors from the directions

Assume initially that $m\ge n$ and that the columns $a_1,\ldots,a_n$ of $A$ are independent: none is a linear combination of the others. Gram–Schmidt produces mutually perpendicular unit vectors $q_i$ and coefficients satisfying

The $j$th input uses only the first $j$ new directions. Put the $q_i$ into columns of $Q$, and their coefficients into $R$:

Here $I_n$ is the identity matrix, with ones on the diagonal and zeros elsewhere. $R$ is upper triangular: entries below the diagonal vanish. This is the reduced QR factorization, with $Q$ of size $m\times n$ and $R$ of size $n\times n$. See LAPACK’s definition and ALAFF’s column-by-column construction.

What each factor does

$R$ converts parameter coordinates into coordinates along the orthonormal model directions. $Q$ places those coordinates in observation space. Because $Q^TQ=I_n$,

Thus $Q$ preserves lengths of vectors entering from $\mathbb R^n$. For rectangular $Q$, this does not say that $Q^T$ preserves every vector in $\mathbb R^m$: perpendicular residual components are discarded. For example, $Q=(1,0)^T$ has $Q^TQ=1$, but $Q^T(0,1)^T=0$.

Shared worked example: the Gram–Schmidt record becomes QR

Supplementary invented example. Continue the Gram–Schmidt worked vectors:

The first column is $\sqrt3q_1$. The second is $\sqrt3q_1+\sqrt2q_2=(0,1,2)^T$. Multiplying verifies $A=QR$, while taking dot products verifies $Q^TQ=I_2$.

The columns of $A$ express a straight-line model $c+dt$ at dimensionless sample positions $t=0,1,2$. Their overlap couples the intercept $c$ and slope $d$. The orthonormal columns instead describe a constant component and a centered trend.

Derive the least-squares solve

Let $y=Q^Tb$ and split the observations into a component in the column space and a perpendicular residual:

The column space contains all possible predictions $Ax$. Its orthogonal projection is $QQ^Tb$. For any parameter vector $x$,

The two terms are perpendicular, so their squared lengths add:

The second term is unavoidable with this model. Independent columns make $R$ nonsingular, so the first term can be made zero by solving

Solve the triangular system by back substitution: use its last row for the last unknown, then work upward. Do not form $R^{-1}$ merely to perform this solve. This is the reduced-form derivation of LAPACK’s QR least-squares method.

Calculate the fit and check the original equations

Use the invented observations $b=(1,2,2)^T$. They do not lie exactly on a line at $t=0,1,2$. First compute

Writing $x=(c,d)^T$ for intercept and slope, back substitution gives

The resulting predictions and residual are

Check both model directions:

The residual is perpendicular to every allowable change in prediction. A least-squares solution need not satisfy $Ax=b$; it makes the remaining error perpendicular to the model space.

Full QR keeps the leftover directions too

Complete the $n$ columns of $Q$ to an orthonormal basis of all $\mathbb R^m$ by adding $Q_\perp$. Then

This is full QR. The $0$ block has $m-n$ rows and $n$ columns. Now $\widetilde Q$ is square and orthogonal: both $\widetilde Q^T\widetilde Q=I_m$ and $\widetilde Q\widetilde Q^T=I_m$ hold. In these full coordinates,

The upper coordinates determine the fit; the lower ones describe the residual. In the worked example, choose $q_3=(1,-2,1)^T/\sqrt6$. Then $q_3^Tb=-1/\sqrt6$, so $r=q_3(-1/\sqrt6)$ and its squared length is $1/6$. The dimensions and residual interpretation follow LAPACK’s full and reduced forms.

Existence, uniqueness, and dependent columns

For $m\ge n$, QR exists even when some columns are dependent. With independent columns, Gram–Schmidt constructs a reduced QR with positive $r_{jj}$. With dependent columns, select additional perpendicular unit directions when a remainder is zero; $R$ can then have zero diagonal entries. The elementary Gram–Schmidt algorithm that divides by every remainder stops unless modified to handle that case.

For full column rank, the convention $r_{jj}>0$ makes the reduced factors unique. To see why, suppose $Q_1R_1=Q_2R_2$. The factors share the same column space, so $Q_2=Q_1S$ with orthogonal $S=Q_1^TQ_2$. Then $S=R_1R_2^{-1}$ is also upper triangular. An upper triangular orthogonal matrix is diagonal with entries $\pm1$; positive diagonals in both $R_i$ force $S=I$. This is an exact-arithmetic proof. The extra columns in full QR are generally not unique.

Without the sign convention, changing a column’s sign in $Q$ and its corresponding row’s sign in $R$ preserves their product. Different software outputs can therefore both be correct. For $m<n$, full $Q$ is $m\times m$ and the factor on the right is $m\times n$, upper trapezoidal, with zeros below its diagonal. See LAPACK’s rectangular cases.

Computing QR: a factorization has several algorithms

Gram–Schmidt explains the factors by constructing columns. Other methods transform the whole matrix while creating zeros. The following are real-arithmetic formulas; complex matrices use conjugate transposes and unitary factors.

Householder reflections: remove a column tail together

A hyperplane is a flat subspace with one fewer dimension than the surrounding space: a line through the origin in two dimensions, or a plane through the origin in three. A Householder reflection reverses the component perpendicular to that hyperplane and leaves components within it unchanged. Its dedicated page derives the geometry, distinguishes reflection from projection, and explains the numerical construction. For nonzero $v$,

satisfies $H^T=H$ and $H^TH=I$. To map a nonzero column segment $z$ to the first coordinate axis, let

where $e_1=(1,0,\ldots,0)^T$ and $\operatorname{sign}(0)=1$ for this construction. Then $Hz=\alpha e_1$. The sign choice avoids subtracting nearly equal leading quantities in $v$. A zero segment needs no reflection.

Invented two-coordinate calculation. For $z=(3,4)^T$, take $\alpha=-5$ and $v=(8,4)^T$:

For a matrix, apply the first reflector to all rows, the second only to rows $2$ onward, and so on, preserving zeros already created. If $H_k\cdots H_1A=[R;0]$, then $\widetilde Q=H_1\cdots H_k$. The opposite order comes from taking the transpose. LAPACK stores reflector vectors and coefficients so it can apply $Q$ or $Q^T$ without explicitly forming a dense $Q$. See its reflector representation and blocked QR algorithm.

Givens rotations: remove one selected entry

A plane rotation changes two coordinates and leaves all other coordinates alone. For a pair $(a,b)$, let $\rho=\sqrt{a^2+b^2}>0$, $\gamma=a/\rho$, and $\sigma=b/\rho$. Then

For $(3,4)$ the coefficients are $3/5$ and $4/5$. If both entries are zero, use the identity rotation. Here $\rho$ is positive; library routines may choose a different sign convention. Robust software scales the calculation to avoid overflow or underflow in $a^2+b^2$. Applying selected rotations can be useful for local zeroing and updating a factorization. See LAPACK’s DLARTG specification.

Choosing the construction

See the Gram–Schmidt rounding example, LAPACK’s reflector implementation, and DLARTG. The identity $A=QR$ alone does not guarantee that a particular implementation is numerically reliable.

The connection to Cholesky, and the accuracy trap

In exact arithmetic, reduced QR gives

With full column rank and positive diagonal $R$, this makes $R^T$ the lower Cholesky factor of $A^TA$. For our worked matrix,

This algebraic relationship does not make the algorithms interchangeable in finite precision. The spectral condition number $\kappa_2(A)=\sigma_{\max}(A)/\sigma_{\min}(A)$ measures the ratio between the largest and smallest stretches of a full-column-rank matrix. Its Gram matrix satisfies

See matrix norms for stretch and conditioning, and ALAFF, §4.2.5 for this identity. Householder QR avoids explicitly forming this squared-condition matrix. It does not remove sensitivity inherent in the original data or guarantee small parameter errors for every ill-conditioned least-squares problem.

Invented counterexample. Take

The columns are independent. If rounding stores $1+\delta^2$ as $1$, the stored Gram matrix becomes singular even though the stored $A$ still contains a nonzero $\delta$. Direct QR keeps that entry in the triangular factor. The threshold depends on the arithmetic; the lesson is about information lost while forming the Gram matrix.

Column pivoting when rank is uncertain

A permutation matrix $P$ reorders coordinates. Column-pivoted QR computes $AP=QR$, choosing column order to help expose independent directions. A small trailing block of $R$ suggests an effective numerical rank at a chosen scale and tolerance. LAPACK’s pivoted factorization routines themselves do not decide the rank. See LAPACK’s rank-sensitive QR discussion.

If columns are dependent, the best prediction remains unique, but multiple parameter vectors can produce it. A basic pivoted-QR solution need not be the parameter vector of smallest Euclidean length. Use an appropriate complete orthogonal factorization or SVD when a minimum-norm solution is required; see LAPACK’s complete orthogonal factorization.

Factorization versus the QR eigenvalue algorithm

A single QR factorization produces factors of one matrix. The QR algorithm for eigenvalues repeatedly factors a square matrix $A_k=Q_kR_k$ and forms $A_{k+1}=R_kQ_k=Q_k^TA_kQ_k$. This is a change of basis, preserving eigenvalues: numbers $\lambda$ for which $Av=\lambda v$ with some nonzero vector $v$. Practical versions use additional structure and shifts. Francis’s 1961 paper concerns that iterative transformation; it should not be treated as the origin of every QR factorization or orthogonalization method.

Interactive example: separate the fit from the residual

This is a modern teaching reconstruction of the original symbolic fitting problem, with the invented matrix and observations used above. The sample positions $0,1,2$ and equal weights stay fixed; edits to $b$ are hypothetical observations with no physical units. The lab displays $Q^Tb$, the triangular solve, predictions, the residual, and $A^Tr$. Change the middle observation to see both the fit and the leftover discrepancy change. Choose observations exactly on a line to make the residual vanish. Restore the labeled baseline to return to $(1,2,2)^T$.

The calculation above is the no-JavaScript baseline. The lab uses the known orthonormal basis of this fixed design matrix; it is an explanation of QR least squares, not a general-purpose numerical solver.

Practice problems

These exercises and numerical values are original teaching examples.

1. Hand calculation: reuse the factors

Using the worked $A,Q,R$, fit the new observations $b=(2,2,5)^T$. Find the intercept, slope, and residual.

2. Interpretation: a rectangular transpose

For reduced $Q\in\mathbb R^{5\times2}$, is $QQ^T=I_5$? What does $QQ^Tb$ represent?

3. Logic: signs and uniqueness

Let $A=QR$ with independent columns and let $D$ be diagonal with entries $\pm1$. Show that $A=(QD)(DR)$ is also a QR factorization. How does a positive diagonal convention remove this freedom?

4. Strategy: distinguish the factor from its computation

A colleague proposes obtaining QR by forming $A^TA$, factoring it by Cholesky, and solving for $Q$. The columns of $A$ are nearly dependent. Why is the exact identity insufficient to justify this choice?

5. Practical application: unequal observation weights

Suppose the three observations have positive weights $w_1,w_2,w_3$. How can QR minimize $\sum_i w_i(b_i-(Ax)_i)^2$? Can the fixed unweighted lab factors be reused unchanged?

Further learning

• Legendre’s original observation model in translation: connect least squares to inconsistent measurement equations; the next page states the squared-error objective.
• LAPACK: QR factorization: compare reduced and full shapes and follow the triangular least-squares solve into software routines.
• LAPACK: representation of orthogonal matrices: understand how software applies an orthogonal factor without storing its full matrix.
• ALAFF: why normal equations can lose accuracy: derive the squared-condition-number identity.

Related pages

• Householder reflection

• Gram–Schmidt orthogonalization
• Matrix factorization and decomposition
• LU decomposition
• Cholesky decomposition
• Matrix norms

References

1. Legendre, A.-M. (1805). Sur la méthode des moindres quarrés, appendix to Nouvelles méthodes pour la détermination des orbites des comètes, pp. 72–75. Translation by H. A. Ruger and H. M. Walker in D. E. Smith, A Source Book in Mathematics, pp. 576–577. The original symbolic model motivates our modern derivation; the line-fitting numbers are invented.
2. Anderson, E., et al. (1999). LAPACK Users’ Guide, QR factorization, column-pivoted QR, and complete orthogonal factorization.
3. Anderson, E., et al. (1999). LAPACK Users’ Guide, representation of orthogonal or unitary matrices and blocked QR.
4. van de Geijn, R. A., and Myers, M. Advanced Linear Algebra: Foundations to Frontiers, classical Gram–Schmidt and §4.2.5, normal-equation conditioning.
5. LAPACK. DLARTG routine specification. Plane rotations and safeguards for their calculation.
6. Francis, J. G. F. (1961). The QR Transformation: A Unitary Analogue to the LR Transformation—Part 1. The Computer Journal, 4(3), 265–271. This concerns the iterative eigenvalue transformation, distinguished here from a single factorization.