# CSE 330 Numerical Analysis Lab 

### Lab 8: LU Decomposition

Let a system of equations be,
\begin{equation}
2\boldsymbol{x}_1 - \boldsymbol{x}_{2}+3\boldsymbol{x}_3 = 4
\end{equation} 
\begin{equation}
4\boldsymbol{x}_1 + 2\boldsymbol{x}_{2}+\boldsymbol{x}_3 = 1
\end{equation} 
\begin{equation}
-6\boldsymbol{x}_1 - \boldsymbol{x}_{2}+2\boldsymbol{x}_3 = 2
\end{equation} 

We can write the whole thing in matrix form,

\begin{equation}
\begin{pmatrix}
2 & -1 & 3 \\
4 & 2 & 1 \\
-6 & -1 & 2 \\
\end{pmatrix}
\begin{pmatrix}
\boldsymbol{x}_1 \\
\boldsymbol{x}_2 \\
\boldsymbol{x}_3 
\end{pmatrix}
=
\begin{pmatrix}
4  \\
1 \\
2
\end{pmatrix}
\end{equation}



In the system of equations, there are 3 equations and 3 unknowns. Therefore, it can be solved for $\boldsymbol{x}_1, \boldsymbol{x}_2, \boldsymbol{x}_3$ unless any 2 of the equations are parallel to each other. 

Gaussian elimination is a commonly used tactics to solve a system of equation. In this method, we create an augmented matrix from a set of equations and then eliminate everything below the diagonal terms. First, let's create the augmented matrix.  

\begin{equation}
\begin{pmatrix}
2 & -1 & 3 &  \qquad   4 \\
4 & 2 & 1  &   \qquad   1\\
-6 & -1 &  2 &   \qquad   2 \\
\end{pmatrix}
\end{equation}

Now, if we multiply row 1 with 2 and subtract it from row 2, we get,

\begin{equation}
\begin{pmatrix}
2 & -1 & 3 &  \qquad   4 \\
0 & 4 & -5  &   \qquad   -7\\
-6 & -1 &  2 &   \qquad   2 \\
\end{pmatrix}
\end{equation}

Additionally, multiply row 1 with -3 and subtract it from row 3,

\begin{equation}
\begin{pmatrix}
2 & -1 & 3 &  \qquad   4 \\
0 & 4 & -5  &   \qquad   -7\\
0 & -4 &  11 &   \qquad   14 \\
\end{pmatrix}
\end{equation}

Finally, multiply row 2 with -1 and subtract it from row 3,

\begin{equation}
\begin{pmatrix}
2 & -1 & 3 &  \qquad   4 \\
0 & 4 & -5  &   \qquad   -7\\
0 & 0 &  6 &   \qquad   7 \\
\end{pmatrix}
\end{equation}

Now, the simplified equations are very easy to solve,

\begin{equation}
6\boldsymbol{x}_3 = 7 \quad or, \quad \boldsymbol{x}_3  = 7/6 
\end{equation}

\begin{equation}
4\boldsymbol{x}_2 - 5\boldsymbol{x}_3  = -7 \quad or, \quad 4\boldsymbol{x}_2 - 35/6 = -7 \quad or, \quad \boldsymbol{x}_2 = -7/24
\end{equation}

\begin{equation}
2\boldsymbol{x}_1 - \boldsymbol{x}_{2}+3\boldsymbol{x}_3 = 4 \quad or, \quad 2\boldsymbol{x}_1 + 7/24 + 21/6 = 4 \quad or, \quad \boldsymbol{x}_1 = 5/48
\end{equation} 


We are solving this for a system $Ax=B$, the complexity of this whole process is somewhere around $O(n^3)$. In many matrix problems you do not only have to solve $Ax=B$, but also have to solve $Ax=C$, $Ax=D$, $Ax=E$ etc., where repeating the same process again and again becomes expensive. 

