Notebook

Topic 17. LU Factorization¶

Author: Spencer Ammermon spencer.ammermon@gmail.com

Engineering Example: Daniel Free, daniel.free816@gmail.com

Introduction¶

Many situations arise where an engineer must solve a linear equation of the form $\textbf{A}\vec{x}=\vec{b}$ . In some cases, matrix $\textbf{A}$ is in a special form (such as a Toeplitz, Vandermonde, or Hankel matrix) and $\vec{x}$ can easily be solved for numerically. When these special forms do not apply to matrix $\textbf{A}$ , one possibility is to factorize $\textbf{A} = \textbf{LU}$ where $\textbf{L}$ is a lower triangular matrix with ones on the diagonal and $\textbf{U}$ is an upper diagonal matrix. The solution to $\textbf{A}\vec{x}=\vec{b}$ can then be found without explicitly computing $\textbf{A}^{-1}\vec{b}$ .

LU Factorization in Detail¶

When finding the factorization of an $m x m$ matrix $\textbf{A}$ , we also compute a permutation matrix $\textbf{P}$ that represents $\textit{pivoting}$ in the factorization. We pivot in the factorization so that we always divide by large numbers to ensure numerical stability. The factorization then becomes $\textbf{PA} = \textbf{LU}$ .

Computing P, L, and U¶

To compute matrices $\textbf{P}$ , $\textbf{L}$ , and $\textbf{U}$ , begin with the original matrix $\textbf{A}$ . The goal is to put matrix $\textbf{A}$ in triangular form while maintaining numerical stability by pivoting to reduce the row with the largest element. Consider Dr. Beard's example matrix:

$\begin{equation} \textbf{A} = \begin{bmatrix} 1 & -2 & 3\\ -4 & 5 & -6\\ 7 & -8 & 9 \end{bmatrix} \end{equation}$ Where we need to swap row 1 and 3 so that the largest element in column 1 is in the 1st position of column 1. To do this, we can compute

$\textbf{P}_{1,3}\textbf{A}$ where

$\textbf{P}_{1,3}$ is found by swapping

$\textit{columns}$ 1 and 3 of an identity matrix.

A function that computes the permutation matrix is as follows:

In [32]:

import numpy as np

A = np.array([[1, -2, 3],
              [-4, 5, -6],
              [7, -8, 9]])

def permutationMatrix(n,col,Q): # n is number of rows and cols, Q is a matrix
    # look down first column, find the largest value in magnitude
    # create a permutation matrix with ones and zeros 
    P = np.identity(n)
    max_index_in_col = np.argmax((np.abs(Q[col:,col:][:,0])))
    # swap the col column with the max column
    P[:,[col, max_index_in_col+col]] = P[:,[max_index_in_col+col, col]]
    return P

P_13 = permutationMatrix(3,0,A)

print('First Permutation Matrix:\n',P_13)

First Permutation Matrix:
 [[0. 0. 1.]
 [0. 1. 0.]
 [1. 0. 0.]]

$\begin{equation} \textbf{P}_{1,3}\textbf{A}= \begin{bmatrix} 0 & 0 & 1\\ 0 & 1 & 0\\ 1 & 0 & 0 \end{bmatrix} \begin{bmatrix} 1 & -2 & 3\\ -4 & 5 & -6\\ 7 & -8 & 9 \end{bmatrix} = \begin{bmatrix} 7 & -8 & 9\\ -4 & 5 & -6\\ 1 & -2 & 3 \end{bmatrix} \end{equation}$

The next step is to zero out rows 2 and 3 in column 1 by dividing by the value in $\textbf{A}[1,1]$ . This is done by finding an elementary matrix that is an identity matrix with a first column that will zero out the first column of $\textbf{P}_{1,3}\textbf{A}$ when multiplied. To find $\textbf{E}_1$ , start with an identity matrix and then copy the first column of $\textbf{P}_{1,3}\textbf{A}$ to the first column of the identity matrix, and then multiply column 1 by $\dfrac{-1}{\textbf{P}_{1,3}\textbf{A}[1,1]}$ . This leads to:

