LU Decomposition for Solving Linear Equations

Learning Objectives

• Understand the use of linear system of equations.
• Use known data to set up linear system of equations for practical problems.
• Describe the factorization \({\bf A} = {\bf LU}\).
• Compare the cost of LU with other operations such as matrix-matrix multiplication.
• Implement an LU decomposition algorithm.
• Given an LU decomposition for \({\bf A}\), solve the system \({\bf Ax} = {\bf b}\).
• Give examples of matrices for which pivoting is needed.
• Implement an LUP decomposition algorithm.
• Manually compute LU and LUP decompositions.
• Compute and use LU decompositions using library functions.

Basic Idea: The "Undo" button for Linear Operations

Matrix-vector multiplication: given the data \({\bf x}\) and the operator \({\bf A}\), we can find \({\bf y}\) such that \({\bf y = Ax}\):

What if we know \({\bf y}\) but not \({\bf x}\)? We need to "undo" the transformation:

Example: Undoing the Transformation

Suppose we have known operator \({\bf A}\), known data \({\bf y}\), and unkown data \({\bf x}\) that satisfies the relationship \({\bf y = Ax}\). The values of \({\bf A}\) and \({\bf y}\) are given below.

How can we solve for \(\textbf{x} = [x_1, x_2]^T\)?

Example: Image Blurring and Recovery

The original image displaying an SSN number is stored as a \(2D\) array of real numbers between \(0\) and \(1\) (\(0\) represents a white pixel, \(1\) represents a black pixel).It has \(40\) rows of pixels and \(100\) columns of pixels. We can flatten the \(2D\) array into a \(1D\) array \(\bf{x}\) containing the \(1D\) data with dimension \(4000\). We can then apply blurring operation to data \(\bf{x}\), i.e.

where \(\bf{A}\) is the blur operator and \(\bf{y}\) is the blurred image:

To "undo" blurring to recover original image, we solve a linear system of equations using the blur operator \(\bf{A}\) and the blurred image \(\bf{y}\). The following graph shows the transformation when \(\bf{y}\)does not have any noise ("clean data"):

It is also possible to recover \(\bf{x}\) with a certain extent of noise:

To answer "How much noise can we add and still be able to recover meaningful information from the original image?" and "At which point this inverse transformation fails?", we need information about the sensitivity of the "undo" operation covered later in this course.

Back Substitution Algorithm for Upper Triangular Systems

To actually solve \({\bf A x} = {\bf b}\), we can start with an "easier" system of equations. Let's consider triangular matrices - the back substitution algorithm solves the linear system \({\bf U x} = {\bf b}\) where \({\bf U}\) is an upper-triangular matrix.

An upper-triangular linear system \({\bf U}{\bf x} = {\bf b}\) can be written in matrix form:

Observe that the upper-triangular system \({\bf U}x = b\) can be written as the set of linear equations:

The back substitution solution works from the bottom up to give:

So the general form of solution is:

Notice that there are \(n\) divisions, \(\frac{n(n-1)}{2}\) subtractions / additions, and \(\frac{n(n-1)}{2}\) multiplications, hence the computational complexity is \(\bf{O(n^2)}\).

Alternatively, we can also write \({\bf U}x = b\) as a linear combination of the columns of \(\bf{U}\):

The properties of the back substitution algorithm are:

1. If any of the diagonal elements \(U_{ii}\) are zero then the system is singular and cannot be solved.
2. If all diagonal elements of \({\bf U}\) are non-zero then the system has a unique solution.
3. The number of operations for the back substitution algorithm is \(O(n^2)\) as \(n \to \infty\).

The code for the back substitution algorithm to solve \({\bf U x} = {\bf b}\) is:

import numpy as np
def back_sub(U, b):
    """x = back_sub(U, b) is the solution to U x = b
       U must be an upper-triangular matrix
       b must be a vector of the same leading dimension as U
    """
    n = U.shape[0]
    x = np.zeros(n)
    for i in range(n-1, -1, -1):
        tmp = b[i]
        for j in range(i+1, n):
            tmp -= U[i,j] * x[j]
        x[i] = tmp / U[i,i]
    return x

Example: Back Substitution for an Upper Triangular System

How can we solve for \(x = [x_1, x_2, x_3, x_4]^T\)?

Forward Substitution Algorithm for Lower Triangular Systems

The forward substitution algorithm solves the linear system \({\bf Lx} = {\bf b}\) where \({\bf L}\) is a lower triangular matrix. It is the reversed version of back substitution.

A lower-triangular linear system \({\bf L}{\bf x} = {\bf b}\) can be written in matrix form:

This can also be written as the set of linear equations:

The forward substitution algorithm solves a lower-triangular linear system by working from the top down and solving each variable in turn. In math this is:

So the general form of solution is:

Notice that there are also \(n\) divisions, \(\frac{n(n-1)}{2}\) subtractions / additions, and \(\frac{n(n-1)}{2}\) multiplications, hence the computational complexity is \(\bf{O(n^2)}\).

The properties of the forward substitution algorithm are:

1. If any of the diagonal elements \(L_{ii}\) are zero then the system is singular and cannot be solved.
2. If all diagonal elements of \({\bf L}\) are non-zero then the system has a unique solution.
3. The number of operations for the forward substitution algorithm is \(O(n^2)\) as \(n \to \infty\).

The code for the forward substitution algorithm to solve \({\bf L x} = {\bf b}\) is:

import numpy as np
def forward_sub(L, b):
    """x = forward_sub(L, b) is the solution to L x = b
       L must be a lower-triangular matrix
       b must be a vector of the same leading dimension as L
    """
    n = L.shape[0]
    x = np.zeros(n)
    for i in range(n):
        tmp = b[i]
        for j in range(i):
            tmp -= L[i,j] * x[j]
        x[i] = tmp / L[i,i]
    return x

Example: Forward-substitution for a Lower Triangular System

How can we solve for \(x = [x_1, x_2, x_3, x_4]^T\)?

LU Decomposition Definition

To solve \({\bf A x} = {\bf b}\) when \(\bf{A}\) is a non-triangular matrix, We can perform LU factorization: given an \(n \times n\) matrix \(\bf{A}\), the LU decomposition of a matrix \({\bf A}\) is the pair of matrices \({\bf L}\) and \({\bf U}\) such that:

1. \({\bf A} = {\bf LU}\)
2. \({\bf L}\) is a lower-triangular matrix with all diagonal entries equal to 1
3. \({\bf U}\) is an upper-triangular matrix.

The properties of the LU decomposition are:

1. The LU decomposition may not exist for a matrix \({\bf A}\).
2. If the LU decomposition exists then it is unique.
3. The LU decomposition provides an efficient means of solving linear equations.
4. The reason that \({\bf L}\) has all diagonal entries set to 1 is that this means the LU decomposition is unique. This choice is somewhat arbitrary (we could have decided that \({\bf U}\) must have 1 on the diagonal) but it is the standard choice.
5. We use the terms decomposition and factorization interchangeably to mean writing a matrix as a product of two or more other matrices, generally with some defined properties (such as lower/upper triangular).

Example: LU Decomposition

Consider a \(3 \times 3\) matrix \(A = \begin{bmatrix} 1 & 2 & 2 \\ 4 & 4 & 2 \\ 4 & 6 & 4 \end{bmatrix} .\)

The LU factorization is \({\bf A} = {\bf LU} = \begin{bmatrix} 1 & 0 & 0 \\ 4 & 1 & 0 \\ 4 & 0.5 & 1 \end{bmatrix} \begin{bmatrix} 1 & 2 & 2 \\ 0 & -4 & -6 \\ 0 & 0 & -1 \end{bmatrix}.\)

Solving LU-Decomposed Linear Systems

Knowing the LU decomposition for a matrix \({\bf A}\) allows us to solve the linear system \({\bf A x} = {\bf b}\) using a combination of forward and back substitution. In equations this is:

where we first evaluate \({\bf L}^{-1} {\bf b}\) using forward substitution and then evaluate \({\bf x} = {\bf U}^{-1} ({\bf L}^{-1} {\bf b})\) using back substitution.

An equivalent way to write this is to introduce a new vector \({\bf y}\) defined by \(\bf{y} = {\bf U x}\). Assuming the LU factorization of matrix \(\bf{A}\) is known, we can solve the general system

By solving two triangular systems:

We have thus replaced \({\bf A x} = {\bf b}\) with two linear systems: \({\bf L y} = {\bf b}\) and \({\bf U x} = {\bf y}\). These two linear systems can then be solved one after the other using forward and back substitution.

The number of operations for the LU solve algorithm is \(O(n^2)\) as \(n \to \infty\).

The LU solve algorithm code for solving the linear system \({\bf L U x} = {\bf b}\) is:

import numpy as np
def lu_solve(L, U, b):
    """x = lu_solve(L, U, b) is the solution to L U x = b
       L must be a lower-triangular matrix
       U must be an upper-triangular matrix of the same size as L
       b must be a vector of the same leading dimension as L
    """
    y = forward_sub(L, b)
    x = back_sub(U, y)
    return x

The LU Decomposition Algorithm

Given a matrix \({\bf A}\) there are many different algorithms to find the matrices \({\bf L}\) and \({\bf U}\) for the LU decomposition. Here we will use the recursive leading-row-column LU algorithm. Let's first consider the simple case when \({\bf A}\) is a \(2\times 2\) matrix:

Use a similar idea to write \({\bf A}\) in block form when \({\bf A}\) is an \(n\times n\) matrix:

Then \({\bf A} = {\bf LU}\) can be re-written as:

In the above block form of the \(n \times n\) matrix \({\bf A}\), the entry \(a_{11}\) is a scalar, \(\boldsymbol{a}_{12}\) is a \(1 \times (n-1)\) row vector, \(\boldsymbol{a}_{12}\) is an \((n-1) \times 1\) column vector, and \({\bf A}_{22}\) is an \((n-1) \times (n-1)\) matrix.

Comparing the left- and right-hand side entries of the above block matrix equation we see that:

These four equations can be rearranged to solve for the components of the \({\bf L}\) and \({\bf U}\) matrices as:

Notice that:

1. The first row of \(\bf{U}\) is the first row of \(\bf{A}\).
2. The first column of \(\bf{L}\) is \(\frac{\text{the first column of }\textbf{A}}{u_{11}}\).
3. \({\bf L}_{22} {\bf U}_{22}\) needs another factorization.

In other words, the first three equations above can be immediately evaluated to give the first row and column of \({\bf L}\) and \({\bf U}\). The last equation can then have its right-hand-side evaluated, which gives the Schur complement \(S_{22}\) of \({\bf A}\). We thus have the equation \({\bf L}_{22} {\bf U}_{22} = {\bf S}_{22}\), which is an \((n-1) \times (n-1)\) LU decomposition problem which we can recursively solve.

The code for the recursive leading-row-column LU algorithm to find \({\bf L}\) and \({\bf U}\) for \({\bf A} = {\bf LU}\) is:

import numpy as np
def lu_decomp(A):
    """(L, U) = lu_decomp(A) is the LU decomposition A = L U
       A is any square matrix
       L will be a lower-triangular matrix with 1 on the diagonal
       U will be an upper-triangular matrix
    """
    n = A.shape[0]
    
    # Initialize L and U as zero matrices of the same shape as A
    L = np.eye(n)  # L is initialized with 1s on the diagonal
    U = np.zeros((n, n))
    
    for i in range(n):
        # Compute the U matrix (upper triangular)
        for j in range(i, n):
            U[i, j] = A[i, j] - np.dot(L[i, :i], U[:i, j])
        
        # Compute the L matrix (lower triangular)
        for j in range(i + 1, n):
            if U[i, i] == 0:
                raise ZeroDivisionError("Zero pivot encountered. Cannot factorize.")
            L[j, i] = (A[j, i] - np.dot(L[j, :i], U[:i, i])) / U[i, i]

    return L, U

Number of divisions: \((n-1)+(n-2)+\ldots+1=\frac{n(n-1)}{2}\)

Number of multiplications: \((n-1)^2+(n-2)^2+\ldots+1^2=\frac{n^3}{3} - \frac{n^2}{2} + \frac{n}{6}\)

Number of subtractions: \((n-1)^2+(n-2)^2+\ldots+1^2=\frac{n^3}{3} - \frac{n^2}{2} + \frac{n}{6}\)

Hence the number of operations for the recursive leading-row-column LU decomposition algorithm is \(\bf{O(n^3)}\) as \(n \to \infty\).

Example: LU decomposition

Consider the matrix \(\bf{A} = \begin{bmatrix} 2 & 8 & 4 & 1 \\ 1 & 2 & 3 & 3 \\ 1 & 2 & 6 & 2 \\ 1 & 3 & 4 & 2 \end{bmatrix} .\)

How can we find an LU Decomposition for this matrix?

Example: Matrix for which LU Decomposition Fails