Therefore, a clever approach is to just drop the right side of the equation at first and decompose the matrix into lower part $\boldsymbol{L}$ and upper part $\boldsymbol{U}$. Where, $\boldsymbol{U}$ is the simplified matrix and $\boldsymbol{L}$ has the multipliers. Then the LU matrices can be repeatedly used to solve multiple systems of equations which is much cheaper. The common structure for these matrices are,

\begin{equation}
L =
\begin{pmatrix}
1 & 0 & 0 \\
L_{21} & 1 & 0 \\
L_{31} & L_{32} & 1 \\
\end{pmatrix}
,
U =
\begin{pmatrix}
U_{11} & U_{12} & U_{13} \\
0 & U_{22} & U_{23} \\
0 & 0 & U_{33} \\
\end{pmatrix}
\end{equation}

Now, let's drop the right side of the original system of equations and isolate the left side,
\begin{pmatrix}
2 & -1 & 3 \\
4 & 2 & 1 \\
-6 & -1 & 2 \\
\end{pmatrix}
Let's calculate the LU decomposition for this.






```
from sympy import Matrix
from sympy import init_printing

A = Matrix([[2,-1,3],
              [4,2,1],
              [-6,-1,2]])
L, U, _ = A.LUdecomposition()

init_printing(use_latex='matplotlib')
L
```




    ⎡1   0   0⎤
    ⎢         ⎥
    ⎢2   1   0⎥
    ⎢         ⎥
    ⎣-3  -1  1⎦




```
U
```




    ⎡2  -1  3 ⎤
    ⎢         ⎥
    ⎢0  4   -5⎥
    ⎢         ⎥
    ⎣0  0   6 ⎦



Using the symbolic python library we can directly compute the LU decomposition for any matrix. Our task for today will be to manually calculate the LU decomposition matrices. LU decompoition has different variants which involve pivoting and partial pivoting in order to avoid 0 on the diagonal and to reduce rounding error. For now, we'll avoid those and calculate in the simplest way possible using gaussian elimination.

An important factor to note: LU decomposition is not a unique value. There might be multiple LU decomposition for the same matrix. Therefore, your manually calculated LU decomposition and the values returned by python libraries may not match.


```
import numpy as np

def LUDecomposition(M):
    #Initializing the Upper matrix that has to be simplified
    U = M
    #Initializing the Lower matrix as an identity matrix (all diagonals are 1)
    L = np.identity(U.shape[0])
    ###Use Gaussian elimination to populate both L and U##### 
    print("L and U of Given Matrix b: \n",L,"\n\n",U,"\n\n")

    pivot =  0
    for i in range(1, len(U)):
        print(f"\n#### Iteration = {i}")
        print("#################################################")
        for row in range(i, len(U)):
            print(f"\n## row = {row}")
            print(f"\npivot = {pivot}, U[{pivot}][{pivot}] = {U[pivot][pivot]}")
            multiplier = U[row][pivot] / U[pivot][pivot]
            L[row][pivot] = multiplier
            print(f"\nmultiplier = {multiplier}\n")
            for col in range(len(U)):
                print(f"col = {col}, U[{row}][{col}] = {U[row][col]}, U[{row-1}][{col}] = {U[row-1][col]}")
                print(f"{U[pivot][col]} - {U[pivot][col]} * {multiplier} = {U[pivot][col] * multiplier}\n")
                U[row][col] -= U[pivot][col] * multiplier
            print(U)
        pivot += 1
            
    return L, U

b = np.array([[2,-1,3],
              [4,2,1],
              [-6,-1,2]])


L, U = LUDecomposition(b)
print("\n\nL and U of Matrix b: \n",L,"\n\n",U,"\n\n")

```

    L and U of Given Matrix b: 
     [[1. 0. 0.]
     [0. 1. 0.]
     [0. 0. 1.]] 
    
     [[ 2 -1  3]
     [ 4  2  1]
     [-6 -1  2]] 
    
    
    
    #### Iteration = 1
    #################################################
    
    ## row = 1
    
    pivot = 0, U[0][0] = 2
    
    multiplier = 2.0
    
    col = 0, U[1][0] = 4, U[0][0] = 2
    2 - 2 * 2.0 = 4.0
    
    col = 1, U[1][1] = 2, U[0][1] = -1
    -1 - -1 * 2.0 = -2.0
    
    col = 2, U[1][2] = 1, U[0][2] = 3
    3 - 3 * 2.0 = 6.0
    
    [[ 2 -1  3]
     [ 0  4 -5]
     [-6 -1  2]]
    
    ## row = 2
    
    pivot = 0, U[0][0] = 2
    
    multiplier = -3.0
    
    col = 0, U[2][0] = -6, U[1][0] = 0
    2 - 2 * -3.0 = -6.0
    
    col = 1, U[2][1] = -1, U[1][1] = 4
    -1 - -1 * -3.0 = 3.0
    
    col = 2, U[2][2] = 2, U[1][2] = -5
    3 - 3 * -3.0 = -9.0
    
    [[ 2 -1  3]
     [ 0  4 -5]
     [ 0 -4 11]]
    
    #### Iteration = 2
    #################################################
    
    ## row = 2
    
    pivot = 1, U[1][1] = 4
    
    multiplier = -1.0
    
    col = 0, U[2][0] = 0, U[1][0] = 0
    0 - 0 * -1.0 = -0.0
    
    col = 1, U[2][1] = -4, U[1][1] = 4
    4 - 4 * -1.0 = -4.0
    
    col = 2, U[2][2] = 11, U[1][2] = -5
    -5 - -5 * -1.0 = 5.0
    
    [[ 2 -1  3]
     [ 0  4 -5]
     [ 0  0  6]]
    
    
    L and U of Matrix b: 
     [[ 1.  0.  0.]
     [ 2.  1.  0.]
     [-3. -1.  1.]] 
    
     [[ 2 -1  3]
     [ 0  4 -5]
     [ 0  0  6]] 
    
    


