Contents

• 1  Introduction
  • 1.1  LU decomposition
• 2  Exercise
  • 2.1  LU decomposition
  • 2.2  Forward substitution Ly = b
  • 2.3  Backward substitution Ux = y
• 3  Other considerations
  • 3.1  Differences between Gauss and LU decomposition
  • 3.2  LU decomposition with elementary matrices

Introduction

LU decomposition

We want to decompose the coefficient matrix $A$ into the product of a lower triangular matrix $L$ by an upper triangular matrix $U$

$$A=LU$$

Not all matrices have an $LU$ decomposition. If $A$ is invertible, then it admits an $LU$ factorization if and only if all its leading principal minors are nonzero.

$LU$ decomposition is not unique. For instance:

• Crout method: the elements of $U$ diagonal are ones.
• Doolittle method: the elements of $L$ diagonal are ones.
• Cholesky method (symmetric positive-definite matrices): diagonal elements are equal in $L$ and $U.$ That is, $u_{ii}=l_{ii}$

We will use Doolittle method and the elements of the diagonal of $L$ will be fixed as ones. In this case:

• $U$ is the upper triangular matrix obtained by Gaussian elimination from $A$.
• The elements of the $L$ matrix (below the main diagonal) are the factors we use in Gaussian row operations.

We consider the system of linear equations $A\mathbf{x}=\mathbf{b}$ where the matrix $A$ admits the $LU$ decomposition to solve the system,

1. Descompose $A=LU$. As $A\mathbf{x}=\mathbf{b}$ is transformed into $LU\mathbf{x}=\mathbf{b},$ if we call $\mathbf{y}$ to $U\mathbf{x}$ we can write the system $L\mathbf{y}=\mathbf{b}.$
2. Solve $L\mathbf{y}=\mathbf{b}$ by forward substitution.
3. Solve $U\mathbf{x}=\mathbf{y}$ by backward substitution.

Exercise

Consider the system $A\mathbf{x}=\mathbf{b}$ with

$$ A=\left(\begin{array}{crc} 1 & 1 & 3\\ 3 & 0 & 1\\ 1 & -2 & 1 \end{array}\right) \quad \mathbf{b}=\left(\begin{array}{c} -2\\ -1\\ -3 \end{array}\right) $$

Find $\mathbf{x}$ using $LU$ decomposition.

---

The steps are

1. Decomposition: $A=LU$
2. Forward substitution: $Ax=b \Longrightarrow LUx=b \Longrightarrow Ly = b$. We solve $Ly=b$.
3. Backward substitution: We solve $Ux=y$.

LU decomposition

In the first step, we make zeros under the element $a_{11}$ by adding the first row multiplied by a real.

$$ \begin{array}{c} r_{1}\\ r_{2}\\ r_{3} \end{array} \left(\begin{array}{rrr} \mathbf{1} & 1 & 3\\ 3 & 0 & 1\\ 1 & -2 & 1 \end{array}\right) \quad \begin{array}{lclcrr} r_{1} & \rightarrow & r_{1} & & & \\ r_{2} & \rightarrow & r_{2} & - & \color{red}{\fbox{3/1}} & r_{1}\\ r_{3} & \rightarrow & r_{3} & - & \color{blue}{\fbox{1/1}} & r_{1} \end{array} $$

The factors, which appear in boxes, are the elements of the $L$ matrix. We insert them into the array, instead of the zeros we have made to save space. The transformed matrix is

$$ \left(\begin{array}{rrr} 1 & 1 & 3\\ \color{red}{\fbox{3}} & \mathbf{-3}& -8\\ \color{blue}{\fbox{1}} & -3 & -2 \end{array}\right)\quad\begin{array}{lclcrr} r_{1} & \rightarrow & r_{1}\\ r_{2} & \rightarrow & r_{2}\\ r_{3} & \rightarrow & r_{3} & - & \color{ForestGreen}{\fbox{$(-3)/(-3)$}} & r_{2} \end{array} $$

We repeat the process creating zeros below $a'_{22}$ and we arrive at a matrix that simultaneously stores $L$ and $U$.

$$ \left(\begin{array}{rrr} 1 & 1 & 3\\ \color{red}{\fbox{3}} & -3 & -8\\ \color{blue}{\fbox{1}} & \color{ForestGreen}{\fbox{1}} & 6 \end{array}\right) $$

And the matrices $L$ and $U$ are:

$$ L=\left(\begin{array}{rrr} 1 & 0 & 0\\ \color{red}{\fbox{3}} & 1 & 0\\ \color{blue}{\fbox{1}} & \color{lime}{\color{ForestGreen}{\fbox{1}}} & 1 \end{array}\right)\quad U=\left(\begin{array}{rrr} 1 & 1 & 3\\ 0 & -3 & -8\\ 0 & 0 & 6 \end{array}\right) \qquad \mathbf{b}=\left(\begin{array}{c} -2\\ -1\\ -3 \end{array}\right) $$

