Cholesky Decomposition

You want to simulate tomorrow's returns for ten correlated stocks. Your computer's random number generator only produces independent numbers — every draw ignorant of every other. How do you turn a stream of independent noise into ten numbers with exactly the correlation structure you specified in a covariance matrix? You need something that plays the role of a "square root" for a matrix, and the Cholesky decomposition is that square root, built specifically for the matrices that show up as covariances.

The analogy: mixing paint to hit a target colour

Suppose you have pure red, pure green, and pure blue paint (independent ingredients) and you want to mix a specific shade of brown (a target with known relationships between its components). You need a recipe: how many parts red, then how many parts green on top, then how many parts blue on top of that. Cholesky is exactly this recipe, written as a matrix. It tells you: take independent standard-normal number one, scale it; take independent number two, mix in a bit of number one, scale the result; take independent number three, mix in bits of both — and so on, until the resulting numbers have exactly the variances and covariances you asked for. Each new ingredient only ever adds to the previous ones, never adjusts something already mixed — that one-directional structure is what makes Cholesky fast and unique.

Writing it down

Take a covariance matrix $\Sigma$ — symmetric, and positive definite (meaning every portfolio built from these assets has strictly positive variance; no combination is riskless or negatively risky). The Cholesky decomposition finds a lower-triangular matrix $L$ (zeros above the diagonal) such that

$$\Sigma = LL^\top .$$

In words: $L$ is a "square root" of $\Sigma$ in matrix form — multiplying it by its own transpose reconstructs the whole covariance matrix. Because $L$ is lower-triangular, its first row only touches the first variable, its second row touches the first two, and so on — that's the "mix in one ingredient at a time" recipe from the paint analogy, made precise.

To generate correlated draws: take a vector $z$ of independent standard-normal numbers (mean 0, variance 1, uncorrelated with each other), and compute

$$x = Lz .$$

In words: multiply your plain independent noise by the recipe matrix $L$, and the output $x$ has exactly covariance $\Sigma$. Why? Because $\text{Cov}(x) = \text{Cov}(Lz) = L\,\text{Cov}(z)\,L^\top = L I L^\top = LL^\top = \Sigma$ — the identity matrix $I$ is the covariance of independent unit-variance noise, and it drops out of the middle.

The entries of $L$ are computed one at a time, top-left to bottom-right, each one from entries already found:

$$L_{ii} = \sqrt{\Sigma_{ii} - \sum_{k<i}L_{ik}^2}, \qquad L_{ji} = \frac{1}{L_{ii}}\left(\Sigma_{ji} - \sum_{k<i}L_{jk}L_{ik}\right) \text{ for } j>i .$$

In words: the diagonal entry is "how much new variance this variable has, after subtracting what's already explained by earlier variables" — under a square root, which is exactly why $\Sigma$ must be positive definite: a negative number under that square root means the matrix wasn't a legitimate covariance matrix at all.

Set the transform to a strong shear (an asymmetric stretch, not a pure rotation) and watch the unit circle become a tilted ellipse. That tilt is correlation: a lower-triangular $L$ applied to independent circular noise produces exactly this kind of tilted, elongated cloud — the shape encodes both the variances (axis lengths) and the correlation (tilt angle).

Worked example 1: two correlated assets, entirely by hand

Two assets with variances $\sigma_1^2 = 0.04$ (20% vol) and $\sigma_2^2=0.09$ (30% vol) and correlation $\rho=0.5$, so $\text{Cov} = \rho\sigma_1\sigma_2 = 0.5(0.2)(0.3)=0.03$. The covariance matrix is $\Sigma = \begin{pmatrix}0.04 & 0.03\\0.03&0.09\end{pmatrix}$.

Step 1, first diagonal: $L_{11}=\sqrt{0.04}=0.20$.

Step 2, below it: $L_{21} = \Sigma_{21}/L_{11} = 0.03/0.20 = 0.15$.

Step 3, second diagonal: $L_{22} = \sqrt{\Sigma_{22} - L_{21}^2} = \sqrt{0.09 - 0.0225} = \sqrt{0.0675} = 0.2598$.

So $L = \begin{pmatrix}0.20 & 0\\0.15 & 0.2598\end{pmatrix}$. Check: $LL^\top = \begin{pmatrix}0.04 & 0.03\\0.03 & 0.0225+0.0675\end{pmatrix} = \begin{pmatrix}0.04&0.03\\0.03&0.09\end{pmatrix}$ — matches $\Sigma$ exactly.

Now generate one correlated draw from independent standard normals $z=(z_1,z_2)=(1.0, -0.5)$ (picked for round arithmetic): $x_1 = L_{11}z_1 = 0.20(1.0)=0.20$; $x_2 = L_{21}z_1 + L_{22}z_2 = 0.15(1.0)+0.2598(-0.5) = 0.15-0.1299=0.0201$. Notice asset 2's draw is built partly from $z_1$ (the shared piece, delivering correlation) and partly from its own independent $z_2$ (the leftover, idiosyncratic piece) — that split is the whole mechanism.

Worked example 2: pricing a spread option via simulation, one path

You want to Monte-Carlo price a spread option paying $\max(S_1-S_2-K,0)$ where $S_1,S_2$ have the $\Sigma$ above, current levels 100 and 100, and $K=5$. Using the correlated draw above as one day's log-return shock scaled by $\sqrt{1/252}$: daily shocks are $x_1/\sqrt{252}=0.0126$ and $x_2/\sqrt{252}=0.00127$, giving $S_1 \approx 100(1.0126) = 101.26$ and $S_2\approx 100(1.00127)=100.13$. The spread is $101.26-100.13=1.13$, below the $5$ strike, so this single path pays $0$. A real pricer repeats this thousands of times, averages the payoffs, and discounts — but every single path starts from exactly this Cholesky step.

What this means in practice

• It's the engine behind every correlated Monte Carlo simulation — risk factor scenarios, VaR simulation, multi-asset option pricing all lean on $\Sigma = LL^\top$.
• It's a fast, reliable way to test if a matrix is a valid covariance matrix. If Cholesky fails (a negative number appears under a square root), $\Sigma$ isn't positive definite — a real problem when covariance matrices are built from noisy, incomplete, or mismatched data.
• It's cheaper than eigen-decomposition. Where the spectral theorem gives a decomposition useful for understanding structure (principal directions), Cholesky is the workhorse for computation — roughly $n^3/3$ operations versus a more expensive full eigendecomposition.