Let's go back to the original matrix form,
\begin{equation}
\begin{pmatrix}
2 & -1 & 3 \\
4 & 2 & 1 \\
-6 & -1 & 2 \\
\end{pmatrix}
\begin{pmatrix}
\boldsymbol{x}_1 \\
\boldsymbol{x}_2 \\
\boldsymbol{x}_3 
\end{pmatrix}
=
\begin{pmatrix}
4  \\
1 \\
2
\end{pmatrix}
\end{equation}

If it is in the form $Ax=B$ then, 
\begin{equation}
A=
\begin{pmatrix}
2 & -1 & 3 \\
4 & 2 & 1 \\
-6 & -1 & 2 \\
\end{pmatrix}
,B=
\begin{pmatrix}
4  \\
1 \\
2
\end{pmatrix}
\end{equation}

We already know that the LU decompositon of A is,
\begin{equation}
L=
\begin{pmatrix}
1 & 0 & 0 \\
2 & 1 & 0 \\
-3 & -1 & 1 \\
\end{pmatrix}
,
U=
\begin{pmatrix}
2 & -1 & 3 \\
0 & 4 & -5 \\
0 & 0 & 6 \\
\end{pmatrix}
\end{equation}

Now, Use $L$, $U$ and $B$ to solve the original system of equations. 


We have to solve y for $Ly=B$
Or, in this case,
\begin{equation}
\begin{pmatrix}
1 & 0 & 0 \\
2 & 1 & 0 \\
-3 & -1 & 1 \\
\end{pmatrix}
\begin{pmatrix}
y_1 \\
y_2 \\
y_3 \\
\end{pmatrix}
=
\begin{pmatrix}
4  \\
1 \\
2
\end{pmatrix}
\end{equation}
It can be solved by forward substitution. It is an iterative process that can be implemented with a nested loop. 
\begin{equation}
\boldsymbol{y}_1 = 4 
\end{equation}

\begin{equation}
2\boldsymbol{y}_1 + \boldsymbol{y}_2  = 1 \quad or, \quad 8 + \boldsymbol{y}_2 = 1 \quad or, \quad \boldsymbol{y}_2 = -7
\end{equation}

