Cholesky Decomposition

Why Cholesky Decomposition Matters

A simulation may need to solve the same stiffness system for many loads. A statistical model may need repeated solves with a covariance matrix and its log determinant. When the matrix is symmetric positive definite, Cholesky decomposition uses that structure: a single lower triangular factor and its transpose replace a general pair of triangular factors.

The standard form is $A=LL^T$, where $L$ has positive diagonal entries. Solving $Ax=b$ then requires a forward solve with $L$ and a backward solve with $L^T$. Dense Cholesky uses about half the leading factorization arithmetic of general dense LU, while storing only one triangular factor.

This lesson also makes pivoting explicit. Ordinary Cholesky does not require numerical pivoting for an exactly positive definite matrix. Our implementation uses a symmetric diagonal-pivoted variant: it reorders rows and columns together and computes $PAP^T=LL^T$. We will trace two swaps, preserve earlier factor entries, solve two right-hand sides, and distinguish a valid factorization from a numerical rejection.

Core Concepts and Terminology

A real matrix $A$ is symmetric when $A=A^T$. It is positive definite when $v^TAv\gt0$ for every nonzero real vector $v$. Positive diagonal entries are necessary but insufficient. For example, $\begin{pmatrix}1&2\\2&1\end{pmatrix}$ has positive diagonal entries, yet $v=(1,-1)^T$ gives $v^TAv=-2$.

For symmetric positive definite $A\in\mathbb{R}^{n\times n}$, there is a unique lower triangular factor with positive diagonal satisfying $A=LL^T$. With a fixed symmetric permutation, the corresponding factor is also unique. Its diagonal is not generally all ones; the positive square roots are essential.

Our zero-based array $p$ specifies the final variable order. Define $P_{i,p_i}=1$ and all other entries zero. Then $(PAP^T)_{ij}=A_{p_i,p_j}$ and $(Pb)_i=b_{p_i}$. To solve, compute $Ly=Pb$, then $L^Tz=y$, then restore the original order by $x_{p_i}=z_i$. Both the right-hand side and the answer require the appropriate permutation operation.

The input contract is deliberately precise: a nonempty finite square matrix that is exactly symmetric in its supplied floating-point entries, plus a finite relative tolerance $0\lt\tau\lt1$. The code rejects nonpositive or too-small residual pivots. It does not return a low-rank factor for a semidefinite matrix or an indefinite factorization.

The LAPACK DPOTRF reference documents standard lower and upper Cholesky forms. Its interface can read only the selected input triangle; our teaching implementation checks the entire supplied matrix for exact symmetry before working on a copy.

How Cholesky Decomposition Works, Step by Step

First consider the unpivoted formulas. At column $k$, match the entries of $A=LL^T$:

$$ L_{kk}=\sqrt{A_{kk}-\sum_{j=0}^{k-1}L_{kj}^2}, \qquad L_{ik}=\frac{A_{ik}-\sum_{j=0}^{k-1}L_{ij}L_{kj}}{L_{kk}} \quad(i\gt k). $$

Earlier columns are already available, so these expressions determine the next column from top to bottom. The diagonal expression must be positive before taking its square root. Our pivoted implementation uses the same formulas after reordering the remaining variables.

1. Validate and scale. Check the square shape, finite entries, exact symmetry, and tolerance. Set $s=\max_{i,j}\lvert A_{ij}\rvert$ and reject $s=0$. Copy $W=A/s$, initialize the normalized factor $C$ to zero, set $p=(0,1,\ldots,n-1)$, and cache $d_i=W_{ii}$.
2. Select a residual pivot. At column $k$, find the largest $d_i$ for $i\ge k$, keeping the first on an exact tie. Reject if that value is at most $\tau$. These are residual diagonals, not the unchanged original diagonal entries.
3. Swap symmetrically. If the selected index is $r\ne k$, swap rows $k,r$ and columns $k,r$ of $W$. Swap $p_k,p_r$ and $d_k,d_r$. Swap only columns before $k$ between rows $k,r$ of $C$, moving the history already computed for those variables.
4. Take the positive square root. Set $C_{kk}=\sqrt{d_k}$. For each $i\gt k$, compute $C_{ik}=(W_{ik}-\sum_{j\lt k}C_{ij}C_{kj})/C_{kk}$.
5. Update the cached diagonals. For every $i\gt k$, set $d_i\leftarrow d_i-C_{ik}^2$. This maintains $d_i=W_{ii}-\sum_{j=0}^{k}C_{ij}^2$ for the next stage, without recomputing all the diagonal sums.
6. Finish and rescale. Repeat until every column is complete. Return $L=\sqrt{s}\,C$ and $p$. Since $P(A/s)P^T=CC^T$, rescaling gives $PAP^T=LL^T$.
7. Solve and restore order. Evaluate

