Why Numerical Solvers Use Matrix Decompositions for Linear Regression – The Pleasure of Finding Things Out: A blog by James Triveri

This post is an attempt to shed light on why the familiar closed-form solution for least-squares coefficients is avoided in statistical software, and how matrix decompositions provide more numerically stable alternatives. We begin by deriving the normal equations for the standard multiple linear regression model, then use that result as a foundation for exploring more reliable computational approaches.

The normal linear regression model specifies that the sampling variability around the mean is i.i.d. from a normal distribution:

\[ \epsilon_{1}, \cdots, \epsilon_{n} \sim \mathrm{i.i.d.}\hspace{.25em}\mathrm{normal}(0, \sigma^{2})\\ y_{i} = \beta^{T}x_{i} + \epsilon_{i}. \]

Therefore the specification for the joint conditional PDF of observed data is given by:

\[ \begin{align*} f(y_{1}, \cdots, y_{n}|x_{1}, \cdots, x_{n}, \beta, \sigma^{2}) &= \prod_{i=1}^{n} f(y_{i}|x_{i}, \beta, \sigma^{2})\\ &= (2\pi \sigma^{2})^{-n/2} \mathrm{exp}\big(-\frac{1}{2\sigma^{2}} \sum_{i=1}^{n}(y_{i} - \beta^{T}x_{i})^{2}\big). \end{align*} \]

Alternatively, the joint PDF can be represented in terms of the multivariate normal distribution. Let \(y\) be the n-dimensional response vector, and \(X\) the \(n \times p\) design matrix whose \(i^{th}\) row is \(x_{i}\). We have:

\[ y|X,\beta, \sigma^{2} \sim \mathrm{MVN}(X\beta, \sigma^{2} \mathrm{I}), \]

where \(I\) represents the \(n \times n\) identity matrix. The density depends on \(\beta\) through the residuals. To find the optimal \(\beta\), we expand the expression for the sum of squared errors:

\[ \begin{align*} \mathrm{SSE} &= \sum_{i=1}^{n} (y_{i} - \beta^{T}x_{i})^2\\ &=(y - X\beta)^{T}(y - X\beta)\\ &=y^{T}y - 2\beta^{T}X^{T}y + \beta^{T}X^{T}X\beta. \end{align*} \]

Computing the first derivative of the last expression above w.r.t. \(\beta\) and setting equal to 0 yields \(-2X^{T}y + 2X^{T}X\beta = 0\), which can be rearranged to obtain the normal equations:

\[ \beta = (X^{T}X)^{-1}X^{T}y. \]

Why isn’t this expression implemented in linear model solvers directly?

The condition number of a matrix is the ratio of maximum-to-minimum singular values (which, for a normal matrix, is the ratio of the maximum-to-minimum absolute value of eigenvalues). Essentially, the condition number tells you how much solving a linear system will magnify noise in the data. It can be thought of as a measure of amplification. The smaller the condition number, the better (the best value being 1).

"""
Demonstrating the equivalence of computing the condition number as the ratio of 
maximum-to-minimum singular values vs. using np.linalg.cond, as well as a 
comparison of the condition numbers of X vs. X^T*X. 
"""
import numpy as np

rng = np.random.default_rng(516)

n = 50

# Random design matrix with first column as intercept.
X = np.column_stack([
    np.ones(n),
    rng.normal(size=(n, 7))
])

# Not used here, will be used in subsequent code cells.
y = rng.normal(scale=5, size=n)

# SVD for X. 
U0, S0, Vt0 = np.linalg.svd(X, full_matrices=True)
c0 = np.linalg.cond(X, p=None)

# SVD for X^T*X. 
U1, S1, Vt1 = np.linalg.svd(X.T @ X, full_matrices=True)
c1 = np.linalg.cond(X.T @ X, p=None)

# S0 and S1 represent the singular values of X and X^T*X.
print(f"S0.max() / S0.min()       : {S0.max() / S0.min():.8f}.")
print(f"Condition number of X     : {c0:.8f}.")
print(f"S1.max() / S1.min()       : {S1.max() / S1.min():.8f}.")
print(f"Condition number of X^T*X : {c1:.8f}.")

S0.max() / S0.min()       : 2.16451700.
Condition number of X     : 2.16451700.
S1.max() / S1.min()       : 4.68513386.
Condition number of X^T*X : 4.68513386.

In terms of numerical precision, forming the Gram matrix \(X^{T}X\) squares the condition number. As an approximation, \(\mathrm{log}_{10}(\mathrm{condition})\) represents the number of digits lost in a given matrix computation. So by merely forming \(X^{T}X\), we’ve doubled the loss of precision in our final result, since

