LU Decomposition

Solving $Ax=b$ by literally computing $A^{-1}$ is the way to think about it, but real numerical code almost never does that — it's wasteful and less stable. Instead, computers split $A$ into two simpler triangular pieces and solve two easy triangular systems in sequence. LU decomposition is that split, and it's the standard machinery quietly running behind almost every "solve this linear system" call in scientific computing.

An analogy: undoing a knot by tracking each step

Solving a general system all at once is like untangling a big knot in one motion — hard to see how. LU decomposition is like methodically recording, one crossing at a time, exactly how the knot was tied (that's the "L," lower-triangular, recording each elimination step) and what shape is left after all the tying (that's the "U," upper-triangular, the simplified end result). Once you have both records, undoing the whole knot for a new starting pull is fast: replay the L record forward, then unwind the U record backward — each in one simple pass, no more searching.

The mathematics, one symbol at a time

LU decomposition factors a square matrix $A$ as

$$A = LU,$$

where $L$ is lower-triangular (all entries above the diagonal are zero, with 1s on the diagonal) and $U$ is upper-triangular (all entries below the diagonal are zero). In words: $L$ records the row-elimination steps used to reduce $A$ to the simplified triangular form $U$ — exactly the steps you'd do by hand doing Gaussian elimination, just kept track of and reusable. To solve $Ax=b$, substitute: $LUx=b$. Let $y=Ux$, first solve $Ly=b$ for $y$ (easy: triangular systems solve top-to-bottom by simple substitution, no elimination needed), then solve $Ux=y$ for $x$ (easy again: triangular, bottom-to-top). Two easy triangular solves replace one hard general solve — and once $L$ and $U$ are computed, solving for many different right-hand sides $b$ is cheap, since the expensive factorization step only needs to be done once.

Worked example 1: factoring a 2×2 matrix by hand

Let $A = \begin{pmatrix} 4 & 3 \\ 6 & 3 \end{pmatrix}$. Eliminate the first column below the pivot: the multiplier is $6/4=1.5$; subtract 1.5 times row 1 from row 2: new row 2 is $(6-1.5(4),\ 3-1.5(3)) = (0, -1.5)$. So $U = \begin{pmatrix}4 & 3\\0&-1.5\end{pmatrix}$ and $L = \begin{pmatrix}1&0\\1.5&1\end{pmatrix}$ (the multiplier recorded below the diagonal). Check: $LU = \begin{pmatrix}1&0\\1.5&1\end{pmatrix}\begin{pmatrix}4&3\\0&-1.5\end{pmatrix} = \begin{pmatrix}4&3\\6&4.5-1.5\end{pmatrix}=\begin{pmatrix}4&3\\6&3\end{pmatrix}$ ✓.

Worked example 2: solving with the factorization, two right-hand sides

Solve $Ax=b_1$ with $b_1=(10,12)$. First $Ly=b_1$: $y_1=10$, $1.5y_1+y_2=12 \Rightarrow y_2=12-15=-3$. Then $Ux=y$: $-1.5x_2=-3 \Rightarrow x_2=2$; $4x_1+3(2)=10 \Rightarrow x_1=(10-6)/4=1$. So $x=(1,2)$. Now a second right-hand side arrives, $b_2=(7,9)$ — no need to refactor $A$ at all: $Ly=b_2$ gives $y_1=7$, $y_2=9-1.5(7)=-1.5$; $Ux=y$ gives $x_2=1$, $x_1=(7-3)/4=1$. Reusing the same $L,U$ made the second solve nearly free — exactly why LU is the standard choice when the same matrix (say, a fixed hedge structure) is solved against many different target vectors over time.

What this means in practice

LU decomposition is what's actually running under the hood of numpy.linalg.solve and virtually every production linear-solve routine, and it's also how determinants are computed efficiently in practice — $\det(A) = \det(L)\det(U)$, and triangular matrix determinants are just the product of their diagonal entries. Any repeated-solve scenario (recomputing hedge weights against a fixed factor structure but changing target exposures daily) benefits from factoring once and reusing $L,U$ repeatedly.