Forward substitution Ly = b

$$ \begin{array}{rrrrrrr} & y_1 & & & & & = &-2\\ & 3 y_1& + &y_2 & & & = & -1\\ & y_1& + & y_2 & + & y_3 & = & -3 \end{array} $$

From the first equation $$ y_1 = -2$$ From the second $$ y_2 = -1 -3 y_1 = -1-3(-2) = 5$$ And from the third $$ y_3 = -3-y_1-y_2=-3-(-2)-5=-6$$ That is $$ y_1 = -2 \quad y_2 = 5 \quad y_3 = -6$$

Backward substitution Ux = y

$$ U=\left(\begin{array}{rrr} 1 & 1 & 3\\ 0 & -3 & -8\\ 0 & 0 & 6 \end{array}\right) \qquad y=\left(\begin{array}{r} -2\\ 5\\ -6 \end{array}\right) $$

$$ \begin{array}{rrrrrrr} & x_1 & + & x_2 &+ & 3x_3 & = &-2\\ & & - &3x_2 & - & 8x_3 & = & 5\\ & & & & & 6x_3 & =& -6& \end{array} $$

From the third equation $$ x_3 = -1$$ From the second one $$ x_2 = \frac{5 +8 x_3}{-3} = \frac{5 +8 (-1)}{-3} = 1$$ From the first one $$ x_1 = -2-x_2-3x_3=-2-1-3(-1)=0$$ And the solution is $$ \fbox{$x_1 = 0 \quad x_2 = 1 \quad x_3 = -1$}$$

---

Other considerations

Differences between Gauss and LU decomposition

$LU$ decomposition method is essentially the Gauss Elimination method. The $U$ matrix is the same upper triangular matrix that we obtain from Gauss. What are the differences?

• In Gauss, we include $\mathbf{b}$ in the row transformations, while in $LU$ decomposition we do not.
• In $LU$ decomposition, the matrix $L$ stores the transformations that we are performing on $\mathbf{b}$ (we calculate $\mathbf{y}$ as $L\mathbf{y}=\mathbf{b}$) to then apply them and obtain $\mathbf{y}$, which would be the vector $\mathbf{b}$ transformed with the Gaussian row operations.

So, if we do the same as with Gauss, why do we not work directly with the augmented matrix? Why do it in two steps?

The $LU$ decomposition method makes sense (compared to Gauss elimination) if we have a series of systems that we have to solve sequentially and where the solution of one determines the independent term of the next $\mathbf{b}$. For example

$$ 1.\quad A\mathbf{x}_1=\mathbf{x}_0,\quad A\mathbf{x}_2=\mathbf{x}_1, \quad A\mathbf{x}_3=\mathbf{x}_2,\dots\\ 2.\quad L\mathbf{y}_1=\mathbf{x}_0,\quad L\mathbf{y}_2=\mathbf{x}_1, \quad L\mathbf{y}_3=\mathbf{x}_2,\dots\\ 3.\quad U\mathbf{x}_1=\mathbf{y}_1,\quad U\mathbf{x}_2=\mathbf{y}_2, \quad U\mathbf{x}_3=\mathbf{y}_y,\dots\\ $$

In this case, with Gauss Elimination we would have to perform it once for each system, changing the $\mathbf{b}$, while with the $LU$ method we would do step 1 only once and repeat steps 2 and 3 for each system.

---

LU decomposition with elementary matrices

The sequence of row operations used to transform the matrix $A$ into an equivalent upper triangular matrix is

$$ A =\left( \begin{array}{ccc} 1 & 1 & 3 \\ 3 & 0 & 1 \\ 1 & -2 & 1 \\ \end{array} \right) $$

$r_2\longleftarrow r_2-(3)r_1$ $$ \left( \begin{array}{rrr} 1 & 1 & 3 \\ 0 & -3 & -8 \\ 1 & -2 & 1 \\ \end{array} \right) $$ $ r_3\longleftarrow r_3-(1)r_1 $ $$ \left( \begin{array}{ccc} 1 & 1 & 3 \\ 0 & -3 & -8 \\ 0 & -3 & -2 \\ \end{array} \right) $$ $ r_3\longleftarrow r_3-(1)r_2 $ $$ \left( \begin{array}{rrr} 1 & 1 & 3 \\ 0 & -3 & -8 \\ 0 & 0 & 6 \\ \end{array} \right) $$

These transformations are achieved left-multiplying $A$ by the elementary matrices

$$ r_2\longleftarrow r_2-(3)r_1 $$

