LU, QR, Cholesky and matrix norms

Problem

The eigen-decomposition and SVD explain what a matrix does, but they are expensive and are rarely how a computer solves anything. Most numerical work reduces to two questions. Given a square A and a vector b, find x with Ax = b. Given a tall A with more equations than unknowns, find the x that makes Ax as close to b as possible. The first appears in Newton steps, Gaussian processes and every solve call; the second is linear regression.

The textbook answers, x = A⁻¹b and x = (A^TA)⁻¹A^Tb, are exactly the ones not to compute. Forming an inverse costs more than solving, and forming A^TA squares the condition number, so a problem that floating point can solve to 10 digits comes back with 3, or not at all. Libraries instead factor A once into triangular and orthogonal pieces, which are cheap to solve with and lose almost no accuracy.

Intuition

A triangular system is easy. If L is lower triangular, the first equation has one unknown, the second has two of which one is already known, and so on down: forward substitution, n²/2 multiplications. An upper triangular system is solved the same way from the bottom up, back substitution. An orthogonal matrix Q is even easier, since Q⁻¹ = Q^T. So the whole game is to rewrite A as a product of triangular and orthogonal matrices.

• LU records Gaussian elimination: PA = LU, where P reorders the rows, L holds the multipliers and U is what elimination leaves. It works for any invertible square A.
• Cholesky is LU for a symmetric positive-definite A, where symmetry lets the two factors be transposes of each other: A = LL^T. Half the work, no pivoting.
• QR makes A's columns orthonormal one at a time: A = QR. It works for tall matrices too, which is why it solves least squares.

Elimination on A = [[2, 1, 1], [4, -6, 0], [-2, 7, 2]] first swaps rows 1 and 2, because 4 is the largest entry in the first column. Each amber number is the multiple of a pivot row that was subtracted to put a zero below it, and U is the matrix left at the end. Multiplying L by U gives back the row-swapped A exactly. Computed by lu_factor, a stdlib copy of what torch.linalg.lu_factor returns unpacked.

Mechanics

LU with partial pivoting

Gaussian elimination clears column k below the diagonal by subtracting ℓ_ik = a_ik/a_kk times row k from each row i > k. Storing each ℓ_ik in L, below a unit diagonal, gives A = LU. The pivot a_kk can be zero, or tiny, and dividing by a tiny pivot makes the multipliers huge and swamps the other entries in rounding error. So before clearing column k, partial pivoting swaps in the row with the largest |a_ik|, which keeps every |ℓ_ik| ≤ 1. The swaps collect in a permutation P, and the factorization is PA = LU. It costs about (2/3)n³ operations.

Solving Ax = b then takes two triangular solves: Ly = Pb by forward substitution, then Ux = y by back substitution, each O(n²). Once A is factored, every new right-hand side is cheap. That is the reason to keep the factors (lu_factor then lu_solve) when the same matrix meets many vectors. torch.linalg.solve does both steps for you, and its docstring notes it is faster and more stable than A.inverse() @ B. PyTorch returns P as LAPACK does, as a 1-indexed vector of row swaps: pivots[i] == j means row i was swapped with row j − 1, in that order.

Cholesky for positive-definite matrices

A covariance matrix, a kernel matrix plus noise, a Hessian at a minimum, and any A^TA of full column rank are symmetric positive-definite: x^TAx > 0 for every x ≠ 0. For these, A = LL^T with L lower triangular and a positive diagonal. Entry by entry,

ℓ_jj = √(a_jj − Σ_k<j ℓ_jk²),   ℓ_ij = (a_ij − Σ_k<j ℓ_ikℓ_jk) / ℓ_jj for i > j.

It costs about n³/3, half of LU, needs no pivoting, and is the fastest way to test positive definiteness: if a square root's argument is not positive, the matrix is not positive-definite. cholesky_ex exists to return that as an info code instead of raising, and its docstring calls this the faster way to check. Cholesky also turns independent noise into correlated noise: if z ~ N(0, I) then Lz has covariance L·I·L^T = Σ. That is how MultivariateNormal samples, and how a Gaussian process draws functions and evaluates its likelihood.

