QR and Cholesky Factorizations in Statistical Algorithms

Why Factorizations Matter in Statistical Computation

Regression, likelihood calculations, and other statistical algorithms often require solving a matrix system. The algebraic solution may be straightforward, but the route used to compute it matters: finite-precision arithmetic can magnify errors, especially when columns of a design matrix are nearly linearly dependent. QR and Cholesky factorization provide structured ways to solve two common problems without treating every system as an unstructured collection of equations.

QR is the general-purpose choice for least-squares problems with a full-column-rank design matrix. Cholesky is efficient when the system matrix is symmetric positive definite. The distinction is important: a design matrix may be rectangular, while Cholesky applies to a square positive-definite matrix. Forming normal equations to turn a least-squares problem into a square system can be less accurate than working with the design matrix directly.

QR Factorization for Least Squares

Definition: For a matrix \(A\in\mathbb{R}^{m\times p}\) with \(m\ge p\) and full column rank, a reduced QR factorization is \(A=QR\), where \(Q\in\mathbb{R}^{m\times p}\) has orthonormal columns (\(Q^TQ=I_p\)) and \(R\in\mathbb{R}^{p\times p}\) is upper triangular with nonzero diagonal entries. A column-pivoted QR factorization is \(A\Pi=QR\), where \(\Pi\) is a permutation matrix that reorders the columns.

Given a response vector \(b\in\mathbb{R}^m\), least squares seeks a vector \(x\in\mathbb{R}^p\) minimizing \(\|b-Ax\|_2\). If \(A=QR\), then orthogonal decomposition gives

\|b-Ax\|_2^2 =\|Q^Tb-Rx\|_2^2+\|(I-QQ^T)b\|_2^2.

The second term does not depend on \(x\). Since \(R\) is nonsingular, the unique minimizer therefore solves \(Rx=Q^Tb\). This is a triangular system, solved by back substitution. The derivation avoids explicitly forming \(A^TA\).

Worked Example: Solving a Least-Squares Problem with QR

Take \(A=\begin{pmatrix}1&1\\1&0\\0&1\end{pmatrix}\) and \(b=\begin{pmatrix}2\\1\\2\end{pmatrix}\). The columns of \(A\) are linearly independent. Applying orthogonalization gives

q_1=\frac{1}{\sqrt{2}}\begin{pmatrix}1\\1\\0\end{pmatrix}, \qquad q_2=\frac{1}{\sqrt{6}}\begin{pmatrix}1\\-1\\2\end{pmatrix}, \qquad R=\begin{pmatrix}\sqrt{2}&1/\sqrt{2}\\0&\sqrt{3/2}\end{pmatrix}.

Indeed, \(q_1^Tq_1=q_2^Tq_2=1\), \(q_1^Tq_2=0\), and multiplying \(Q=(q_1,q_2)\) by \(R\) recovers the two columns of \(A\). Next, \(Q^Tb=(3/\sqrt{2},5/\sqrt{6})^T\). Back substitution gives

\sqrt{3/2}\,x_2=\frac{5}{\sqrt{6}} \quad\Longrightarrow\quad x_2=\frac{5}{3}, \qquad \sqrt{2}\,x_1+\frac{x_2}{\sqrt{2}}=\frac{3}{\sqrt{2}} \quad\Longrightarrow\quad x_1=\frac{2}{3}.

The fitted vector is \(Ax=(7/3,2/3,5/3)^T\), so the residual is \((-1/3,1/3,1/3)^T\), with squared norm \(1/3\). It is orthogonal to both columns of \(A\), as a least-squares residual must be.

Theorem: If \(A\) has full column rank and \(A=QR\) is a reduced QR factorization, the unique least-squares solution is the solution to \(Rx=Q^Tb\).

Proof. The columns of \(Q\) are orthonormal, so \(QQ^T\) is the orthogonal projection onto their span. Decompose \(b-Ax=(I-QQ^T)b+Q(Q^Tb-Rx)\). The two terms are orthogonal: the first is orthogonal to every column of \(Q\), and the second belongs to their span. The Pythagorean identity yields the displayed sum of squared norms. The first term is fixed as \(x\) varies, and the second is minimized to zero by solving \(Rx=Q^Tb\). Full column rank makes \(R\) nonsingular, so this solution is unique. \(\square\)