An example of a matrix which has no LU decomposition is

Why?

Solving Linear Systems Using LU Decomposition

We can put the above sections together to produce an algorithm for solving the system \({\bf A x} = {\bf b}\), where we first compute the LU decomposition of \({\bf A}\) and then use forward and back substitution to solve for \({\bf x}\).

The properties of this algorithm are:

1. The algorithm may fail, even if \({\bf A}\) is invertible.
2. The number of operations in the algorithm is \(\mathcal{O}(n^3)\) as \(n \to \infty\).

The code for the linear solver using LU decomposition is:

import numpy as np
def linear_solve_without_pivoting(A, b):
    """x = linear_solve_without_pivoting(A, b) is the solution to A x = b (computed without pivoting)
       A is any matrix
       b is a vector of the same leading dimension as A
       x will be a vector of the same leading dimension as A
    """
    (L, U) = lu_decomp(A)
    x = lu_solve(L, U, b)
    return x

Pivoting

The LU decomposition can fail when the top-left entry in the matrix \({\bf A}\) is zero or very small compared to other entries. Pivoting is a strategy to mitigate this problem by rearranging the rows and/or columns of \({\bf A}\) to put a larger element in the top-left position.

There are many different pivoting algorithms. The most common of these are full pivoting, partial pivoting, and scaled partial pivoting. We will only discuss partial pivoting in detail.

1) Partial pivoting only rearranges the rows of \({\bf A}\) and leaves the columns fixed.

2) Full pivoting rearranges both rows and columns.

3) Scaled partial pivoting approximates full pivoting without actually rearranging columns.

LU Decomposition with Partial Pivoting

The LU decomposition with partial pivoting (LUP) of an \(n \times n\) matrix \({\bf A}\) is the triple of matrices \({\bf L}\), \({\bf U}\), and \({\bf P}\) such that:

1. \({\bf P A} = {\bf LU} \)
2. \({\bf L}\) is an \(n \times n\) lower-triangular matrix with all diagonal entries equal to 1.
3. \({\bf U}\) is an \(n \times n\) upper-triangular matrix.
4. \({\bf P}\) is an \(n \times n\) permutation matrix.

The properties of the LUP decomposition are:

1. The permutation matrix \({\bf P}\) acts to permute the rows of \({\bf A}\). This attempts to put large entries in the top-left position of \({\bf A}\) and each sub-matrix in the recursion, to avoid needing to divide by a small or zero element.
2. The LUP decomposition always exists for a matrix \({\bf A}\).
3. The LUP decomposition of a matrix \({\bf A}\) is not unique.
4. The LUP decomposition provides a more robust method of solving linear systems than LU decomposition without pivoting, and it is approximately the same cost.

LU factorization with partial pivoting can be completed for any matrix A: Suppose you are at stage k and there is no non-zero entry on or below the diagonal in column k. At this point, there is nothing else you can do, so the algorithm leaves a zero in the diagonal entry of U. Note that the matrix U is singular, and so is the matrix A. Subsequent back substitutions using U will fail, but the LU factorization itself is still completed.

Solving LUP decomposition linear systems

Knowing the LUP decomposition for a matrix \({\bf A}\) allows us to solve the linear system \({\bf A x} = {\bf b}\) by first applying \({\bf P}\) and then using the LU solver. In equations we start by taking \({\bf A x} = {\bf b}\) and multiplying both sides by \({\bf P}\), giving

The code for the LUP solve algorithm to solve the linear system \({\bf L U x} = {\bf P b}\) is:

import numpy as np
def lup_solve(L, U, P, b):
    """x = lup_solve(L, U, P, b) is the solution to L U x = P b
       L must be a lower-triangular matrix
       U must be an upper-triangular matrix of the same shape as L
       P must be a permutation matrix of the same shape as L
       b must be a vector of the same leading dimension as L
    """
    z = np.dot(P, b)
    x = lu_solve(L, U, z)
    return x

The number of operations for the LUP solve algorithm is \(\mathcal{O}(n^2)\) as \(n \to \infty\).

The LUP decomposition algorithm

Just as there are different LU decomposition algorithms, there are also different algorithms to find a LUP decomposition. Here we use the recursive leading-row-column LUP algorithm.

This algorithm is a recursive method for finding \({\bf L}\), \({\bf U}\), and \({\bf P}\) so that \({\bf P A} = {\bf L U}\). It consists of the following steps.