$$ y_i=\frac{b_{p_i}-\sum_{j=0}^{i-1}L_{ij}y_j}{L_{ii}} \quad\text{from top to bottom}, \qquad z_i=\frac{y_i-\sum_{j=i+1}^{n-1}L_{ji}z_j}{L_{ii}} \quad\text{from bottom to top}. $$

Finally assign $x_{p_i}=z_i$. Empty sums are zero. The second pass uses the transpose of the same factor; neither an explicit inverse nor a separately stored upper factor is needed.

flowchart TD
    A["Validate shape, symmetry, finite data, and tolerance"] --> B["Scale A; initialize factor, order, residual diagonals"]
    B --> C{"More columns?"}
    C -->|Yes| D["Select largest active residual diagonal"]
    D --> E{"Pivot exceeds tolerance?"}
    E -->|No| F["Report rejected positive pivot"]
    E -->|Yes| G["Swap rows and columns, order, diagonals, factor prefix"]
    G --> H["Take positive square root; compute lower column"]
    H --> I["Update residual diagonals; advance column"]
    I --> C
    C -->|No| J["Rescale factor; forward solve with permuted b"]
    J --> K["Back solve transpose; restore variable order"]
    K --> M["Check factor product and original-system residual"]

Worked Example with a Full Trace

Take

$$ A=\begin{pmatrix}10&5&4\\5&9&8\\4&8&16\end{pmatrix}, \qquad b=\begin{pmatrix}32\\47\\68\end{pmatrix}. $$

The code uses $s=16$, stores $C=L/4$ during the computation, and caches residual diagonals divided by $16$. Below, factors and residuals are shown in original units so the integers remain visible. The final product will also independently establish that this symmetric matrix is positive definite.

Initialize. The variable order is $p=(0,1,2)$, the factor is zero, and the original-unit diagonal candidates are $(10,9,16)$.

Stage $k=0$: first symmetric swap. Choose index $2$, whose candidate is $16$. Swap indices $0$ and $2$ in both rows and columns. Now $p=(2,1,0)$ and

$$ B_1=\begin{pmatrix}16&8&4\\8&9&5\\4&5&10\end{pmatrix}. $$

The pivot is $16$, so $L_{00}=4$. Compute $L_{10}=8/4=2$ and $L_{20}=4/4=1$. The next residual diagonal candidates are $9-2^2=5$ and $10-1^2=9$. Thus the partial factor is

$$ L_{\mathrm{partial}}=\begin{pmatrix}4&0&0\\2&0&0\\1&0&0\end{pmatrix}, \qquad d_{\mathrm{active}}=(5,9). $$

Stage $k=1$: move a previously computed entry. Choose residual $9$, at current index $2$. Swap indices $1$ and $2$ symmetrically. Update $p=(2,0,1)$ and swap the completed factor entries $L_{10}=2$ and $L_{20}=1$. Immediately after the swap,

$$ B_2=\begin{pmatrix}16&4&8\\4&10&5\\8&5&9\end{pmatrix}, \qquad L_{\mathrm{partial}}=\begin{pmatrix}4&0&0\\1&0&0\\2&0&0\end{pmatrix}. $$

Take $L_{11}=\sqrt9=3$. Then $L_{21}=(5-L_{20}L_{10})/3=(5-2\cdot1)/3=1$. Update the final residual diagonal to $5-1^2=4$; equivalently, it is $9-2^2-1^2=4$.

Stage $k=2$: last positive pivot. The remaining pivot is $4$, so $L_{22}=2$. All normalized pivot candidates, $1$, $9/16$, and $1/4$, exceed $10^{-12}$. The final factor and its transpose are

$$ L=\begin{pmatrix}4&0&0\\1&3&0\\2&1&2\end{pmatrix}, \qquad L^T=\begin{pmatrix}4&1&2\\0&3&1\\0&0&2\end{pmatrix}. $$

Verify the product. Multiplying gives

$$ LL^T=\begin{pmatrix}16&4&8\\4&10&5\\8&5&9\end{pmatrix} =PAP^T, \qquad P=\begin{pmatrix}0&0&1\\1&0&0\\0&1&0\end{pmatrix}. $$

Since $L$ is invertible, $v^TAv=\lVert L^TPv\rVert_2^2\gt0$ for nonzero $v$. This verifies positive definiteness independently of the input's appearance.

Forward solve. Gather $Pb=(68,32,47)^T$. Compute $y_0=68/4=17$, $y_1=(32-1\cdot17)/3=5$, and $y_2=(47-2\cdot17-1\cdot5)/2=4$.

