4.4. LU Decomposition#

4.4.1. Learning Objectives#

After studying this notebook, completing the activities, and asking questions in class, you should be able to:

  • Explain to a classmate how Gaussian elimination and LU factorization are similar and different

  • In Python, use LU factorization to solve a linear system of equations

  • For LU factorization with partial pivoting, explain the following (same as Gaussian elimination last class):

    • Why is it important? What are the benefits?

    • How does the algorithm work? What are the basic steps?

    • What are the limitations? What cases will partial pivoting not work?

import numpy as np

4.4.2. Motivation: Solving the “Same” System Multiple Times#

In practice, you often want to solve the same linear system multiple times with different right-hand side vectors \(\vec{b}\). For example, imagine resolving the same mass balance problem with different inlet and outlet compositions.

From the previous class, we know how to use Gaussian elimination to solve a linear system. But with Gaussian elimination, we need to redo the entire algorithm separately for each new instance of \(\vec{b}\). Is there an algorithm that lets us reuse some of our work? Spoiler alert: yes, and it is called LU factorization.

4.4.3. LU Factorization#

LU factorization writes a matrix \(\mathbf{A}\) as the product of a lower triangular matrix (that is a matrix with nonzero elements on the diagonal and below) times an upper triangular matrix (that is a matrix with nonzero elements above the diagonal). In other words, \(\mathbf{A} = \mathbf{LU},\) where

\[\begin{split} \mathbf{L} = \begin{pmatrix} * & 0 & \dots\\ * & * & 0 & \dots\\ * & * & * &0 & \dots\\ \vdots & \vdots & \vdots & \ddots\\ * & * & * & \dots & * \end{pmatrix} \end{split}\]

and

\[\begin{split} \mathbf{U} = \begin{pmatrix} * & * & * & \dots & * \\ 0 & * & * & \dots & * \\ 0 & 0 & \ddots & \dots & *\\ \vdots & \vdots & \vdots & \ddots\\ 0 & \dots & &0 & *\end{pmatrix} \end{split}\]

where the \(*\) denote a potentially nonzero matrix element.

The question remains; how we find the entries in the matrices in the LU factorization? To see this, let’s look at a matrix we used in last class:

\[\begin{split} \begin{pmatrix} 3 & 2 & 1\\ -1 & 4 & 5\\ 2& -8 & 10\end{pmatrix} \end{split}\]

we want to find the matrix LU such that

\[\begin{split} \begin{pmatrix} l_{11} & 0 & 0\\ l_{21} & l_{22} & 0\\ l_{31}& l_{32} & l_{33}\end{pmatrix} \begin{pmatrix} u_{11} & u_{12} & u_{13}\\ 0 & u_{22} & u_{23}\\ 0& 0 & u_{33}\end{pmatrix}=\begin{pmatrix} 3 & 2 & 1\\ -1 & 4 & 5\\ 2& -8 & 10\end{pmatrix}. \end{split}\]

Home Activity

Please read through the explanation below. Think of at least one question to ask in class.

It’s not exactly clear right now, but we actually have a few choices at this point. The first choice we have is which diagonal, either \(\mathbf{L}\) or \(\mathbf{U}\), we want to be all ones. We will choose to have all the \(l_{ii} = 1\), that is the \(\mathbf{L}\) matrix has a diagonal of all 1’s. (This means we are choosing Doolittle’s method for LU factorization.) Performing the matrix multiplication we get 9 equations and 9 unknowns for the factorization. The product of the first column of \(\mathbf{U}\) by the first column of \(\mathbf{L}\) gives

\[\begin{split} u_{11} = 3\\ l_{21} u_{11} = -1\\ l_{31} u_{11} = 2.\end{split}\]

Which gives us

\[\begin{split}u_{11} = 3,\\ l_{21} = -\frac{1}{3},\\ l_{31} = \frac{2}{3}.\end{split}\]

Notice that the solutions below the diagaonal, the \(l_{21}\) and \(l_{31}\), are the factors we used in our Gaussian elimination algorithm to remove the first column in the second and third equations. And the \(u_{11}\) value is that found in the same position in our Gaussian-eliminated matrix.