Pivoting and Nearly Dependent Columns

Column pivoting reorders columns so that the algorithm selects a column with a large remaining component at each step. In practice, pivoted QR is useful when columns have very different scales or some are nearly linear combinations of others. A small diagonal entry of \(R\) can flag a direction that contributes little new information relative to earlier selected columns. Pivoting can make this structure easier to see; it does not make an intrinsically ill-conditioned problem well-conditioned.

Householder transformations are commonly used to compute QR: each applies an orthogonal transformation to zero entries below a selected diagonal position. Because orthogonal transformations preserve Euclidean lengths, Householder QR is generally preferred to explicitly forming normal equations when numerical accuracy matters. Column-pivoted Householder QR adds the reordering described above.

Worked Example: A Nearly Dependent Pair of Columns

Consider \(A=\begin{pmatrix}1&1\\0&0.01\end{pmatrix}\). Its columns are \(a_1=(1,0)^T\) and \(a_2=(1,0.01)^T\). Their norms are \(1\) and \(\sqrt{1.0001}\), so a pivoting step can place \(a_2\) first. After projecting \(a_1\) onto \(a_2\), the remaining component has length

\left\|a_1-\frac{a_2^Ta_1}{a_2^Ta_2}a_2\right\|_2 =\sqrt{\|a_1\|_2^2-\frac{(a_2^Ta_1)^2}{\|a_2\|_2^2}} =\frac{0.01}{\sqrt{1.0001}}.

Thus, after this ordering, the second diagonal entry of \(R\) is small relative to the first: the second direction adds little beyond the first. For a concrete least-squares solution, take \(b=(2,0.0101)^T\). The equations \(x_1+x_2=2\) and \(0.01x_2=0.0101\) give \(x_2=1.01\) and \(x_1=0.99\). The small second direction means that small changes in the data can have a substantial effect on the fitted coefficients, even though the system has an exact solution.

Cholesky Factorization for Positive-Definite Systems

Definition: A real symmetric matrix \(G\in\mathbb{R}^{p\times p}\) is positive definite if \(z^TGz>0\) for every nonzero \(z\in\mathbb{R}^p\). Its Cholesky factorization is \(G=LL^T\), where \(L\) is lower triangular with strictly positive diagonal entries.

For a system \(Gx=c\), first solve \(Ly=c\) by forward substitution, then solve \(L^Tx=y\) by backward substitution. The factorization is available when \(G\) is symmetric positive definite; symmetry alone is not enough.

Theorem: Every real symmetric positive-definite matrix \(G\) has a unique Cholesky factorization \(G=LL^T\) with \(L\) lower triangular and positive diagonal.

Proof. We use induction on the matrix dimension. In dimension one, positive definiteness means \(G=(g)\) with \(g>0\), and the factor is \(L=(\sqrt{g})\), uniquely. For a larger dimension, partition \(G=\begin{pmatrix}a&r^T\\r&C\end{pmatrix}\). Positive definiteness implies \(a>0\), so set \(L=\begin{pmatrix}\sqrt a&0\\r/\sqrt a&L_2\end{pmatrix}\). The Schur complement \(C-rr^T/a\) is positive definite: for every nonzero \(z\), evaluate the quadratic form of \(G\) at \((-r^Tz/a,z^T)^T\) to obtain \(z^T(C-rr^T/a)z>0\). By induction it has a unique factorization \(L_2L_2^T\) with positive diagonal. Multiplication then gives \(LL^T=G\). The first diagonal entry and first column of any such factor are forced by \(a\) and \(r\); induction forces the remaining factor, proving uniqueness. \(\square\)

Worked Example: Solving a Positive-Definite System with Cholesky

Let \(G=\begin{pmatrix}4&2&2\\2&5&1\\2&1&3\end{pmatrix}\) and \(c=(6,8,5)^T\). The factor