Target Σ has var(x₁) = 4, var(x₂) = 10 and covariance 2; 4000 samples from random.Random(7) land within a few percent of each. Every sample cost one triangular matrix-vector product.

QR: Gram–Schmidt and Householder

Gram–Schmidt builds Q one column at a time: take column j of A, subtract its projection onto each earlier q_k (the coefficient r_kj = q_k·a_j), and normalize what is left (its length is r_jj). The coefficients form an upper triangular R with A = QR. For an m × n matrix with m ≥ n, the reduced QR has Q of shape m × n with orthonormal columns and R of shape n × n, which is torch.linalg.qr's default mode='reduced'.

Written the obvious way, classical Gram–Schmidt projects the original column onto every q_k, and rounding errors leave the q's less and less orthogonal on ill-conditioned input. Modified Gram–Schmidt projects the running remainder instead, the same arithmetic in a different order, and loses far less. On the 6 × 6 Hilbert matrix (condition number 1.5e+07), ‖Q^TQ − I‖_F is 7.1e-05 for classical and 2.7e-10 for modified. Libraries go further and use Householder reflections: each step reflects the remaining part of a column onto a coordinate axis with an exactly orthogonal matrix I − 2vv^T/v^Tv, so Q is orthogonal to machine precision whatever the conditioning. LAPACK's geqrf, which PyTorch calls, stores Q as those reflectors rather than as a matrix.

Least squares: QR versus the normal equations

With more equations than unknowns, minimize ‖Ax − b‖₂. Setting the gradient to zero gives the normal equations A^TAx = A^Tb, solvable by Cholesky. They are fast and correct in exact arithmetic, but κ(A^TA) = κ(A)², so they lose twice as many digits as the problem requires. QR avoids forming A^TA: since Q has orthonormal columns, ‖Ax − b‖ is minimized by solving Rx = Q^Tb, a single back substitution, with error proportional to κ(A) for small-residual problems.

At ε = 10⁻⁷, κ(A) = 1.4e+07 and κ(A^TA) ≈ 2·10¹⁴: the normal equations are off by 1.1e-03. At ε = 10⁻⁸, 1 + ε² rounds to 1, A^TA becomes [[1, 1], [1, 1]] and there is nothing left to solve. QR gets [1, −1] at every ε.

Matrix norms and the condition number

A matrix norm measures the size of a matrix. torch.linalg.matrix_norm supports the ones that matter, shown here for the worked example's A = [[2, 1, 1], [4, -6, 0], [-2, 7, 2]]:

The spectral norm is the one most results are stated in: ‖Ax‖ ≤ ‖A‖₂‖x‖, and ‖A‖₂ ≤ ‖A‖_F ≤ ‖A‖_nuc. The condition number κ(A) = ‖A‖‖A⁻¹‖, which for the 2-norm is σ_max/σ_min = 20.83 here, bounds how much a relative error in b can grow in x. A rule of thumb: with float64's 16 digits, expect to lose about log₁₀κ of them. The same norms appear elsewhere in ML: weight decay penalizes the squared Frobenius norm, spectral normalization divides a weight by its spectral norm, the nuclear norm is the convex stand-in for rank in matrix completion, and Muon sets every singular value of its update to about 1.

Where they appear in ML

• Linear and ridge regression: least squares with QR or SVD, or Cholesky on A^TA + λI, where the ridge term λ also fixes the conditioning.
• Gaussian processes and Gaussians: Cholesky of the kernel matrix gives both samples and log det Σ = 2Σ log ℓ_ii, needed in the likelihood.
• Orthogonal initialization: nn.init.orthogonal_ takes the QR of a Gaussian matrix and multiplies Q's columns by the signs of R's diagonal, which makes Q uniformly distributed over orthogonal matrices rather than biased by the QR convention.
• Second-order methods and natural gradients: each step solves a linear system in a Hessian or Fisher matrix, usually by Cholesky with damping.
• Orthogonalized updates: Muon's Newton–Schulz iteration approximates the orthogonal polar factor of the gradient using only matrix products, because QR and SVD are slow on GPUs in bfloat16.

Worked example

Solving a 3 × 3 system with LU

Solve Ax = b with A = [[2, 1, 1], [4, -6, 0], [-2, 7, 2]] and b = [5, -2, 9].

1. Column 1: the largest |entry| is 4, in row 2, so swap rows 1 and 2 (pivot 2). Subtract 0.5 × row 1 from row 2 and -0.5 × row 1 from row 3.
2. Column 2: the remaining entries are 4 and 4; the first is kept (pivot 2, no swap). Subtract 1 × row 2 from row 3.
3. That leaves U = [[4, -6, 0], [0, 4, 1], [0, 0, 1]], with L = [[1, 0, 0], [0.5, 1, 0], [-0.5, 1, 1]] and pivots [2, 2, 3] in PyTorch's 1-indexed form.
4. Permute b to Pb = [-2, 5, 9]. Forward substitution in Ly = Pb gives y = [−2, 5 − 0.5·(−2), 9 + 0.5·(−2) − 6] = [−2, 6, 2].
5. Back substitution in Ux = y gives x₃ = 2, x₂ = (6 − 2)/4 = 1, x₁ = (−2 + 6)/4 = 1, so x = [1, 1, 2].

Cholesky by hand

Factor Σ = [[4, 2, 2], [2, 10, 7], [2, 7, 21]]. Column 1: ℓ₁₁ = √4 = 2, ℓ₂₁ = 2/2 = 1, ℓ₃₁ = 2/2 = 1. Column 2: ℓ₂₂ = √(10 − 1²) = 3, ℓ₃₂ = (7 − 1·1)/3 = 2. Column 3: ℓ₃₃ = √(21 − 1² − 2²) = 4. So L = [[2, 0, 0], [1, 3, 0], [1, 2, 4]]. To solve Σx = [10, 26, 51], forward-substitute Ly = b to get y = [5, 7, 8], then back-substitute L^Tx = y to get x = [1, 1, 2]. Its determinant is (2·3·4)² = 576, and its log-determinant 2(ln 2 + ln 3 + ln 4) = 2 ln 24 ≈ 6.356, read straight off the diagonal of L.

Fitting a line with QR

Fit y ≈ c + mx to the points x = 0, 1, 2, 3, y = [1, 2, 2, 4]. The design matrix has a column of ones and a column of x. Gram–Schmidt gives the first column of Q as [0.5, 0.5, 0.5, 0.5] with r₁₁ = 2 and r₁₂ = 3, and the second as [-0.6708, -0.2236, 0.2236, 0.6708] with r₂₂ = √5 ≈ 2.236. Then Q^Ty = [4.5, 2.012], and back substitution in Rx = Q^Ty gives c = 0.9 and m = 0.9. The normal equations give the same line here, because this matrix is well conditioned.

Implementation

import torch

A = torch.tensor([[2., 1., 1.], [4., -6., 0.], [-2., 7., 2.]])
b = torch.tensor([5., -2., 9.])
torch.linalg.solve(A, b)                       # tensor([1., 1., 2.])

LU, pivots = torch.linalg.lu_factor(A)         # pivots = tensor([2, 2, 3], dtype=torch.int32)
torch.linalg.lu_solve(LU, pivots, b[:, None])  # reuse the factors for any right-hand side
P, L, U = torch.linalg.lu(A)                   # A == P @ L @ U

S = torch.tensor([[4., 2., 2.], [2., 10., 7.], [2., 7., 21.]])
L = torch.linalg.cholesky(S)                   # [[2, 0, 0], [1, 3, 0], [1, 2, 4]]
L, info = torch.linalg.cholesky_ex(-S)         # info > 0: not positive-definite, no exception
z = torch.randn(3, 10_000)
x = torch.linalg.cholesky(S) @ z               # samples with covariance S
torch.cov(x)

X = torch.tensor([[1., 0.], [1., 1.], [1., 2.], [1., 3.]])
y = torch.tensor([[1.], [2.], [2.], [4.]])
Q, R = torch.linalg.qr(X)                      # reduced: Q [4, 2], R [2, 2]; signs may differ
torch.linalg.solve_triangular(R, Q.T @ y, upper=True)   # [[0.9], [0.9]]
torch.linalg.lstsq(X, y).solution              # same, through LAPACK's QR

torch.linalg.matrix_norm(A), torch.linalg.matrix_norm(A, ord=2), torch.linalg.cond(A)

The same pieces from scratch, in plain Python:

import math

def lu_factor(a):                          # partial pivoting; returns order, L, U with A[order] = L U
    n = len(a); u = [list(map(float, r)) for r in a]
    lo = [[0.0] * n for _ in range(n)]; order = list(range(n))
    for k in range(n):
        p = max(range(k, n), key=lambda r: abs(u[r][k]))
        u[k], u[p] = u[p], u[k]; lo[k], lo[p] = lo[p], lo[k]; order[k], order[p] = order[p], order[k]
        for r in range(k + 1, n):
            lo[r][k] = u[r][k] / u[k][k]
            u[r] = [x - lo[r][k] * y for x, y in zip(u[r], u[k])]
    for i in range(n):
        lo[i][i] = 1.0
    return order, lo, u

def forward(lo, b):
    y = []
    for i, row in enumerate(lo):
        y.append((b[i] - sum(row[j] * y[j] for j in range(i))) / row[i])
    return y

def back(u, y):
    n = len(u); x = [0.0] * n
    for i in reversed(range(n)):
        x[i] = (y[i] - sum(u[i][j] * x[j] for j in range(i + 1, n))) / u[i][i]
    return x

def cholesky(a):
    n = len(a); lo = [[0.0] * n for _ in range(n)]
    for i in range(n):
        for j in range(i + 1):
            s = a[i][j] - sum(lo[i][k] * lo[j][k] for k in range(j))
            lo[i][j] = math.sqrt(s) if i == j else s / lo[j][j]   # sqrt of a negative: not PD
    return lo

order, lo, u = lu_factor([[2, 1, 1], [4, -6, 0], [-2, 7, 2]])
b = [5, -2, 9]
back(u, forward(lo, [b[i] for i in order]))   # [1.0, 1.0, 2.0]

Tradeoffs

• Never invert to solve. inv(A) @ b does the LU anyway, then extra work, and is less accurate. Keep the factors when the matrix is reused.
• Normal equations versus QR versus SVD. Normal equations with Cholesky cost about mn² + n³/3 and are fine when A is well conditioned; QR costs about 2mn² and loses only κ(A); the SVD (gelsd) costs several times more but handles rank deficiency by dropping singular values below rcond. On CUDA lstsq only offers 'gels', which assumes full rank.
• Pivoting is not optional. Without it, elimination on [[ε, 1], [1, 1]] divides by ε and wipes out the answer in rounding. pivot=False exists only on CUDA, for matrices known to be safe, such as diagonally dominant ones.
• Positive-definite in theory, not in floating point. A covariance or kernel matrix can have eigenvalues that round below zero. Adding a small jitter εI before Cholesky is standard in Gaussian-process code; check cholesky_ex's info instead of catching exceptions.
• GPU and precision. These factorizations are sequential column by column and run in float32 or float64, not bfloat16. On GPUs, for many small matrices batch them; for one big orthogonalization in low precision, matrix-product iterations such as Newton–Schulz can be faster than QR.
• Non-uniqueness. LU depends on pivot choices and QR on sign conventions, so different devices can return different valid factors. Compare products or solutions, not factors.

Check yourself

Connections

References

• Trefethen and Bau. Numerical Linear Algebra. SIAM, 1997: lectures 7–11 on QR and least squares, 20–23 on LU and Cholesky.
• Golub and Van Loan. Matrix Computations, 4th ed. Johns Hopkins University Press, 2013: chapters 3–5.
• Higham. Accuracy and Stability of Numerical Algorithms, 2nd ed. SIAM, 2002: chapter 19 on QR, 20 on least squares.
• Läuchli. Jordan-Elimination und Ausgleichung nach kleinsten Quadraten. Numerische Mathematik 3, 226–240, 1961.
• Mezzadri. How to generate random matrices from the classical compact groups. 2006. arXiv:math-ph/0609050 (cited by orthogonal_).
• Rasmussen and Williams. Gaussian Processes for Machine Learning. MIT Press, 2006: Algorithm 2.1, Cholesky-based GP regression.
• PyTorch v2.15.0: linalg.lu_factor, linalg.cholesky, linalg.qr, linalg.lstsq, linalg.matrix_norm and nn.init.orthogonal_.