Transpose solve. Compute $z_2=4/2=2$, $z_1=(5-1\cdot2)/3=1$, and $z_0=(17-1\cdot1-2\cdot2)/4=3$. Thus $z=(3,1,2)^T$ in permuted coordinates. Restore $x_2=z_0=3$, $x_0=z_1=1$, and $x_1=z_2=2$, giving $x=(1,2,3)^T$. Direct multiplication yields $Ax=(32,47,68)^T$.

Reuse and determinant. For $b^{(2)}=(10,5,4)^T$, the same factors give $Pb^{(2)}=(4,10,5)^T$, $y^{(2)}=(1,3,0)^T$, $z^{(2)}=(0,1,0)^T$, and $x^{(2)}=(1,0,0)^T$. Also $\det(A)=(4\cdot3\cdot2)^2=576$. Symmetric permutations do not introduce a determinant sign change because the row and column permutation signs multiply to $1$.

Implementation in Python, C++, and Java

Each panel is a complete program with identical sample data, pivot selection, scaling, and output. Run Python 3 directly, use C++14 or newer, or save the Java program as Main.java. The factor routine preserves $A$ and returns $L$ and a final order array $p$; the solve routine preserves its factors and $b$ and expects an unchanged successful factorization result.

The default tolerance is $\tau=10^{-12}$. It applies to residual diagonals of $A/s$, not to the square-root diagonal of $L$. A positive pivot smaller than the threshold is rejected even when the exact matrix is positive definite. For example, $\operatorname{diag}(1,10^{-14})$ is SPD but is rejected at this default tolerance. The threshold is a numerical policy, not a condition-number estimate.

The exact-symmetry check intentionally avoids silently replacing the user's matrix with a different one. If an application supplies approximately symmetric data, choose and document an appropriate symmetry test or symmetrization policy before calling this routine. A nonsymmetric matrix must not be accepted simply because one of its triangles produces plausible factor entries.

The factor diagnostic is $\max_{i,j}\lvert(LL^T)_{ij}-A_{p_i,p_j}\rvert/s$, while each solve diagnostic is $\lVert Ax-b\rVert_\infty$ in original coordinates. The sample checks both against $10^{-10}$ because its magnitudes are modest. General applications need scale-aware solve criteria. Scaling the factorization reduces some overflow hazards, but the implementation uses ordinary floating-point arithmetic, and the simple diagnostic dot products are intended for moderate inputs.

import math


def cholesky_factor(a, tol=1e-12):
    if not math.isfinite(tol) or not 0.0 < tol < 1.0:
        raise ValueError("tol must be finite and between 0 and 1")
    n = len(a)
    if n == 0 or any(len(row) != n for row in a):
        raise ValueError("A must be nonempty and square")
    if any(not math.isfinite(v) for row in a for v in row):
        raise ValueError("A must contain finite values")
    if any(a[i][j] != a[j][i] for i in range(n) for j in range(i)):
        raise ValueError("A must be exactly symmetric")
    scale = max(abs(v) for row in a for v in row)
    if scale == 0.0:
        raise ValueError("zero matrix is not positive definite")
    w = [[v / scale for v in row] for row in a]
    l = [[0.0] * n for _ in range(n)]
    p = list(range(n))
    diagonal = [w[i][i] for i in range(n)]
    for k in range(n):
        pivot, best = k, -math.inf
        for i in range(k, n):
            d = diagonal[i]
            if not math.isfinite(d):
                raise ValueError("nonfinite residual diagonal")
            if d > best:
                pivot, best = i, d
        if best <= tol:
            raise ValueError("positive pivot rejected at the chosen tolerance")
        if pivot != k:
            w[k], w[pivot] = w[pivot], w[k]
            for row in w:
                row[k], row[pivot] = row[pivot], row[k]
            p[k], p[pivot] = p[pivot], p[k]
            diagonal[k], diagonal[pivot] = diagonal[pivot], diagonal[k]
            for j in range(k):
                l[k][j], l[pivot][j] = l[pivot][j], l[k][j]
        l[k][k] = math.sqrt(best)
        for i in range(k + 1, n):
            value = w[i][k]
            for j in range(k):
                value -= l[i][j] * l[k][j]
            l[i][k] = value / l[k][k]
            if not math.isfinite(l[i][k]):
                raise ValueError("nonfinite factor entry")
            diagonal[i] -= l[i][k] * l[i][k]
            if not math.isfinite(diagonal[i]):
                raise ValueError("nonfinite residual diagonal")
    root_scale = math.sqrt(scale)
    for i in range(n):
        for j in range(i + 1):
            l[i][j] *= root_scale
            if not math.isfinite(l[i][j]):
                raise ValueError("factor is not representable")
        if l[i][i] == 0.0:
            raise ValueError("factor diagonal underflowed")
    return l, p