Continuing on, the equations from the second column are

\[\begin{split} u_{12} = 2,\\ u_{12}l_{21} + u_{22} = 4,\\ u_{12} l_{31} + u_{22} l_{32} = -8,\end{split}\]

the solutions to these equations are

\[\begin{split}u_{12} = 2,\\ u_{22} = 4\frac{2}{3}, \\ l_{32} = -2.\end{split}\]

You probably don’t remember, because I didn’t, but the value of \(l_{32}\) is the value of the factor used to eliminate the second column below the diagonal in the Gaussian elimination example from the last notebook. Finally, the last three equations are

\[\begin{split} u_{13} = 1,\\u_{13} l_{21} + u_{23} = 5,\\ u_{13} l_{31} + u_{23}l_{32} + u_{33} = 10.\end{split}\]

The solutions to these equations are

\[\begin{split}u_{13} = 1, \\u_{23} = 5\frac{1}{3},\\ u_{33} = 20.\end{split}\]

These are precisely the values of the third column in or Gaussian elimination example.

In other words our factorization looks like:

\[\begin{split}\underbrace{\begin{pmatrix} 1 & 0 & 0\\ -\frac{1}{3} & 1& 0\\ \frac{2}{3} & -2 & 1\end{pmatrix}}_{\mathbf{L}} \underbrace{\begin{pmatrix} 3 & 2 & 1\\ 0 & 4\frac{2}{3} & 5\frac{1}{3}\\ 0& 0 & 20\end{pmatrix}}_{\mathbf{U}}=\underbrace{\begin{pmatrix} 3 & 2 & 1\\ -1 & 4 & 5\\ 2& -8 & 10\end{pmatrix}}_{\mathbf{A}}.\end{split}\]

Home Activity

Let’s verify the above equation is correct. Create numpy arrays L and U for the example above. Then multiple matrices L and U and store the results in A.

# Add your solution here
# Removed autograder test. You may delete this cell.

Some observations:

  1. The upper triangular matrix is the same as the matrix we received after doing Gaussian elimination

  2. The non-zero and non-diagonal elements of the lower triangular matrix are the factors we used to arrive at our Gaussian matrix.

This suggests that we can reformulate our Gaussian elimination example to give us an LU factorization.

Home Activity

Write down questions you can ask in class about the example above.

Record your question(s) here.

Class Activity

Walk through the LU factorization example together.

4.4.4. Our First Function for LU Factorization#

We will now generalize the LU factorization algorithm into a Python function.

import numpy as np
def LU_factor1(A,LOUD=True):
    '''Compute LU (Doolitte) decomposition of A
    
    Args:
        A: N by N array
    
    Returns:
        L: N by N array
        U: N by N array
    
    Side Effect:
        Does not modify A
    '''
    
    # Check dimensions
    [Nrow, Ncol] = A.shape
    assert Nrow == Ncol
    N = Nrow
    
    # Allocate L as identify matrix
    L = np.eye(N)
    
    # Copy A into U
    U = A.copy()
    
    # Loop over columns
    for c in range(0,N):
    
        # Loop over rows
        for r in range(c+1,N):
        
            # Extract row to modify
            mod_row = U[r,:]
            
            # Calculate factor (just like Gaussian Elimination)
            factor = mod_row[c]/U[c,c]
            
            # Save factor into L
            L[r,c] = factor
            
            # Perform elimination, save into U
            mod_row -= factor*U[c]
            
            #mod_row = U[r-1,:]
            
            if LOUD:
                print("After working on column",c,"and row",r)
                print("L = \n",L)
                print("U = \n",U)
                print("\n")
    return L,U
            
#let's try it on a 3 x 3 to start
A = np.array([(3.0,2,1),(-1,4,5),(2,-8,10)])
print("The original matrix is\n",A)
L,U = LU_factor1(A)
#print("The LU factored A is\n",A) 
    