1. First choose \(i\) so that row \(i\) in \({\bf A}\) has the largest absolute first entry. That is, \(\vert A_{i1}\vert \ge \vert A_{j1}\vert\) for all \(j\). Let \({\bf P}_1\) be the permutation matrix that pivots (shifts) row \(i\) to the first row, and leaves all other rows in order. We can explicitly write \({\bf P}_1\) as

   \[{\bf P}_1 = \begin{bmatrix} 0_{1(i-1)} & 1 & 0_{1(n-i)} \\ I_{(i-1)(i-1)} & 0 & 0_{(i-1)(n-i)} \\ 0_{(n-i)(i-1)} & 0 & I_{(n-i)(n-i)} \end{bmatrix} = \begin{bmatrix} 0 & \ldots & 0 & 1 & 0 & \ldots & 0 \\ 1 & \ldots & 0 & 0 & 0 & \ldots & 0 \\ \vdots & \ddots & \vdots & \vdots & \vdots & \ldots & 0 \\ 0 & \ldots & 1 & 0 & 0 & \ldots & 0 \\ 0 & \ldots & 0 & 0 & 1 & \ldots & 0 \\ \vdots & \ddots & \vdots & \vdots & \vdots & \ldots & 0 \\ 0 & \ldots & 0 & 0 & 0 & \ldots & 1 \end{bmatrix}.\]

2. Write \(\bar{ {\bf A} }\) to denote the pivoted \({\bf A}\) matrix, so \(\bar{ {\bf A} } = {\bf P}_1 {\bf A}\).

3. Let \({\bf P}_2\) be a permutation matrix that leaves the first row where it is, but permutes all other rows. We can write \({\bf P}_2\) as \({\bf P}_2 = \begin{bmatrix} 1 & \boldsymbol{0} \\ \boldsymbol{0} & P_{22} \end{bmatrix},\) where \({\bf P}_{22}\) is an \((n-1) \times (n-1)\) permutation matrix.

4. Factorize the (unknown) full permutation matrix \({\bf P}\) as the product of \({\bf P}_2\) and \({\bf P}_1\), so \({\bf P} = {\bf P}_2 {\bf P}_1\). This means that \({\bf P} A = {\bf P}_2 {\bf P}_1 A = {\bf P}_2 \bar{ {\bf A} }\), which first shifts row \(i\) of \({\bf A}\) to the top, and then permutes the remaining rows. This is a completely general permutation matrix \({\bf P}\), but this factorization is key to enabling a recursive algorithm.

5. Using the factorization \({\bf P} = {\bf P}_2 {\bf P}_1\), now write the LUP factorization in block form as

   \[\begin{aligned} {\bf P A} &= {\bf L U} \\ {\bf P_2} \bar{ {\bf A} } &= {\bf L U} \\ \begin{bmatrix} 1 & \boldsymbol{0} \\ \boldsymbol{0} & {\bf P}_{22} \end{bmatrix} \begin{bmatrix} \bar{a}_{11} & \bar{\boldsymbol{a}}_{12} \\ \bar{\boldsymbol{a}}_{21} & \bar{ {\bf A} }_{22} \end{bmatrix} &= \begin{bmatrix} 1 & \boldsymbol{0} \\ \boldsymbol{\ell}_{21} & {\bf L}_{22} \end{bmatrix} \begin{bmatrix} u_{11} & \boldsymbol{u}_{12} \\ \boldsymbol{0} & {\bf U}_{22} \end{bmatrix} \\ \begin{bmatrix} \bar{a}_{11} & \bar{\boldsymbol{a}}_{12} \\ {\bf P}_{22} \bar{\boldsymbol{a}}_{21} & {\bf P}_{22} \bar{ {\bf A} }_{22} \end{bmatrix} &= \begin{bmatrix} u_{11} & \boldsymbol{u}_{12} \\ u_{11} \boldsymbol{\ell}_{21} & (\boldsymbol{\ell}_{21} \boldsymbol{u}_{12} + {\bf L}_{22} {\bf U}_{22}) \end{bmatrix} \end{aligned}\]

6. Equating the entries in the above matrices gives the equations

   \[\begin{aligned} \bar{a}_{11} &= u_{11} \\ \bar{\boldsymbol{a}}_{12} &= \boldsymbol{u}_{12} \\ {\bf P}_{22} \bar{\boldsymbol{a}}_{21} &= u_{11} \boldsymbol{\ell}_{21} \\ {\bf P}_{22} \bar{A}_{22} &= \boldsymbol{\ell}_{21} \boldsymbol{u}_{12} + {\bf L}_{22} {\bf U}_{22}. \end{aligned}\]