def cholesky_solve(factors, b):
    # factors must be an unchanged result of cholesky_factor.
    l, p = factors
    n = len(p)
    if len(b) != n or any(not math.isfinite(v) for v in b):
        raise ValueError("b must have n finite entries")
    y = [0.0] * n
    for i in range(n):
        value = b[p[i]]
        for j in range(i):
            value -= l[i][j] * y[j]
        y[i] = value / l[i][i]
        if not math.isfinite(y[i]):
            raise ValueError("nonfinite forward solve")
    z = [0.0] * n
    for i in range(n - 1, -1, -1):
        value = y[i]
        for j in range(i + 1, n):
            value -= l[j][i] * z[j]
        z[i] = value / l[i][i]
        if not math.isfinite(z[i]):
            raise ValueError("nonfinite transpose solve")
    x = [0.0] * n
    for i in range(n):
        x[p[i]] = z[i]
    return x


def factor_error(a, factors):
    l, p = factors
    n = len(p)
    scale = max(abs(v) for row in a for v in row)
    return max(abs(sum(l[i][k] * l[j][k] for k in range(n)) - a[p[i]][p[j]])
               for i in range(n) for j in range(n)) / scale


def solve_error(a, x, b):
    return max(abs(sum(row[j] * x[j] for j in range(len(x))) - b[i])
               for i, row in enumerate(a))


def vector_text(v):
    return "[" + ", ".join(f"{(0.0 if value == 0.0 else value):.6f}" for value in v) + "]"


def matrix_text(a):
    return "[" + ", ".join(vector_text(row) for row in a) + "]"


def main():
    a = [[10.0, 5.0, 4.0], [5.0, 9.0, 8.0], [4.0, 8.0, 16.0]]
    b = [32.0, 47.0, 68.0]
    b2 = [10.0, 5.0, 4.0]
    factors = cholesky_factor(a)
    l, p = factors
    x, x2 = cholesky_solve(factors, b), cholesky_solve(factors, b2)
    fe, se, re = factor_error(a, factors), solve_error(a, x, b), solve_error(a, x2, b2)
    print("permutation = " + str(p))
    print("L = " + matrix_text(l))
    print("x = " + vector_text(x))
    print("reused x = " + vector_text(x2))
    print(f"relative factor error = {fe:.3e}")
    print(f"solve residual = {se:.3e}")
    print(f"reuse residual = {re:.3e}")
    print("check = " + ("PASS" if max(fe, se, re) < 1e-10 else "FAIL"))


if __name__ == "__main__":
    main()

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <limits>
#include <numeric>
#include <sstream>
#include <stdexcept>
#include <string>
#include <vector>

using Matrix = std::vector<std::vector<double>>;
struct Factors { Matrix l; std::vector<std::size_t> p; };

Factors cholesky_factor(const Matrix& a, double tol = 1e-12) {
    if (!std::isfinite(tol) || !(tol > 0.0 && tol < 1.0))
        throw std::invalid_argument("tol must be finite and between 0 and 1");
    const std::size_t n = a.size();
    if (n == 0) throw std::invalid_argument("A must be nonempty");
    double scale = 0.0;
    for (const auto& row : a) {
        if (row.size() != n) throw std::invalid_argument("A must be square");
        for (double v : row) {
            if (!std::isfinite(v)) throw std::invalid_argument("A must be finite");
            scale = std::max(scale, std::abs(v));
        }
    }
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t j = 0; j < i; ++j)
            if (a[i][j] != a[j][i]) throw std::invalid_argument("A must be exactly symmetric");
    if (scale == 0.0) throw std::runtime_error("zero matrix is not positive definite");
    Matrix w = a, l(n, std::vector<double>(n, 0.0));
    for (auto& row : w) for (auto& v : row) v /= scale;
    std::vector<std::size_t> p(n);
    std::iota(p.begin(), p.end(), 0);
    std::vector<double> diagonal(n);
    for (std::size_t i = 0; i < n; ++i) diagonal[i] = w[i][i];
    for (std::size_t k = 0; k < n; ++k) {
        std::size_t pivot = k;
        double best = -std::numeric_limits<double>::infinity();
        for (std::size_t i = k; i < n; ++i) {
            const double d = diagonal[i];
            if (!std::isfinite(d)) throw std::runtime_error("nonfinite residual diagonal");
            if (d > best) { pivot = i; best = d; }
        }
        if (best <= tol) throw std::runtime_error("positive pivot rejected at the chosen tolerance");
        if (pivot != k) {
            std::swap(w[k], w[pivot]);
            for (auto& row : w) std::swap(row[k], row[pivot]);
            std::swap(p[k], p[pivot]);
            std::swap(diagonal[k], diagonal[pivot]);
            for (std::size_t j = 0; j < k; ++j) std::swap(l[k][j], l[pivot][j]);
        }
        l[k][k] = std::sqrt(best);
        for (std::size_t i = k + 1; i < n; ++i) {
            double value = w[i][k];
            for (std::size_t j = 0; j < k; ++j) value -= l[i][j] * l[k][j];
            l[i][k] = value / l[k][k];
            if (!std::isfinite(l[i][k])) throw std::runtime_error("nonfinite factor entry");
            diagonal[i] -= l[i][k] * l[i][k];
            if (!std::isfinite(diagonal[i])) throw std::runtime_error("nonfinite residual diagonal");
        }
    }
    const double root_scale = std::sqrt(scale);
    for (std::size_t i = 0; i < n; ++i) {
        for (std::size_t j = 0; j <= i; ++j) {
            l[i][j] *= root_scale;
            if (!std::isfinite(l[i][j])) throw std::runtime_error("factor is not representable");
        }
        if (l[i][i] == 0.0) throw std::runtime_error("factor diagonal underflowed");
    }
    return {l, p};
}