print("L * U = \n")
print(L,"\n*\n",U,"\n=\n",L.dot(U))
The original matrix is
 [[ 3.  2.  1.]
 [-1.  4.  5.]
 [ 2. -8. 10.]]
After working on column 0 and row 1
L = 
 [[ 1.          0.          0.        ]
 [-0.33333333  1.          0.        ]
 [ 0.          0.          1.        ]]
U = 
 [[ 3.          2.          1.        ]
 [ 0.          4.66666667  5.33333333]
 [ 2.         -8.         10.        ]]


After working on column 0 and row 2
L = 
 [[ 1.          0.          0.        ]
 [-0.33333333  1.          0.        ]
 [ 0.66666667  0.          1.        ]]
U = 
 [[ 3.          2.          1.        ]
 [ 0.          4.66666667  5.33333333]
 [ 0.         -9.33333333  9.33333333]]


After working on column 1 and row 2
L = 
 [[ 1.          0.          0.        ]
 [-0.33333333  1.          0.        ]
 [ 0.66666667 -2.          1.        ]]
U = 
 [[ 3.          2.          1.        ]
 [ 0.          4.66666667  5.33333333]
 [ 0.          0.         20.        ]]


L * U = 

[[ 1.          0.          0.        ]
 [-0.33333333  1.          0.        ]
 [ 0.66666667 -2.          1.        ]] 
*
 [[ 3.          2.          1.        ]
 [ 0.          4.66666667  5.33333333]
 [ 0.          0.         20.        ]] 
=
 [[ 3.  2.  1.]
 [-1.  4.  5.]
 [ 2. -8. 10.]]

We see that our function gives the same L and U as before. At the start of class, we will work through the LU decomposition example on the board together for this example. We will then revist the function above.

Home Activity

Perform LU decomposition for matrix A1, which is defined below. Store the results in variables L1 and U1.

A1 = np.array([[1., 2, 3],[-1, 2, 1], [1, 2, 0]])
print("A1 = \n",A1,"\n")

# Add your solution here
# Removed autograder test. You may delete this cell.

Home Activity

Spend 5 honest minutes reading through the function LU_factor1. Write down at least two questions you have right now.

My questions while inspecting the code below: 1. 2.

Class Activity

Revist your list of questions about LU_factor1 with your neighbor. We’ll then regroup and discuss as a class.

4.4.5. A More Memory Efficient Version#

Our function LU_factor1 makes a copy of A to initialize U. This is nice because it does not overwrite the original A. But imagine if A was a 10 million by 10 million matrix. Making a copy would be expensive and require a lot of memory.

Instead, the function below also performs LU factorization. But instead of returning L and U, it overwrites A with the factorization result. This avoids copying large matrices.

import numpy as np
def LU_factor(A):
    """Factor in place A into L*U=A. 
    The lower triangular parts of A
    are the L matrix.  
    The L has implied ones on the diagonal.
    Args:
        A: N by N array
    Side Effects:
        A is factored in place.
    """
    [Nrow, Ncol] = A.shape
    assert Nrow == Ncol
    N = Nrow
    for column in range(0,N):
        for row in range(column+1,N):
            mod_row = A[row]
            factor = mod_row[column]/A[column,column]
            mod_row = mod_row - factor*A[column,:]
            #put the factor in the correct place in the modified row
            mod_row[column] = factor
            #only take the part of the modified row we need
            mod_row = mod_row[column:N]
            A[row,column:N] = mod_row
    return

#let's try it on a 3 x 3 to start
A = np.array([(3.0,2,1),(-1,4,5),(2,-8,10)])
print("The original matrix is\n",A)
LU_factor(A)
print("The LU factored A is\n",A)
The original matrix is
 [[ 3.  2.  1.]
 [-1.  4.  5.]
 [ 2. -8. 10.]]
The LU factored A is
 [[ 3.          2.          1.        ]
 [-0.33333333  4.66666667  5.33333333]
 [ 0.66666667 -2.         20.        ]]

Notice how our matrix A has changed:

print(A)
[[ 3.          2.          1.        ]
 [-0.33333333  4.66666667  5.33333333]
 [ 0.66666667 -2.         20.        ]]

