LU Decomposition

An $N\times N$ square matrix ${\bf A}$ can be decomposed into a product of a lower triangular matrix ${\bf L}$ and an upper triangular matrix ${\bf U}$,

$${\bf A}=\left[ \begin{array}{cccc} a_{11} & a_{12} & \cdots & a_{... ...ddots & \vdots \\ 0 & 0 & \cdots & u_{NN} \end{array} \right] ={\bf L} {\bf U}$$ (186)

and the equation system ${\bf A}{\bf x}={\bf b}$ can be trivially solved by two rounds of back substitution:

• Rewrite the equation as below and solve for ${\bf y}={\bf U}{\bf x}$ by back substitution (top-down):

  $${\bf A}{\bf x}={\bf L}{\bf U}{\bf x}={\bf L}{\bf y}={\bf b}$$ (187)

• Solve the equation below for ${\bf x}$ by back substitution (bottom-up):

  $${\bf U}{\bf x}={\bf y}$$ (188)

We now consider how to find ${\bf L}$ and ${\bf U}$. The ijth element of the equation ${\bf A}={\bf LU}$ can be written as

$$a_{ij}=\sum_{k=1}^N l_{ik} u_{kj}=\sum_{k=1}^{min(i,j)} l_{ik} u_{kj}, \;\;\;\;\;\;\;\;(i,j=1,\cdots,N)$$ (189)

Note that there are $N^2$ equations (one for each element) but $N^2+N$ unknowns in both ${\bf L}$ and ${\bf U}$. We can therefore arbitrarily set the all diagonal elements of ${\bf L}$ to 1: $l_{11}=\cdots=l_{NN}=1$.

Moreover, as ${\bf L}$ and ${\bf U}$ are both triangular, we do not have to solve an $N^2$ nonlinear equations, as they can be disentangled if they are solved in a specific order. Actually these equations can be solved in $N$ iterations, each carried out in two steps for the jth column ( $j=1,\cdots,N$):

• $j=1$, find $u_{11}$ in step 1 and $l_{21},\cdots,l_{N1}$ in step 2;
• $j=2$, find $u_{12}$ and $u_{22}$ in step 1 and $l_{31},\cdots,l_{N1}$ in step 2;
• $\cdots, \cdots, \cdots, \cdots, \cdots, \cdots$
• $j=N-1$, find $u_{1\;N-1},\cdots,u_{N\;N-1}$ in step 1 and $l_{N\;N-1}$ in step 2;
• $j=N$, find $u_{1N},\cdots,u_{N\;N}$ in step 1.

Note that in each iteration for the jth column, we find all $u_{ij}$ ($i\le j$) in the upper triangle of ${\bf U}$, and all $l_{ij}$ ($i>j$) in the lower triangle of ${\bf L}$:

• For all $i\le j$, i.e., $i=1,\cdots,j-1$ (upper triangle), find all $u_{ij}$ by

  $$a_{ij}=\sum_{k=1}^i l_{ik} u_{kj}$$ (190)

  If $i=1$ (when $j=1,\cdots,N$)

  $$a_{1j}=l_{11}u_{1j}=u_{1j},\;\;\;\;\;u_{1j}=a_{1j}$$ (191)

  We see that the first row of ${\bf U}$ is the same as the first two of ${\bf A}$.

  If $i=2$ (when $j=2,\cdots,N$),

  $$a_{2j}=l_{21}u_{1j}+l_{22}u_{2j},\;\;\;\;\;u_{2j}=a_{2j}-l_{21}u_{1j}$$ (192)

  where $l_{21}$ has been obtained in previous iteration for $j=1$.

  In general, for all $i=1,\cdots,j$, we have

  $$a_{ij}=\sum_{k=1}^i l_{ik}u_{kj},\;\;\;\;\;u_{ij}=a_{ij}-\sum_{k=1}^{i-1}l_{ik}u_{kj}$$ (193)

  All $l_{i1},\cdots,l_{ii}$ in the ith row of ${\bf L}$ have already been obtained from previous iterations of $j$.

• For all $i>j$, i.e., $i=j+1,\cdots,N$ (lower triangle), find all $l_{ij}$ by

  $$a_{ij}=\sum_{k=1}^j l_{ik} u_{kj}$$ (194)

  When $j=1$ and $i=2,\cdots,N$,

  $$a_{i1}=l_{i1}u_{11},\;\;\;\;\;\;l_{i1}=a_{i1}/u_{11}$$ (195)

  where $u_{11}$ has been found in step 1.

  When $j=2$ and $i=3,\cdots,N$,

  $$a_{i2}=l_{i1}u_{12}+l_{i2}u_{22},\;\;\;\;\;l_{i2} =\frac{1}{u_{22}}\left(a_{i2}-l_{i1}u_{12}\right)$$ (196)

  where $u_{12}$ and $u_{22}$ have been found in step 1.

  In general, when $j=1,\cdots,N$ and $i>j$, we have

  $$a_{ij}=\sum_{k=1}^{j-1} l_{ik}u_{kj},\;\;\;\;\;\;l_{i2} =\frac{1}{u_{jj}}\left(a_{ij}-\sum_{k=1}^{j-1} l_{ik}u_{kj}\right)$$ (197)

  where $u_{1j},\cdots,u_{jj}$ in the jth column of ${\bf U}$ have been found in step 1.

Example

$${\bf A}=\left[\begin{array}{lll} a_{11} & a_{12} & a_{13}\\ a_{2... ...rray}{ccc}u_{11}&u_{12}&u_{13}\\ 0&u_{22}&u_{23}\\ 0&0&u_{33}\end{array}\right]$$

• $j=1$ (for column 1)

  $$\begin{array}{ll} i=1: & u_{11}=a_{11}=1\\ i=2: & l_{21}=a_{... ...1}=a_{21}=2\\ i=3: & l_{31}=a_{31}/u_{11}=a_{31}=3 \end{array}$$

• $j=2$ (for column 2)

  $$\begin{array}{ll} i=1: & u_{12}=a_{12}=2\\ i=2: & u_{22}=a_{... ...}=(a_{32}-l_{31}u_{12})/u_{22}=(1-3\times 2)/(-1)=5 \end{array}$$

• $j=3$ (for column 3)

  $$\begin{array}{ll} i=1: & u_{13}=a_{13}=3\\ i=2: & u_{23}=a_{... ...{31}u_{13}-l_{32}u_{23}=2-3\times 3-5\times (-5)=18 \end{array}$$

Now the LU decomposition of ${\bf A}$ is:

$$\left[\begin{array}{lll} 1 & 2 & 3\\ 2 & 3 & 1\\ 3 & 1 & 2\end{ar... ...\left[\begin{array}{rrr} 1 & 2 & 3\\ 0 & -1 & -5\\ 0 & 0 & 18\end{array}\right]$$

The LU decomposition can be used to solve linear systems.

$${\bf A}{\bf x}={\bf L}{\bf U}{\bf x}={\bf L}{\bf y}={\bf b}, \;\;\;\;\;\; \;\;\;\;\;\;{\bf y}={\bf U}{\bf x}$$ (198)

We can first solve ${\bf L}{\bf y}={\bf b}$ for ${\bf y}$ by forward substitution, and then solve ${\bf y}={\bf U}{\bf x}$ for ${\bf x}$ by backward substitution:

$$\left[ \begin{array}{cccc} l_{11} & 0 & 0 & 0 \\ l_{21} & l_{22}... ...rray}\right]= \left[\begin{array}{c} b_1\\ b_2\\ \vdots\\ b_N\end{array}\right]$$ (199)

which can be easily solved for ${\bf y}$ 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), \;\;\;\;\;\;(i=2,\cdots,N)$$ (200)

And then ${\bf x}$ can be found by solving ${\bf U}{\bf x}={\bf y}$

$$\left[ \begin{array}{cccc} u_{11} & u_{12} & \cdots & u_{1N} \\ ... ...ay}\right] \left[\begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_N\end{array}\right]$$ (201)

which can be easily solved by backward 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), \;\;\;\;\;\;\;(i=N-1,\cdots,1)$$ (202)

