We consider solving a system of linear equations by using Gaussian elimination. Given a nonsingular square matrix $A \in \mathbb{R}^{n \times n}$, and a column vector $\vec{b} = \mathbb{R}^{n}$, find $\vec{x} \in \mathbb{R}^{n}$ such that 
$$A\vec{\mathstrut x} = \vec{\mathstrut b}$$

# Backward and forward substitution
Assume $A$ is an upper triangular matrix, denoted as $U$, consider to solve:

$$U\vec{x} = \left[ 
\begin{array}{cccc}
u_{11} & u_{12} & \cdots & u_{1n} \\
       & u_{22} & \cdots & u_{2n} \\
       &        & \ddots & \vdots \\
       &        &        & u_{nn}
\end{array}\right] \cdot 
\left[ 
\begin{array}{c}
x_1 \\
x_2 \\
\vdots \\
x_n
\end{array} \right] = \left[ 
\begin{array}{c}
b_1 \\
b_2 \\
\vdots \\
b_n
\end{array} \right] = \vec{b}$$

$Algorithm$: **Backward substitution** 

    for j = n to 1  
      if u_{jj} = 0, then display('Error: Singular Matrix');  
      x_j = b_j /u_{jj};  
      for i = 1 to j − 1  
        b_i = b_i − u_{ij} x_j ;  
      end  
    end

And for lower triangular matrix, denoted as $L$, consider to solve:

$$L\vec{x} = \left[ 
\begin{array}{cccc}
l_{11} &        &        &        \\
l_{21} & l_{22} &        &        \\
\vdots & \vdots & \ddots &        \\
l_{n1} & l_{n2} & \cdots & l_{nn}
\end{array}\right] \cdot 
\left[ 
\begin{array}{c}
x_1 \\
x_2 \\
\vdots \\
x_n
\end{array} \right] = \left[ 
\begin{array}{c}
b_1 \\
b_2 \\
\vdots \\
b_n
\end{array} \right] = \vec{b}$$

$Algorithm$: **Forward substitution** 

    for j = 1 to n  
      if u_{jj} = 0, then display('Error: Singular Matrix');  
      x_j = b_j /u_{jj};  
      for i = j + 1 to n  
        b_i = b_i − u_{ij} x_j ;  
      end  
    end

# Simple examples of Gaussian elimination
Let's look at an example

$$A = \left[ 
\begin{array}{ccc}
2 & 2 & 2 \\
4 & 3 & 2 \\
4 & 6 & 4
\end{array}\right]$$

We apply Gaussian elimination by multiplying elementary matrices to the left of $A$, we have

$$\begin{align}
M_1 A = \left[ 
\begin{array}{ccc}
1 & 0 & 0 \\
-2 & 1 & 0 \\
-2 & 0 & 1
\end{array}\right] \cdot \left[ 
\begin{array}{ccc}
2 & 2 & 2 \\
4 & 3 & 2 \\
4 & 6 & 4
\end{array}\right] &= \left[ 
\begin{array}{ccc}
2 & 2 & 2 \\
0 & -1 & -2 \\
0 & 2 & 0
\end{array}\right] \\ 
M_2 A = \left[ 
\begin{array}{ccc}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 2 & 1
\end{array}\right] \cdot \left[ 
\begin{array}{ccc}
2 & 2 & 2 \\
0 & -1 & -2 \\
0 & 2 & 0
\end{array}\right] &= \left[ 
\begin{array}{ccc}
2 & 2 & 2 \\
0 & -1 & -2 \\
0 & 0 & -4
\end{array}\right]
\end{align}$$

which is upper triangular, so that now we have:

$$M_2 M_1 A = \left[ 
\begin{array}{ccc}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 2 & 1
\end{array}\right] \cdot \left[
\begin{array}{ccc}
1 & 0 & 0 \\
-2 & 1 & 0 \\
-2 & 0 & 1
\end{array}\right] \cdot \left[ 
\begin{array}{ccc}
2 & 2 & 2 \\
4 & 3 & 2 \\
4 & 6 & 4
\end{array}\right] = \left[ 
\begin{array}{ccc}
2 & 2 & 2 \\
0 & -1 & -2 \\
0 & 0 & -4
\end{array}\right] = U$$