\[ \mathrm{log}_{10}(\mathrm{condition}(X^{T}X)) \approx 2 \times \mathrm{log}_{10}(\mathrm{condition}(X)). \]

If the condition number of \(X\) is small, forming the Gram matrix and solving the system via \(\beta = (X^{T}X)^{-1}X^{T}y\) should be fine. But as the condition number grows, solving the normal equations becomes increasingly unstable, and solutions become increasingly sensitive to numerical error. Statistical software needs to accommodate design matrices with widely varying numerical properties, so they generally rely on more stable methods rather than explicitly solving the normal equations. In the examples that follow, we demonstrate how the least-squares coefficients can be estimated using the Cholesky, QR, and Singular Value decompositions, then highlight the pros and cons of each approach.

The Cholesky Decomposition

The Cholesky decomposition factors a symmetric positive-definite matrix \(X\) as

\[ X=LL^T, \]

where \(L\) is lower triangular and \(L^T\) is upper triangular. In the least-squares setting, Unlike QR and SVD, the Cholesky decomposition is not applied directly to the design matrix \(X\). It is applied to the symmetric \(p\times p\) Gram matrix \(X^TX\):

\[ X^TX = LL^T. \]

If \(X\) has full column rank, then \(X^TX\) is symmetric positive-definite and therefore has a unique Cholesky decomposition:

\[ \begin{aligned} X^TX\beta &= X^Ty,\\ LL^T\beta &= X^Ty. \end{aligned} \]

Instead of computing \((X^TX)^{-1}\), the coefficients can be obtained by solving two triangular systems. First, define an intermediate vector \(z\) and use forward substitution to solve:

\[ Lz = X^Ty. \]

Then use back substitution to solve:

\[ L^T\beta = z \]

Combining the two systems gives:

\[ LL^T\beta = X^Ty, \]

which is equivalent to the normal equations.

Solving triangular systems is computationally efficient, therefore the Cholesky decomposition provides a fast way to solve the normal equations without explicitly computing a matrix inverse.

The next cell demonstrates the use of the Cholesky Decomposition to obtain least squares coefficients.

"""
Obtaining least-squares coefficients via Cholesky decomposition.
"""

from scipy.linalg import solve_triangular


# Normal equations solved directly.
B0 = np.linalg.solve(X.T @ X, X.T @ y)

# Cholesky decomposition.
L = np.linalg.cholesky(X.T @ X)
z = solve_triangular(L, X.T @ y, lower=True)
B_cd = solve_triangular(L.T, z, lower=False)

print(f"B0: {B0}")
print(f"B_cd: {B_cd}")
print(f"np.allclose(B0, B_cd): {np.allclose(B0, B_cd)}")

B0: [ 0.21363461  0.55927513  1.12519467 -0.16427514 -0.21705275  0.43228936
  1.04517888 -1.72699424]
B_cd: [ 0.21363461  0.55927513  1.12519467 -0.16427514 -0.21705275  0.43228936
  1.04517888 -1.72699424]
np.allclose(B0, B_cd): True

Cholesky eliminates the need to invert \(X^TX\), but it does not eliminate the need to form it. It is faster than QR or SVD for well-conditioned, full-rank problems, but comes at the cost of reduced numerical stability. If \(X\) is rank deficient, \(X^TX\) will not be positive definite and the decomposition will fail.

The QR Decomposition

The QR decomposition improves upon Cholesky since it operates directly on \(X\), and avoids the loss of numerical precision as a result of forming \(X^TX\) and squaring its condition number. It factors a matrix \(X\) into the product of an orthonormal matrix \(Q\) and an upper triangular matrix \(R\):

\[ X = QR. \]

The columns of \(Q\) are normalized, mutually orthogonal vectors. The entries of \(R\) describe how the original columns of \(X\) can be reconstructed from the orthonormal basis:

The diagonal elements \(r_{11}, r_{22}, \ldots\) represent the lengths of the residual vectors used to construct and normalize each successive column of \(Q\).
The off-diagonal elements, such as \(r_{12}\), represent the amount of \(\mathbf{q}_1\) contained in the original vector \(\mathbf{x}_2\). \(r_{12} = \mathbf{q}_1^T\mathbf{x}_2\) is the projection scalar removed from \(\mathbf{x}_2\) during orthogonalization.

\(Q\) represents the orthogonal coordinate system, and \(R\) specifies how to combine the orthonormal directions to reconstruct \(X\).

Starting with the normal equations and substituting \(X=QR\), we obtain:

\[ \begin{align*} X^{T}X \beta &= X^{T}y\\ (QR)^{T}(QR) \beta &= (QR)^{T}y\\ R^{T}(Q^{T}Q)R \beta &= R^{T} Q^{T} y\\ R^{T}R \beta &= R^{T} Q^{T} y\\ (R^{T})^{-1} R^{T} R \beta &= (R^{T})^{-1} R^{T} Q^{T} y\\ R \beta &= Q^{T} y. \end{align*} \]

