The forward substitution algorithm solves the linear system where is a lower triangular matrix.
A lower-triangular linear system 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:
The properties of the forward substitution algorithm are:
The code for the forward substitution algorithm to solve 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
The back substitution algorithm solves the linear system where is an upper-triangular matrix. It is the backwards version of forward substitution.
The upper-triangular system can be written as the set of linear equations:
The back substitution solution works from the bottom up to give:
The properties of the back substitution algorithm are:
The code for the back substitution algorithm to solve 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
The LU decomposition of a matrix is the pair of matrices and such that:
The properties of the LU decomposition are:
Consider the matrix
The LU factorization is
An example of a matrix which has no LU decomposition is
If we try and find the LU decomposition of this matrix then we get
Equating the individual entries gives us four equations to solve. The top-left and bottom-left entries give the two equations:
These equations have no solution, so does not have an LU decomposition.
Knowing the LU decomposition for a matrix allows us to solve the linear system using a combination of forward and back substitution. In equations this is:
where we first evaluate using forward substitution and then evaluate using back substitution.
An equivalent way to write this is to introduce a new vector defined by . This means we can rewrite as:
We have thus replaced with two linear systems: and . These two linear systems can then be solved one after the other using forward and back substitution.
The LU solve algorithm for solving the linear system written as code 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 number of operations for the LU solve algorithm is as .
Given a matrix there are many different algorithms to find the matrices and for the LU decomposition. Here we will use the recursive leading-row-column LU algorithm. This algorithm is based on writing in block form as:
In the above block form of the matrix , the entry is a scalar, is a row vector, is an column vector, and is an 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 and matrices as:
The first three equations above can be immediately evaluated to give the first row and column of and . The last equation can then have its right-hand-side evaluated, which gives the Schur complement of . We thus have the equation , which is an LU decomposition problem which we can recursively solve.
The code for the recursive leading-row-column LU algorithm to find and for is:
import numpy as np
def lu_decomp(A):
"""(L, U) = lu_decomp(A) is the LU decomposition 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
"""
n = A.shape[0]
if n == 1:
L = np.array([[1]])
U = A.copy()
return (L, U)
A11 = A[0,0]
A12 = A[0,1:]
A21 = A[1:,0]
A22 = A[1:,1:]
L11 = 1
U11 = A11
L12 = np.zeros(n-1)
U12 = A12.copy()
L21 = A21.copy() / U11
U21 = np.zeros(n-1)
S22 = A22 - np.outer(L21, U12)
(L22, U22) = lu_decomp(S22)
L = np.block([[L11, L12], [L21, L22]])
U = np.block([[U11, U12], [U21, U22]])
return (L, U)
The number of operations for the recursive leading-row-column LU decomposition algorithm is as .
We can put the above sections together to produce an algorithm for solving the system , where we first compute the LU decomposition of and then use forward and backward substitution to solve for .
The properties of this algorithm are:
The code for the linear solver using LU decomposition is: import numpy as np
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
The LU decomposition can fail when the top-left entry in the matrix 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 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 and leaves the columns fixed.
2) Full pivoting rearranges both rows and columns.
3) Scaled partial pivoting approximates full pivoting without actually rearranging columns.
The LU decomposition with partial pivoting (LUP) of an matrix is the triple of matrices , , and such that:
The properties of the LUP decomposition are:
Knowing the LUP decomposition for a matrix allows us to solve the linear system by first applying and then using the LU solver. In equations we start by taking and multiplying both sides by , 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 as .
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 , , and so that . It consists of the following steps.
1) First choose so that row in has the largest absolute first entry. That is, for all . Let be the permutation matrix that pivots (shifts) row to the first row, and leaves all other rows in order. We can explicitly write as
2) Write to denote the pivoted matrix, so .
3) Let be a permutation matrix that leaves the first row where it is, but permutes all other rows. We can write as where is an permutation matrix.
4) Factorize the (unknown) full permutation matrix as the product of and , so . This means that , which first shifts row of to the top, and then permutes the remaining rows. This is a completely general permutation matrix , but this factorization is key to enabling a recursive algorithm.
5) Using the factorization , now write the LUP factorization in block form as
6) Equating the entries in the above matrices gives the equations
7) Substituting the first three equations above into the last one and rearranging gives
8) Recurse to find the LUP decomposition of , resulting in , , and that satisfy the above equation.
9) Solve for the first rows and columns of and with the above equations to give
10) Finally, reconstruct the full matrices , , and from the component parts.
In code the recursive leading-row-column LUP algorithm for finding the LU decomposition of 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
U 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:
The computational complexity (number of operations) of the algorithm is as .
The last step in the code that computes does not do so by constructing and multiplying and . This is because this would be an step, making the whole algorithm . Instead we take advantage of the special structure of and to compute with work.
Just as with the plain LU decomposition, we can use LUP decomposition to solve the linear system . This is the linear solver using LUP decomposition algorithm.
The properties of this algorithm are:
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
Recall our example of a matrix which has no LU decomposition:
To find the LUP decomposition of , we first write the permutation matrix that shifts the second row to the top, so that the top-left entry has the largest possible magnitude. This gives