std::vector<double> cholesky_solve(const Factors& f, const std::vector<double>& b) {
    // f must be an unchanged result of cholesky_factor.
    const std::size_t n = f.p.size();
    if (b.size() != n) throw std::invalid_argument("b must have n entries");
    for (double v : b) if (!std::isfinite(v)) throw std::invalid_argument("b must be finite");
    std::vector<double> y(n), z(n), x(n);
    for (std::size_t i = 0; i < n; ++i) {
        double value = b[f.p[i]];
        for (std::size_t j = 0; j < i; ++j) value -= f.l[i][j] * y[j];
        y[i] = value / f.l[i][i];
        if (!std::isfinite(y[i])) throw std::runtime_error("nonfinite forward solve");
    }
    for (std::size_t i = n; i-- > 0;) {
        double value = y[i];
        for (std::size_t j = i + 1; j < n; ++j) value -= f.l[j][i] * z[j];
        z[i] = value / f.l[i][i];
        if (!std::isfinite(z[i])) throw std::runtime_error("nonfinite transpose solve");
    }
    for (std::size_t i = 0; i < n; ++i) x[f.p[i]] = z[i];
    return x;
}

double factor_error(const Matrix& a, const Factors& f) {
    double error = 0.0, scale = 0.0;
    for (const auto& row : a) for (double v : row) scale = std::max(scale, std::abs(v));
    const std::size_t n = a.size();
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t j = 0; j < n; ++j) {
            double value = 0.0;
            for (std::size_t k = 0; k < n; ++k) value += f.l[i][k] * f.l[j][k];
            error = std::max(error, std::abs(value - a[f.p[i]][f.p[j]]));
        }
    return error / scale;
}

double solve_error(const Matrix& a, const std::vector<double>& x, const std::vector<double>& b) {
    double error = 0.0;
    for (std::size_t i = 0; i < a.size(); ++i) {
        double value = 0.0;
        for (std::size_t j = 0; j < x.size(); ++j) value += a[i][j] * x[j];
        error = std::max(error, std::abs(value - b[i]));
    }
    return error;
}

std::string vector_text(const std::vector<double>& v) {
    std::ostringstream s;
    s << std::fixed << std::setprecision(6) << '[';
    for (std::size_t i = 0; i < v.size(); ++i) {
        if (i) s << ", ";
        s << (v[i] == 0.0 ? 0.0 : v[i]);
    }
    s << ']';
    return s.str();
}

std::string matrix_text(const Matrix& a) {
    std::string s = "[";
    for (std::size_t i = 0; i < a.size(); ++i) {
        if (i) s += ", ";
        s += vector_text(a[i]);
    }
    return s + "]";
}

int main() {
    const Matrix a = { {10.0, 5.0, 4.0}, {5.0, 9.0, 8.0}, {4.0, 8.0, 16.0} };
    const std::vector<double> b = {32.0, 47.0, 68.0}, b2 = {10.0, 5.0, 4.0};
    const auto f = cholesky_factor(a);
    const auto x = cholesky_solve(f, b), x2 = cholesky_solve(f, b2);
    const double fe = factor_error(a, f), se = solve_error(a, x, b), re = solve_error(a, x2, b2);
    std::cout << "permutation = [";
    for (std::size_t i = 0; i < f.p.size(); ++i) {
        if (i) std::cout << ", ";
        std::cout << f.p[i];
    }
    std::cout << "]\nL = " << matrix_text(f.l) << "\nx = " << vector_text(x)
              << "\nreused x = " << vector_text(x2) << std::scientific << std::setprecision(3)
              << "\nrelative factor error = " << fe << "\nsolve residual = " << se
              << "\nreuse residual = " << re
              << "\ncheck = " << (std::max({fe, se, re}) < 1e-10 ? "PASS" : "FAIL") << '\n';
}

