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.

Algorithm (Cholesky-Banachiewicz):

For $j = 1, \ldots, n$: $L_{jj} = \sqrt{A_{jj} - \sum_{k=1}^{j-1} L_{jk}^2}$ $L_{ij} = \frac{1}{L_{jj}} \left( A_{ij} - \sum_{k=1}^{j-1} L_{ik} L_{jk} \right) \quad \text{for } i = j+1, \ldots, n.$

The square root $\sqrt{A_{jj} - \sum_{k=1}^{j-1} L_{jk}^2}$ is real and positive iff $A$ is SPD — this is both the algorithm and the constructive proof of existence.

Proof sketch. ($\Rightarrow$) If $A = LL^\top$ and $L$ has positive diagonal, then for any $x \neq 0$: $x^\top A x = x^\top L L^\top x = \|L^\top x\|^2 > 0$ (since $L$ invertible $\Rightarrow L^\top x \neq 0$). ($\Leftarrow$) If $A \succ 0$, the algorithm runs to completion without encountering a square root of a non-positive number — this follows from the Schur complement characterisation of positive definiteness. $\square$

3. Application: Correlated Monte Carlo Sampling

Given a covariance matrix $\Sigma \in \mathbb{R}^{n \times n}$ (SPD), to generate $m$ samples from $\mathcal{N}(0, \Sigma)$:

1. Compute $\Sigma = LL^\top$ via Cholesky.
2. Draw $Z \in \mathbb{R}^{n \times m}$ with i.i.d. $\mathcal{N}(0,1)$ entries.
3. Set $X = LZ$.

Then $\mathbb{E}[XX^\top / m] \approx L \cdot I \cdot L^\top = \Sigma$ for large $m$. Each column of $X$ is a sample from $\mathcal{N}(0, \Sigma)$.

4. Solving $Ax = b$ Using Cholesky

For SPD systems, Cholesky gives the same two-pass strategy as LU:

Given $A = LL^\top$:

1. Forward substitution: solve $Ly = b$ — $O(n^2)$.
2. Backward substitution: solve $L^\top x = y$ — $O(n^2)$.

This is the fastest general method for SPD systems, used in all production-grade libraries (LAPACK dpotrs, Eigen's LLT solver, Armadillo's solve).

Validation

The companion notebook verifies:

1. LU correctness — factorises a $4 \times 4$ matrix, checks $PA = LU$ to machine precision, solves $Ax = b$ and verifies the residual.
2. Cholesky correctness — factorises a $3 \times 3$ SPD matrix by the explicit formula, verifies $LL^\top = A$.
3. Correlated sampling — generates 100,000 bivariate normal samples from Example 3.1, verifies the empirical covariance converges to $\Sigma$.
4. Failure detection — confirms that a near-singular (rank-deficient) matrix causes Cholesky to fail, and shows the smallest eigenvalue is the signal.

(Answer: $L_{11} = 2$, $L_{21} = 1$, $L_{22} = \sqrt{3 - 1} = \sqrt{2}$. So $L = \begin{pmatrix} 2 & 0 \\ 1 & \sqrt{2} \end{pmatrix}$ and $x = (2, 1 - \sqrt{2})^\top$.)

Limitations

Scope limitations:

• Cholesky assumes exact SPD. For SPSD matrices (rank-deficient covariance), use the modified Cholesky or the pivoted Cholesky factorisation (which stops when the diagonal would become zero, producing a low-rank factor).
• These direct methods are $O(n^3)$ and impractical for $n > 10^4$. For large sparse systems (e.g., PDE grids), iterative methods are required (Module 5).

Interview Angle