---
# Section 1.7: Gaussian Elimination and the LU Decomposition
---

Let $A \in \mathbb{R}^{n \times n}$ be **nonsingular**, and let $b \in \mathbb{R}^n$.

We will now look at solving

$$
Ax = b
$$

by **Gaussian Elimination**:

> Transform $Ax = b$ into an **equivalent** system
>
> $$Ux = y,$$
>
> where $U$ is **upper-triangular**, by using the following **elementary row operations**:
>
> 1. Add a multiple of one equation to another equation.
> 2. Interchange two equations.
> 3. Multiply an equation by a nonzero constant.

This can be done using **augmented matrices**:

$$
\left[\begin{array}{c|c}
A & b
\end{array}\right]
\rightarrow \cdots \rightarrow
\left[\begin{array}{c|c}
U & y
\end{array}\right].
$$

We know that $Ax = b$ and $Ux = y$ will be equivalent (i.e., have exactly the same set of solutions) by the following fact.

> ### Proposition 1:
Let $\hat{A}x = \hat{b}$ be obtained by an elementary row operation on $Ax = b$. Then $Ax = b$ and $\hat{A}x = \hat{b}$ are equivalent .
>
> Moreover, 
> $$
\text{$A$ is nonsingular} 
\quad \Longleftrightarrow \quad
\text{$\hat{A}$ is nonsingular}.
$$

---

## Example 1

Use Gaussian elimination to solve $Ax = b$, where

$$
A = 
\begin{bmatrix}
2 & 1 & 1 \\
2 & 2 & -1 \\
4 & -1 & 6
\end{bmatrix},
\qquad
b = 
\begin{bmatrix}
7 \\ 3 \\ 20
\end{bmatrix}.
$$

## Solution

$$
\begin{array}{rcl}
\left[\begin{array}{c|c}
A & b
\end{array}\right] = 
\left[\begin{array}{ccc|c}
2 & 1 & 1 & 7 \\
2 & 2 & -1 & 3 \\
4 & -1 & 6 & 20
\end{array}\right]
& \xrightarrow{(r_2 \gets r_2 - r_1)} &
\left[\begin{array}{ccc|c}
2 & 1 & 1 & 7 \\
0 & 1 & -2 & -4 \\
4 & -1 & 6 & 20
\end{array}\right] \\
& \xrightarrow{(r_3 \gets r_3 - 2r_1)} &
\left[\begin{array}{ccc|c}
2 & 1 & 1 & 7 \\
0 & 1 & -2 & -4 \\
0 & -3 & 4 & 6
\end{array}\right] \\
& \xrightarrow{(r_3 \gets r_3 - (-3)r_2)} &
\left[\begin{array}{ccc|c}
2 & 1 & 1 & 7 \\
0 & 1 & -2 & -4 \\
0 & 0 & -2 & -6
\end{array}\right] = 
\left[\begin{array}{c|c}
U & y
\end{array}\right].
\end{array}
$$

Then, by backward substitution, we have

$$
x = \begin{bmatrix} 1 \\ 2 \\ 3 \end{bmatrix}.
$$

---

We solved the example above using only the first type of row operation, but this is not always possible.

---

## Example 2

Use Gaussian elimination to solve $Ax = b$, where

$$
A = 
\begin{bmatrix}
2 & 2 & 1 \\
2 & 2 & -1 \\
4 & -1 & 6
\end{bmatrix},
\qquad
b = 
\begin{bmatrix}
5 \\ 3 \\ 14
\end{bmatrix}.
$$

## Solution

$$
\begin{array}{rcl}
\left[\begin{array}{c|c}
A & b
\end{array}\right] = 
\left[\begin{array}{ccc|c}
2 & 2 & 1 & 5 \\
2 & 2 & -1 & 3 \\
4 & -1 & 6 & 14
\end{array}\right]
& \xrightarrow{(r_2 \gets r_2 - r_1)} &
\left[\begin{array}{ccc|c}
2 & 1 & 1 & 5 \\
0 & 0 & -2 & -2 \\
4 & -1 & 6 & 14
\end{array}\right] \\
& \xrightarrow{(r_3 \gets r_3 - 2r_1)} &
\left[\begin{array}{ccc|c}
2 & 1 & 1 & 5 \\
0 & 0 & -2 & -2 \\
0 & -3 & 4 & 4
\end{array}\right] \\
& \xrightarrow{(r_2 \leftrightarrow r_3)} &
\left[\begin{array}{ccc|c}
2 & 1 & 1 & 5 \\
0 & -3 & 4 & 4 \\
0 & 0 & -2 & -2 \\
\end{array}\right] = 
\left[\begin{array}{c|c}
U & y
\end{array}\right].
\end{array}
$$

