LU and Cholesky Factorisation

Setup

The central operation: solving $Ax = b$

Almost every numerical routine in quantitative finance reduces to solving a linear system. Yield curve bootstrapping solves $Ax = b$ for discount factors. Greeks computed by bump-and-reval solve for a perturbed price. Calibration via Gauss-Newton solves the normal equations $J^\top J \delta = -J^\top r$. Even generating correlated random variables for Monte Carlo requires factorising a covariance matrix.

The naïve approach — computing $A^{-1}$ explicitly and then forming $A^{-1}b$ — is both slower and less numerically stable than direct factorisation. The standard alternatives are:

• LU factorisation with partial pivoting: general square $A$, $O(n^3/3)$ flops.
• Cholesky factorisation: symmetric positive definite $A$, $O(n^3/6)$ flops — twice as fast as LU, always stable without pivoting.

Assumptions

• $A \in \mathbb{R}^{n \times n}$ is square. For LU: assume $A$ is non-singular (otherwise partial pivoting detects the singularity). For Cholesky: assume $A$ is symmetric positive definite (SPD).
• All arithmetic is in exact real arithmetic unless otherwise stated; floating-point effects are discussed in Limitations.
• Conventions: LU factorisation produces $PA = LU$ where $P$ is a permutation matrix, $L$ is unit lower triangular (diagonal entries all 1), and $U$ is upper triangular. Cholesky produces $A = LL^\top$ where $L$ is lower triangular with positive diagonal.

Theory

1. LU Factorisation

The idea is systematic Gaussian elimination: reduce $A$ to upper triangular form $U$ by applying elementary row operations, then record those operations as a lower triangular matrix $L$.

Algorithm (Gaussian elimination with partial pivoting):

At step $k$ (for $k = 1, \ldots, n-1$):

1. Find the row $p \geq k$ with the largest absolute value in column $k$: $p = \arg\max_{j \geq k} |A_{jk}|$ (partial pivot).
2. Swap rows $k$ and $p$ in $A$.
3. For each row $i > k$: compute the multiplier $\ell_{ik} = A_{ik} / A_{kk}$ and subtract $\ell_{ik}$ times row $k$ from row $i$.
4. Store $\ell_{ik}$ in the lower triangular part of $L$.

After $n-1$ steps, $A$ has been transformed to $U$ and the multipliers fill $L$.

Solving $Ax = b$ using LU:

Given $PA = LU$, the system $Ax = b$ becomes $LUx = Pb$. Solve in two triangular passes:

1. Forward substitution: solve $Ly = Pb$ for $y$ — $O(n^2)$ operations.
2. Backward substitution: solve $Ux = y$ for $x$ — $O(n^2)$ operations.

Total cost: $O(n^3/3)$ for the factorisation (done once), then $O(n^2)$ per right-hand side. For bootstrapping with multiple right-hand sides (e.g., different rate scenarios), reusing the factorisation is a significant efficiency gain.

2. Cholesky Factorisation

When $A$ is symmetric and positive definite, a more efficient and always-stable factorisation exists.