Home Activity

Use LU_factor for the example below. Your matrix A1 should change.

A1 = np.array([[1., 2, 3],[-1, 2, 1], [1, 2, 0]])
print("A1 (original) = \n",A1,"\n")

# Add your solution here

print("A1 (factored) = \n",A1,"\n")
# Removed autograder test. You may delete this cell.

4.4.6. Solving Linear Systems#

This algorithm gives us the same thing our onerous, by-hand procedure did. The question we have not answered is how we can take an LU factored matrix and get the solution to the system \(\mathbf{A}\vec{x} = \vec{b}\).

Let’s start by substituting \(\mathbf{L U} = \mathbf{A}\) into our original equation:

\[\mathbf{A}\vec{x} = \mathbf{LU}\vec{x} = \vec{b}\]

Now let’s create a new intermediate variable \(\vec{y} = \mathbf{U} \vec{x}\). If we substitute \(\vec{y}\) into the \(\mathbf{LU}\vec{x} = \vec{b}\), we get:

\[\mathbf{L}\vec{y} = \vec{b}\]

Home Activity

Compute \(\vec{y}\) for the example below by hand or with Python. Store the results in a numpy array y.

\[\begin{split}\mathbf{L} = \begin{pmatrix} 1 & 0 & 0\\ -\frac{1}{3} & 1& 0\\ \frac{2}{3} & -2 & 1\end{pmatrix}, \quad \vec{b} = \begin{pmatrix}6 \\ 8 \\ 4 \end{pmatrix} \end{split}\]
# Add your solution here
# Removed autograder test. You may delete this cell.

You may have noticed in the example, we can easily compute \(\vec{y}\) because \(\mathbf{L}\) is a lower triangular matrix. We start at the top to calculate the first element of \(\vec{y}\). We then move to the next row of \(\mathbf{L}\), substitute and compute the next element of \(\vec{y}\). We call this procedure forward substitution because we start at the top.

We now have \(\vec{y}\). We are half way there! We just need to solve the linear system \(\mathbf{U}\vec{x} = \vec{y}\).

Home Activity

Compute \(\vec{x}\) for the example below either by hand or with Python. Store your answer in the numpy array x.

\[\begin{split}\mathbf{U} = \begin{pmatrix} 3 & 2 & 1\\ 0 & 4\frac{2}{3} & 5\frac{1}{3}\\ 0& 0 & 20\end{pmatrix}, \quad \vec{y} = \begin{pmatrix}6 \\ 10 \\ 20 \end{pmatrix} \end{split}\]
# Add your solution here
# Removed autograder test. You may delete this cell.

Notice again this calculation is easy because \(\mathbf{U}\) is upper triangular. We start at the bottom and calculate the last element of \(\vec{x}\). We then move up a row, substitute, and calculate the second to last element of \(\vec{x}\). This is called backward substitution and is the same procedure we used with Gaussian elimination last class.

4.4.7. Our First Function for Forward and Backward Solves#

Let’s return to our first factorization function LU_factor which returned a distinct L and U. We will create Python code to generalize the forward and backward substitution steps to solve the linear system.

def LU_solve1(L,U,b,LOUD=True):
    '''Solve L*U*x = b for x where L and U
    are lower and upper triangular matrices and
    b is RHS vector.
    
    Args:
        L: N by N lower triangular matrix with 1 on diagonal
        U: N by N upper triangular matrix
        b: N by 1 vector
        
    Returns
        x: N by 1 solution to linear system
    '''
    
    # Check dimensions
    [Nrow_L, Ncol_L] = L.shape
    assert Nrow_L == Ncol_L
    [Nrow_U, Ncol_U] = U.shape
    assert Nrow_U == Ncol_U
    assert Nrow_U == Nrow_L
    N = Nrow_U
    assert N == b.size
    
    # forward solve
    y = np.zeros(N)
    if LOUD:
        print("\nForward solve...")
    for r in range(0,N): # first for loop
        RHS = b[r]
        
        # loop over some of the columns to evaluated summation
        for c in range(0,r):
            RHS -= y[c]*L[r,c]
        
        if LOUD:
            print("y[",r,"]=",RHS)
        
        y[r] = RHS
    
    # backward solve
    x = np.zeros(N)
    if LOUD:
        print("\nBackward solve...")
    for r in range(N-1,-1,-1): # second for loop
        RHS = y[r]
        
        # loop over some of the columns to evaluated summation
        for c in range(r+1,N):
            RHS -= x[c]*U[r,c]
        
        if LOUD:
            print("x[",r,"]=",RHS/U[r,r])
        x[r] = RHS/U[r,r]
    
    return x