import java.util.Arrays;
import java.util.Locale;

public class Main {
    static class Factors {
        final double[][] l;
        final int[] p;
        Factors(double[][] l, int[] p) { this.l = l; this.p = p; }
    }

    static Factors choleskyFactor(double[][] a, double tol) {
        if (!Double.isFinite(tol) || !(tol > 0.0 && tol < 1.0))
            throw new IllegalArgumentException("tol must be finite and between 0 and 1");
        int n = a.length;
        if (n == 0) throw new IllegalArgumentException("A must be nonempty");
        double scale = 0.0;
        for (double[] row : a) {
            if (row == null || row.length != n) throw new IllegalArgumentException("A must be square");
            for (double v : row) {
                if (!Double.isFinite(v)) throw new IllegalArgumentException("A must be finite");
                scale = Math.max(scale, Math.abs(v));
            }
        }
        for (int i = 0; i < n; ++i)
            for (int j = 0; j < i; ++j)
                if (a[i][j] != a[j][i]) throw new IllegalArgumentException("A must be exactly symmetric");
        if (scale == 0.0) throw new IllegalArgumentException("zero matrix is not positive definite");
        double[][] w = new double[n][n], l = new double[n][n];
        int[] p = new int[n];
        double[] diagonal = new double[n];
        for (int i = 0; i < n; ++i) {
            p[i] = i;
            for (int j = 0; j < n; ++j) w[i][j] = a[i][j] / scale;
        }
        for (int i = 0; i < n; ++i) diagonal[i] = w[i][i];
        for (int k = 0; k < n; ++k) {
            int pivot = k;
            double best = Double.NEGATIVE_INFINITY;
            for (int i = k; i < n; ++i) {
                double d = diagonal[i];
                if (!Double.isFinite(d)) throw new IllegalArgumentException("nonfinite residual diagonal");
                if (d > best) { pivot = i; best = d; }
            }
            if (best <= tol) throw new IllegalArgumentException("positive pivot rejected at the chosen tolerance");
            if (pivot != k) {
                double[] row = w[k]; w[k] = w[pivot]; w[pivot] = row;
                for (double[] values : w) {
                    double value = values[k]; values[k] = values[pivot]; values[pivot] = value;
                }
                int index = p[k]; p[k] = p[pivot]; p[pivot] = index;
                double d = diagonal[k]; diagonal[k] = diagonal[pivot]; diagonal[pivot] = d;
                for (int j = 0; j < k; ++j) {
                    double value = l[k][j]; l[k][j] = l[pivot][j]; l[pivot][j] = value;
                }
            }
            l[k][k] = Math.sqrt(best);
            for (int i = k + 1; i < n; ++i) {
                double value = w[i][k];
                for (int j = 0; j < k; ++j) value -= l[i][j] * l[k][j];
                l[i][k] = value / l[k][k];
                if (!Double.isFinite(l[i][k])) throw new IllegalArgumentException("nonfinite factor entry");
                diagonal[i] -= l[i][k] * l[i][k];
                if (!Double.isFinite(diagonal[i])) throw new IllegalArgumentException("nonfinite residual diagonal");
            }
        }
        double rootScale = Math.sqrt(scale);
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j <= i; ++j) {
                l[i][j] *= rootScale;
                if (!Double.isFinite(l[i][j])) throw new IllegalArgumentException("factor is not representable");
            }
            if (l[i][i] == 0.0) throw new IllegalArgumentException("factor diagonal underflowed");
        }
        return new Factors(l, p);
    }

    static double[] choleskySolve(Factors f, double[] b) {
        // f must be an unchanged result of choleskyFactor.
        int n = f.p.length;
        if (b.length != n) throw new IllegalArgumentException("b must have n entries");
        for (double v : b) if (!Double.isFinite(v)) throw new IllegalArgumentException("b must be finite");
        double[] y = new double[n], z = new double[n], x = new double[n];
        for (int i = 0; i < n; ++i) {
            double value = b[f.p[i]];
            for (int j = 0; j < i; ++j) value -= f.l[i][j] * y[j];
            y[i] = value / f.l[i][i];
            if (!Double.isFinite(y[i])) throw new IllegalArgumentException("nonfinite forward solve");
        }
        for (int i = n - 1; i >= 0; --i) {
            double value = y[i];
            for (int j = i + 1; j < n; ++j) value -= f.l[j][i] * z[j];
            z[i] = value / f.l[i][i];
            if (!Double.isFinite(z[i])) throw new IllegalArgumentException("nonfinite transpose solve");
        }
        for (int i = 0; i < n; ++i) x[f.p[i]] = z[i];
        return x;
    }

    static double factorError(double[][] a, Factors f) {
        double error = 0.0, scale = 0.0;
        int n = a.length;
        for (double[] row : a) for (double v : row) scale = Math.max(scale, Math.abs(v));
        for (int i = 0; i < n; ++i)
            for (int j = 0; j < n; ++j) {
                double value = 0.0;
                for (int k = 0; k < n; ++k) value += f.l[i][k] * f.l[j][k];
                error = Math.max(error, Math.abs(value - a[f.p[i]][f.p[j]]));
            }
        return error / scale;
    }

    static double solveError(double[][] a, double[] x, double[] b) {
        double error = 0.0;
        for (int i = 0; i < a.length; ++i) {
            double value = 0.0;
            for (int j = 0; j < x.length; ++j) value += a[i][j] * x[j];
            error = Math.max(error, Math.abs(value - b[i]));
        }
        return error;
    }

    static String vectorText(double[] v) {
        StringBuilder s = new StringBuilder("[");
        for (int i = 0; i < v.length; ++i) {
            if (i > 0) s.append(", ");
            s.append(String.format(Locale.ROOT, "%.6f", v[i] == 0.0 ? 0.0 : v[i]));
        }
        return s.append(']').toString();
    }

    static String matrixText(double[][] a) {
        StringBuilder s = new StringBuilder("[");
        for (int i = 0; i < a.length; ++i) {
            if (i > 0) s.append(", ");
            s.append(vectorText(a[i]));
        }
        return s.append(']').toString();
    }

    public static void main(String[] args) {
        double[][] a = { {10.0, 5.0, 4.0}, {5.0, 9.0, 8.0}, {4.0, 8.0, 16.0} };
        double[] b = {32.0, 47.0, 68.0}, b2 = {10.0, 5.0, 4.0};
        Factors f = choleskyFactor(a, 1e-12);
        double[] x = choleskySolve(f, b), x2 = choleskySolve(f, b2);
        double fe = factorError(a, f), se = solveError(a, x, b), re = solveError(a, x2, b2);
        System.out.println("permutation = " + Arrays.toString(f.p));
        System.out.println("L = " + matrixText(f.l));
        System.out.println("x = " + vectorText(x));
        System.out.println("reused x = " + vectorText(x2));
        System.out.printf(Locale.ROOT, "relative factor error = %.3e%nsolve residual = %.3e%nreuse residual = %.3e%n", fe, se, re);
        System.out.println("check = " + (Math.max(fe, Math.max(se, re)) < 1e-10 ? "PASS" : "FAIL"));
    }
}

All three programs print these eight lines:

• permutation = [2, 0, 1]
• L = [[4.000000, 0.000000, 0.000000], [1.000000, 3.000000, 0.000000], [2.000000, 1.000000, 2.000000]]
• x = [1.000000, 2.000000, 3.000000]
• reused x = [1.000000, 0.000000, 0.000000]
• relative factor error = 0.000e+00
• solve residual = 0.000e+00
• reuse residual = 0.000e+00
• check = PASS

Complexity, Correctness, and Invariants

For a successful dense $n\times n$ input, these implementations take $\Theta(n^3)$ factorization time in the best, average, and worst cases. They execute the dot-product loops even when entries are zero. Pivot searches, symmetric swaps, prefix swaps, and cached diagonal updates add $O(n^2)$ work overall. The leading factorization arithmetic is about $n^3/3$ operations when multiplication and subtraction are counted separately, compared with about $2n^3/3$ for general dense LU.

The code uses $\Theta(n^2)$ storage for the input copy and explicit square factor arrays, plus $O(n)$ permutation and diagonal arrays. Although only $n(n+1)/2$ factor entries are mathematically needed, the teaching code allocates full nested arrays. Each solve costs $\Theta(n^2)$ time and $\Theta(n)$ extra workspace. Factoring once and solving $q$ right-hand sides costs $\Theta(n^3+qn^2)$ total time. Rejection may stop early, but valid square input validation and copying already require quadratic work.

To justify correctness in exact arithmetic, consider the currently permuted, normalized matrix $B=P(A/s)P^T$. After $k$ columns, let $C_k$ contain the completed columns of the normalized factor. The residual relation has the block form

$$ B=C_kC_k^T+\begin{pmatrix}0&0\\0&S_k\end{pmatrix}, $$

where $S_k$ is the remaining Schur complement. The cached diagonal satisfies $d_i=(S_k)_{ii}$ for active indices. Initially there are no completed columns and $S_0=A/s$.

A symmetric swap of active variables permutes both rows and columns of the residual block and the corresponding rows of $C_k$. Swapping the completed factor prefixes therefore preserves the relation. A positive pivot $d$ creates the new column with first active entry $\sqrt d$ and lower entries equal to the residual column divided by $\sqrt d$. Subtracting that column's outer product clears its row and column in the residual, leaving the next Schur complement. The dot products and diagonal updates in the code implement exactly this relation without explicitly storing the entire Schur complement.