$i.e.$

$$\begin{align}
A = M_{1}^{-1} M_{2}^{-1} U &= \left[
\begin{array}{ccc}
1 & 0 & 0 \\
-2 & 1 & 0 \\
-2 & 0 & 1
\end{array}\right]^{-1} \cdot \left[ 
\begin{array}{ccc}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 2 & 1
\end{array}\right]^{-1} \cdot \left[ 
\begin{array}{ccc}
2 & 2 & 2 \\
0 & -1 & -2 \\
0 & 0 & -4
\end{array}\right] \\
&= \left[
\begin{array}{ccc}
1 & 0 & 0 \\
2 & 1 & 0 \\
2 & 0 & 1
\end{array}\right] \cdot \left[ 
\begin{array}{ccc}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & -2 & 1
\end{array}\right] \cdot \left[ 
\begin{array}{ccc}
2 & 2 & 2 \\
0 & -1 & -2 \\
0 & 0 & -4
\end{array}\right] \\
&= \left[
\begin{array}{ccc}
1 & 0 & 0 \\
2 & 1 & 0 \\
2 & -2 & 1
\end{array}\right] \cdot \left[ 
\begin{array}{ccc}
2 & 2 & 2 \\
0 & -1 & -2 \\
0 & 0 & -4
\end{array}\right] = LU
\end{align}$$

where $L$ is *lower triangular* and $U$ is *upper triangular*. We call $LU$ the **LU factorization** of $A$.

# Elementary matrices
This in a different elementary matrics with the following structure

$$M_{k} = \left[ 
\begin{array}{ccccccc}
1 & \cdots & 0 & 0 & 0 & \cdots & 0 \\
\vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\
0 & \cdots & 1 & 0 & 0 & \cdots & 0 \\
0 & \cdots & 0 & 1 & 0 & \cdots & 0 \\
0 & \cdots & 0 & m_{k+1, k} & 1 & \cdots & 0 \\
\vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\
0 & \cdots & 0 & m_{n,k} & 0 & \cdots & 1
\end{array}\right] $$

where the only nontrivial entries are located below the $(k, k)$ diagonal on column $k$.

$Conclusion$
1. Gaussian elimination
$$M_{k} \vec{a} = \left[ 
\begin{array}{ccccccc}
1 & \cdots & 0 & 0 & 0 & \cdots & 0 \\
\vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\
0 & \cdots & 1 & 0 & 0 & \cdots & 0 \\
0 & \cdots & 0 & 1 & 0 & \cdots & 0 \\
0 & \cdots & 0 & -a_{k+1}/a_{k} & 1 & \cdots & 0 \\
\vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\
0 & \cdots & 0 & -a_{n}/a_{k} & 0 & \cdots & 1
\end{array}\right] \cdot \left[ 
\begin{array}{c}
a_1 \\
a_2 \\
\vdots \\
a_{k} \\
a_{k+1} \\
\vdots \\
a_n
\end{array} \right] = \left[ 
\begin{array}{c}
a_1 \\
a_2 \\
\vdots \\
a_{k} \\
0 \\
\vdots \\
0
\end{array} \right] $$

2. Inverse
$$M_{k}^{-1} = \left[ 
\begin{array}{ccccccc}
1 & \cdots & 0 & 0 & 0 & \cdots & 0 \\
\vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\
0 & \cdots & 1 & 0 & 0 & \cdots & 0 \\
0 & \cdots & 0 & 1 & 0 & \cdots & 0 \\
0 & \cdots & 0 & m_{k+1, k} & 1 & \cdots & 0 \\
\vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\
0 & \cdots & 0 & m_{n,k} & 0 & \cdots & 1
\end{array}\right]^{-1} =\;\; \left[ 
\begin{array}{ccccccc}
1 & \cdots & 0 & 0 & 0 & \cdots & 0 \\
\vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\
0 & \cdots & 1 & 0 & 0 & \cdots & 0 \\
0 & \cdots & 0 & 1 & 0 & \cdots & 0 \\
0 & \cdots & 0 & -m_{k+1, k} & 1 & \cdots & 0 \\
\vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\
0 & \cdots & 0 & -m_{n,k} & 0 & \cdots & 1
\end{array}\right]$$