L=\begin{pmatrix}2&0&0\\1&2&0\\1&0&\sqrt{2}\end{pmatrix}

satisfies \(LL^T=G\): its diagonal entries give \(4,5,3\), and its off-diagonal entries give \(2,2,1\). Solve \(Ly=c\) first. The equations give \(y_1=3\), \(y_2=(8-3)/2=5/2\), and \(y_3=(5-3)/\sqrt{2}=\sqrt{2}\). Now solve \(L^Tx=y\): \(x_3=1\), \(x_2=(5/2-0\cdot x_3)/2=5/4\), and \(2x_1+x_2+x_3=3\), so \(x_1=3/8\). Direct substitution verifies the answer: \(G(3/8,5/4,1)^T=(6,8,5)^T\).

Why Normal Equations Can Lose Accuracy

The normal equations for least squares are \(A^TAx=A^Tb\). In exact arithmetic, they encode the same minimizer as QR when \(A\) has full column rank. But they replace \(A\) by a cross-product matrix whose condition number is squared. To see why, let the singular values of full-column-rank \(A\) be \(\sigma_1\ge\cdots\ge\sigma_p>0\). The eigenvalues of \(A^TA\) are \(\sigma_i^2\). Therefore, in the Euclidean norm,

\kappa_2(A^TA) =\frac{\sigma_1^2}{\sigma_p^2} =\left(\frac{\sigma_1}{\sigma_p}\right)^2 =\kappa_2(A)^2.

This identity is the key numerical warning: if the original problem is ill-conditioned, forming \(A^TA\) makes its conditioning worse. In floating-point computation, rounding errors introduced while forming and solving the normal equations can then have a larger effect on the answer. QR avoids this extra squaring, although it cannot remove sensitivity that is already present in the underlying least-squares problem.

Worked Example: Comparing the Condition Numbers

For the nearly dependent matrix in the pivoting example, \(A^TA=\begin{pmatrix}1&1\\1&1.0001\end{pmatrix}\). Its eigenvalues are

\lambda_{\pm} =\frac{2.0001\pm\sqrt{(0.0001)^2+4}}{2} \approx 2.00005,\;0.00005.

Consequently, \(\kappa_2(A^TA)\approx 2.00005/0.00005\approx 40002\), while \(\kappa_2(A)\approx\sqrt{40002}\approx200\). These values agree with the squared-condition-number identity. The example does not imply that every computed solution from normal equations must fail; it shows that this route is more exposed to rounding error when the design is poorly conditioned.

Choosing the Factorization

Problem	Appropriate method	Important condition
Least squares with a rectangular design matrix	QR, often Householder QR	Full column rank for a unique coefficient solution
Least squares where columns may be nearly dependent	Column-pivoted QR can help reveal the structure	Pivoting diagnoses ordering and small directions; it does not eliminate ill-conditioning
Square system \(Gx=c\)	Cholesky	\(G\) must be symmetric positive definite

A common mistake is to assume that every least-squares problem should be solved by forming \(A^TA\), because the normal equations are short to write. The correction is to distinguish algebra from computation: the equations are equivalent in exact arithmetic, but their numerical conditioning is not. Another mistake is to apply Cholesky merely because a matrix is symmetric. Cholesky requires positive definiteness; if that condition fails, the positive-diagonal factor described here does not exist.

Key takeaway: Use QR to solve full-rank least-squares problems without squaring the design matrix's condition number. Use Cholesky for symmetric positive-definite systems, and use column pivoting when the column structure needs to be exposed.

Check Your Understanding

Use the definitions and derivations above to answer these questions.

Why does the least-squares solution from a reduced QR factorization require solving an upper-triangular system?
What does a small later diagonal entry in a column-pivoted \(R\) suggest about the columns of the design matrix?
What two systems are solved after a Cholesky factorization \(G=LL^T\) is obtained?
If \(\kappa_2(A)=50\) and \(A\) has full column rank, what is \(\kappa_2(A^TA)\)?
Why is symmetry alone insufficient to guarantee that the Cholesky factorization described here exists?