\begin{equation}
-3\boldsymbol{y}_1 - \boldsymbol{y}_{2}+\boldsymbol{y}_3 = 2 \quad or, \quad -12 + 7 + \boldsymbol{y}_3 = 2 \quad or, \quad \boldsymbol{y}_3 = 7
\end{equation} 


After solving $y$, in order to calculate $x$ we need to solve, 
$Ux = y$ or,
\begin{equation}
\begin{pmatrix}
2 & -1 & 3 \\
0 & 4 & -5 \\
0 & 0 & 6 \\
\end{pmatrix}
\begin{pmatrix}
x_1 \\
x_2 \\
x_3 \\
\end{pmatrix}
=
\begin{pmatrix}
y_1 \\
y_2 \\
y_3 \\
\end{pmatrix}
		\Longrightarrow
  \begin{pmatrix}
2 & -1 & 3 \\
0 & 4 & -5 \\
0 & 0 & 6 \\
\end{pmatrix}
\begin{pmatrix}
x_1 \\
x_2 \\
x_3 \\
\end{pmatrix}
=
\begin{pmatrix}
4 \\
-7 \\
7 \\
\end{pmatrix}
\end{equation}
Which can be done through backward substitution.

If you find it hard to iterate backwards, you can just flip $U$, $x$ and $y$, then it results into,
\begin{equation}
\begin{pmatrix}
6 & 0 & 0 \\
-5 & 4 & 0 \\
3 & -1 & 2 \\
\end{pmatrix}
\begin{pmatrix}
x_3 \\
x_2 \\
x_1 \\
\end{pmatrix}
=
\begin{pmatrix}
7 \\
-7 \\
4 \\
\end{pmatrix}
\end{equation}

Solving which would be the same as the forward substitution!

The forward substitution part is done for you. Complete the backward substitution part for correct output.


```
import numpy as np
L = np.array([[1,0,0],
              [2,1,0],
              [-3,-1,1]])

U = np.array([[2,-1,3],
              [0,4,-5],
              [0,0,6]])

B = np.array([4,1,2])


#Forward Substitution process
B_L = np.zeros(B.shape[0]) 

for i in range(L.shape[0]):
    summation=0
    for j in range(L.shape[0]):
        if i == j:
            B_L[j] = B[j] - summation
            B_L[j] = B_L[j]/L[i,j]
        break
    else:
        summation = summation + L[i,j]*B_L[j]

# Backward Substitution, your task is to complete this task
#Flip the U and B_L matrices if necessary using np.flip(array, axis) method
U = np.flip(U, 0)
U = np.flip(U, 1)
B_L = np.flip(B)

B_LU = np.zeros(B.shape[0]) 
#Use U and B_L to populate B_LU, just like how L and B was used to populate B_L in forward substitution

#Place your code here (You may need nested loop)

for i in range(L.shape[0]):
    summation=0
    for j in range(L.shape[0]):
        if i == j:
            B_LU[j] = B_L[j] - summation
            B_LU[j] = B_LU[j]/L[i,j]
        break
    else:
        summation = summation + U[i,j]*B_LU[j]

################################################################


final_result = np.flip(B_LU)

print(final_result)

```

    [0. 0. 0.]


Both the forward and the backward substitution method can accomplish the task with $O(n^2)$ complexity!

b39046a7312bd7d1922200cda8c623c3d0bff989
ipynb
Jupyter Notebook
LU Decomposition.ipynb
sheikhmishar/Numerical-Analysis-Python
03a737ba38b372fb52ad773f52cd029f7da2b307
MIT
LU Decomposition.ipynb
sheikhmishar/Numerical-Analysis-Python
03a737ba38b372fb52ad773f52cd029f7da2b307
MIT
LU Decomposition.ipynb
sheikhmishar/Numerical-Analysis-Python
03a737ba38b372fb52ad773f52cd029f7da2b307
MIT
Qwen/Qwen-72B
1. YES
2. YES