3. Production
$$\left[ 
\begin{array}{ccccccc}
1 & \cdots & 0 & 0 & 0 & \cdots & 0 \\
\vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\
0 & \cdots & 1 & 0 & 0 & \cdots & 0 \\
0 & \cdots & 0 & 1 & 0 & \cdots & 0 \\
0 & \cdots & 0 & m_{k+1, k} & 1 & \cdots & 0 \\
\vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\
0 & \cdots & 0 & m_{n,k} & 0 & \cdots & 1
\end{array}\right] \cdot \left[ 
\begin{array}{ccccccc}
1 & \cdots & 0 & 0 & 0 & \cdots & 0 \\
\vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\
0 & \cdots & 1 & 0 & 0 & \cdots & 0 \\
0 & \cdots & 0 & 1 & 0 & \cdots & 0 \\
0 & \cdots & 0 & m_{j+1, j} & 1 & \cdots & 0 \\
\vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\
0 & \cdots & 0 & m_{n,j} & 0 & \cdots & 1
\end{array}\right] \\ = \left[ 
\begin{array}{ccccccccccc}
1 & \cdots & 0 & 0 & 0 & \cdots & 0 & 0 & 0 & \cdots & 0 \\
\vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\
0 & \cdots & 1 & 0 & 0 & \cdots & 0 & 0 & 0 & \cdots & 0 \\
0 & \cdots & 0 & 1 & 0 & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\
0 & \cdots & 0 & m_{k+1, k} & 1 & \cdots & 0 & 0 & 0 & \cdots & 0 \\
\vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\
0 & \cdots & 0 & m_{j-1, k} & 0 & \cdots & 1 & 0 & 0 & \cdots & 0 \\
0 & \cdots & 0 & m_{j, k} & 0 & \cdots & 0 & 1 & 0 & \cdots & 0 \\
0 & \cdots & 0 & m_{j+1, k} & 0 & \ddots & 0 & m_{j+1, j} & 1 & \cdots & 0 \\
\vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\
0 & \cdots & 0 & m_{j+1, n} & 0 & \ddots & 0 & m_{j+1, n} & 0 & \cdots & 0 \\
\end{array}\right]$$

# LU factorization
For a general invertible matrix $A = [a_{i,j}]_{n \times n}$, 

$$\begin{align}
M_1A = \left[ 
\begin{array}{ccccc}
1 & & & & \\
-a_{21}/a_{11} & 1 & & &  \\
-a_{31}/a_{11} & 0 & 1 & & \\
\vdots & \vdots & \vdots & \ddots & \\
-a_{n1}/a_{11} & 0 & 0 & \cdots & 1
\end{array}\right] \cdot \left[ 
\begin{array}{ccccc}
a_{11} & a_{12} & a_{13} & \cdots & a_{1n} \\
a_{21} & a_{22} & a_{23} & \cdots & a_{2n} \\
a_{31} & a_{32} & a_{33} & \cdots & a_{3n} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
a_{n1} & a_{n2} & a_{n3} & \cdots & a_{nn} \\
\end{array}\right] = \left[ 
\begin{array}{ccccc}
a_{11} & a_{12} & a_{13} & \cdots & a_{1n} \\
0 & \tilde{a}_{22} & \tilde{a}_{23} & \cdots & \tilde{a}_{2n} \\
0 & \tilde{a}_{32} & \tilde{a}_{33} & \cdots & \tilde{a}_{3n} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0 & \tilde{a}_{n2} & \tilde{a}_{n3} & \cdots & \tilde{a}_{nn} \\
\end{array}\right]
\end{align}$$

Keep multiply these $M$ Lower matrix, finally we can see that

$$M_{n-1}M_{n-2}\cdots M_2M_1A = \left[ 
\begin{array}{cccc}
u_{11} & u_{12} & \cdots & u_{1n} \\
       & u_{22} & \cdots & u_{2n} \\
       &        & \ddots & \vdots \\
       &        &        & u_{nn}
\end{array}\right] = U$$

$i.e.$