$\begin{equation} \textbf{E}_1\textbf{P}_{1,3}\textbf{A} = \begin{bmatrix} 1 & 0 & 0\\ \frac{4}{7} & 1 & 0\\ \frac{-1}{7} & 0 & 1 \end{bmatrix} \begin{bmatrix} 0 & 0 & 1\\ 0 & 1 & 0\\ 1 & 0 & 0 \end{bmatrix} \begin{bmatrix} 1 & -2 & 3\\ -4 & 5 & -6\\ 7 & -8 & 9 \end{bmatrix} = \begin{bmatrix} 7 & -8 & 9\\ 0 & 0.4286 & -0.8571\\ 0 & -0.8571 & 1.7143 \end{bmatrix} \end{equation}$

The elementary matrices can be found with the function "zeroOutColumn(n,col,Q)" found below:

In [21]:

def zeroOutColumn(n,col,Q):
    E = np.identity(n)          # start out with an identity matrix
    divisor = -1*Q[col][col]    # divisor is the first element in the matrix
    size = len(Q[col:,col:])    
    for i in range(size-1):
        # E[col:,col:] is a square submatrix of E
        # Q[col:,col:] is a square submatrix of the input matrix Q
        # Go through the first column and divide all elements after the first 
        #  by the divisor
        E[col:,col:][i+1][0] = Q[col:,col:][i+1][0]/divisor
    return E

P_13 = permutationMatrix(3,0,A)
E_1 = zeroOutColumn(3,0,P_13 @ A)
res_1 = E_1 @ P_13 @ A
print('E_1:\n',np.round(E_1,4))
print('E_1@P_13@A:\n',np.round(res_1,4))

E_1:
 [[ 1.      0.      0.    ]
 [ 0.5714  1.      0.    ]
 [-0.1429  0.      1.    ]]
E_1@P_13@A:
 [[ 7.     -8.      9.    ]
 [-0.      0.4286 -0.8571]
 [ 0.     -0.8571  1.7143]]

Utilizing the above code, I chose to solve for the remaining permutation and elementary matrices by looking at submatrices of $\textbf{A}$ . Now that the first column contains a pivot, we can focus on the second column by looking at the submatrix starting at $\textbf{E}_1\textbf{P}_{1,3}\textbf{A}[2,2]$ :

$\begin{equation} \textbf{E}_1\textbf{P}_{1,3}\textbf{A}_{[2,2]}= \begin{bmatrix} 0.4286 & -0.8571\\ -0.8571 & 1.7143 \end{bmatrix} \end{equation}$ and then computing the permutation and elementary matrix in the same manner as

$\textbf{P}_{1,3}$ and

$\textbf{E}_1$ . The remaining permutation and elementary matrices will remain

$m x m$ , not the size of the submatrix. The result should be:

$\begin{equation} \textbf{E}_2\textbf{P}_{2,3}\textbf{E}_1\textbf{P}_{1,3}\textbf{A} = \begin{bmatrix} 7 & -8 & 9\\ 0 & -0.8571 & 1.7143\\ 0 & 0 & 0 \end{bmatrix} = \textbf{U} \end{equation}$

$\textbf{L}$ can be computed with an intermediate step by first computing:

$\begin{equation} \textbf{V} = \textbf{P}_{1,3}\textbf{E}^{-1}_1\textbf{P}_{2,3}\textbf{E}^{-1}_{2,3} \end{equation}$

$\begin{equation} \textbf{P}_{2,3}\textbf{P}_{1,3}\textbf{V} = \begin{bmatrix} 1 & 0 & 0\\ 0.1429 & 1 & 0\\ -0.5714 & -0.5 & 1 \end{bmatrix} = \textbf{L} \end{equation}$