Example

$$\left\{ \begin{array}{rrrr} x_1+&2x_2+&3x_3&=14\\ 2x_1+&3x_2+&x_3&=11\\ 3x_1+&x_2+&2x_3&=11\end{array} \right.$$

By LU decomposition, this system can be written as:

$${\bf A}{\bf x}= \left[\begin{array}{lll} 1 & 2 & 3\\ 2 & 3 & 1\\ ... ...end{array}\right]={\bf b}= \left[\begin{array}{l}14\\ 11\\ 11\end{array}\right]$$

We first solve the system

$${\bf L}{\bf y}=\left[\begin{array}{lll} 1 & 0 & 0\\ 2 & 1 & 0\\ ... ...end{array}\right]={\bf b}= \left[\begin{array}{l}14\\ 11\\ 11\end{array}\right]$$

to get

$${\bf y}=[14, -17, 54]^T$$ (203)

We then solve the system

$${\bf U}{\bf x}= \left[\begin{array}{rrr} 1 & 2 & 3\\ 0 & -1 & -5\... ...nd{array}\right]={\bf y}= \left[\begin{array}{r}14\\ -17\\ 54\end{array}\right]$$

to get ${\bf x}=[1, 2, 3]^T$.

The LU decomposition ${\bf A}={\bf L}{\bf U}$ can also be used to find the inverse ${\bf A}^{-1}$ that satisfies

$${\bf A}{\bf A}^{-1}={\bf I}=diag(1,\cdots,1)=\left[{\bf e}_1,\cdots,{\bf e}_n \right]$$ (204)

where ${\bf e}_i=[0,\cdots,0,1,0,\cdots,0]^T$ is the ith column vector of the identity matrix ${\bf I}$ with all element equal to 0 except the ith one which is 1. The ith column ${\bf x}_i$ of ${\bf A}^{-1}=[{\bf x}_1\cdots {\bf x}_n]$ can be found by solving the following

$${\bf A}{\bf x}_i={\bf e}_i,\;\;\;\;\;\;\;(i=1,\cdots,n)$$ (205)

Example

Find the inverse of

$${\bf A}= \left[\begin{array}{lll} 1 & 2 & 3\\ 2 & 3 & 1\\ 3 & 1 & 2\end{array}\right]$$

Note that

$${\bf A}{\bf A}^{-1}={\bf A}[{\bf x}_1\;{\bf x}_2\;{\bf x}_3]={\bf I} =[{\bf e}_1\;{\bf e}_2\;{\bf e}_3]$$

where ${\bf x}_1$, ${\bf x}_1$, and ${\bf x}_1$ are the three columns of ${\bf A}$, which can be found by solving three equation systems:

$${\bf A}{\bf x}_1={\bf e}_1=\left[\begin{array}{c}1\\ 0\\ 0\end{ar... ...\; {\bf A}{\bf x}_3={\bf e}_3=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]$$

We use LU decomposition to get

$${\bf A}[{\bf x}_1\;{\bf x}_2\;{\bf x}_3] ={\bf L}{\bf U}[{\bf x}_1\;{\bf x}_2\;{\bf x}_3]={\bf I}$$

i.e.,

$$\left[\begin{array}{lll} 1 & 0 & 0\\ 2 & 1 & 0\\ 3 & 5 & 1\end{ar... ...] =\left[\begin{array}{rrr} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\end{array}\right]$$

Solving these three equation systems we get

$${\bf x}_1=\frac{1}{18}\left[\begin{array}{r}-5\\ 1\\ 7\end{array}... ...\;\;\; {\bf x}_3=\frac{1}{18}\left[\begin{array}{r}7\\ -5\\ 1\end{array}\right]$$

i.e.,

$${\bf A}^{-1}=[{\bf x}_1, {\bf x}_2, {\bf x}_3]=\frac{1}{18} \left[\begin{array}{rrr} -5 & 1 & 7\\ 1 & 7 & -5\\ 7 & -5 & 1\end{array}\right]$$