For SPD input, every such Schur complement is SPD in exact arithmetic, so its diagonal pivots are positive and the process can continue. The code imposes the additional tolerance requirement. At termination, no residual block remains: $P(A/s)P^T=CC^T$, hence $PAP^T=LL^T$ after rescaling.

The solves then give $LL^Tz=Pb$. Since $z=Px$ and $P$ is invertible, this implies $Ax=b$. Positive factor diagonals make both triangular solves unique. In floating-point arithmetic these identities are approximate; failure on a nearly singular SPD matrix can reflect rounding and the chosen tolerance, while successful reconstruction does not by itself establish a small forward solution error for an ill-conditioned system.

Variants, Pitfalls, and Optimizations

Positive diagonal is not enough. The matrix $\begin{pmatrix}1&2\\2&1\end{pmatrix}$ has a first Cholesky entry $1$ and a lower entry $2$, but its next residual pivot is $1-2^2=-3$. Stop before taking a real square root. Replacing the pivot by its absolute value would create $L_{\mathrm{fake}}=\begin{pmatrix}1&0\\2&\sqrt3\end{pmatrix}$, whose product has bottom-right entry $7$ instead of $1$: it factors a different matrix. A zero pivot likewise needs a separate semidefinite or rank-deficient policy.

Preserve symmetry and order. Swapping only rows destroys the symmetric problem. Swap rows and columns together, move earlier factor entries, gather $b_{p_i}$ for the forward solve, and scatter $x_{p_i}=z_i$ afterward. In the transpose solve, use $L_{ji}$ for an upper-triangular entry. A factor may look triangular while a missing prefix swap makes its product wrong.

Pivoted, unpivoted, and indefinite variants. Standard SPD Cholesky usually runs without numerical pivoting, possibly after a structural fill-reducing permutation. The LAPACK DPSTRF reference describes a complete-pivoting variant for semidefinite problems; its convention is $P_{\mathrm{lib}}^TAP_{\mathrm{lib}}=LL^T$, corresponding to our $P=P_{\mathrm{lib}}^T$. Its rank-return and stopping interface differ from our full-factor-or-error routine. For symmetric indefinite input, a suitably pivoted $LDL^T$ factorization is the appropriate different tool, rather than forcing Cholesky to continue.

Blocked, sparse, and approximate factors. Blocked Cholesky groups trailing updates into matrix operations for better cache use and optimized BLAS execution. Sparse Cholesky adds symbolic analysis and ordering to reduce fill-in; the largest-diagonal ordering used here is not a general fill-reduction strategy. Incomplete Cholesky drops selected contributions to build a preconditioner, with separate breakdown and stability considerations. Banded and tridiagonal SPD matrices admit specialized routines with much lower costs than dense storage and loops.

Updates and modeling changes. A rank-one update replaces $A$ by $A+uu^T$ and can update a dense factor in $O(n^2)$ work. A downdate $A-uu^T$ needs a positive-definiteness check and can fail. Adding a diagonal term $\lambda I$ can be an intentional regularization choice, but it changes the system; it must not be an unexplained repair for a failed pivot. Approximately symmetric data likewise require a deliberate policy.

Other uses of the factor. Compute $\log\det(A)=2\sum_i\log L_{ii}$ to avoid multiplying many diagonal values. For covariance $\Sigma$ with $P\Sigma P^T=LL^T$, an independent standard normal vector $g$ gives a correlated draw $P^TLg$. Full-rank covariance is needed for this invertible factor; degenerate covariance can require a semidefinite method. For complex Hermitian positive definite matrices, replace the transpose by the conjugate transpose. Avoid forming $A^{-1}$ merely to solve a system, and remember that forming normal equations can square the condition number of an underlying least-squares matrix.

Practice Problems and Real-World Use

Begin with a known triangular matrix $C$ whose diagonal is positive, form $A=CC^T$, and permute its variables. This guarantees an SPD test matrix in exact arithmetic and lets you check the factor product independently of the solver. To test solves, choose a known $x$, generate $b=Ax$, factor once, and recover several vectors from the same factors. Different pivot orders can yield different valid factors, so always compare with the appropriate symmetrically permuted matrix.

Include one-by-one inputs, a second-stage swap, zero and negative pivots, exact asymmetry, very small accepted pivots under a changed tolerance, and uniformly scaled systems. Verify lower-triangular shape, positive diagonals, input preservation, $LL^T\approx PAP^T$, and the original-system residual. Use a scale-aware residual when the matrix and vectors vary greatly in magnitude.