Finally, $\textbf{P}$ can be computed as:

$\begin{equation} \textbf{P}_{2,3}\textbf{P}_{1,3} = \begin{bmatrix} 0 & 0 & 1\\ 1 & 0 & 0\\ 0 & 1 & 0 \end{bmatrix} = \textbf{P} \end{equation}$

Homebrew LU Factorization Code¶

We now have all the components necessary to solve an equation of the form $\textbf{Ax} = \textbf{b}$ using LU factorization. To put all the code together and create our own LU Factorization:

In [14]:

"""
Created on Tue Oct 27 10:28:34 2020

@author: sma

LU Factorization 
"""
import numpy as np

A = np.array([[1, -2, 3],
              [-4, 5, -6],
              [7, -8, 9]])

def permutationMatrix(n,col,Q): # n is number of rows and cols, Q is a matrix
    # look down first column, find the largest value in magnitude
    # create a permutation matrix with ones and zeros 
    P = np.identity(n)
    max_index_in_col = np.argmax((np.abs(Q[col:,col:][:,0])))
    # flip the col column with the max column
    P[:,[col, max_index_in_col+col]] = P[:,[max_index_in_col+col, col]]
    return P

def zeroOutColumn(n,col,Q):
    E = np.identity(n)          # start out with an identity matrix
    divisor = -1*Q[col][col]    # divisor is the first element in the matrix
    size = len(Q[col:,col:])    
    for i in range(size-1):
        # E[col:,col:] is a square submatrix of E
        # Q[col:,col:] is a square submatrix of the input matrix Q
        # Go through the first column and divide all elements after the first 
        #  by the divisor
        E[col:,col:][i+1][0] = Q[col:,col:][i+1][0]/divisor
    return E


def myLUFactorization(A):
    n = len(A) # all matrices for LU must be n x n
    P = []
    E = []
    U = np.zeros((n,n))
    
    for i in range(n-1):
        P.append(permutationMatrix(n,i,A))      # Save the permutation matrices
        E.append(zeroOutColumn(n,i,P[i] @ A))   # Save the zeroing matrices 
        U = E[i] @ P[i] @ A 
        A = U
      
    V_prev = P[0]@np.linalg.inv(E[0])
    V = np.zeros((n,n))
    
    for j in range(len(P)-1): 
        V = V_prev @ P[j+1] @ np.linalg.inv(E[j+1]) 
        V_prev = V
    
    permutation_prev = P[0]
    for k in range(len(P)-1): #counting down
        permutation = P[k+1] @ permutation_prev 
        permutation_prev = permutation

    L = permutation @ V

    return permutation, L, U #P,L,U

myP, myL, myU = myLUFactorization(A)
print('myP:\n',np.round(myP,4))
print('myL:\n',np.round(myL,4))
print('myU:\n',np.round(myU,4))

myP:
 [[0. 0. 1.]
 [1. 0. 0.]
 [0. 1. 0.]]
myL:
 [[ 1.      0.      0.    ]
 [ 0.1429  1.      0.    ]
 [-0.5714 -0.5     1.    ]]
myU:
 [[ 7.     -8.      9.    ]
 [ 0.     -0.8571  1.7143]
 [-0.      0.      0.    ]]

Simple Numerical Examples¶

We should first check that $\textbf{PA} = \textbf{LU}$

In [27]:

print(myP@A,)
print('\n')
print(myL@myU)

[[ 7. -8.  9.]
 [ 1. -2.  3.]
 [-4.  5. -6.]]


[[ 7. -8.  9.]
 [ 1. -2.  3.]
 [-4.  5. -6.]]

Success! We can now solve an equation in the form $\textbf{A}\vec{x}=\vec{b}$ . Suppose

$\begin{equation} \vec{b} = \begin{bmatrix} 10\\ 11\\ 12 \end{bmatrix} \end{equation}$

We then write

$\begin{equation} \textbf{A}\vec{x}=\vec{b} \end{equation}$ as