.

$$ \left( \begin{array}{rrr} 1 & 0 & 0 \\ -3 & 1 & 0 \\ 0 & 0 & 1 \\ \end{array} \right) \left( \begin{array}{rrr} 1 & 1 & 3 \\ 3 & 0 & 1 \\ 1 & -2 & 1 \\ \end{array} \right) = \left( \begin{array}{rrr} 1 & 1 & 3 \\ 0 & -3 & -8 \\ 1 & -2 & 1 \\ \end{array} \right) $$

.

$$ r_3\longleftarrow r_3-(1)r_1 $$

.

$$ \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -1 & 0 & 1 \\ \end{array} \right) \left( \begin{array}{ccc} 1 & 0 & 0 \\ -3 & 1 & 0 \\ 0 & 0 & 1 \\ \end{array} \right) \left( \begin{array}{ccc} 1 & 1 & 3 \\ 3 & 0 & 1 \\ 1 & -2 & 1 \\ \end{array} \right) = \left( \begin{array}{ccc} 1 & 1 & 3 \\ 0 & -3 & -8 \\ 0 & -3 & -2 \\ \end{array} \right) $$

.

$$ r_3\longleftarrow r_3-(1)r_2 $$

.

$$ \left( \begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -1 & 1 \\ \end{array} \right) \left( \begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -1 & 0 & 1 \\ \end{array} \right) \left( \begin{array}{rrr} 1 & 0 & 0 \\ -3 & 1 & 0 \\ 0 & 0 & 1 \\ \end{array} \right) \left( \begin{array}{rrr} 1 & 1 & 3 \\ 3 & 0 & 1 \\ 1 & -2 & 1 \\ \end{array} \right) = \left( \begin{array}{rrr} 1 & 1 & 3 \\ 0 & -3 & -8 \\ 0 & 0 & 6 \\ \end{array} \right) $$

Each elementary matrix $G_1$, $G_2$ y $G_3$ corresponds with a row operation

$$G_3G_2G_1A=U$$

and therefore

$$A=(G_3G_2G_1)^{-1}U=G_1^{-1}G_2^{-1}G_3^{-1}U$$

If we call $L$ to the following product

$$L=G_1^{-1}G_2^{-1}G_3^{-1}$$

${\color{white}.}$ $$A=\text{LU}$$

where $L$ is a lower triangular matrix. Let us calculate this $L$ matrix

$$ G_1^{-1}= \left( \begin{array}{rrr} 1 & 0 & 0 \\ -3 & 1 & 0 \\ 0 & 0 & 1 \\ \end{array} \right)^{-1} = \left( \begin{array}{rrr} 1 & 0 & 0 \\ 3 & 1 & 0 \\ 0 & 0 & 1 \\ \end{array} \right) $$

.

$$ G_2^{-1}= \left( \begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -1 & 0 & 1 \\ \end{array} \right)^{-1} = \left( \begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 1 & 0 & 1 \\ \end{array} \right) $$

.

$$ G_3^{-1}=\left( \begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -1 & 1 \\ \end{array} \right)^{-1} = \left( \begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 1 & 1 \\ \end{array} \right) $$

And then

$$ G_1^{-1}G_2^{-1}= \left( \begin{array}{rrr} 1 & 0 & 0 \\ 3 & 1 & 0 \\ 0 & 0 & 1 \\ \end{array} \right) \left( \begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 1 & 0 & 1 \\ \end{array} \right) = \left( \begin{array}{rrr} 1 & 0 & 0 \\ 3 & 1 & 0 \\ 1 & 0 & 1 \\ \end{array} \right) $$

.

$$ G_1^{-1}G_2^{-1}G_3^{-1}= \left( \begin{array}{rrr} 1 & 0 & 0 \\ 3 & 1 & 0 \\ 0 & 0 & 1 \\ \end{array} \right) \left( \begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 1 & 0 & 1 \\ \end{array} \right) \left( \begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 1 & 1 \\ \end{array} \right) = \left( \begin{array}{rrr} 1 & 0 & 0 \\ 3 & 1 & 0 \\ 1 & 1 & 1 \\ \end{array} \right) $$

And we have already decomposed $A$ as $LU$

$$ A =\left( \begin{array}{ccc} 1 & 1 & 3 \\ 3 & 0 & 1 \\ 1 & -2 & 1 \\ \end{array} \right) = \left( \begin{array}{rrr} 1 & 0 & 0 \\ 3 & 1 & 0 \\ 1 & 1 & 1 \\ \end{array} \right) \left( \begin{array}{rrr} 1 & 1 & 3 \\ 0 & -3 & -8 \\ 0 & 0 & 6 \\ \end{array} \right) = LU $$