$$A = \left(\prod _{i = 1}^{n-1} M_{i}^{-1} \right) U = LU = \left[ 
\begin{array}{cccc}
   1   &        &        &        \\
l_{21} &    1   &        &        \\
\vdots & \vdots & \ddots &        \\
l_{n1} & l_{n2} & \cdots &   1   
\end{array}\right] \cdot \left[ 
\begin{array}{cccc}
u_{11} & u_{12} & \cdots & u_{1n} \\
       & u_{22} & \cdots & u_{2n} \\
       &        & \ddots & \vdots \\
       &        &        & u_{nn}
\end{array}\right] $$

where $L$ is a lower triangular matrix with $1$ on its diagonals, and $U$ is an upper normal triangular matrix.

$Algorithm$: **LU factorization**

    for k = 1, . . . , n − 1         % loop over the first n − 1 columns
      if a_{kk} = 0 then stop
        for i = k + 1, . . . , n     % loop over Row k + 1 through Row n
          a_ {ik} ← a_ {ik}/a_{kk}   % store the multiplier, which gives L, NOT included in the flops
          for j = k + 1, . . . , n   % update other entries of that row
            a_{ij} ← a_{ij} − a_ {ik} a_{kj}
          end
        end
    end

To evaluate the c*omputational expense* of an algorithm we can analyze the number of floating point operation (flops). One operation of addition, subtraction, multiplication, division, etc, is called a **floating point operation**. For LU factorization, we have

$$\begin{align}
\text{flops} &= \sum_{k = 1} ^{n - 1} \sum_{i = k + 1} ^{n} \sum_{j = k + 1} ^{n} 2 \\
&= 2 \sum_{k = 1} ^{n} (n - k)^2 = 2\sum_{l = 1} ^{n - 1} l^2 \\
&= 2\left[n^3 - 1 - \frac{3n(n - 1)} {2} - (n - 1) \right]\\
&\approx \frac{2} {3} n^3
\end{align}$$

It is a big number, when we are dealing with a huge matrix. But with certain special structures we can modify the algorithm to make it more effective. Especially the **banded matrix**, where the nonzero entries of a matrix $A_{n \times n}$ only appear in the main $m$ off diagonals.

$Algorithm$: **LU factorization for banded matrix**

    for k = 1, . . . , n − 1                     % loop over the first n − 1 columns
      if a_{kk} = 0 then `stop`
        for i = k + 1, . . . , min{k + m, n}     % loop over Row k + 1 to the last nonzero entry
          a_ {ik} ← a_ {ik}/a_{kk}               % store the multiplier, which gives L, NOT included in the flops
          for j = k + 1, . . . , min{k + m, n}   % update other entries of that row
            a_{ij} ← a_{ij} − a_ {ik} a_{kj}
          end
        end
    end

$$\begin{align}
\text{flops} &= \sum_{k = 1} ^{n - 1} \sum_{i = k + 1} ^{k + m} \sum_{j = k + 1} ^{k + m} 2 \\
&= 2 \sum_{k = 1} ^{n} m^2 \\
&= 2m^2 (n - 1)
\end{align}$$
***
So how do we use this factorization? Actually we can use it for solving linear equations

$$A\vec{x} = \vec{b} \Leftrightarrow 
\begin{cases}
L\vec{\mathstrut y} = \vec{\mathstrut b} \\
U\vec{\mathstrut x} = \vec{\mathstrut y} 
\end{cases}
$$

