LU factorization

LU factorization is a method for decomposing an $N\times N$ matrix $\mathbf {A}$ into the product of a lower triangular matrix $\mathbf {L}$ and an upper triangular matrix $\mathbf {U}$,

$\mathbf {L} \,\mathbf {U} =\mathbf {A}$

For a $3\times 3$ matrix, this becomes:

${\begin{bmatrix}l_{11}&0&0\\l_{12}&l_{22}&0\\l_{13}&l_{23}&l_{33}\\\end{bmatrix}}\cdot {\begin{bmatrix}u_{11}&u_{12}&u_{13}\\0&u_{22}&u_{23}\\0&0&u_{33}\\\end{bmatrix}}={\begin{bmatrix}a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33}\\\end{bmatrix}}$

This gives three types of equations

$i<j\qquad l_{i1}u_{1j}+l_{i2}u_{2j}+\dots +l_{ii}u_{ij}=a_{ij}$

$i=j\qquad l_{i1}u_{1j}+l_{i2}u_{2j}+\dots +l_{ii}u_{jj}=a_{ij}$

$i>j\qquad l_{i1}u_{1j}+l_{i2}u_{2j}+\dots +l_{ii}u_{jj}=a_{ij}$

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

$\mathbf {A} \mathbf {x} =(\mathbf {L} \mathbf {U} )\mathbf {x} =\mathbf {L} (\mathbf {U} \mathbf {x} )=\mathbf {b} ,$

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

$y_{1}={\frac {b_{1}}{l_{11}}}$

$y_{i}={\frac {1}{l_{ii}}}\left(b_{i}-\sum _{j=1}^{i-1}l_{ij}y_{j}\right)$

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

$x_{N}={\frac {y_{N}}{u_{NN}}}$

$x_{i}={\frac {1}{u_{ii}}}\left(y_{i}-\sum _{j=i+1}^{N}u_{ij}x_{j}\right)$

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