Then, by backward substitution, we have

$$
x = \begin{bmatrix} 2 \\ 0 \\ 1 \end{bmatrix}.
$$

---

To do Gaussian elimination **without row interchanges**, the matrix $A$ needs the following property:

$$
\begin{equation} \tag{1}
\fbox{$A_k$ is nonsingular, for $k = 1,\ldots,n$.}
\end{equation}
$$

Here $A_k$ is the **$\mathbf{k^{th}}$ leading principal submatrix** (i.e., $A_k =$ `A[1:k, 1:k]`).

---

## Example 1 (revisited)

$$
A = 
\begin{bmatrix}
2 & 1 & 1 \\
2 & 2 & -1 \\
4 & -1 & 6
\end{bmatrix}
$$


```julia
using LinearAlgebra
```


```julia
A = [2 1 1; 2 2 -1; 4 -1 6]
```


```julia
A1 = A[1:1, 1:1]
```


```julia
det(A1)
```


```julia
A2 = A[1:2, 1:2]
```


```julia
det(A2)
```


```julia
A3 = A[1:3, 1:3]
```


```julia
det(A3)
```

Therefore, the matrix $A$ satisfies property (1).

---

## Example 2 (revisited)

$$
A = 
\begin{bmatrix}
2 & 2 & 1 \\
2 & 2 & -1 \\
4 & -1 & 6
\end{bmatrix}
$$


```julia
A = [2 2 1; 2 2 -1; 4 -1 6]
```


```julia
A1 = A[1:1, 1:1]
```


```julia
det(A1)
```


```julia
A2 = A[1:2, 1:2]
```


```julia
det(A2)
```

Therefore, the matrix $A$ does not satisfy property (1).

---

>### Proposition 2:
Let $A \in \mathbb{R}^{n \times n}$ be nonsingular.
>
> If we can solve $Ax = b$ using Gaussian elimination **without row interchanges**, then property (1) holds.

### Proof.

Suppose we do the Gaussian elimination 

$$
\left[\begin{array}{c|c}
A & b
\end{array}\right]
\rightarrow \cdots \rightarrow
\left[\begin{array}{c|c}
U & y
\end{array}\right]
$$

without row interchanges. Then $U$ is an upper-triangular matrix which is nonsingular since $A$ is nonsingular.

Since $U$ is nonsingular and triangular, the diagonal entries of $U$ are nonzero. Therefore, $U_k$ is nonsingular, for $k = 1,\ldots,n$.

Since no row interchanges occurred, we have that $U_k$ was obtained from $A_k$ by elementary row operations. Since $U_k$ is nonsingular, we have that $A_k$ is nonsingular.

Therefore, property (1) holds. $\blacksquare$

---

### The converse of Proposition 2

We will now prove that the converse of Proposition 2 holds:

> If property (1) holds for the matrix $A$, then we can do Gaussian elimination on $Ax = b$ without row interchanges.

### Proof.

First we note that if property (1) holds, then $a_{11} \neq 0$.

We will prove this statement by using **block Gaussian elimination**:

$$
A = 
\begin{bmatrix}
a_{11} & b^T \\
c & \hat{A}
\end{bmatrix}
\quad \xrightarrow{\left(r_2 \gets r_2 - \frac{1}{a_{11}} c r_1\right)} \quad
\begin{bmatrix}
a_{11} & b^T \\
0 & \hat{A} - \frac{1}{a_{11}} cb^T
\end{bmatrix}.
$$

We let 

$$
\tilde{A} = \hat{A} - \frac{1}{a_{11}} cb^T,
$$

which is known as the **Schur complement** of $a_{11}$ in $A$.

The next step of Gaussian elimination will repeat the above row operation on $\tilde{A}$.

Thus, if we can show that $\tilde{A}$ inherits property (1) from $A$, then we know that we can continue the process of Gaussian elimination on $Ax = b$ without row interchanges.

So we need to show that $\tilde{A}_k$ is nonsingular for $k = 1,\ldots,n-1$. Suppose, for the sake of contradiction, that $\tilde{A}_k$ is singular, for some $k$.

Then, $\exists \tilde{y} \in \mathbb{R}^k$ nonzero such that 

$$
\tilde{A}_k \tilde{y} = 0.
$$