The derivation is valid only when \(X\) has full column rank, so that \(R\) is non-singular. If \(X\) is rank deficient, ordinary back substitution on \(R\beta=Q^Ty\) does not produce a unique solution.

We’ve taken advantage of how transpose distributes over matrix products (\((AB)^{T} = B^{T}A^{T}\)), and the fact that since \(Q\) is orthonormal, \(Q^{T}Q = I\). Then because \(R\) is upper triangular, \(\beta\) can be solved for using back substitution.

In the next cell, the call to np.linalg.qr sets mode to “reduced”. Doing so keeps only the \(n\) orthogonal vectors in \(Q\) that directly span the column space of our original data matrix. We use reduced mode because it gives us the exact same solution as mode = "complete" but avoids redundant math.

"""
Obtaining least-squares coefficients via QR decomposition.
"""

from scipy.linalg import solve_triangular

# QR decomposition.
Q, R = np.linalg.qr(X, mode="reduced")
B_qr = solve_triangular(R, Q.T @ y)

print(f"B0: {B0}")
print(f"B_qr: {B_qr}")
print(f"np.allclose(B0, B_qr): {np.allclose(B0, B_qr)}")

B0: [ 0.21363461  0.55927513  1.12519467 -0.16427514 -0.21705275  0.43228936
  1.04517888 -1.72699424]
B_qr: [ 0.21363461  0.55927513  1.12519467 -0.16427514 -0.21705275  0.43228936
  1.04517888 -1.72699424]
np.allclose(B0, B_qr): True

Using the QR decomposition, we no longer need to explicitly form the Gram matrix \(X^TX\) and avoid squaring the condition number of the original design matrix. We also no longer need to invert a matrix: After decomposing \(X\) into \(QR\), the coefficients can be obtained by solving the triangular system \(R\beta = Q^Ty\) via substitution.

QR is more stable than Cholesky, but the standard approach assumes \(X\) has full column rank, so that \(R\) is non-singular and coefficient estimates are unique. When \(X\) is rank-deficient or nearly rank-deficient, some diagonal elements of \(R\) will be zero or very small, making it difficult or impossible to obtain a unique coefficient vector.

The Singular Value Decomposition provides a more direct and robust way to identify rank deficiency and compute a minimum-norm least-squares solution.

The Singular Value Decomposition

The Singular Value Decomposition (SVD) improves upon QR by handling rank-deficient and severely ill-conditioned design matrices, though at greater computational cost. It is a generalization of the eigendecomposition of a square matrix to any \(n \times p\) matrix. The SVD decomposes a matrix \(X\) into 3 matrices \(X = U \Sigma V^{T}\):

\(U\) is a \(n \times p\) orthogonal matrix (assuming \(X\) is real); columns represent left singular vectors.
\(\Sigma\) is a \(p \times p\) diagonal matrix with diagonal entries representing the singular values of \(X\).
\(V^{T}\) is a \(p \times p\) orthogonal matrix (assuming \(X\) is real); rows represent right singular vectors.

Starting with the normal equations and assuming \(X\) has full column rank (e.g., \(\Sigma\) is invertible), replace \(X\) with \(U \Sigma V^{T}\) and solve for \(\beta\):

\[ \begin{align*} X^{T}X \beta &= X^{T}y\\ (U \Sigma V^{T})^{T}U \Sigma V^{T}\beta &= (U \Sigma V^{T})^{T}y\\ V \Sigma^{T} U^{T} U \Sigma V^{T} \beta &= V \Sigma U^{T} y\\ V \Sigma^{T} \Sigma V^{T} \beta &=V \Sigma U^{T} y\\ V V^{T} \beta &= V \Sigma^{-1} U^{T} y\\ \beta &= V \Sigma^{-1} U^{T} y \end{align*} \]

Since \(\Sigma\) is diagonal, its pseudoinverse \(\Sigma^{-1}\) is obtained by taking the reciprocal of each nonzero singular value. Because \(V\) is orthogonal, \(VV^{T}=I\). This gives the SVD least-squares solution:

\[ \beta = V \Sigma^{-1} U^{T} y. \]

If \(X\) is rank deficient, one or more singular values will be zero. In that case, those values are excluded from the division. This allows SVD to produce the coefficient vector with the smallest Euclidean norm among all possible least-squares solutions.

The regression coefficients can be estimated using SVD as follows:

"""
Obtaining least-squares coefficients via SVD decomposition.
"""

# SVD.
U, S, Vt = np.linalg.svd(X, full_matrices=False)
B_svd = Vt.T @ np.diag(1 / S) @ U.T @ y

print(f"B0: {B0}")
print(f"B_svd: {B_svd}")
print(f"np.allclose(B0, B_svd): {np.allclose(B0, B_svd)}")

B0: [ 0.21363461  0.55927513  1.12519467 -0.16427514 -0.21705275  0.43228936
  1.04517888 -1.72699424]
B_svd: [ 0.21363461  0.55927513  1.12519467 -0.16427514 -0.21705275  0.43228936
  1.04517888 -1.72699424]
np.allclose(B0, B_svd): True

SVD is the most robust approach for ill-conditioned or rank-deficient matrices, but comes at a higher computational cost than Cholesky or QR. We’ll examine this trade-off more closely in the closing section.

An Ill-Conditioned Design Matrix

Let’s estimate the least-squares coefficients using an ill-conditioned design matrix and compare the performance of each method.

"""
Create a matrix with nearly linearly-dependent columns.
"""

rng = np.random.default_rng(516)

n = 100
x = np.linspace(0, 1, n)

# Columns 2 and 3 are nearly identical.
X = np.column_stack([
    np.ones(n),
    x,
    x + 1e-8 * rng.normal(size=n)
])

B_true = np.array([1.0, 2.0, -2.0])
y = X @ B_true

print(f"cond(X):     {np.linalg.cond(X):.2e}")
print(f"cond(X.T@X): {np.linalg.cond(X.T @ X):.2e}")

cond(X):     1.59e+08
cond(X.T@X): 6.17e+16

\(X\) has full column rank, but it is ill-conditioned because two of its columns are nearly linearly dependent. A condition number of \(10^{16}\) is pretty close to the limit of double-precision arithmetic. Lets compare the solutions:


# Explicit normal-equations solution.
try:
    B_ne = np.linalg.inv(X.T @ X) @ X.T @ y
except np.linalg.LinAlgError as e:
    print(f"Normal equations solution failed: {e}")
    B_ne = None

# Cholesky decomposition.
try:
    L = np.linalg.cholesky(X.T @ X)
    B_chol = solve_triangular(
        L.T, 
        solve_triangular(L, X.T @ y, lower=True)
    )
except np.linalg.LinAlgError as e:
    print(f"Cholesky decomposition failed: {e}")
    B_chol = None

# QR decomposition.
try:
    Q, R = np.linalg.qr(X, mode="reduced")
    B_qr = solve_triangular(R, Q.T @ y)
except np.linalg.LinAlgError as e:
    print(f"QR decomposition failed: {e}")
    B_qr = None

try:
    B_svd = np.linalg.lstsq(X, y, rcond=None)[0]
except np.linalg.LinAlgError as e:
    print(f"SVD failed: {e}")
    B_svd = None

print("True coefficients:   ", B_true)
print("Normal inverse:      ", B_ne)
print("Cholesky solution:   ", B_chol)
print("QR solution:         ", B_qr)
print("SVD solution:        ", B_svd)

Normal equations solution failed: Singular matrix
Cholesky decomposition failed: Matrix is not positive definite
True coefficients:    [ 1.  2. -2.]
Normal inverse:       None
Cholesky solution:    None
QR solution:          [ 1.  2. -2.]
SVD solution:         [ 1.  2. -2.]

For SVD, we used np.linalg.lstsq, which used xGELSD from LAPACK. xGELSD solves least-squares problems using a divide-and-conquer SVD strategy. Here is the source.

Forming \(X^TX\) magnified the numerical instability to the point that the normal equations failed with a singular matrix error and Cholesky failed because the computed Gram matrix was not positive definite. QR and SVD operated directly on \(X\), and were able to recover the true coefficients exactly. This highlights why factorization methods that do not form \(X^TX\) are generally preferred for poorly conditioned least-squares problems.

Computational Considerations

The cost of an SVD-based least-squares solver depends on the specific algorithm and on whether singular vectors are accumulated [Gu, et. al 1994]. The particular estimate of \(4mn^2+8n^3\) is a conventional approximate count for a dense \(m \times n\) matrix, with \(m\geq n\), including the singular vectors.

Method	Approximate flop count
Cholesky via normal equations	\(mn^2+\frac{1}{3}n^3\)
Householder QR	\(2mn^2-\frac{2}{3}n^3\)
Full SVD	Approximately \(4mn^2+8n^3\)

Overall, Cholesky is the fastest but least numerically stable. SVD is the most robust but also the most computationally expensive. QR offers a good balance between these extremes, which explains its frequent use in least-squares solvers.

To learn more about how the QR decomposition is used in R specifically, check out A Deep Dive Into How R Fits a Linear Model, an excellent post by Matthew Drury that I found to be very informative.