3. Linear Systems#
In this chapter we explore the implementation of a both direct and iterative methods for solving systems of linear equations. The corresponding theory is covered in Chapter 3 of the lecture notes.
3.1. Forward substitution#
import numpy as np
def forward_substitution(L,b):
'''
Implements the forward substitution algorithm to solve a lower triangular system.
- L: an n x n real lower triangular matrix (numpy array)
- b: an n dimensional real vector (numpy array)
'''
# Identifies the dimension of the matrix
n = L.shape[0]
# Initializes the solution as a vector of zeros
x = np.zeros(n)
for i in range(n):
# Initializes the numerator of the expression of x_i
numerator = b[i]
for j in range(i):
# Updates the numerator
numerator = numerator - L[i,j] * x[j]
# Divides by \ell_ii
x[i] = numerator / L[i,i]
# Returns the solution
return x
3.1.1. Example 1#
L = np.matrix([[2,0, 0], [1,3,0],[2,2,1]])
b = np.array([1, 2, 3])
forward_substitution(L,b)
array([0.5, 0.5, 1. ])
3.2. Backward substitution#
import numpy as np
def backward_substitution(U,b):
'''
Implements the backward substitution algorithm to solve an upper triangular system.
- U: an n x n real upper triangular matrix (numpy array)
- b: an n dimensional real vector (numpy array)
'''
# Identifies the dimension of the matrix U
n = U.shape[0]
# Initializes x with all zeros.
x = np.zeros(n)
# For i = n-1,...,0 with a step of -1
for i in range(n-1,-1,-1):
# Initializes the fraction that gives x_i
numerator = b[i]
# For j = i+1,...,n-1
for j in range(i+1,n):
numerator = numerator - U[i,j] * x[j]
x[i] = numerator / U[i,i]
# Returns x
return x
3.2.1. Example 1#
U = np.matrix([[2,2,1],[0,2,1],[0,0, 2]])
b = np.array([1, 2, 3])
backward_substitution(U,b)
array([-0.5 , 0.25, 1.5 ])
3.3. Swapping rows#
Swapping rows will be useful in the methods introduced later. In particular, it will be used to prevent division by zero.
def swap_rows(A,b,row1,row2):
'''
Swaps row1 and row2 of matrix A and vector b.
It assumes that the rows are numbered from 1 to n.
So it actually swaps the rows at positions row1-1 and row2-1.
'''
# Swaps the rows of A
temp = A[[row1-1],:]
A[[row1-1],:] = A[[row2-1],:]
A[[row2-1],:] = temp
# Swap row k and row max_index of b
temp = b[row1-1]
b[row1-1] = b[row2-1]
b[row2-1] = temp
return A,b
3.3.1. Example#
A = np.array([[1.0, 1.0, 1.0],
[0.0, 2.0, 5.0],
[2.0, 5.0, -1.0]])
b = np.array([6.0, -4.0, 27.0])
swap_rows(A,b,2,3)
(array([[ 1., 1., 1.],
[ 2., 5., -1.],
[ 0., 2., 5.]]),
array([ 6., 27., -4.]))
3.4. Gaussian elimination#
We provide an implementation of Gaussian elimination in terms of elementary operations.
import numpy as np
def gaussian_elimination_elementary(A, b):
'''
Performs Gaussian elimination on the matrix A with the vector b.
- A: n x n numpy matrix
- b: n numpy vector
- The number of steps to perform.
'''
# Determines the size of A
n = A.shape[0]
print(f"Size {n}")
# Sets the iteration number
k = 1
# Loops until k = n-1
while k <= n-1:
print(f"Step {k}")
# Checks if the pivot element is null
if A[k-1, k-1] == 0: # Recall that we count from k=1 but elements are indexed from 0 to n-1.
# If the pivot element is 0 we stop
print("The pivot element is 0")
# Looks for the index of the largest element below the pivot
max_element = float('-inf')
max_index = k
for l in range(k+1,n+1): # We loop up to n+1 because we need to include row n too.
if abs(A[l-1,k-1] > max_element):
max_element = A[l-1,k-1]
max_index = l
print(f"Swapping rows {k} and {max_index}")
if max_index == k:
print("The matrix is singular")
return
else:
# Swap row k and row max_index
swap_rows(A,b,k,max_index)
# Iterates over the rows, from k+1 to n
for i in range(k+1,n+1):
# Computes the pivot element
m = A[i-1, k-1] / A[k-1, k-1]
# Iterates of the columns, from k to n
for j in range(k,n+1):
# Performs the elimination
if j == k:
A[i-1, j-1] = 0
else:
A[i-1, j-1] = A[i-1, j-1] - (m * A[k-1, j-1])
# Updates the right-hand side
b[i-1] = b[i-1] - m * b[k-1]
# Prints the matrix A and the vector b after the k-th iteration.
print(f"A^({k}): {A}")
print(f"b^({k}): {b}")
# Increments the iteration counter
k = k + 1
return A, b
3.4.1. Example 1#
A = np.array([[1.0, 1.0, 1.0],
[0.0, 2.0, 5.0],
[2.0, 5.0, -1.0]])
b = np.array([6.0, -4.0, 27.0])
U, y = gaussian_elimination_elementary(A,b)
x = backward_substitution(U,y)
print(f"Solution x = {x}")
Size 3
Step 1
A^(1): [[ 1. 1. 1.]
[ 0. 2. 5.]
[ 0. 3. -3.]]
b^(1): [ 6. -4. 15.]
Step 2
A^(2): [[ 1. 1. 1. ]
[ 0. 2. 5. ]
[ 0. 0. -10.5]]
b^(2): [ 6. -4. 21.]
Solution x = [ 5. 3. -2.]
3.4.2. Example 2#
Here is an example that requires swapping rows.
A = np.array([[1.0, 1.0, 1.0],
[2.0, 2.0, 5.0],
[2.0, 5.0, -1.0]])
b = np.array([1.0,2.0,3.0])
U, y = gaussian_elimination_elementary(A,b)
x = backward_substitution(U,y)
print(f"Solution x = {x}")
Size 3
Step 1
A^(1): [[ 1. 1. 1.]
[ 0. 0. 3.]
[ 0. 3. -3.]]
b^(1): [1. 0. 1.]
Step 2
The pivot element is 0
Swapping rows 2 and 3
A^(2): [[ 1. 1. 1.]
[ 0. 3. -3.]
[ 0. 0. 3.]]
b^(2): [1. 1. 0.]
Solution x = [0.66666667 0.33333333 0. ]
3.5. LU decomposition#
We provide an implementation of Doolittle’s factorization.
import numpy as np
def LUfactorization(A):
'''
Implements the LU factorization of a matrix A.
It uses a recursive algorithm.
- A: n x n numpy matrix
'''
# Get the length of the matrix
n=len(A)
# Initialize L and U with zeros
L = np.zeros(shape=(n,n))
U = np.zeros(shape=(n,n))
# The operations on L_11 and U_11
# are always performed, independently of the size of the matrix.
L[0,0] = 1.0
U[0,0] = A[0,0]
# Base case
if n==1:
# If n = 1 we are done
return L,U
else:
L[0,0] = 1.0
U[0,0] = A[0,0]
# Populates row 1 (index 0) of U, from column 2 (index 1) to n (index n-1)
U[[0],1:n] = A[[0],1:n]/L[0,0]
# Populates column 1 (index 0) of L, from row 2 (index 1) to n (index n-1)
L[1:n,[0]] = A[1:n,[0]]/U[0,0]
# Recursion step on the matrix M = A22-\ell_1u_1
M = A[1:n,1:n] - np.matmul(L[1:n,[0]],U[[0],1:n])
(L22,U22) = LUfactorization(M)
# Puts the decomposed submatrices in their respective positions
L[1:n,1:n] = L22
U[1:n,1:n] = U22
return L,U
3.5.1. Example 1#
A = np.array([[1.0, 1.0, 1.0],
[0.0, 2.0, 5.0],
[2.0, 5.0, -1.0]])
L,U = LUfactorization(A)
print(f"L={L}")
print(f"U={U}")
L=[[1. 0. 0. ]
[0. 1. 0. ]
[2. 1.5 1. ]]
U=[[ 1. 1. 1. ]
[ 0. 2. 5. ]
[ 0. 0. -10.5]]
Note how the \(U\) matrix is the same as the matrix produced by Gaussian elimination in Exercise 1.
3.6. Cholesky decomposition#
Here we provide a basic implementation of Cholesky decomposition.
import numpy as np
def cholesky_decomposition(A):
'''
Performs Cholesky decomposition on a given symmetric and positive-definite matrix.
- A (numpy.ndarray): A symmetric positive-definite matrix.
'''
# Here we introduce a safety check.
# We ensure that the matrix is a NumPy array.
A = np.array(A, dtype=float)
# Check if the matrix is square
if A.shape[0] != A.shape[1]:
# This is how you signal an error.
raise ValueError("The matrix A must be square")
# Gets the dimension of the matrix
n = A.shape[0]
# Initializes L with zeros
L=np.zeros(shape=(n,n))
# k iterates over the columns
for k in range(n):
# Calculates the diagonal element (k,k)
l_kk = np.sqrt(A[k, k] - np.sum(L[k, :k]**2)) # np.sum sums the square of the elements in row k from position 1 (index 0) to position k (index k-1).
# Inserts the element into L at position (k,k).
L[k, k] = l_kk
# Calculates the remaining elements of column k.
# i iterates over the rows with index larger than k.
for i in range(k+1, n):
l_ik = (A[i, k] - np.sum(L[i, :k] * L[k, :k])) / l_kk # np.sum sums the products of the elements in row i up to position k and row k up to poisition k
# Inserts the element into L row i column k
L[i, k] = l_ik
return L
3.6.1. Example 1#
A = np.array([[4, 3, 2,1],
[3, 3, 2,1],
[2, 2, 2, 1],
[1,1,1,1]])
L = cholesky_decomposition(A)
print("\nLower triangular matrix L:")
print(L)
Lower triangular matrix L:
[[2. 0. 0. 0. ]
[1.5 0.8660254 0. 0. ]
[1. 0.57735027 0.81649658 0. ]
[0.5 0.28867513 0.40824829 0.70710678]]
Let us verify that \(LL^\top\) gives \(A\).
print("Product LL^T:")
np.matmul(L,L.T)
Product LL^T:
array([[4., 3., 2., 1.],
[3., 3., 2., 1.],
[2., 2., 2., 1.],
[1., 1., 1., 1.]])
3.7. Jacobi Iteration Method#
Here we provide a basic implementation of the Jacobi iteration method.
We begin with an implementation of the \(\ell_\infty\) norm which we will use to compute the error.
import numpy as np
def norm_infinity(x):
'''
Computes the infinity norm of a vector x.
- x: numpy vector
'''
return np.max(np.abs(x))
Here is a basic Jacobi iteration method. We arbitrary choose to compute the error as \(\frac{\lvert\lvert \mathbf{x}_{k+1}-\mathbf{x}_{k}\rvert\rvert}{\lvert\lvert\mathbf{x}_{k+1}\rvert\rvert}\). Observe that we also pass a maximum number of iterations. This is often a good idea in iterative methods in case convergence to a solution with the desider error is too slow.
import numpy as np
def jacobi_iteration_method(A, b, x0, epsilon, max_iterations):
'''
Solves a the linear system Ax = b using the Jacobi iteration method.
- A: n x n numpy matrix
- b: n numpy vector
- x0: initial guess for the solution
- max_iterations: maximum number of iterations that can be performed
'''
# Extracts n
n = A.shape[0]
# Sets the current x.
currentx = x0
# Sets the iteration counter
k = 1
# Initializes the error
err = float('inf')
# Performs the main iterations
while (err >= epsilon) and (k <= max_iterations):
# Initializes the new approximate solution as a vector of zeros
newx = np.zeros(n)
for i in range(n):
# Populates component i (observe that here i=0,...,n-1)
newx[i] = (b[i] - sum([A[i,j]*currentx[j] if j != i else 0 for j in range(n)])) / A[i, i]
# Computes the error.
err = norm_infinity(newx - currentx)/norm_infinity(newx)
# Logs solution information
print(f"k={k:d} x=[{', '.join([f'{x:.6f}' for x in newx])}] err={err:.6f}")
# Updates the current
currentx = newx
# Updates the iteration counter
k = k + 1
return currentx
3.7.1. Example 1#
A = np.array([[3, 2, 1], [1, 6, 2], [1, 2, 4]])
b = np.array([1, 2, 1])
# Initial guess
x0 = [0.5, 2, 1]
xJ = jacobi_iteration_method(A, b, x0, 10**-4, 50)
print("Approximate solution: ", xJ)
k=1 x=[-1.333333, -0.083333, -0.875000] err=1.562500
k=2 x=[0.680556, 0.847222, 0.625000] err=2.377049
k=3 x=[-0.439815, 0.011574, -0.343750] err=2.547368
k=4 x=[0.440201, 0.521219, 0.354167] err=1.688379
k=5 x=[-0.132202, 0.141911, -0.120660] err=4.033530
k=6 x=[0.278946, 0.395587, 0.212095] err=1.039336
k=7 x=[-0.001090, 0.216144, -0.017530] err=1.295596
k=8 x=[0.195081, 0.339358, 0.142200] err=0.578062
k=9 x=[0.059694, 0.253420, 0.031551] err=0.534237
k=10 x=[0.153870, 0.312867, 0.108367] err=0.301008
k=11 x=[0.088633, 0.271566, 0.055099] err=0.240225
k=12 x=[0.133923, 0.300195, 0.092059] err=0.150869
k=13 x=[0.102517, 0.280327, 0.066422] err=0.112033
k=14 x=[0.124308, 0.294107, 0.084207] err=0.074093
k=15 x=[0.109193, 0.284546, 0.071870] err=0.053120
k=16 x=[0.119679, 0.291178, 0.080429] err=0.036013
k=17 x=[0.112405, 0.286577, 0.074491] err=0.025383
k=18 x=[0.117451, 0.289769, 0.078610] err=0.017415
k=19 x=[0.113951, 0.287555, 0.075753] err=0.012174
k=20 x=[0.116379, 0.289091, 0.077735] err=0.008400
k=21 x=[0.114695, 0.288025, 0.076360] err=0.005849
k=22 x=[0.115863, 0.288764, 0.077314] err=0.004047
k=23 x=[0.115053, 0.288252, 0.076652] err=0.002813
k=24 x=[0.115615, 0.288607, 0.077111] err=0.001949
k=25 x=[0.115225, 0.288360, 0.076793] err=0.001353
k=26 x=[0.115495, 0.288532, 0.077014] err=0.000938
k=27 x=[0.115308, 0.288413, 0.076860] err=0.000651
k=28 x=[0.115438, 0.288495, 0.076967] err=0.000451
k=29 x=[0.115348, 0.288438, 0.076893] err=0.000313
k=30 x=[0.115410, 0.288478, 0.076944] err=0.000217
k=31 x=[0.115367, 0.288450, 0.076909] err=0.000151
k=32 x=[0.115397, 0.288469, 0.076933] err=0.000105
k=33 x=[0.115376, 0.288456, 0.076916] err=0.000073
Approximate solution: [0.11537604 0.28845612 0.07691608]
Let us compare it with the true solution.
L,U = LUfactorization(A)
y = forward_substitution(L,b)
x = backward_substitution(U,y)
x
array([0.11538462, 0.28846154, 0.07692308])
3.8. Exercises#
Exercise 3.1
Implement the Gaussian elimination algorithm in terms of row operations.
Exercise 3.2
Implement the Gaussian elimination algorithm in terms of matrix operations.
Exercise 3.3
Observe how the LUfactorization method fails if A[0,0]=0. Adjust the method by swapping rows whenever the pivot element is zero.
Exercise 3.4
Implement Crout’s factorization.
Exercise 3.5
Implement the Gauss-Seidel iteration method and solve the same example solved with the Jacobi iteration method.