Letting $M = \begin{bmatrix} a_{11} & b^T \\ 0 & \tilde{A} \end{bmatrix}$, we have

$$
M_{k+1} = 
\begin{bmatrix} a_{11} & b_k^T \\ 0 & \tilde{A}_k \end{bmatrix},
$$

where $b_k = b[1 \colon k]$.

Since $M_{k+1}$ was obtained from $A_{k+1}$ by elementary row operations, and since $A_{k+1}$ is nonsingular, we have that $M_{k+1}$ is nonsingular.

Now, let's get a contradiction by finding a vector $y \neq 0$ such that $M_{k+1} y = 0$.

Let 

$$
y = 
\begin{bmatrix}
-\frac{1}{a_{11}}b_k^T \tilde{y} \\
\tilde{y}
\end{bmatrix}.
$$

Then $y \neq 0$ since $\tilde{y} \neq 0$, and

$$
M_{k+1} y = 
\begin{bmatrix} a_{11} & b_k^T \\ 0 & \tilde{A}_k \end{bmatrix}
\begin{bmatrix}
-\frac{1}{a_{11}}b_k^T \tilde{y} \\
\tilde{y}
\end{bmatrix} =
\begin{bmatrix}
-b_k^T \tilde{y} + b_k^T \tilde{y} \\
\tilde{A}_k \tilde{y}
\end{bmatrix} =
\begin{bmatrix}
0 \\ 0
\end{bmatrix},
$$

which gives us our contradiction to the fact $M_{k+1}$ should be nonsingular.

Therefore, $\tilde{A}_k$ is nonsingular for $k = 1,\ldots,n-1$, which implies we can continue to do Gaussian elimination on $Ax = b$ without row interchanges. $\blacksquare$

---

Thus, we have the following result:

> ### Proposition 3:
Let $A \in \mathbb{R}^{n \times n}$ be nonsingular. Then we can solve $Ax = b$ using Gaussian elimination without row interchanges if and only if all leading principal submatrices of $A$ are nonsingular.

---

## LU Decomposition

We will now see that Gaussian elimination without row interchanges gives us a decomposition of the matrix $A$ into a lower-triangular matrix $L$ times an upper-triangular matrix $U$:

$$
A = LU.
$$

Recall that in Example 1, we had

$$
\begin{array}{rcl}
\left[\begin{array}{c|c}
A & b
\end{array}\right] 
= 
\left[\begin{array}{ccc|c}
2 & 1 & 1 & 7 \\
2 & 2 & -1 & 3 \\
4 & -1 & 6 & 20
\end{array}\right]
& \xrightarrow{(r_2 \gets r_2 - r_1)} &
\left[\begin{array}{ccc|c}
2 & 1 & 1 & 7 \\
0 & 1 & -2 & -4 \\
4 & -1 & 6 & 20
\end{array}\right] \\
& \xrightarrow{(r_3 \gets r_3 - 2r_1)} &
\left[\begin{array}{ccc|c}
2 & 1 & 1 & 7 \\
0 & 1 & -2 & -4 \\
0 & -3 & 4 & 6
\end{array}\right] \\
& \xrightarrow{(r_3 \gets r_3 - (-3)r_2)} &
\left[\begin{array}{ccc|c}
2 & 1 & 1 & 7 \\
0 & 1 & -2 & -4 \\
0 & 0 & -2 & -6
\end{array}\right] 
= 
\left[\begin{array}{c|c}
U & y
\end{array}\right].
\end{array}
$$

Let's look carefully at what is happening to the right-hand-side.

We can see that $y$ has been computed from $b$ as follows:

\begin{array}{lllllr}
y_1 &=&  b_1 &=& 7 & = & 7, \\
y_2 &=&  b_2 - y_1 &=& 3 - 7 &=& -4,\\
y_3 &=&  b_3 - 2y_1 - (-3)y_2 &=& 20 - 2 \cdot 7 + 3 \cdot (-4) &=& -6.
\end{array}

Putting all the $y$ terms on the left-hand-side, we obtain the following linear system.

\begin{array}{rrrrrrr}
y_1 & &       & &     &=&  b_1 \\
y_1 &+& y_2   & &     &=&  b_2 \\
2y_1 &-& 3y_2 &+& y_3 &=&  b_3
\end{array}

We can write these equations in matrix form as

$$
\begin{bmatrix}
1 & 0 & 0 \\
1 & 1 & 0 \\
2 & -3 & 1 \\
\end{bmatrix}
\begin{bmatrix}
y_1 \\ y_2 \\ y_3
\end{bmatrix}
= 
\begin{bmatrix}
b_1 \\ b_2 \\ b_3
\end{bmatrix}
\qquad
(Ly = b),
$$

