M.1.3 Linear Systems & Matrix Decompositions (LU, Cholesky)

Basic Prerequisites: 02 · Vectors, Spaces & Matrices.

---

1. Intuition & Practical Objective

Nearly every computation in the toolbox reduces to solving $Ax=b$ - find the weight/parameter vector $b$‘s pre-image. Regression’s normal equations $X^{\prime }X\beta =X^{\prime }y$, portfolio optimization’s quadratic-program KKT system, and Monte Carlo’s correlated-normal draw $X=\mu +LZ$ are all linear solves. The practical objective of this page: learn not to solve by brute-force inversion, but to factor once and solve twice. A factorization (LU, Cholesky) splits $A$ into easy triangular pieces, so solving $Ax=b$ becomes two $O(n^2)$ back/forward substitutions - and, critically, the same factorization solves a hundred different right-hand sides $b$ for free.

The single most finance-relevant special case is the Cholesky factorization of a covariance matrix (Glasserman §2.2):

$$\Sigma =L{L}^{\prime },X=\mu +LZ,Z\sim N(0,I),$$

which is the standard way to simulate correlated asset returns: draw independent normals $Z$, then rotate-and-scale by $L$ to get the correct covariance. If your covariance is not positive definite, Cholesky throws an error - which is exactly the signal that the covariance is broken (06).

---

2. Mathematical Ground Truth & Derivations

Gaussian elimination / LU. Any invertible square matrix can be reduced to an upper-triangular $U$ by row operations encoded in a unit lower-triangular $L$ (with partial pivoting, $PA=LU$). Then $Ax=b$ becomes

$$Ly=Pb,Ux=y,$$

two triangular solves - each trivially $O(n^2)$. Flop cost of the factorization is $O(n^3)$.

Cholesky (Glasserman eqs. 2.29–2.31). For a symmetric positive definite $\Sigma$, there is a unique lower-triangular $L$ with positive diagonal such that $\Sigma =L{L}^{\prime }$. The entries come from the recursion

$$L_{ii}=\sqrt{{\Sigma }_{ii}-\sum _{k<i}{L}_{ik}^2},L_{ij}=\frac{{\Sigma }_{ij}-{\sum }_{k<j}{L}_{ik}{L}_{jk}}{L_{jj}}(j<i).$$

For the $2\times 2$ case this is the textbook result $L=\left(\begin{array}{cc}{\sigma }_1& 0\\ \rho {\sigma }_2& \sqrt{1-{\rho }^2}{\sigma }_2\end{array}\right)$, so $X_2={\mu }_2+{\sigma }_2\rho Z_1+{\sigma }_2\sqrt{1-{\rho }^2}{Z}_2$. Cholesky costs about $\frac{1}{3}{n}^3$ flops - half of LU, because it exploits symmetry. If any ${\Sigma }_{ii}-{\sum }_{k<i}{L}_{ik}^2$ goes non-positive, the matrix is not PD and the recursion hits a square root of a negative number: the first-principles failure of Cholesky on a non-PSD covariance.

The eigen/PC alternative (Glasserman eq. 2.32). Instead of $L$, use $\Sigma =V\Lambda V^{\prime }$ and take $A=V{\Lambda }^{1/2}$; $A{A}^{\prime }=\Sigma$ too. This is the PCA route to the same simulation and connects directly to 04 and 05.

---

3. Computational Implementation - Cholesky from scratch, then use it

Implementing the recursion above and then sampling correlated normals with it demonstrates the whole loop: factor, verify $L{L}^{\prime }=\Sigma$, draw, and confirm the empirical covariance matches the target. Stdlib only.

The factorization reproduces $L{L}^{\prime }=\Sigma$ exactly (it is an algebraic identity, not an approximation), and the Monte Carlo covariance converges to the target - this is the machinery inside every multi-asset path generator (Glasserman Ch 3). The same $L$ is then reused for thousands of draws at negligible marginal cost, which is why you factor once and sample many times.

---

4. Failure Modes & First-Principles Breakdowns

1. The non-PSD crash is a diagnosis, not a bug. Cholesky of a non-PSD covariance hits $\sqrt{\text{negative}}$ - the routine throws (Matrix is not positive definite). That exception is telling you the covariance violates ${\lambda }_i\ge 0$, which happens with $N>T$, pairwise-assembled correlations, or missing-data covariance estimates. Fix the covariance (06, shrinkage), not the solver.
2. Naive inversion is wasteful and unstable. Computing $A^{-1}$ to solve $Ax=b$ costs $O(n^3)$ and squares the round-off; the factor-and-substitute route is cheaper and more accurate. Never form $X^{\prime }X$ explicitly when $X$ is wide or ill-conditioned (05 has the SVD alternative).
3. No pivoting $\Rightarrow$ division-by-near-zero. LU without partial pivoting can divide by an exactly/tiny pivot and blow up. Cholesky avoids pivoting only because PD guarantees positive diagonal pivots - another reason PD-ness matters.

---

5. References

• Glasserman, Monte Carlo Methods in Financial Engineering
• Strang, Introduction to Linear Algebra (5th ed.)
• Horn & Johnson, Matrix Analysis
• Tsay, Analysis of Financial Time Series

---

6. Connected Graph Bridges

• Back: 02 · Vectors & Matrices · Index Hub
• Continue: 04 · Eigenvalues & Covariance · 05 · SVD, PCA & Regression
• Forward: Stochastic Calculus (Cholesky drives correlated-GBM/multi-asset SDE simulation) · Backtesting & Simulation Engines (production correlated-path generation) · Numerical Methods (05 · Numerical Linear Algebra)