Let’s take a minute to look at the first and second for loops in LU_solve1.

N = 5

# loop 1
print("Loop 1...")
for i in range(0,N):
    print(i)
    
# loop 2
print("\nLoop 2...")
for i in range(N-1,-1,-1):
    print(i)
Loop 1...
0
1
2
3
4

Loop 2...
4
3
2
1
0

Home Activity

Why does the first for loop iterate over range(0,N) and the second for loop iterate over range(N-1,-1,-1) in LU_solve1? Write a sentence below.

Answer:

Class Activity

Discuss the code above with a partner. We’ll then regroup as a class.

And let’s test with the example above.

A = np.array([(3.0,2,1),(-1,4,5),(2,-8,10)])
b = np.array([6,8,4])
print("The original matrix is\n",A)

L,U = LU_factor1(A,False)
print("L = \n",L)
print("\nU= \n",U)

x = LU_solve1(L,U,b)
The original matrix is
 [[ 3.  2.  1.]
 [-1.  4.  5.]
 [ 2. -8. 10.]]
L = 
 [[ 1.          0.          0.        ]
 [-0.33333333  1.          0.        ]
 [ 0.66666667 -2.          1.        ]]

U= 
 [[ 3.          2.          1.        ]
 [ 0.          4.66666667  5.33333333]
 [ 0.          0.         20.        ]]

Forward solve...
y[ 0 ]= 6
y[ 1 ]= 10.0
y[ 2 ]= 20.0

Backward solve...
x[ 2 ]= 1.0
x[ 1 ]= 1.0
x[ 0 ]= 1.0

And now verify the solution.

print("A*x =",A.dot(x))
print("b =",b)
A*x = [6. 8. 4.]
b = [6 8 4]

4.4.8. Another Implementation (Optional)#

Now let’s consider the case where L and U are stored in the same matrix. This complements our more memory efficient function LU_factor.

def LU_solve(A,b):
    """Take a LU factorized matrix and 
    solve it for RHS b
    Args:
        A: N by N array that has been LU factored with
        assumed 1's on the diagonal of the L matrix
        b: N by 1 array of righthand side
    Returns:
        x: N by 1 array of solutions
    """
    [Nrow, Ncol] = A.shape
    assert Nrow == Ncol
    N = Nrow
    
    x = np.zeros(N)
    #temporary vector for L^-1 b
    y = np.zeros(N)
    #forward solve
    for row in range(N):
        RHS = b[row]
        for column in range(0,row):
            RHS -= y[column]*A[row,column]
        y[row] = RHS
    #back solve
    for row in range(N-1,-1,-1):
        RHS = y[row]
        for column in range(row+1,N):
            RHS -= x[column]*A[row,column]
        x[row] = RHS/A[row,row]
    return x

Let’s test on the same system as above.

#let's try it on a 3 x 3 to start
A = np.array([(3.0,2,1),(-1,4,5),(2,-8,10)])
Aorig = A.copy()
print("Originally A =\n",Aorig)
LU_factor(A)
print("\nAfter factorization A =\n",A)

b = np.array([6,8,4])
x = LU_solve(A,b)
print("\nThe solution to the system is",x)
Originally A =
 [[ 3.  2.  1.]
 [-1.  4.  5.]
 [ 2. -8. 10.]]

After factorization A =
 [[ 3.          2.          1.        ]
 [-0.33333333  4.66666667  5.33333333]
 [ 0.66666667 -2.         20.        ]]