$\begin{equation} \textbf{LU}\vec{x} = \textbf{P}\vec{b} \end{equation}$ which we can then let

$\textbf{U}\vec{x} = \vec{y}$ , which leads to the system

$\begin{equation} \textbf{L}\vec{y} = \textbf{P}\vec{b} = \vec{c} \end{equation}$

Let's start with $\textbf{L}\vec{y} = \textbf{P}\vec{b}$ by solving for $\vec{y}$ with

$\begin{equation} \vec{y}=\textbf{L}^{-1}\textbf{P}\vec{b} \end{equation}$ Which is easy to solve numerically because

$\textbf{L}$ is lower triangular

In [40]:

import scipy.linalg
b = np.array([[1],
             [2],
             [3]])
y = np.linalg.pinv(myL) @ myP @ b 
print(y)

[[3.        ]
 [0.57142857]
 [4.        ]]

Now we can solve $\textbf{U}\vec{x}=\vec{y}$ for $\vec{x}$ numerically by using

$\begin{equation} \vec{x}=\textbf{U}^{-1}\vec{y} \end{equation}$

In [46]:

x = np.linalg.pinv(myU) @ y
print('x =\n',x)

x =
 [[-0.05555556]
 [-0.11111111]
 [ 0.27777778]]

An Engineering Application¶

Suppose you are sending data, in this case an image, but it has been encrypted or scrambled in some way. Through other techniques, you have determined the matrix that can help unscramble the image, given in the form $A\mathbf{x}=\mathbf{b}$ . However, this matrix applies only to a 50x50 grid of pixels at a time (i.e. $\mathbf{b}$ is the 50x50 matrix of pixel values reshaped into a single vector of length 2500). This means that the matrix $A$ is 2500x2500 and an explicit inverse could be costly. Instead, lets do the LU factorization and use forward and backward substitution to resolve the image. The code below loads in the mixed data and the determined filter, computes the LU factorization, then resolves each 50x50 grid of pixels for each color. Finally, the image is converted to the right data type and displayed.

In [1]:

# Import libraries necessary to run the code
import matplotlib.pyplot as plt
import numpy as np
import scipy.linalg as alg
import scipy.io as sio

# Plot the images inline
%matplotlib inline

# Load the data. Comes with an A matrix and a matrix mx with size (300,150,3) representing rgb values. This is the scrambled image
data = sio.loadmat('mixed_img_andFilter.mat')
A = data['A']
mx = data['mx']

# Plot the scrambled image
fig1 = plt.imshow(mx)
plt.show()

# Compute the LU factorization. the scipy.linalg library computes it as A = P*L*U as opposed to Matlab's A = P'*L*U
P,L,U = alg.lu(A)
# Transpose the permutation matrix P for easier readability later
Pt = np.transpose(P)

# Loop through the scambled image doing 50x50 grids at a time 
for i in range(3):
    for j in range(0,251,50):
        for k in range(0,101,50):
            # Extract the current data
            temp = mx[j:j+50,k:k+50,i]
            tt = np.transpose(temp)
            # Reshape the matrix into a vector
            b = np.reshape(tt,2500)
            Pb = np.dot(Pt,b)
            # Solve the first triangular system (Ly = Pb)
            y = alg.solve_triangular(L,Pb,lower=True)
            # Solve the second triangular system (Ux = y)
            x = alg.solve_triangular(U,y,lower = False)
            # Reshape the output to be a 50x50 matrix again
            mx[j:j+50,k:k+50,i] = np.transpose(np.reshape(x,(50,50)))

# Convert the matrix back into unsigned integers for image plotting
final = np.uintc(mx)

# Show the final image. Merry Christmas!
fig2 = plt.imshow(final)
plt.show()

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

Challenge Problem¶

Add a homework assignment that might take 10 minutes to complete. Make sure you can work the problem yourself, but you do not need to submit a solution to the problem.

In [ ]: