LU Decomposition

A procedure for decomposing an $N\times N$ matrix ${\hbox{\sf A}}$ into a product of a lower Triangular Matrix ${\hbox{\sf L}}$ and an upper Triangular Matrix ${\hbox{\sf U}}$,

$${\hbox{\sf L}}{\hbox{\sf U}}={\hbox{\sf A}}.$$ (1)

Written explicitly for a $3\times 3$ Matrix, the decomposition is

$$\left[{\matrix{l_{11} & 0 & 0\cr l_{21} & l_{22} & 0\cr l_{3... ...{21} & a_{22} & a_{23}\cr a_{31} & a_{32} & a_{33}\cr}}\right]$$ (2)

$$\left[{\matrix{l_{11}u_{11} & l_{11}u_{12} & l_{11}u_{13}\cr... ...21} & a_{22} & a_{23}\cr a_{31} & a_{32} & a_{33}\cr}}\right].$$ (3)

This gives three types of equations

$i<j \qquad l_{i1}u_{1j}+l_{i2}u_{2j}+\ldots+l_{ii}u_{ij}=a_{ij}$	(4)
$i=j \qquad l_{i1}u_{1j}+l_{i2}u_{2j}+\ldots+l_{ii}u_{jj}=a_{ij}$	(5)
$i>j \qquad l_{i1}u_{1j}+l_{i2}u_{2j}+\ldots+l_{ij}u_{jj}=a_{ij}.$	(6)

This gives $N^2$ equations for $N^2+N$ unknowns (the decomposition is not unique), and can be solved using Crout's Method. To solve the Matrix equation

$${\hbox{\sf A}}{\bf x}=({\hbox{\sf L}}{\hbox{\sf U}}){\bf x} = {\hbox{\sf L}}({\hbox{\sf U}}{\bf x})={\bf b},$$ (7)

first solve ${\hbox{\sf L}}{\bf y}={\bf b}$ for ${\bf y}$. This can be done by forward substitution

$$y_1 = {b_1\over l_{11}}$$ (8)

$$y_i = {1\over l_{ii}} \left({b_i-\sum_{j=1}^{i-1} l_{ij}y_j}\right)$$ (9)

for $i=2$, ..., $N$. Then solve ${\hbox{\sf U}}{\bf x}={\bf y}$ for ${\bf x}$. This can be done by back substitution

$$x_N = {y_N\over u_{NN}}$$ (10)

$$x_i = {1\over u_{ii}}\left({y_i-\sum_{j=i+1}^N u_{ij}x_j}\right)$$ (11)

for $i=N-1$, ..., $1$.

See also Cholesky Decomposition, QR Decomposition, Triangular Matrix

References

Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; and Vetterling, W. T. ``LU Decomposition and Its Applications.'' §2.3 in Numerical Recipes in FORTRAN: The Art of Scientific Computing, 2nd ed. Cambridge, England: Cambridge University Press, pp. 34-42, 1992.