where $L$ is lower-triangular with ones on the diagonal. We say that $L$ is a **unit lower-triangular matrix**.

Thus, we can solve $Ax = b$ using the following algorithm:

1. Solve $Ly = b$ for $y$.
2. Solve $Ux = y$ for $x$.

This implies that $LUx = b$, and in fact

$$
LU = 
\begin{bmatrix}
1 &  &  \\
1 & 1 &  \\
2 & -3 & 1 \\
\end{bmatrix}
\begin{bmatrix}
2 & 1 & 1 \\
 & 1 & -2 \\
 &  & -2
\end{bmatrix}
=
\begin{bmatrix}
2 & 1 & 1  \\
2 & 2 & -1 \\
4 & -1 & 6
\end{bmatrix}
= A
$$

Therefore, Gaussian elimination without row interchanges gives us the $LU$-factorization

$$
A = LU,
$$

where $L$ is **unit lower-triangular** and $U$ is **upper-triangular**.

Since $L$ is known to have ones along its diagonal, we do not need to store the diagonal of $L$. This allows us to store the $LU$-factorization of $A$ in a single matrix:

$$
"LU" \quad = \quad 
\begin{bmatrix}
2 & 1 & 1 \\
1 & 1 & -2 \\
2 & -3 & -2
\end{bmatrix}
$$

In fact, we can overwrite the matrix $A$ with its $L$ and $U$ factors.

---

## `lu`


```julia
?lu
```


```julia
A = Float64[2 1 1; 2 2 -1; 4 -1 6]
b = Float64[7; 3; 20]
pivot = Val(false)

F = lu(A, pivot)
```


```julia
dump(F)
```


```julia
F.factors
```


```julia
L = F.L
```


```julia
U = F.U
```


```julia
L*U - A
```


```julia
x = F\b
```


```julia
x = U\(L\b)
```

---

As mentioned above, if you don't need $A$ after doing the $LU$-factorization, and you want to save some memory, you can have $A$ **overwritten** by the $L$ and $U$ factors.

This is done using the `lu!` function.


```julia
A = Float64[2 1 1; 2 2 -1; 4 -1 6]
```


```julia
pivot = Val(false)
lu!(A, pivot)
```


```julia
A
```

---


```julia
A = Float64[2 2 1; 2 2 -1; 4 -1 6]
```


```julia
L, U, p = lu(A, Val(false))
```


```julia
L, U, p = lu(A)
```


```julia
L*U
```


```julia
A
```


```julia
p
```


```julia
A[p,:]
```


```julia
L*U - A[p,:]
```

---

## A larger numerical example


```julia
n = 3000

A = randn(n, n) + n*I
x = rand(n)
b = A*x;
```


```julia
@time F = lu!(A);
```


```julia
xhat = copy(b)
@time ldiv!(F, xhat);
```


```julia
x ≈ xhat
```

---

## Outer-product form of the $LU$-algorithm

We partition $A = LU$ as 

$$
\begin{bmatrix}
a_{11} & b^T \\
c & \hat{A} \\
\end{bmatrix}
=
\begin{bmatrix}
1 & 0 \\
l & \hat{L} \\
\end{bmatrix}
\begin{bmatrix}
u_{11} & u^T \\
0 & \hat{U} \\
\end{bmatrix}
=
\begin{bmatrix}
u_{11} & u^T \\
u_{11}l & lu^T + \hat{L}\hat{U} \\
\end{bmatrix}.
$$

This gives us the following algorithm.

### Algorithm:

1. Let $u_{11} = a_{11}$ and $u^T = b^T$.
2. Let $l = c/u_{11}$.
3. Let $\tilde{A} = \hat{A} - lu^T$.
4. Factor $\tilde{A} = \hat{L}\hat{U}$.

### Flop count of the $LU$-algorithm

We can use the above recursive algorithm to count the number of flops.

Let $f_n$ be the number of flops required to compute the $LU$-decomposition of an $n \times n$ matrix $A$.

- Step 1 requires 0 flops.
- Step 2 requires $n-1$ flops.
- Step 3 requires $2(n-1)^2$ flops.
- Step 4 requires $f_{n-1}$ flops.

Therefore,

$$
f_n = (n-1) + 2(n-1)^2 + f_{n-1}, \qquad n = 2, 3, 4, \ldots.
$$