7. Substituting the first three equations above into the last one and rearranging gives

   \[{\bf P}_{22} \underbrace{\Bigl(\bar{A}_{22} - \bar{\boldsymbol{a}}_{21} (\bar{a}_{11})^{-1} \bar{\boldsymbol{a}}_{12}\Bigr)}_{\text{Schur complement } {\bf S}_{22}} = {\bf L}_{22} {\bf U}_{22}.\]

8. Recurse to find the LUP decomposition of \(S_{22}\), resulting in \({\bf L}_{22}\), \({\bf U}_{22}\), and \({\bf P}_{22}\) that satisfy the above equation.

9. Solve for the first rows and columns of \({\bf L}\) and \({\bf U}\) with the above equations to give

   \[\begin{aligned} u_{11} &= \bar{a}_{11} \\ \boldsymbol{u}_{12} &= \bar{\boldsymbol{a}}_{12} \\ \boldsymbol{\ell}_{21} &= \frac{1}{\bar{a}_{11}} {\bf P}_{22} \bar{\boldsymbol{a}}_{21}. \end{aligned}\]

10. Finally, reconstruct the full matrices \({\bf L}\), \({\bf U}\), and \({\bf P}\) from the component parts.

In code the recursive leading-row-column LUP algorithm for finding the LU decomposition of \({\bf A}\) with partial pivoting is:

import numpy as np
def lup_decomp(A):
    """(L, U, P) = lup_decomp(A) is the LUP decomposition P A = L U
       A is any matrix
       L will be a lower-triangular matrix with 1 on the diagonal, the same shape as A
       U will be an upper-triangular matrix, the same shape as A
       P will be a permutation matrix, the same shape as A
    """
    n = A.shape[0]
    if n == 1:
        L = np.array([[1]])
        U = A.copy()
        P = np.array([[1]])
        return (L, U, P)

    i = np.argmax(A[:,0])
    A_bar = np.vstack([A[i,:], A[:i,:], A[(i+1):,:]])

    A_bar11 = A_bar[0,0]
    A_bar12 = A_bar[0,1:]
    A_bar21 = A_bar[1:,0]
    A_bar22 = A_bar[1:,1:]

    S22 = A_bar22 - np.dot(A_bar21, A_bar12) / A_bar11

    (L22, U22, P22) = lup_decomp(S22)

    L11 = 1
    U11 = A_bar11

    L12 = np.zeros(n-1)
    U12 = A_bar12.copy()

    L21 = np.dot(P22, A_bar21) / A_bar11
    U21 = np.zeros(n-1)

    L = np.block([[L11, L12], [L21, L22]])
    U = np.block([[U11, U12], [U21, U22]])
    P = np.block([
        [np.zeros((1, i-1)), 1,                  np.zeros((1, n-i))],
        [P22[:,:(i-1)],      np.zeros((n-1, 1)), P22[:,i:]]
    ])
    return (L, U, P)

The properties of the recursive leading-row-column LUP decomposition algorithm are:

1. The computational complexity (number of operations) of the algorithm is \(\mathcal{O}(n^3)\) as \(n \to \infty\).

2. The last step in the code that computes \({\bf P}\) does not do so by constructing and multiplying \({\bf P}_2\) and \({\bf P}_1\). This is because this would be an \(\mathcal{O}(n^3)\) step, making the whole algorithm \(\mathcal{O}(n^4)\). Instead we take advantage of the special structure of \({\bf P}_2\) and \({\bf P}_1\) to compute \({\bf P}\) with \(\mathcal{O}(n^2)\) work.

Solving General Linear Systems using LUP Decomposition

Just as with the plain LU decomposition, we can use LUP decomposition to solve the linear system \({\bf A x }= {\bf b}\). This is the linear solver using LUP decomposition algorithm.

The properties of this algorithm are:

1. The algorithm may fail. In particular if \({\bf A}\) is singular (or singular in finite precision), U will have a zero on it's diagonal.
2. As LU factorization has \(\mathcal{O}(n^3)\) opeations, and solving \(Ly = b\) and \(Ux = y\) both have \(\mathcal{O}(n^2)\) opeations, the number of total operations in this algorithm is \(\bf{\mathcal{O}(n^3)}\) as \(n \to \infty\).
3. We decouple the factorization step from the actual solve, since in enigeering applications, it might be possible to reuse the LU factorization of the same matrix in solving multiple linear systems.

The code for the linear solver using LUP decomposition is:

import numpy as np
def linear_solve(A, b):
    """x = linear_solve(A, b) is the solution to A x = b (computed with partial pivoting)
       A is any matrix
       b is a vector of the same leading dimension as A
       x will be a vector of the same leading dimension as A
    """
    (L, U, P) = lup_decomp(A)
    x = lup_solve(L, U, P, b)
    return x

Example: LUP Decomposition

Consider the matrix \(\bf{A} = \begin{bmatrix} 2 & 1 & 1 & 0 \\ 4 & 3 & 3 & 1 \\ 8 & 7 & 9 & 5 \\ 6 & 7 & 9 & 8 \end{bmatrix} .\)

How can we find an LUP Decomposition for this matrix?

Example: Matrix for which LUP Decomposition Succeeds but LU Decomposition Fails

Consider a matrix which has no LU decomposition:

How can we find an LUP Decomposition for this matrix?

Review Questions

1. Given a factorization \({\bf P A} = {\bf LU}\), how would you solve the system \({\bf A}\mathbf{x} = \mathbf{b}\)?
2. Understand the process of solving a triangular system. Solve an example triangular system.
3. Recognize and understand Python code implementing forward substitution, back substitution, and LU factorization.
4. When does an LU factorization exist?
5. When does an LUP factorization exist?
6. What special properties do \({\bf P}\), \({\bf L}\) and \({\bf U}\) have?
7. Can we find an LUP factorization of a singular matrix?
8. What happens if we try to solve a system \({\bf A}\mathbf{x} = \mathbf{b}\) with a singular matrix \({\bf A}\)?
9. Compute the LU factorization of a small matrix by hand.
10. Why do we use pivoting when solving linear systems?
11. How do we choose a pivot element?
12. What effect does a given permutation matrix have when multiplied by another matrix?
13. What is the cost of matrix-matrix multiplication?
14. What is the cost of computing an LU or LUP factorization?
15. What is the cost of forward or back substitution?
16. What is the cost of solving \({\bf A}\mathbf{x} = \mathbf{b}\) for a general matrix?
17. What is the cost of solving \({\bf A}\mathbf{x} = \mathbf{b}\) for a triangular matrix?
18. What is the cost of solving \({\bf A}\mathbf{x} = \mathbf{b_i}\) with the same matrix \({\bf A}\) and several right-hand side vectors \(\mathbf{b}_i\)?
19. Given a process that takes time \(\mathcal{O}(n^k)\), what happens to the runtime if we double the input size (i.e. double \(n\))? What if we triple the input size?

Changelog

• October 2nd, 2024: Dev Singh (dsingh14) — fix lu_decomp
• September 30th, 2024: Kaiyao Ke (kaiyaok2) — minor bug fix
• February 13th, 2024: Kaiyao Ke (kaiyaok2) — aligned notes with slides, added examples and refactored existing notes

• View Remaining Entries

  • October 20th, 2023: Yuxuan Chen (yuxuan19) — update lu_decomp() code
  • February 28th, 2018: Erin Carrier (ecarrie2) — fix error in ludecomp() code
  • February 22nd, 2018: Erin Carrier (ecarrie2) — update properties for solving using LUP
  • January 14th, 2018: Erin Carrier (ecarrie2) — removes demo links
  • November 2nd, 2017: Arun Lakshmanan (lakshma2) — minor fix in lup_solve(), add changelog
  • November 2nd, 2017: John Doherty (jjdoher2) — fixed typo in back substitution
  • October 25th, 2017: Nathan Bowman (nlbowma2) — added review questions
  • October 23rd, 2017: Erin Carrier (ecarrie2) — fix links
  • October 20th, 2017: Matthew West (mwest) — minor fix in back_sub()
  • October 19th, 2017: Nathan Bowman (nlbowma2) — minor existence of LUP
  • October 17th, 2017: Erin Carrier (ecarrie2) — fixes
  • October 17th, 2017: Luke Olson (lukeo) — update links
  • October 16th, 2017: Matthew West (mwest) — first draft complete

Authors

• CS 357 Course Staff