>**e.g.** Solve
>
>$$A\vec{x} = \left[ 
\begin{array}{ccc}
2 & 2 & 2 \\
4 & 3 & 2 \\
4 & 6 & 4
\end{array}\right] \cdot \left[ 
\begin{array}{c}
x_1 \\
x_2 \\
x_3
\end{array}\right] = \left[ 
\begin{array}{c}
2 \\
3 \\
2
\end{array}\right]$$
>
>First we factorization matrix $A$ like
>
>$$\begin{align}
A\vec{x} &= LU\vec{x} = \vec{b} \\
&= \left( \begin{bmatrix}
1 & 0 & 0 \\
2 & 1 & 0 \\
2 & -2 & 1
\end{bmatrix} \cdot
\begin{bmatrix}
2 & 2 & 2 \\
0 & -1 & -2 \\
0 & 0 & -4
\end{bmatrix}\right ) \cdot \left[ 
\begin{array}{c}
x_1 \\
x_2 \\
x_3
\end{array}\right] = \left[ 
\begin{array}{c}
2 \\
3 \\
2
\end{array}\right]
\end{align}$$
>
>$i.e.$
>
>$$L\vec{y} = \begin{bmatrix}
1 & 0 & 0 \\
2 & 1 & 0 \\
2 & -2 & 1
\end{bmatrix} \cdot \left[ 
\begin{array}{c}
y_1 \\
y_2 \\
y_3
\end{array}\right] = \left[ 
\begin{array}{c}
2 \\
3 \\
2
\end{array}\right] \\ \Downarrow \\
\vec{y} = \left[ 
\begin{array}{c}
y_1 \\
y_2 \\
y_3
\end{array}\right] = \left[ 
\begin{array}{c}
2 \\
-1 \\
-4
\end{array}\right]$$
>
>Then we solve
>$$U\vec{x} = \begin{bmatrix}
2 & 2 & 2 \\
0 & -1 & -2 \\
0 & 0 & -4
\end{bmatrix} \cdot \left[ 
\begin{array}{c}
x_1 \\
x_2 \\
x_3
\end{array}\right] = \left[ 
\begin{array}{c}
2 \\
-1 \\
-4
\end{array}\right] \\ \Downarrow \\
\vec{x} = \left[ 
\begin{array}{c}
x_1 \\
x_2 \\
x_3
\end{array}\right] = \left[ 
\begin{array}{c}
1 \\
-1 \\
1
\end{array}\right]$$

# Gaussian elimination with pivoting
**Pivot**, the $a_{kk}$. And the elimination without pivoting may fail when meeting situation like the following:

$$A = \begin{bmatrix}
0 & 1 \\ 
1 & 1
\end{bmatrix}$$

First we define the **permutation matrix**, obtained by exchanging any two rows of the identity matrix. An example can show how it works.

$$PMP = \begin{bmatrix}
1 & & & & \\ 
 & & & 1 & \\
 & & 1 & & \\
 & 1 & & & \\
 & & & & 1 
\end{bmatrix}\begin{bmatrix}
1 & & & & \\ 
m_{21} & 1 & & & \\
m_{31} & & 1 & & \\
m_{41} & & & 1 & \\
m_{51} & & & & 1 
\end{bmatrix}\begin{bmatrix}
1 & & & & \\ 
 & & & 1 & \\
 & & 1 & & \\
 & 1 & & & \\
 & & & & 1 
\end{bmatrix} = \begin{bmatrix}
1 & & & & \\ 
m_{41} & 1 & & & \\
m_{31} & & 1 & & \\
m_{21} & & & 1 & \\
m_{51} & & & & 1 
\end{bmatrix}$$

Another property of the permutation matrix is that for any permutation matrix $P$, it holds that $P^2 = I$, $i.e.$, $P
^{−1} = P$.

Now we can get down to the Gaussian elimination with pivoting. Before we eliminate column $k$, we first need to exchange the largest entry among all $a_{ik}, i = k, k+1, \dots, n$ to the entry $a_{kk}$, meaning that change the row with the largest entry with row $k$.

So that in general, this process can be represented as

$$M_{n-1}P_{n-1}M_{n-2}P_{n-2}\cdots M_{2}P_{2}M_{1}P_{1}A = U$$

$Conclusion$

For any nonsingular matrix $A$, there exists a permutation matrix $P = P_{n-1}P_{n-2}\dots P_2P_1$ such that $PA = LU$.

$Proof$

We've already got that $A = P_{1}^{-1}M_{1}^{-1}P_{2}^{-1}M_{2}^{-1}\dots P_{n-2}^{-1}M_{n-2}^{-1}P_{n - 1}^{-1}M_{n-1}^{-1} = P_{1}L_{1}P_{2}L_{2}\dots P_{n-2}L_{n-2}P_{n-1}L_{n-1}$. The $L$ is the inverse of $M$, also a lower triangular matrix. So that now we have

$\begin{align}
A &= P_{1}L_{1}P_{2}L_{2}\dots P_{n-2}L_{n-2}P_{n-1}L_{n-1}\\
P_{n-1}P_{n-2}\dots P_2P_1A &= (P_{n-1}P_{n-2}\dots P_2) L_{1}P_{2}L_{2}\dots P_{n-2}L_{n-2}P_{n-1}L_{n-1} \\
&= (P_{n-1}P_{n-2}\dots P_2 L_{1}P_2P_3\dots P_{n-2}P_{n-1}) (P_{n-1}P_{n-2}\dots P_2)P_2 L_2 \dots P_{n-2}L_{n-2}P_{n-1}L_{n-1} \\
& = \tilde{L}_1 (P_{n-1}P_{n-2}\dots P_3 L_{2}P_3\dots P_{n-2}P_{n-1}) (P_{n-1}P_{n-2}\dots P_3)P_3 L_3 \dots P_{n-2}L_{n-2}P_{n-1}L_{n-1} \\
&\cdots \\
&= \tilde{L}_1\tilde{L}_2\cdots \tilde{L}_{n-1} L_n U
\end{align}$

To summary, 

$$\begin{cases}
\tilde{L}_1 &= P_{n-1}P_{n-2}\dots P_{2}L_{1}P_{2}\dots P_{n-2}P_{n-1} \\
\tilde{L}_2 &= P_{n-1}P_{n-2}\dots P_{3}L_{2}P_{3}\dots P_{n-2}P_{n-1} \\
\vdots &\vdots \\
\tilde{L}_{n-2} &=  P_{n-1}L_{n-2}P_{n-1} \\
\tilde{L}_{n-1} &= L_{n-1} \\
P &= P_{n-1} P_{n-2} \dots P_2 P_1 \\
L &= \tilde{L}_1 \tilde{L}_2 \cdots \tilde{L}_{n-2} \tilde{L}_{n-1}
\end{cases}$$

$Algorithm$: **LU factorization using pivoting**

    for k = 1, . . . , n − 1
      find the row number i^∗, where |a_{i^∗k}| = max_{k≤i≤n} |a_{ik}|
      exchange the two rows A(k,:) ↔ A(i^∗,:)
      p(k) = i^∗
      for i = k + 1, . . . , n
        a_{ik} ← a_{ik}/a{kk}
        for j = k + 1, . . . , n
          a_{ij} ← a_{ij} − a_{ik} a_{kj}
        end
      end
    end

So now to solve a linear system $A\vec{x} = \vec{b}$, we can solve $PA\vec{x} = LU\vec{x} = P\vec{b}$, so actually what's needed is permutation matrix $P$, which can be obtained during the process of pivoting.

# Stability of Gaussian elimination
For matrix

$$A = \begin{bmatrix}
10^{-20} & 1 \\ 
1 & 1
\end{bmatrix}$$

Easily we can find that $L = \begin{bmatrix}
1 & 0 \\ 
10^{20} & 1 
\end{bmatrix}$ and $U = \begin{bmatrix}
10^{-20} & 1 \\ 
0 & 1 - 10^{20}
\end{bmatrix}$. But without pivoting, in computer, we will see that $\tilde{U} = \begin{bmatrix}
10^{-20} & 1 \\ 
0 & - 10^{20}
\end{bmatrix}$. And the truth is that

$$\tilde{L}\tilde{U} =  \begin{bmatrix}
1 & 0 \\ 
10^{20} & 1 
\end{bmatrix} \cdot \begin{bmatrix}
10^{-20} & 1 \\ 
0 & - 10^{20}
\end{bmatrix} = \begin{bmatrix}
10^{-20} & 1 \\ 
1 & 0
\end{bmatrix}$$

So the first process should be 

$$P_1A = \begin{bmatrix}
0 & 1 \\ 
1 & 0
\end{bmatrix}\begin{bmatrix}
10^{-20} & 1 \\ 
1 & 1
\end{bmatrix} = \begin{bmatrix}
1 & 1 \\ 
10^{-20} & 1
\end{bmatrix}$$

Then

$$M_1P_1A = \begin{bmatrix}
1 & 0 \\ 
-10^{-20} & 1
\end{bmatrix}\begin{bmatrix}
0 & 1 \\ 
1 & 0
\end{bmatrix}\begin{bmatrix}
10^{-20} & 1 \\ 
1 & 1
\end{bmatrix} = \begin{bmatrix}
1 & 1 \\ 
0 & 1-10^{-20}
\end{bmatrix} = \begin{bmatrix}
1 & 1 \\ 
0 & 1
\end{bmatrix}$$

Now, $\tilde{U} = \begin{bmatrix}
1 & 1 \\ 
0 & 1
\end{bmatrix}$, $\tilde{L} = \begin{bmatrix}
1 & 0 \\ 
10^{-20} & 1
\end{bmatrix}$, $P = \begin{bmatrix}
0 & 1 \\ 
1 & 0
\end{bmatrix}$, and we can find that

$$\tilde{L}\tilde{U} = \begin{bmatrix}
1 & 1 \\ 
10^{-20} & 1+10^{-20}
\end{bmatrix}$$

which is the true *LU factorization* of the nearby matrix to $PA$. And since been pivoted, generally the lower trianuglar entries in the $L$ matrix in the Gaussian elimination with pivoting are always at most $1$ in absolute value.

$Theorem$ 1

Let $A = LU$ be computed by Gaussian elimination, with or without partial pivoting, on a computer with machine epsilon $\epsilon_{\text{machine}}$ . Denote the computed $L$ and $U$ matrices by $\tilde{L}$ and $\tilde{U}$, respectively. Then we have

$$\tilde{L}\tilde{U} = A + \delta A$$

for a certain $\delta A$ which satisfies $\displaystyle \frac{\left\| \delta A \right\|} {\left\| L \right\| \left\| U \right\|} = O(\epsilon_{\text{machine}})$

This help us to estimate the relative difference between $\tilde{L}\tilde{U}$ and $A$:

$$\frac{\left\| \delta A \right\|} {\left\| A \right\|} = \frac{\left\| L \right\| \left\| U \right\|} {\left\| A \right\|}O(\epsilon_{\text{machine}})$$

$i.e.$

The stability ($\tilde{f}(x) = f(\tilde{x})$) of the Gaussian elimination depends on the norms of $L$ and $U$. After pivoting, all entries in $L$ is equal or less than $1$, so that the norm of $L$ is not large, then the stability will be determined by $\displaystyle \frac{\left\| U \right\|} {\left\| A \right\|}$

Now we define $\rho = \displaystyle \frac{max_{ij} u_{ij}} {max_{ij} a_{ij}}$, the **growth factor** of the LU factorization. Because $\left\| U \right\| = O( \rho \left\| A \right\| )$.

$Theorem$ 2

Let $PA = LU$ be the LU factorization computed by Gaussian elimination **with partial pivoting** on a **computer** with machine epsilon $\epsilon_{\text{machine}}$. Denote the computed $P$, $L$ and $U$ matrices by $\tilde{P}$, $\tilde{L}$ and $\tilde{U}$, respectively. Then we have

$$\tilde{L}\tilde{U} = \tilde{P}A + \delta A$$

for a certain $\delta A$, which now satisfies $\displaystyle \frac{\left\|U \right\|} {\left\|A\right\|} = O(\rho \epsilon_{\text{machine}})$

A typical size for the growth factor $\rho$ for an $n \times n$ matrix is $O(\sqrt{n})$. And the worst case is that $\rho = 2^{n-1}$

# Cholesky Factorization
If an $n$ by $n$ matrix $A$ is symmetric positive definite, $i.e.$, for all $n\text{-dimensional}$ column vector $\vec{x} \neq \vec{0}$, $x^{\mathrm{T}} Ax > 0$, then we can write the matrix $A$ as

$$\begin{align}
A &= \begin{bmatrix}
a_{11} & \vec{v}_1^{\mathrm{T}} \\ 
\vec{v}_1 & K_1
\end{bmatrix} \\
&= \begin{bmatrix}
\sqrt{a_{11}} & \vec{0} \\ 
\frac{\vec{v}_1}{\sqrt{a_{11}}} & I
\end{bmatrix}
\begin{bmatrix}
1 & \vec{0} \\ 
\vec{0} & K_1 - \frac{\vec{v}_1\vec{v}_1^{\mathrm{T}}}{a_{11}}
\end{bmatrix}
\begin{bmatrix}
\sqrt{a_{11}} & \frac{\vec{v}_1^{\mathrm{T}}}{\sqrt{a_{11}}} \\ 
\vec{0} & I
\end{bmatrix} \\
&= L_1 A_1 L_1^{\mathrm{T}}
\end{align}$$

where $L_1$ is lower triangular and $A_1$ is also symmetric positive definite. Continuing the process, we have a series of $A_i$ that $A_i = L_{i+1}A_{i+1}L_{i+1}^{\mathrm{T}}$, where $L_i$ are all lower triangular and $A_i$ are all symmetric positive definite, with $A_{n+1} = 1$. Therefore

$$A = L_1 A_1 L_1^{\mathrm{T}} = L_1 L_2 A_2 L_2^{\mathrm{T}} L_1^{\mathrm{T}} = \cdots = LL^{\mathrm{T}}$$

here
$$\begin{align}L &= L_1 L_2 \cdots L_n \\ 
&= \begin{bmatrix}
\sqrt{a_{11}} & \vec{0} \\ 
\frac{\vec{v}_1}{\sqrt{a_{11}}} & I
\end{bmatrix}
\left[
\begin{array}{cc}
1 & \vec{0}\\
\vec{0} & \begin{array}{cc}
 {\sqrt{a_{22}}} & \vec{0} \\
 \frac{\vec{v}_2} {\sqrt{a_{22}}} & I
\end{array}
\end{array}
\right] \cdots
\begin{bmatrix}
1 & & & \\
 & \ddots & & \\
 & & 1 & \\
 & & & \sqrt{a_{22}}
\end{bmatrix} \\
&= \left[
\begin{array}{cc}
\sqrt{a_{11}} & \vec{0} \\
\frac{\vec{v}_1} {\sqrt{a_{11}}} & \begin{array}{cc}
\sqrt{a_{22}} & \vec{0} \\
\frac{\vec{v}_2} {\sqrt{a_{22}}} & \begin{array}{cc}
\ddots & \vdots \\
\cdots & \sqrt{a_{nn}}
\end{array}
\end{array}
\end{array}
\right]
\end{align}$$

which is still a lower triangular. We call this the **Cholesky factorization** of the matrix $A$.

$Algorithm$: **Cholesky Factorization**

```
for k = 1, . . . , n                         % loop over columns
    for i = k + 1, . . . , n                 % update each row behind row k
        A(i, k + 1 : n) = A(i, k + 1 : n) − a_{ik}/a_{kk} A(k, k + 1 : n)
    end
    A(k : n, k) = A(k : n, k)/√(a_{kk})
end
```

With this algorithm, $L$ is saved as the lower triangular part of the result matrix $A$. And actually the update of the upper part is not necessary, so here's the efficient version

$Algorithm$: **Cholesky Factorization Efficient Version**

```
for k = 1, . . . , n                         % loop over columns
    for i = k + 1, . . . , n                 % update each row behind row k
        A(i, k + 1 : i) = A(i, k + 1 : i) − a_{ik}/a_{kk} (A(k + 1 : i, k))^T
    end
    A(k : n, k) = A(k : n, k)/√(a_{kk})
end
```

At last, we claim that The Cholesky factorization of symmetric positive definite matrices is **always stable**.

deeda4cf64f85b32879a0d2c783ac9774e6c902e
ipynb
Jupyter Notebook
Computational/Intro to Numerical Computing/Note_Chap05_Gaussian elimination, LU factorization, Cholesky factorization.ipynb
XavierOwen/Notes
d262a9103b29ee043aa198b475654aabd7a2818d
MIT
2018-11-27T10:31:08.000Z
2019-01-20T03:11:58.000Z
Computational/Intro to Numerical Computing/Note_Chap05_Gaussian elimination, LU factorization, Cholesky factorization.ipynb
XavierOwen/Notes
d262a9103b29ee043aa198b475654aabd7a2818d
MIT
Computational/Intro to Numerical Computing/Note_Chap05_Gaussian elimination, LU factorization, Cholesky factorization.ipynb
XavierOwen/Notes
d262a9103b29ee043aa198b475654aabd7a2818d
MIT
2020-07-14T19:57:23.000Z
2020-07-14T19:57:23.000Z
Qwen/Qwen-72B
1. YES
2. YES