Also, $f_1 = 0$.

Letting $g(n) = 2n^2 + n$, we have

$$
\begin{align}
f_n 
&= g(n-1) + f_{n-1} \\
&= g(n-1) + g(n-2) + f_{n-2} \\
&= g(n-1) + g(n-2) + g(n-3) + f_{n-3} \\
&\quad\vdots \\
&= g(n-1) + g(n-2) + g(n-3) + \cdots + g(1) + f_1 \\
&= \sum_{i=1}^{n-1} g(i) + 0 \\
&= \sum_{i=1}^{n-1} \left( 2i^2 + i \right) \\
&= \int_0^n \left(2x^2 + x\right)dx + O(n^2) \\
&= \fbox{$\frac23 n^3 + O(n^2)$}. \\
\end{align}
$$

So, we conclude that Gaussian elimination requires about $\frac23 n^3$ flops (i.e., roughly twice as many flops as Cholesky's algorithm).

---

## The LU Decomposition Theorem

> ### Theorem: (LU Decomposition)
>
> Let $A$ be an $n \times n$ matrix whose leading principal submatrices are nonsingular.
>
> Then there is a **unique** factorization
>
> $$A = LU$$
>
> where $L$ is unit lower-triangular and $U$ is upper-triangular.

### Proof.

We will prove this theorem by induction on $n$.

If $n = 1$, then $A = [a_{11}]$. If $A = LU$, where $L$ is unit lower-triangular and $U$ is upper-triangular, then $L = [1]$ and $U = [a_{11}]$. Therefore, $A$ has a unique $LU$-decomposition.

Now suppose that all $(n-1) \times (n-1)$ matrices with nonsingular leading principal submatrices have a unique $LU$-decomposition. This is our *induction hypothesis*.



Let $A$ be an $n \times n$ matrix with nonsingular leading principal submatrices and consider the following "bordered" partition of $A = LU$.

$$
\tag{2}
\begin{bmatrix}
\hat{A} & b \\
c^T & a_{nn}
\end{bmatrix}
=
\begin{bmatrix}
\hat{L} & 0 \\
l^T & 1
\end{bmatrix}
\begin{bmatrix}
\hat{U} & u \\
0 & u_{nn}
\end{bmatrix}
$$

We want to show that there is a unique solution to equation (2), where $\hat{L}$ is unit lower-triangular and $\hat{U}$ is upper-triangular.

We can rewrite equation (2) as:

$$
\hat{A} = \hat{L}\hat{U}, \quad b = \hat{L}u, \quad c^T = l^T\hat{U}, \quad a_{nn} = l^Tu + u_{nn}.
$$

- Since $A$ has nonsingular leading principal submatrices, so does $\hat{A}$. And since $\hat{A}$ is $(n-1) \times (n-1)$, we know that $\hat{A}$ has a unique $LU$-decomposition, by our induction hypothesis. Thus, there exist a **unique** $\hat{L}$ that is unit lower-triangular and a **unique** $\hat{U}$ that is upper-triangular such that $\hat{A} = \hat{L}\hat{U}$.

- Since $\hat{L}$ is nonsingular, there exists a **unique** solution $u$ to the equation $\hat{L}u = b$.

- Also, $\hat{U}$ is nonsingular because $\hat{A}$ is nonsingular. Therefore, there exists a **unique** solution $l$ to the equation $\hat{U}^Tl = c$.

- Finally, $u_{nn} = a_{nn} - l^Tu$, so $u_{nn}$ is also **uniquely determined**.

Therefore, $A$ has a **unique** $LU$-decomposition. $\blacksquare$

---

e38a8074a1736ade0497304beebf6cef052c7de3
ipynb
Jupyter Notebook
Section 1.7 - Gaussian Elimination and the LU Decomposition.ipynb
jesus-lua-8/fall2021math434
e6c6931676922193a6718b51c92232a3de7fa650
MIT
2021-08-31T21:01:22.000Z
2021-08-31T21:01:22.000Z
Section 1.7 - Gaussian Elimination and the LU Decomposition.ipynb
jesus-lua-8/fall2021math434
e6c6931676922193a6718b51c92232a3de7fa650
MIT
Section 1.7 - Gaussian Elimination and the LU Decomposition.ipynb
jesus-lua-8/fall2021math434
e6c6931676922193a6718b51c92232a3de7fa650
MIT
2021-11-16T19:28:56.000Z
2021-11-16T19:28:56.000Z
Qwen/Qwen-72B
1. YES
2. YES