Appendix O — Cholesky Decomposition

O.1 Cholesky Decomposition

The Cholesky decomposition is a method for factorizing a positive definite matrix into the product of a lower triangular matrix and its transpose. It is particularly useful for solving systems of linear equations and for generating samples from multivariate normal distributions.

Definition O.1 (Definition of the Cholesky Decomposition) André-Louis Cholesky (1875–1918) was a cartographer in the French army, who introduced a method for decomposing a symmetric positive definite matrix into the product of a lower triangular matrix and its transpose. This decomposition is known as the Cholesky decomposition.

Let A be a symmetric positive definite matrix. The Cholesky decomposition of A is a factorization of the form: A = LL^\top \tag{O.1} where L is a lower triangular matrix with positive diagonal entries.

The Cholesky decomposition is unique for a given positive definite matrix, and it can be computed efficiently using algorithms such as the Doolittle algorithm or the Crout algorithm.

O.1.1 Doolittle Algorithm

Algorithm O.1

\begin{algorithm} \caption{Doolittle's algorithm for LU decomposition} \begin{algorithmic} \Require $A$ is an $n \times n$ matrix \Ensure $L$ is a unit lower triangular matrix and $U$ is an upper triangular matrix such that $A = LU$ \State Initialize $L$ as an $n \times n$ identity matrix and $U$ as an $n \times n$ zero matrix \For{$i = 0$ to $n-1$} \For{$k = i$ to $n-1$} \State $U[i][k] = A[i][k] - \sum_{j=0}^{i-1} L[i][j] U[j][k]$ \EndFor \For{$k = i+1$ to $n-1$} \State $L[k][i] = \frac{A[k][i] - \sum_{j=0}^{i-1} L[k][j] U[j][i]}{U[i][i]}$ \EndFor \EndFor \end{algorithmic} \end{algorithm}

def doolittle(A):
    """Performs Doolittle LU decomposition: A = LU, with L unit lower triangular and U upper triangular."""
    n = len(A)
    L = [[0.0]*n for _ in range(n)]
    U = [[0.0]*n for _ in range(n)]
    
    for i in range(n):
        # Upper Triangular
        for k in range(i, n):
            U[i][k] = A[i][k] - sum(L[i][j]*U[j][k] for j in range(i))
        
        # Lower Triangular
        L[i][i] = 1.0
        for k in range(i+1, n):
            if U[i][i] == 0:
                raise ZeroDivisionError("Zero pivot encountered.")
            L[k][i] = (A[k][i] - sum(L[k][j]*U[j][i] for j in range(i))) / U[i][i]
    
    return L, U


def print_matrix(M, name="Matrix"):
    print(f"{name} =")
    for row in M:
        print("  [" + "  ".join(f"{val:8.3f}" for val in row) + "]")
    print()

A = [
    [2, 3, 1],
    [4, 7, 7],
    [6, 18, 22]
]
print_matrix(A, name="A")

L, U = doolittle(A)
print_matrix(L, name="L")
print_matrix(U, name="U")

A =
  [   2.000     3.000     1.000]
  [   4.000     7.000     7.000]
  [   6.000    18.000    22.000]

L =
  [   1.000     0.000     0.000]
  [   2.000     1.000     0.000]
  [   3.000     9.000     1.000]

U =
  [   2.000     3.000     1.000]
  [   0.000     1.000     5.000]
  [   0.000     0.000   -26.000]

O.1.2 Doolittle’s Algorithm with Partial Pivoting

When performing LU decomposition, it is often necessary to use partial pivoting to ensure numerical stability and to handle cases where the matrix may be singular or nearly singular. Partial pivoting involves swapping rows of the matrix to place the largest absolute value in the pivot position.

Adding partial pivoting is algebraically equivalent to multiplying the original matrix by a permutation matrix P, such that PA = LU, where P is a permutation matrix, L is a lower triangular matrix, and U is an upper triangular matrix.

def doolittle_partial_pivoting(A):
    """Performs LU decomposition with partial pivoting: PA = LU."""
    n = len(A)
    # Deep copy of A
    A = [row[:] for row in A]
    P = list(range(n))
    L = [[0.0]*n for _ in range(n)]
    U = [[0.0]*n for _ in range(n)]

    for k in range(n):
        # Partial pivoting: find row with max abs value in column k
        pivot_row = max(range(k, n), key=lambda i: abs(A[i][k]))
        if A[pivot_row][k] == 0:
            raise ZeroDivisionError("Matrix is singular.")

        # Swap rows in A and record permutation
        A[k], A[pivot_row] = A[pivot_row], A[k]
        P[k], P[pivot_row] = P[pivot_row], P[k]
        for i in range(k):
            L[k][i], L[pivot_row][i] = L[pivot_row][i], L[k][i]

        # Compute U[k][k:] and L[k+1:][k]
        L[k][k] = 1.0
        for j in range(k, n):
            U[k][j] = A[k][j] - sum(L[k][s]*U[s][j] for s in range(k))
        for i in range(k+1, n):
            L[i][k] = (A[i][k] - sum(L[i][s]*U[s][k] for s in range(k))) / U[k][k]

    # Permutation matrix P as a 2D matrix
    P_matrix = [[1 if j == P[i] else 0 for j in range(n)] for i in range(n)]
    return P_matrix, L, U

Demo for Doolittle’s algorithm with partial pivoting:

A = [
    [0, 3, 1],
    [4, 7, 7],
    [6, 18, 22]
]

print_matrix(A, name="A")


P, L, U = doolittle_partial_pivoting(A)
print_matrix(P, name="P")
print_matrix(L, name="L")
print_matrix(U, name="U")

A =
  [   0.000     3.000     1.000]
  [   4.000     7.000     7.000]
  [   6.000    18.000    22.000]

P =
  [   0.000     0.000     1.000]
  [   0.000     1.000     0.000]
  [   1.000     0.000     0.000]

L =
  [   1.000     0.000     0.000]
  [   0.667     1.000     0.000]
  [   0.000    -0.600     1.000]

U =
  [   6.000    18.000    22.000]
  [   0.000    -5.000    -7.667]
  [   0.000     0.000    -3.600]

A = [
    [2, 1, 1, 3, 2],
    [1, 2, 2, 1, 1],
    [3, 2, 3, 2, 1],
    [2, 1, 2, 2, 1],
    [1, 1, 1, 1, 1]
]

P, L, U = doolittle_partial_pivoting(A)

print_matrix(P, "P")
print_matrix(L, "L")
print_matrix(U, "U")

def matmul(A, B):
    return [[sum(A[i][k] * B[k][j] for k in range(len(B)))
             for j in range(len(B[0]))]
            for i in range(len(A))]

def permute(A, P):
    """P is a permutation matrix; return PA."""
    return matmul(P, A)

PA = permute(A, P)
LU = matmul(L, U)

# Print comparison
print_matrix(PA, "PA")
print_matrix(LU, "LU")

P =
  [   0.000     0.000     1.000     0.000     0.000]
  [   0.000     1.000     0.000     0.000     0.000]
  [   0.000     0.000     0.000     1.000     0.000]
  [   1.000     0.000     0.000     0.000     0.000]
  [   0.000     0.000     0.000     0.000     1.000]

L =
  [   1.000     0.000     0.000     0.000     0.000]
  [   0.333     1.000     0.000     0.000     0.000]
  [   0.667    -0.250     1.000     0.000     0.000]
  [   0.667    -0.250    -3.000     1.000     0.000]
  [   0.333     0.250    -1.000     0.250     1.000]

U =
  [   3.000     2.000     3.000     2.000     1.000]
  [   0.000     1.333     1.000     0.333     0.667]
  [   0.000     0.000     0.250     0.750     0.500]
  [   0.000     0.000     0.000     4.000     3.000]
  [   0.000     0.000     0.000     0.000     0.250]

PA =
  [   3.000     2.000     3.000     2.000     1.000]
  [   1.000     2.000     2.000     1.000     1.000]
  [   2.000     1.000     2.000     2.000     1.000]
  [   2.000     1.000     1.000     3.000     2.000]
  [   1.000     1.000     1.000     1.000     1.000]

LU =
  [   3.000     2.000     3.000     2.000     1.000]
  [   1.000     2.000     2.000     1.000     1.000]
  [   2.000     1.000     2.000     2.000     1.000]
  [   2.000     1.000     1.000     3.000     2.000]
  [   1.000     1.000     1.000     1.000     1.000]

O.1.3 Vectorization of the doolittle algorithm

import numpy as np

def doolittle_numpy(A):
    """LU decomposition with partial pivoting using NumPy. Returns P, L, U such that PA = LU."""
    A = np.array(A, dtype=float)
    n = A.shape[0]
    P = np.eye(n)
    L = np.zeros((n, n))
    U = A.copy()

    for k in range(n):
        # Partial pivoting
        pivot = np.argmax(abs(U[k:, k])) + k
        if U[pivot, k] == 0:
            raise ZeroDivisionError("Matrix is singular.")
        if pivot != k:
            U[[k, pivot]] = U[[pivot, k]]
            P[[k, pivot]] = P[[pivot, k]]
            L[[k, pivot], :k] = L[[pivot, k], :k]

        L[k, k] = 1.0
        L[k+1:, k] = U[k+1:, k] / U[k, k]
        U[k+1:] -= np.outer(L[k+1:, k], U[k])

    return P, L, U

def random_sign_matrix(n, seed=None):
    """Generate an n×n matrix with random entries in {-1, 0, 1}."""
    rng = np.random.default_rng(seed)
    return rng.choice([-1, 0, 1], size=(n, n))

A = [
    [2, 3, 1],
    [4, 7, 7],
    [6, 18, 22]
]
A = random_sign_matrix(16, seed=42)

print_matrix(A, name="A")

P, L, U = doolittle_numpy(A)

print_matrix(P, name="P")

print_matrix(L, name="L")

print_matrix(U, name="U")

A =
  [  -1.000     1.000     0.000     0.000     0.000     1.000    -1.000     1.000    -1.000    -1.000     0.000     1.000     1.000     1.000     1.000     1.000]
  [   0.000    -1.000     1.000     0.000     0.000     0.000    -1.000     1.000     1.000     0.000     0.000     1.000     0.000     0.000     0.000    -1.000]
  [  -1.000     0.000     1.000    -1.000     1.000     1.000    -1.000     0.000    -1.000     1.000     1.000     0.000    -1.000     1.000     0.000     1.000]
  [   1.000     1.000     1.000    -1.000     0.000     0.000     0.000    -1.000     0.000    -1.000     1.000     1.000     1.000     1.000     0.000     1.000]
  [   0.000    -1.000     1.000     0.000    -1.000     0.000     1.000    -1.000     0.000    -1.000     1.000     0.000    -1.000    -1.000     0.000     1.000]
  [   1.000     0.000    -1.000     1.000     0.000     1.000    -1.000    -1.000     1.000     1.000     0.000     1.000     1.000     0.000     1.000    -1.000]
  [  -1.000     1.000     0.000    -1.000     1.000    -1.000     1.000    -1.000     1.000     1.000     1.000     0.000     0.000     1.000    -1.000     1.000]
  [   0.000     0.000     0.000     0.000    -1.000    -1.000    -1.000    -1.000     0.000     1.000     0.000     0.000     1.000     0.000    -1.000     1.000]
  [   0.000     0.000     0.000     0.000    -1.000     0.000     1.000    -1.000     0.000    -1.000     0.000     0.000     1.000    -1.000    -1.000     0.000]
  [   1.000     1.000    -1.000    -1.000     1.000    -1.000     1.000    -1.000     1.000    -1.000     0.000     0.000    -1.000     0.000     0.000     1.000]
  [   1.000     0.000     0.000     0.000     0.000     1.000    -1.000    -1.000     0.000    -1.000    -1.000    -1.000     1.000     1.000     1.000     0.000]
  [   1.000    -1.000     1.000     0.000     1.000    -1.000     0.000     1.000     0.000     0.000    -1.000     0.000    -1.000    -1.000     1.000     0.000]
  [   0.000     0.000     1.000    -1.000     0.000    -1.000     0.000     1.000     0.000     1.000     0.000     1.000     0.000    -1.000     1.000     1.000]
  [  -1.000     1.000    -1.000     1.000     1.000     0.000     1.000    -1.000    -1.000    -1.000     0.000     1.000    -1.000     0.000     1.000    -1.000]
  [   1.000    -1.000     1.000     0.000     0.000    -1.000     0.000     1.000    -1.000     1.000     0.000     1.000     0.000     0.000    -1.000     0.000]
  [  -1.000     0.000    -1.000     0.000     0.000    -1.000     1.000     0.000     1.000    -1.000    -1.000     0.000    -1.000    -1.000     0.000    -1.000]

P =
  [   1.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000]
  [   0.000     0.000     0.000     1.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000]
  [   0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     1.000     0.000     0.000     0.000     0.000     0.000     0.000]
  [   0.000     0.000     0.000     0.000     0.000     1.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000]
  [   0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     1.000     0.000     0.000     0.000     0.000]
  [   0.000     0.000     1.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000]
  [   0.000     0.000     0.000     0.000     0.000     0.000     0.000     1.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000]
  [   0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     1.000     0.000     0.000]
  [   0.000     0.000     0.000     0.000     0.000     0.000     1.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000]
  [   0.000     1.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000]
  [   0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     1.000     0.000     0.000     0.000     0.000     0.000]
  [   0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     1.000     0.000     0.000     0.000]
  [   0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     1.000     0.000]
  [   0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     1.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000]
  [   0.000     0.000     0.000     0.000     1.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000]
  [   0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     1.000]

L =
  [   1.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000]
  [  -1.000     1.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000]
  [  -1.000     1.000     1.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000]
  [  -1.000     0.500     0.750     1.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000]
  [  -1.000     0.000    -0.500     0.000     1.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000]
  [   1.000    -0.500    -0.750    -1.000     0.667     1.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000]
  [  -0.000     0.000    -0.000     0.000    -0.667    -0.571     1.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000]
  [   1.000     0.000     0.500     0.667     0.667    -0.714    -0.938     1.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000]
  [   1.000     0.000    -0.000    -0.667     0.333    -0.143    -0.187     0.537     1.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000]
  [  -0.000    -0.500    -0.750    -0.333     0.333     0.286     0.375    -0.232     0.670     1.000     0.000     0.000     0.000     0.000     0.000     0.000]
  [  -1.000     0.500     0.250     0.333     0.000     0.429     0.125    -0.189    -0.393     0.619     1.000     0.000     0.000     0.000     0.000     0.000]
  [  -0.000     0.000    -0.500    -0.667     0.000     0.000     0.438    -0.242     0.214    -0.636    -0.177     1.000     0.000     0.000     0.000     0.000]
  [  -1.000     0.000    -0.500     0.000     0.333    -0.143     0.250    -0.211    -0.625    -0.996    -0.786     0.408     1.000     0.000     0.000     0.000]
  [  -0.000     0.000    -0.000     0.000    -0.667    -0.143    -0.188     0.032    -0.161     0.360     0.480     0.771     0.733     1.000     0.000     0.000]
  [  -0.000    -0.500    -0.750    -0.333    -0.333     0.143    -0.250     0.042    -0.045     0.933     0.070     0.108     0.143     0.431     1.000     0.000]
  [   1.000    -0.500     0.250    -0.333    -0.333    -0.286     0.062     0.158     0.554     0.895     0.568     0.802     0.025     0.140    -0.221     1.000]

U =
  [  -1.000     1.000     0.000     0.000     0.000     1.000    -1.000     1.000    -1.000    -1.000     0.000     1.000     1.000     1.000     1.000     1.000]
  [   0.000     2.000     1.000    -1.000     0.000     1.000    -1.000     0.000    -1.000    -2.000     1.000     2.000     2.000     2.000     1.000     2.000]
  [   0.000     0.000    -2.000     0.000     1.000    -1.000     1.000     0.000     1.000     0.000    -1.000    -1.000    -2.000    -1.000     0.000     0.000]
  [   0.000     0.000     0.000     1.500    -0.750     2.250    -2.250     0.000    -0.250     1.000     0.250     1.750     2.500     0.750     1.500    -1.000]
  [   0.000     0.000     0.000     0.000     1.500    -0.500    -0.500     2.000    -0.500    -1.000    -1.500     0.500    -1.000    -0.500     2.000     1.000]
  [   0.000     0.000     0.000     0.000     0.000     2.333    -1.667    -2.333     0.333     2.667     2.000     0.667     0.667     1.333    -0.333    -0.667]
  [   0.000     0.000     0.000     0.000     0.000     0.000    -2.286    -1.000    -0.143     1.857     0.143     0.714     0.714     0.429     0.143     1.286]
  [   0.000     0.000     0.000     0.000     0.000     0.000     0.000    -5.938     0.104     3.646     2.896     0.146    -0.854     0.688    -2.437    -1.271]
  [   0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     1.965     1.772     0.425     0.151     1.688     0.568    -0.379    -0.172]
  [   0.000     0.000     0.000     0.000     0.000    -0.000     0.000     0.000     0.000    -2.134     0.095     1.141    -1.120    -0.096     0.064    -1.137]
  [   0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000    -1.551    -2.328     1.486     0.788     0.474     0.854]
  [   0.000     0.000     0.000     0.000     0.000    -0.000     0.000     0.000     0.000     0.000     0.000     1.671    -0.663    -1.065     1.553    -1.072]
  [   0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     1.449     2.208    -1.697    -0.149]
  [   0.000     0.000     0.000     0.000     0.000    -0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000    -2.134     0.123     1.760]
  [   0.000     0.000     0.000     0.000     0.000    -0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     1.765     2.843]
  [   0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.802]

O.1.4 Crout’s Algorithm

Algorithm O.2

\begin{algorithm} \caption{Crout's algorithm for LU decomposition} \begin{algorithmic} \Require $A$ is an $n \times n$ matrix \Ensure $L$ is a lower triangular matrix and $U$ is an upper triangular matrix such that $A = LU$ \State Initialize $L$ as an $n \times n$ zero matrix and $U$ as an $n \times n$ identity matrix \For{$j = 0$ to $n-1$} \For{$i = j$ to $n-1$} \State $L[i][j] = A[i][j] - \sum_{k=0}^{j-1} L[i][k] U[k][j]$ \EndFor \For{$i = j+1$ to $n-1$} \State $U[j][i] = \frac{A[j][i] - \sum_{k=0}^{j-1} L[j][k] U[k][i]}{L[j][j]}$ \EndFor \EndFor \end{algorithmic} \end{algorithm}

Crout’s algorithm is another method for performing LU decomposition, similar to Doolittle’s algorithm. It constructs the lower triangular matrix L directly, while the upper triangular matrix U is obtained from the original matrix A.

the key parts are the formula for computing the entries of L:

L[i][j] = A[i][j] - \sum_{k=0}^{j-1} L[i][k] U[k][j]

and using this to computing the entries of U:

U[j][i] = \frac{A[j][i] - \sum_{k=0}^{j-1} L[j][k] U[k][i]}{L[j][j]}

import numpy as np

def crout_lu(A, *, pivot=True, rtol=1e-9, atol=None):
    """
    Robust Crout LU with partial pivoting.

    Improvements
    ------------
    1. Pivot is selected from the *updated* column, eliminating the
       “false zero-pivot” issue.
    2. Test `abs(pivot) < max(atol, rtol*‖A‖∞)` so tolerance scales with data.

    Returns L, U, P such that  P @ A = L @ U.
    """
    A = np.asarray_chkfinite(A, dtype=float)
    n  = A.shape[0]
    if A.ndim != 2 or n != A.shape[1]:
        raise ValueError("square matrix required")

    L = np.zeros_like(A)
    U = np.eye(n, dtype=A.dtype)
    P = np.eye(n, dtype=A.dtype)
    rows = np.arange(n)

    if atol is None:
        atol = np.finfo(A.dtype).eps * np.linalg.norm(A, np.inf) * 10

    for k in range(n):
        # current residual column
        col = A[rows[k:], k] - L[k:, :k] @ U[:k, k]

        if pivot:
            j = k + np.argmax(np.abs(col))        # best row
            if np.abs(col[j - k]) < max(atol, rtol*np.abs(col).max()):
                raise np.linalg.LinAlgError("matrix is numerically singular")
            if j != k:                            # swap logical rows
                rows[[k, j]] = rows[[j, k]]
                L[[k, j], :k] = L[[j, k], :k]
                P[[k, j]] = P[[j, k]]
                col[[0, j - k]] = col[[j - k, 0]]

        L[k:, k] = col
        if np.abs(L[k, k]) < max(atol, rtol*np.abs(col).max()):
            raise np.linalg.LinAlgError("zero pivot encountered")

        # row k of U (unit diagonal)
        U[k, k+1:] = (A[rows[k], k+1:] - L[k, :k] @ U[:k, k+1:]) / L[k, k]

    return L, U, P

Test the Crout LU decomposition with a random full rank sign matrix:

def random_full_rank_sign_matrix(n, seed=42):
    rng = np.random.default_rng(seed)
    while True:
        A = rng.choice([-1,0, 1], size=(n, n))
        if np.linalg.matrix_rank(A) == n:
            return A

A = random_full_rank_sign_matrix(12, seed=42)
print_matrix(A, name="A")
A_orig = A.copy()  # Keep original for verification
L, U, P = crout_lu(A)

# Check correctness
assert np.allclose(P @ A_orig, L @ U, atol=1e-8)

A =
  [  -1.000     1.000     0.000     0.000     0.000     1.000    -1.000     1.000    -1.000    -1.000     0.000     1.000]
  [   1.000     1.000     1.000     1.000     0.000    -1.000     1.000     0.000     0.000     0.000    -1.000     1.000]
  [   1.000     0.000     0.000     1.000     0.000     0.000     0.000    -1.000    -1.000     0.000     1.000    -1.000]
  [   1.000     1.000    -1.000     0.000    -1.000     1.000     1.000     0.000    -1.000     1.000     0.000     1.000]
  [   1.000     1.000     1.000    -1.000     0.000     0.000     0.000    -1.000     0.000    -1.000     1.000     1.000]
  [   1.000     1.000     0.000     1.000     0.000    -1.000     1.000     0.000    -1.000     0.000     1.000    -1.000]
  [   0.000    -1.000     1.000     0.000    -1.000    -1.000     0.000     1.000     1.000     0.000    -1.000     1.000]
  [   0.000     1.000    -1.000    -1.000     1.000     1.000     0.000     1.000     1.000     0.000     1.000    -1.000]
  [  -1.000     1.000     0.000    -1.000     1.000    -1.000     1.000    -1.000     1.000     1.000     1.000     0.000]
  [   0.000     1.000    -1.000     1.000     0.000     0.000     0.000     0.000    -1.000    -1.000    -1.000    -1.000]
  [   0.000     1.000     0.000     0.000     1.000     0.000    -1.000     1.000     0.000     0.000     0.000     0.000]
  [  -1.000     0.000     1.000    -1.000     0.000    -1.000     0.000     0.000     1.000    -1.000    -1.000     0.000]

The Cholesky decomposition has several important properties:

If A is positive definite, then L is unique and has positive diagonal entries.
The Cholesky decomposition can be used to solve linear systems of equations of the form Ax = b by first solving Ly = b for y, and then solving L^\top x = y.
The Cholesky decomposition can be used to generate samples from a multivariate normal distribution by transforming samples from a standard normal distribution using the Cholesky factor L.

The Cholesky decomposition is widely used in various applications, including numerical optimization, Bayesian inference, and machine learning. It is particularly useful for solving linear systems efficiently and for generating samples from multivariate normal distributions.

--- title: "Cholesky Decomposition" subtitle: "Appendix" description: "This appendix explains the Cholesky decomposition, which is a method for solving linear equations efficiently using linear algebra." categories: - "numerical linear algebra" keywords: - "Moore-Penrose inversion" - "Cholesky decomposition" - "linear equations" - "matrix factorization" - "numerical stability" - "linear algebra" --- ## Cholesky Decomposition {#sec-cholesky} The Cholesky decomposition is a method for factorizing a positive definite matrix into the product of a lower triangular matrix and its transpose. It is particularly useful for solving systems of linear equations and for generating samples from multivariate normal distributions. ::: {#def-cholesky-decomposition} ## Definition of the Cholesky Decomposition André-Louis Cholesky (1875–1918) was a cartographer in the French army, who introduced a method for decomposing a symmetric positive definite matrix into the product of a lower triangular matrix and its transpose. This decomposition is known as the Cholesky decomposition. Let $A$ be a symmetric positive definite matrix. The Cholesky decomposition of $A$ is a factorization of the form: $$ A = LL^\top $$ {#eq-cholesky-decomposition} where $L$ is a lower triangular matrix with positive diagonal entries. ::: The Cholesky decomposition is unique for a given positive definite matrix, and it can be computed efficiently using algorithms such as the [Doolittle algorithm]() or the [Crout algorithm](). ### Doolittle Algorithm {#sec-doolittle-algorithm} ::: {#alg-doolittle-algorithm} ```pseudocode #| html-indent-size: "1.2em" #| html-comment-delimiter: "//" #| html-line-number: true #| html-line-number-punc: ":" #| html-no-end: false #| pdf-placement: "htb!" #| pdf-line-number: true #| pdf-comment-delimiter: "//" \begin{algorithm} \caption{Doolittle's algorithm for LU decomposition} \begin{algorithmic} \Require $A$ is an $n \times n$ matrix \Ensure $L$ is a unit lower triangular matrix and $U$ is an upper triangular matrix such that $A = LU$ \State Initialize $L$ as an $n \times n$ identity matrix and $U$ as an $n \times n$ zero matrix \For{$i = 0$ to $n-1$} \For{$k = i$ to $n-1$} \State $U[i][k] = A[i][k] - \sum_{j=0}^{i-1} L[i][j] U[j][k]$ \EndFor \For{$k = i+1$ to $n-1$} \State $L[k][i] = \frac{A[k][i] - \sum_{j=0}^{i-1} L[k][j] U[j][i]}{U[i][i]}$ \EndFor \EndFor \end{algorithmic} \end{algorithm} ``` ::: ```{python} def doolittle(A): """Performs Doolittle LU decomposition: A = LU, with L unit lower triangular and U upper triangular.""" n = len(A) L = [[0.0]*n for _ in range(n)] U = [[0.0]*n for _ in range(n)] for i in range(n): # Upper Triangular for k in range(i, n): U[i][k] = A[i][k] - sum(L[i][j]*U[j][k] for j in range(i)) # Lower Triangular L[i][i] = 1.0 for k in range(i+1, n): if U[i][i] == 0: raise ZeroDivisionError("Zero pivot encountered.") L[k][i] = (A[k][i] - sum(L[k][j]*U[j][i] for j in range(i))) / U[i][i] return L, U def print_matrix(M, name="Matrix"): print(f"{name} =") for row in M: print(" [" + " ".join(f"{val:8.3f}" for val in row) + "]") print() ``` ```{python} A = [ [2, 3, 1], [4, 7, 7], [6, 18, 22] ] print_matrix(A, name="A") L, U = doolittle(A) print_matrix(L, name="L") print_matrix(U, name="U") ``` ### Doolittle's Algorithm with Partial Pivoting {#sec-doolittle-partial-pivoting} When performing LU decomposition, it is often necessary to use partial pivoting to ensure numerical stability and to handle cases where the matrix may be singular or nearly singular. Partial pivoting involves swapping rows of the matrix to place the largest absolute value in the pivot position. Adding partial pivoting is algebraically equivalent to multiplying the original matrix by a permutation matrix $P$, such that $PA = LU$, where $P$ is a permutation matrix, $L$ is a lower triangular matrix, and $U$ is an upper triangular matrix. ```{python} def doolittle_partial_pivoting(A): """Performs LU decomposition with partial pivoting: PA = LU.""" n = len(A) # Deep copy of A A = [row[:] for row in A] P = list(range(n)) L = [[0.0]*n for _ in range(n)] U = [[0.0]*n for _ in range(n)] for k in range(n): # Partial pivoting: find row with max abs value in column k pivot_row = max(range(k, n), key=lambda i: abs(A[i][k])) if A[pivot_row][k] == 0: raise ZeroDivisionError("Matrix is singular.") # Swap rows in A and record permutation A[k], A[pivot_row] = A[pivot_row], A[k] P[k], P[pivot_row] = P[pivot_row], P[k] for i in range(k): L[k][i], L[pivot_row][i] = L[pivot_row][i], L[k][i] # Compute U[k][k:] and L[k+1:][k] L[k][k] = 1.0 for j in range(k, n): U[k][j] = A[k][j] - sum(L[k][s]*U[s][j] for s in range(k)) for i in range(k+1, n): L[i][k] = (A[i][k] - sum(L[i][s]*U[s][k] for s in range(k))) / U[k][k] # Permutation matrix P as a 2D matrix P_matrix = [[1 if j == P[i] else 0 for j in range(n)] for i in range(n)] return P_matrix, L, U ``` Demo for Doolittle's algorithm with partial pivoting: ````{python} A = [ [0, 3, 1], [4, 7, 7], [6, 18, 22] ] print_matrix(A, name="A") P, L, U = doolittle_partial_pivoting(A) print_matrix(P, name="P") print_matrix(L, name="L") print_matrix(U, name="U") ```` ```{python} A = [ [2, 1, 1, 3, 2], [1, 2, 2, 1, 1], [3, 2, 3, 2, 1], [2, 1, 2, 2, 1], [1, 1, 1, 1, 1] ] P, L, U = doolittle_partial_pivoting(A) print_matrix(P, "P") print_matrix(L, "L") print_matrix(U, "U") def matmul(A, B): return [[sum(A[i][k] * B[k][j] for k in range(len(B))) for j in range(len(B[0]))] for i in range(len(A))] def permute(A, P): """P is a permutation matrix; return PA.""" return matmul(P, A) PA = permute(A, P) LU = matmul(L, U) # Print comparison print_matrix(PA, "PA") print_matrix(LU, "LU") ``` ### Vectorization of the doolittle algorithm {#sec-doolittle-vectorization} ```{python} import numpy as np def doolittle_numpy(A): """LU decomposition with partial pivoting using NumPy. Returns P, L, U such that PA = LU.""" A = np.array(A, dtype=float) n = A.shape[0] P = np.eye(n) L = np.zeros((n, n)) U = A.copy() for k in range(n): # Partial pivoting pivot = np.argmax(abs(U[k:, k])) + k if U[pivot, k] == 0: raise ZeroDivisionError("Matrix is singular.") if pivot != k: U[[k, pivot]] = U[[pivot, k]] P[[k, pivot]] = P[[pivot, k]] L[[k, pivot], :k] = L[[pivot, k], :k] L[k, k] = 1.0 L[k+1:, k] = U[k+1:, k] / U[k, k] U[k+1:] -= np.outer(L[k+1:, k], U[k]) return P, L, U ``` ```{python} def random_sign_matrix(n, seed=None): """Generate an n×n matrix with random entries in {-1, 0, 1}.""" rng = np.random.default_rng(seed) return rng.choice([-1, 0, 1], size=(n, n)) A = [ [2, 3, 1], [4, 7, 7], [6, 18, 22] ] A = random_sign_matrix(16, seed=42) print_matrix(A, name="A") P, L, U = doolittle_numpy(A) print_matrix(P, name="P") print_matrix(L, name="L") print_matrix(U, name="U") ``` ### Crout's Algorithm {#sec-crout-algorithm} ::: {#alg-crout-algorithm} ```pseudocode #| html-indent-size: "1.2em" #| html-comment-delimiter: "//" #| html-line-number: true #| html-line-number-punc: ":" #| html-no-end: false #| pdf-placement: "htb!" #| pdf-line-number: true #| pdf-comment-delimiter: "//" \begin{algorithm} \caption{Crout's algorithm for LU decomposition} \begin{algorithmic} \Require $A$ is an $n \times n$ matrix \Ensure $L$ is a lower triangular matrix and $U$ is an upper triangular matrix such that $A = LU$ \State Initialize $L$ as an $n \times n$ zero matrix and $U$ as an $n \times n$ identity matrix \For{$j = 0$ to $n-1$} \For{$i = j$ to $n-1$} \State $L[i][j] = A[i][j] - \sum_{k=0}^{j-1} L[i][k] U[k][j]$ \EndFor \For{$i = j+1$ to $n-1$} \State $U[j][i] = \frac{A[j][i] - \sum_{k=0}^{j-1} L[j][k] U[k][i]}{L[j][j]}$ \EndFor \EndFor \end{algorithmic} \end{algorithm} ``` ::: Crout's algorithm is another method for performing LU decomposition, similar to Doolittle's algorithm. It constructs the lower triangular matrix $L$ directly, while the upper triangular matrix $U$ is obtained from the original matrix $A$. the key parts are the formula for computing the entries of $L$: $$ L[i][j] = A[i][j] - \sum_{k=0}^{j-1} L[i][k] U[k][j] $$ and using this to computing the entries of $U$: $$ U[j][i] = \frac{A[j][i] - \sum_{k=0}^{j-1} L[j][k] U[k][i]}{L[j][j]} $$ ```{python} import numpy as np def crout_lu(A, *, pivot=True, rtol=1e-9, atol=None): """ Robust Crout LU with partial pivoting. Improvements ------------ 1. Pivot is selected from the *updated* column, eliminating the “false zero-pivot” issue. 2. Test `abs(pivot) < max(atol, rtol*‖A‖∞)` so tolerance scales with data. Returns L, U, P such that P @ A = L @ U. """ A = np.asarray_chkfinite(A, dtype=float) n = A.shape[0] if A.ndim != 2 or n != A.shape[1]: raise ValueError("square matrix required") L = np.zeros_like(A) U = np.eye(n, dtype=A.dtype) P = np.eye(n, dtype=A.dtype) rows = np.arange(n) if atol is None: atol = np.finfo(A.dtype).eps * np.linalg.norm(A, np.inf) * 10 for k in range(n): # current residual column col = A[rows[k:], k] - L[k:, :k] @ U[:k, k] if pivot: j = k + np.argmax(np.abs(col)) # best row if np.abs(col[j - k]) < max(atol, rtol*np.abs(col).max()): raise np.linalg.LinAlgError("matrix is numerically singular") if j != k: # swap logical rows rows[[k, j]] = rows[[j, k]] L[[k, j], :k] = L[[j, k], :k] P[[k, j]] = P[[j, k]] col[[0, j - k]] = col[[j - k, 0]] L[k:, k] = col if np.abs(L[k, k]) < max(atol, rtol*np.abs(col).max()): raise np.linalg.LinAlgError("zero pivot encountered") # row k of U (unit diagonal) U[k, k+1:] = (A[rows[k], k+1:] - L[k, :k] @ U[:k, k+1:]) / L[k, k] return L, U, P ``` Test the Crout LU decomposition with a random full rank sign matrix: ```{python} def random_full_rank_sign_matrix(n, seed=42): rng = np.random.default_rng(seed) while True: A = rng.choice([-1,0, 1], size=(n, n)) if np.linalg.matrix_rank(A) == n: return A A = random_full_rank_sign_matrix(12, seed=42) print_matrix(A, name="A") A_orig = A.copy() # Keep original for verification L, U, P = crout_lu(A) # Check correctness assert np.allclose(P @ A_orig, L @ U, atol=1e-8) ``` The Cholesky decomposition has several important properties: 1. If $A$ is positive definite, then $L$ is unique and has positive diagonal entries. 2. The Cholesky decomposition can be used to solve linear systems of equations of the form $Ax = b$ by first solving $Ly = b$ for $y$, and then solving $L^\top x = y$. 3. The Cholesky decomposition can be used to generate samples from a multivariate normal distribution by transforming samples from a standard normal distribution using the Cholesky factor $L$. The Cholesky decomposition is widely used in various applications, including numerical optimization, Bayesian inference, and machine learning. It is particularly useful for solving linear systems efficiently and for generating samples from multivariate normal distributions.