Notebook

Topic 17. LU Factorization¶

Author: Daniel Free, daniel.free816@gmail.com

Introduction¶

The LU factorization method provides a way to break down a square matrix into two triangular matrices, one lower and one upper (hence L and U). Triangular matrices have very desirable properties. For example, systems of equations with a triangular matrix are easily solved with forward or backward substitution. Solving a system in this manner is numerically superior to computing the inverse of a matrix explicitly.

Explanation of the theory¶

Definition¶

A square matrix $A$ can be factored as $PA=LU$ where $P$ is a permutation matrix, $L$ is a lower triangular matrix with ones along the diagonal, and $U$ is an upper triangular matrix. Permutation matrices are orthogonal (i.e. $P^T P=PP^T=I$ ), so the factorization can also be written as $A=P^T LU$ .

Implementation¶

The main idea for computing the $LU$ factorization is to use Gauss elimination or row reduction on the matrix $A$ and keep track of each row-reducing step. To aid in this, we define an elementary matrix $E$ as an identify matrix with a few elements in a single column off the diagonal that are non-zero. This represents the row reduction step. Doing these reduction steps to progressively make an upper triangular matrix results in the matrix $U$ . For example, a $3\times3$ matrix will have 2 $E$ matrices, each one zeroing out the lower elements of a column to make the upper triangular matrix, as shown below.

$A = \begin{bmatrix} \times & \times & \times \\ \times & \times & \times \\ \times & \times & \times \end{bmatrix}$

$E_1A = \begin{bmatrix} 1 & 0 & 0 \\ \times & 1 & 0 \\ \times & 0 & 1 \end{bmatrix} \begin{bmatrix} \times & \times & \times \\ \times & \times & \times \\ \times & \times & \times \end{bmatrix} = \begin{bmatrix} \times & \times & \times \\ 0 & \times & \times \\ 0 & \times & \times \end{bmatrix}$

$E_2E_1A = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & \times & 1 \end{bmatrix} \begin{bmatrix} \times & \times & \times \\ 0 & \times & \times \\ 0 & \times & \times \end{bmatrix} = \begin{bmatrix} \times & \times & \times \\ 0 & \times & \times \\ 0 & 0 & \times \end{bmatrix} = U$

where $\times$ is a non-zero element.

To find the off diagonal elements of $E$ , simply divide the current element by the negative of the diagonal element. For example, if the first column of a $3\times 3$ matrix was $\begin{bmatrix} 4 \\ 9 \\ -1 \end{bmatrix}$ , then the first elementary matrix would be $\begin{bmatrix} 1 & 0 & 0 \\ \frac{-9}{4} & 1 & 0 \\ \frac{1}{4} & 0 & 1 \end{bmatrix}$ . However, since finding the elements of the $E$ matrices involves division, there is a possibility of accruing numerical error. To mititgate the risk of dividing by a large number by a small number and causing numerical issues, we use a permutation matrix $P_{ij}$ to reorder the rows of the current column so that the largest absolute value in the lower portion of the column is on the diagonal, where $P_{ij}$ is the identity matrix with rows $i$ and $j$ swapped (e.g. $P_{13} = \begin{bmatrix} 0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{bmatrix}, P_{23} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \end{bmatrix}$ ). Following the same example just given, we would create a permutation matrix $P_{12}$ and then multiply the resulting matrix by $E = \begin{bmatrix} 1 & 0 & 0 \\ \frac{-4}{9} & 1 & 0 \\ \frac{1}{9} & 0 & 1 \end{bmatrix}$ . In the general case, this would yield

$E_2P_{cd}E_1P_{ab}A = \begin{bmatrix} \times & \times & \times \\ 0 & \times & \times \\ 0 & 0 & \times \end{bmatrix} = U$ without fear of numerical instabilites. Since the elementary matrices and permutation matrices are all invertible, we can solve for

$A$ as shown below. Additionally, the inverse of a permutation matrix is itself, and the inverse of an elementary matrix is achieved by negating the off diagonal elements.

$A = (E_2P_{cd}E_1P_{ab})^{-1}U = P_{ab}^{-1}E_1^{-1}P_{cd}^{-1}E_2^{-1}U = VU$

The matrix $V$ is not neccessarily lower triangular. However, by applying the same permutation matrices in the original order, we can obtain the lower triangular matrix $L$ .

$P_{cd}P_{ab}A = PA = P_{cd}P_{ab}VU = LU$

where $P = P_{cd}P_{ab}$ .

Applications¶

Suppose you have a system of equations to solve in the form of $A\mathbf{x} = \mathbf{b}$ with $A$ being a square $n \times n$ matrix. Multiplying both sides of the equation by $P$ and substituting in the $LU$ factorization, we obtain $LU\mathbf{x} = P\mathbf{b}$ . This can be rewritten as $L\mathbf{y} = \mathbf{c}$ , where $\mathbf{y} = U\mathbf{x}$ and $\mathbf{c} = P\mathbf{b}$ . Written out, this looks like

$\begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ l_{21} & 1 & 0 & \cdots & 0 \\ l_{31} & l_{32} & 1 & \cdots & 0 \\ \vdots & \vdots & l_{k(n-1)} & \ddots & \vdots \\ l_{n1} & l_{n2} & \cdots & l_{n(n-1)} & 1 \end{bmatrix} \begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_n \end{bmatrix} = \begin{bmatrix} c_1 \\ c_2 \\ c_3 \\ \vdots \\ c_n \end{bmatrix}$

From this matrix equation, we can see that we will end up with $n$ linear equations. However, since $L$ is lower triangular, solving for $\mathbf{y}$ can be done using forward substitution. For example, the first equation will be

$y_1=c_1$ The next equation is

$l_{21} y_1+y_2=c_2$ Solving for

$y_2$ , we obtain

$y_2=c_2-l_{21} y_1$ which only has one unknown, since we just solved for

$y_1$ . Proceeding forward, the

$j^{th}$ equation is solved by

$y_j=c_j-\sum_{i=1}^{j-1}l_{ji} y_i$

Now that we know $\mathbf{y}$ we can solve the equation $\mathbf{y}=U\mathbf{x}$ for $\mathbf{x}$ . Similar to the lower triangular case, solving a matrix equation with an upper triangular matrix makes use of backwards substitution. The equation for the $j^{th}$ element of $\mathbf{x}$ can be solved by

$x_j=\frac{1}{u_{jj}}(y_j-\sum_{k=j+1}^n u_{jk} x_k)$ where the

$u$ terms refer to the upper triangular matrix elements.

In addition to solving matrix equations, the LU factorization can be used to compute the determinate of a square matrix. Substituting in the definition of the factorization, we get

$det⁡(A)=det⁡(P^T LU)$ By properties of a determinate, this can be written as

$det(A) = det⁡(P)det⁡(L)det⁡(U)$ The determinate of a permutation matrix

$P$ is

$\pm 1$ , where an even number of permutations yields the positive value and an odd number of permutations yields the negative value. The determinate of a triangular matrix is just the product of its diagonal values, so

$det(L) = \prod_{i=1}^n l_{ii} = \prod_{i=1}^n 1 = 1$ . Therefore,

$det⁡(A)= \pm\prod_{i=1}^n u_{ii}$ where the sign again is determined by the number of rows that have been switched in

$P$ .

Simple Numerical Examples¶

Suppose we have a matrix $A$ and we want to find the LU factorization of it. The example below shows the step by step process.

In [1]:

import numpy as np
A = np.array([[9,4,8],[1,2,7],[-6,3,5]])
print("A =\n", A)
E1 = np.array([[1.,0.,0.],[-1/9,1.,0.],[6/9,0.,1.]])
print("E1 =\n", E1)
E1A = np.dot(E1,A)
print("E1*A =\n",E1A)
P23 = np.array([[1,0,0],[0,0,1],[0,1,0]])
P23E1A = np.dot(P23,E1A)
print("P23*E1*A = \n",P23E1A)
E2 = np.array([[1,0,0],[0,1,0],[0,-P23E1A[2,1]/P23E1A[1,1],1]])
print("E2 =\n",E2)
U = np.dot(E2,P23E1A)
print("E2*P23*E1*A = U =\n", U)
V = np.dot(np.dot((2*np.identity(3) - E1),P23),(2*np.identity(3)-E2))
print("V = \n",V)
L = np.dot(P23,V)
print("L = \n",L)

A =
 [[ 9  4  8]
 [ 1  2  7]
 [-6  3  5]]
E1 =
 [[ 1.          0.          0.        ]
 [-0.11111111  1.          0.        ]
 [ 0.66666667  0.          1.        ]]
E1*A =
 [[ 9.          4.          8.        ]
 [ 0.          1.55555556  6.11111111]
 [ 0.          5.66666667 10.33333333]]
P23*E1*A = 
 [[ 9.          4.          8.        ]
 [ 0.          5.66666667 10.33333333]
 [ 0.          1.55555556  6.11111111]]
E2 =
 [[ 1.         0.         0.       ]
 [ 0.         1.         0.       ]
 [ 0.        -0.2745098  1.       ]]
E2*P23*E1*A = U =
 [[ 9.          4.          8.        ]
 [ 0.          5.66666667 10.33333333]
 [ 0.          0.          3.2745098 ]]
V = 
 [[ 1.          0.          0.        ]
 [ 0.11111111  0.2745098   1.        ]
 [-0.66666667  1.          0.        ]]
L = 
 [[ 1.          0.          0.        ]
 [-0.66666667  1.          0.        ]
 [ 0.11111111  0.2745098   1.        ]]

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 [ ]: