LU and Cholesky factorization

Deep dive

How it works

In: square matrix A (SPD for Cholesky), one or more right-hand sides b · Out: x with Ax = b (and the factors P, L, U or L) · State: the factors, reusable for further right-hand sides · Cost: LU 2n³/3 flops, Cholesky n³/3, each solve 2n²

LU (right-looking, column k at a time):

1. Find p = argmax_{i≥k} |A[i][k]|; swap rows k and p, record perm[k] = p.
2. For i > k: L[i][k] = A[i][k] / A[k][k] (stored in place of A[i][k]).
3. Update trailing block: A[i][j] −= L[i][k]·A[k][j] for i, j > k (a rank-1 update).
4. Solve: apply the permutation to b, forward substitution with L (unit diagonal), back substitution with U.

Cholesky (row by row, lower triangle): L[j][j] = √(A[j][j] − Σ_{k<j} L[j][k]²), L[i][j] = (A[i][j] − Σ_{k<j} L[i][k]·L[j][k]) / L[j][j] for i > j.

Fast implementations block the rank-1 updates into rank-b updates (b ≈ 64-256) so almost all flops (all but O(b/n)) run in a matrix multiply (tiled matmul), which reaches near-peak throughput. Unblocked code is memory-bound.

Code

/// Row-major dense matrix.
#[derive(Clone)]
pub struct Mat { pub n: usize, pub a: Vec<f64> }
impl Mat {
    pub fn at(&self, i: usize, j: usize) -> f64 { self.a[i * self.n + j] }
}

/// In-place LU with partial pivoting: A = P L U. Returns row permutation, or None if singular.
pub fn lu_factor(m: &mut Mat) -> Option<Vec<usize>> {
    let n = m.n;
    let a = &mut m.a;
    let mut perm: Vec<usize> = (0..n).collect();
    for k in 0..n {
        let p = (k..n).max_by(|&i, &j| a[i * n + k].abs().total_cmp(&a[j * n + k].abs())).unwrap();
        if a[p * n + k] == 0.0 { return None; }
        if p != k {
            for j in 0..n { a.swap(k * n + j, p * n + j); }
            perm.swap(k, p);
        }
        let pivot = a[k * n + k];
        let (top, rest) = a.split_at_mut((k + 1) * n);
        let row_k = &top[k * n..];
        for row in rest.chunks_exact_mut(n) {
            let l = row[k] / pivot;
            row[k] = l;
            if l != 0.0 {
                for j in k + 1..n { row[j] -= l * row_k[j]; } // contiguous: vectorizes
            }
        }
    }
    Some(perm)
}

pub fn lu_solve(lu: &Mat, perm: &[usize], b: &[f64]) -> Vec<f64> {
    let n = lu.n;
    let mut x: Vec<f64> = perm.iter().map(|&p| b[p]).collect();
    for i in 0..n { for j in 0..i { x[i] -= lu.at(i, j) * x[j]; } }            // L y = P b
    for i in (0..n).rev() {                                                     // U x = y
        for j in i + 1..n { x[i] -= lu.at(i, j) * x[j]; }
        x[i] /= lu.at(i, i);
    }
    x
}

/// In-place Cholesky, lower triangle holds L. Err(j) if the j-th pivot is not positive.
pub fn cholesky(m: &mut Mat) -> Result<(), usize> {
    let n = m.n;
    let a = &mut m.a;
    for j in 0..n {
        let d = a[j * n + j] - (0..j).map(|k| a[j * n + k] * a[j * n + k]).sum::<f64>();
        if !(d > 0.0) { return Err(j); } // also catches NaN
        let d = d.sqrt();
        a[j * n + j] = d;
        for i in j + 1..n {
            let s: f64 = (0..j).map(|k| a[i * n + k] * a[j * n + k]).sum();
            a[i * n + j] = (a[i * n + j] - s) / d;
        }
        for i in 0..j { a[i * n + j] = 0.0; } // clear the upper triangle
    }
    Ok(())
}

pub fn cholesky_solve(l: &Mat, b: &[f64]) -> Vec<f64> {
    let n = l.n;
    let mut x = b.to_vec();
    for i in 0..n { for k in 0..i { x[i] -= l.at(i, k) * x[k]; } x[i] /= l.at(i, i); }      // L y = b
    for i in (0..n).rev() { for k in i + 1..n { x[i] -= l.at(k, i) * x[k]; } x[i] /= l.at(i, i); } // L^T x = y
    x
}

/// One step of iterative refinement: r = b - A x (ideally in higher precision), x += A^{-1} r.
pub fn refine(a: &Mat, lu: &Mat, perm: &[usize], b: &[f64], x: &mut [f64]) {
    let n = a.n;
    let r: Vec<f64> = (0..n).map(|i| b[i] - (0..n).map(|j| a.at(i, j) * x[j]).sum::<f64>()).collect();
    let d = lu_solve(lu, perm, &r);
    for i in 0..n { x[i] += d[i]; }
}

Complexity and constants