The solution to the system is [1. 1. 1.]

Let’s resolve for \(\vec{x} = [2, 12, 2019]^T\)

solution2 = np.array([2,12,2019.])
b = np.dot(Aorig,solution2)
x = LU_solve(A,b)
print("The solution to the system is",x)
The solution to the system is [2.000e+00 1.200e+01 2.019e+03]

4.4.9. LU with Pivoting#

Pivoting with LU is needed in the same cases as it is in Gaussian elimination. However, when we switch equations around, we need to make sure that we keep track of where we’ve switched rows, so we can do the same to the right-hand side when we do our forward and back solves. This complicates things considerably, but only because we have to add that bookkeeping to our algorithm.

First, we’ll need to be able to swap rows in the matrix:

def swap_rows(A, a, b):
    """Rows to rows in a matrix, switch row a with row b
    args: 
    A: matrix to perform row swaps on
    a: row index of matrix
    b: row index of matrix
    
    returns: 
    nothing
    
    side effects:
    changes A to have rows a and b swapped
    """
    assert (a>=0) and (b>=0)
    N = A.shape[0] #number of rows
    assert (a<N) and (b<N) #less than because 0-based indexing
    temp = A[a,:].copy()
    A[a,:] = A[b,:].copy()
    A[b,:] = temp.copy()

Next, we’ll change our algorithm to handle pivoting and keep track of the pivots:

def LU_factor_pivot(A,LOUD=True):
    """Factor in place A in L*U=A. 
    The lower triangular parts of A
    are the L matrix.  
    The L has implied ones on the diagonal.

    Args:
        A: N by N array
    Returns:
        a vector holding the order of the rows, 
        relative to the original order
    Side Effects:
        A is factored in place.
    """
    
    # Check dimensions
    [Nrow, Ncol] = A.shape
    assert Nrow == Ncol
    N = Nrow
    
    # Create scale factors 
    s = np.zeros(N)
    count = 0
    row_order = np.arange(N)
    for row in A:
        # Save largest absolute element in each row
        s[count] = np.max(np.fabs(row))
        count += 1
    if LOUD:
        print("s =",s)
    if LOUD:
        print("Original Matrix is\n",A)
    
    # loop over columns
    for column in range(0,N):
        #swap rows if needed
        largest_pos = np.argmax(np.fabs(A[column:N,column]/
                                        s[column])) + column
        if (largest_pos != column):
            if (LOUD):
                print("Swapping row",column,"with row",
                      largest_pos)
                print("Pre swap\n",A)
            swap_rows(A,column,largest_pos)
            #keep track of changes to RHS
            tmp = row_order[column]
            row_order[column] = row_order[largest_pos]
            row_order[largest_pos] = tmp
            #re-order s
            tmp = s[column]
            s[column] = s[largest_pos]
            s[largest_pos] = tmp
            if (LOUD):
                print("A =\n",A)
        for row in range(column+1,N):
            mod_row = A[row]
            factor = mod_row[column]/A[column,column]
            mod_row = mod_row - factor*A[column,:]
            #put the factor in the correct place in the modified row
            mod_row[column] = factor
            #only take the part of the modified row we need
            mod_row = mod_row[column:N]
            A[row,column:N] = mod_row
    return row_order

Now let’s verify with a 4x4 system

#let's try it on a 4 x 4 to start
A = np.array([(3.0,2,1,-2),(-1,4,5,4),(2,-8,10,3),(-2,-8,10,0.1)])
print("A =\n",A,"\n")
answer = np.ones(4)
print("True solution. x =",answer,"\n")
b = np.dot(A,answer)
print("Thus b =",b,"\n")

# Apply LU with pivoting
row_order = LU_factor_pivot(A,False)

# Print A and the row order
print("After factorization A =\n",A,"\n")
print("The new row order is",row_order)
A =
 [[ 3.   2.   1.  -2. ]
 [-1.   4.   5.   4. ]
 [ 2.  -8.  10.   3. ]
 [-2.  -8.  10.   0.1]] 

True solution. x = [1. 1. 1. 1.] 

Thus b = [ 4.  12.   7.   0.1] 

After factorization A =
 [[ 3.          2.          1.         -2.        ]
 [ 0.66666667 -9.33333333  9.33333333  4.33333333]
 [-0.33333333 -0.5        10.          5.5       ]
 [-0.66666667  0.71428571  0.4        -6.52857143]] 

The new row order is [0 2 1 3]

We need to change the LU_solve function to take advantage of the swapped rows.

def LU_solve_pivot(A,b,row_order):
    """Take a LU factorized matrix and solve it for RHS b

    Args:
        A: N by N array that has been LU factored with
        assumed 1's on the diagonal of the L matrix
        b: N by 1 array of righthand side
        row_order:  list giving the re-ordered equations
        from the the LU factorization with pivoting
    Returns:
        x: N by 1 array of solutions
        
    Side Effect:
        Reorders b using row_order
    """
    [Nrow, Ncol] = A.shape
    assert Nrow == Ncol
    assert b.size == Ncol
    assert row_order.max() == Ncol-1
    N = Nrow
    
    #reorder the equations
    tmp = b.copy()
    for row in range(N):
        b[row_order[row]] = tmp[row]
    #b now has the correct order
    
    # We can either reuse our LU_solve code from above now that b has
    # been reorded or copy and paste
    return LU_solve(A,b)
    '''
    x = np.zeros(N)
    #temporary vector for L^-1 b
    y = np.zeros(N)
    #forward solve
    for row in range(N):
        RHS = b[row]
        for column in range(0,row):
            RHS -= y[column]*A[row,column]
        y[row] = RHS
    #back solve
    for row in range(N-1,-1,-1):
        RHS = y[row]
        for column in range(row+1,N):
            RHS -= x[column]*A[row,column]
        x[row] = RHS/A[row,row]
    return x
    '''


print("A =\n",A,"\n\nRow order =",row_order,"\n")
x = LU_solve_pivot(A,b,row_order)
print("The solution to the system is",x)
A =
 [[ 3.          2.          1.         -2.        ]
 [ 0.66666667 -9.33333333  9.33333333  4.33333333]
 [-0.33333333 -0.5        10.          5.5       ]
 [-0.66666667  0.71428571  0.4        -6.52857143]] 

Row order = [0 2 1 3] 

The solution to the system is [1. 1. 1. 1.]

Class Activity

Verify that our LU decomposition with pivoting code works for other test cases.

epsilon = 1e-14
A_eps = np.array([(epsilon,-1,1),(-1,2,-1),(2,-1,0)])
print(A_eps)
b_eps = np.array([0,0,1.0])
print(b_eps)

Asave = A_eps.copy()
bsave = b_eps.copy()
[[ 1.e-14 -1.e+00  1.e+00]
 [-1.e+00  2.e+00 -1.e+00]
 [ 2.e+00 -1.e+00  0.e+00]]
[0. 0. 1.]
row_order = LU_factor_pivot(A_eps,False)
x_eps = LU_solve_pivot(A_eps,b_eps,row_order)

print("The solution is",x_eps)
print("The residual is",bsave - np.dot(Asave,x_eps))
The solution is [1. 1. 1.]
The residual is [0.00000000e+00 1.11022302e-16 0.00000000e+00]

4.4.10. Take Home Messages#

LU factorization is very similar to Gaussian elimination:

  • \(\mathbf{L}\) matrix saves the scaling factors from Gaussian elimination

  • \(\mathbf{U}\) matrix is the same as the reduced matrix

  • \(\mathbf{L}\) and \(\mathbf{U}\) save all of the steps from matrix reduction, making LU factorization efficient to resolve \(\mathbf{A} \vec{x} = \vec{b}\) for many different \(\vec{b}\)

  • Partial pivoting works the same way for both LU factorization and Gaussian elimination

Solving a linear system with LU factorization is easy:

  • First, perform forward substitution, which exploits that \(\mathbf{L}\) is lower triangular

  • Second, perform backward substitution, which exploits that \(\mathbf{U}\) is upper triangular