• Flops: LU 2n³/3, Cholesky n³/3, LDLᵀ (symmetric indefinite, Bunch-Kaufman) n³/3, QR (Householder) 4n³/3.
• Measured with the code above (unblocked, one Apple-silicon core, -O): LU at n = 1000 takes 0.057 s, n = 2000 0.45 s (about 12 GFLOP/s); the dot-product Cholesky reaches only about 4-5 GFLOP/s. Residuals on random matrices: 2·10⁻¹³ at n = 200, 2.5·10⁻¹⁴ after one refinement step. A blocked, multithreaded LAPACK is far faster on the same machine because GEMM reuses cache.
• Memory: n² doubles; n = 20,000 is 3.2 GB.
• Backward error: computed x solves (A + δA)x = b with ‖δA‖ ≈ growth·n·ε·‖A‖; forward error ≈ κ(A)·ε, so with κ = 10¹⁰ expect only about 6 correct digits in f64.

Efficiency

• Cost: LU 2n³/3 flops, Cholesky n³/3, a solve 2n² each; memory n² doubles (n = 20,000 is 3.2 GB). Measured with the code above (unblocked, one Apple-silicon core): LU at n = 1000 takes 0.057 s, n = 2000 0.45 s (about 12 GFLOP/s); the dot-product Cholesky reaches only 4-5 GFLOP/s. Bound: memory bandwidth when unblocked, compute (GEMM peak) when blocked.
• Levers:
  • blocking into rank-b updates (b ≈ 64-256) — nearly all flops run in matrix multiply at near-peak throughput (Golub-Van Loan 2013) — dense n of a few hundred and up
  • Cholesky instead of LU — half the flops and no pivoting — SPD matrices
  • factor once, solve many — each extra right-hand side costs 2n² instead of n³ — repeated solves, Newton with a frozen Jacobian
  • mixed-precision factorization with f64 refinement — up to 2× faster on CPUs, more on tensor cores (page above) — κ·ε_low < 1
  • fill-reducing ordering (AMD, nested dissection) — keeps sparse factors small (page above) — sparse A
• Naive → standard → state of the art:
  • explicit A⁻¹·b · about 3× the work of a solve, larger error
  • unblocked LU with partial pivoting · 12 GFLOP/s at n = 2000 (measured above)
  • blocked, multithreaded LAPACK or GPU (cuSOLVER, MAGMA), mixed precision · far faster on the same machine

Implementation notes

• Row-major with row operations (as above) keeps the inner loop contiguous; column-major LAPACK does the same with columns.
• Blocking: factor a panel of b columns unblocked, then update the trailing matrix with one GEMM (dgetrf, dpotrf); recursive LU is a simple alternative.
• Determinant: product of U's diagonal times (−1)^(number of swaps); use log-sum to avoid overflow.
• Mixed precision: factor in f32 (or f16 on GPUs), refine in f64; converges to f64 accuracy when κ·ε_low < 1 and is up to 2× faster on CPUs, more on tensor cores.
• Cholesky as SPD test: cheaper than eigenvalues; add a small jitter (A + δI) for Gaussian-process kernels that are SPD only in exact arithmetic.
• Sparse direct: reorder (AMD, nested dissection via METIS) to limit fill-in, symbolic factorization, then numeric (supernodal or multifrontal).

Pitfalls

• No pivoting in LU: [[1e-20, 1], [1, 1]] gives a completely wrong answer; always pivot unless the matrix is SPD or diagonally dominant.
• Testing for singularity with == 0.0: nearly singular matrices pass and return garbage; check the condition estimate (dgecon) or pivot ratio.
• Computing A⁻¹·b explicitly: 3× the work and larger error.
• Normal equations AᵀA x = Aᵀb for least squares: squares κ; fine for κ(A) < 10⁴ or so in f64, otherwise QR.
• Calling Cholesky on a matrix that is symmetric only up to rounding (e.g. assembled covariance): symmetrize first, (A + Aᵀ)/2.

Variants and alternatives

In production

Libraries: LAPACK (dgetrf/dgetrs, dpotrf/dpotrs) via OpenBLAS, Intel MKL, Apple Accelerate; cuSOLVER and MAGMA on GPUs; Eigen (PartialPivLU, LLT), nalgebra and faer in Rust; NumPy/SciPy linalg.solve and cho_factor. Benchmarking: HPL (the TOP500 benchmark) is a blocked LU. Machine learning: Gaussian-process regression factors the kernel matrix with Cholesky (Rasmussen and Williams, Algorithm 2.1). Finance: Cholesky factors turn independent normals into correlated asset draws in Monte Carlo (Glasserman).

Implementations

Sources

• Gene H. Golub, Charles F. Van Loan, Matrix Computations, 4th ed., Johns Hopkins University Press, 2013. book
• Nicholas J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002. book
• Lloyd N. Trefethen, David Bau III, Numerical Linear Algebra, SIAM, 1997. book
• E. Anderson et al., LAPACK Users' Guide, 3rd ed., SIAM, 1999. book · code
• Carl E. Rasmussen, Christopher K. I. Williams, Gaussian Processes for Machine Learning, MIT Press, 2006. book
• Paul Glasserman, Monte Carlo Methods in Financial Engineering, Springer, 2004. book