WS05: Sparse systems and pivoting#
These exercises are indented to give you practice at using the material on numerical approximation and are intended to reinforce the material that was covered in lectures.
Please attempt the worksheet before your tutorial. Support is available in your tutorial or in the Class Team.
Partial solutions are available to check your understanding. Please create Issues and Pull requests with your solutions or share your answers for peer feedback in the class team Worksheet solutions channel
Part a (pen and paper warm up)#
1. Sparse matrix storage#
Write down the three arrays: array of floating point numbers (say A_real) to store the non-zero entries, array of integers (say I_row) to store the row number and array of integers (say I_col) to store the column number.
\[\begin{split}
\begin{pmatrix}
3 & 0 & 1.5 & 0 \\ 2 & 4 & 2.2 & 0 \\ 0 & 2 & 6 & 0 \\ 0 & 0 & 4 & -9
\end{pmatrix}
\end{split}\]
2. Gaussian elimination with pivoting#
Using GE with pivoting to solve the following linear system of equations
\[\begin{split}
\begin{pmatrix}
0 & 1 & 0 \\
1 & 1.001 & 1 \\
100 & 100 & 0 \\
\end{pmatrix}
\begin{pmatrix}
x_1 \\
x_2 \\
x_3
\end{pmatrix} =
\begin{pmatrix}
0 \\ 1002 \\ 200
\end{pmatrix}.
\end{split}\]
3. Sparse matrix storage and access#
To implement a direct solver such as Gaussian elimination, is it convenient to store a matrix as sparse matrix? What’s about iterative solver: Jacobi iteration or Gauss-Seidel iteration?
Part b (code implementations and testing)#
4. Sparse matrix multiplication#
import numpy as np
A_real = np.array([3, 1.5, 2,4 ,2.2, 2, 6, 4, -9], dtype=np.float64)
I_row = np.array([0, 0, 1, 1, 1, 2, 2, 3, 3], dtype=np.int32)
I_col = np.array([0, 2, 0, 1, 2, 1, 2, 2, 3], dtype=np.int32)
nonzero = len(A_real)
dim = 4
y = np.zeros(dim) + 1.
z = np.zeros(dim)
for k in range(nonzero):
# insert code to compute z=A*y
print(z)
5. implementation of Gaussian elimination with pivoting#
import numpy as np
# copy from ws03
def upper_triangular_solve(A, b):
"""
Solve the system A x = b where A is assumed to be lower triangular,
i.e. A(i,j) = 0 for j > i, and the diagonal is assumed to be nonzero,
i.e. A(i,i) != 0.
The code checks that A is lower triangular and converts A and b to
double precision before computing.
ARGUMENTS: A lower triangular n x n array
b right hand side column n-vector
RETURNS: x column n-vector solution
"""
"""
Copy your previous implementation here
"""
def Gaussian_elimination_pivoting(A, b):
# To ensure that arrays are stored in double precision.
A = A.astype(np.float64)
b = b.astype(np.float64)
# size of solution vector
n = len(b)
for i in range(n):
# find the index of the maximal value in array in column i below A[i,i]
maximum = abs(A[i, i])
for j in range(i + 1, n):
if abs(A[j, i]) > maximum:
maximum = abs(A[j, i])
"""
swap row i and row max_index
"""
if np.abs(A[i, i]) < 1.0e-15:
print("A is singular!")
return
"""
Apply Gaussian elimination to both matrix A and right-hand side vector b
A becomes an upper triangular matrix at the end of this loop
"""
# solve using previous worksheet
return upper_triangular_solve(A, b)
6. Test your implementation#
A = np.array(
[[-10, 2, 0, 67], [-2, 50, -77, 1.0e-5], [1, 7, 30, 8], [-10, -7, 0.001, 80]]
)
b = np.array([1, 2, 9, 0])
# numpy linear solvers
x0 = np.linalg.solve(A, b)
print("Solution by numpy solver:", x0)
x = Gaussian_elimination_pivoting(A, b)
print("Gaussian elimination: ", x)
print("Residual: ", np.matmul(A, x) - b)
n = 20
A = np.random.rand(n, n)
b = np.random.rand(n)
# numpy linear solvers
x0 = np.linalg.solve(A, b)
print("Solution by numpy solver:", x0)
x = Gaussian_elimination_pivoting(A, b)
print("Gaussian elimination: ", x)
print("Residual: ", np.matmul(A, x) - b)
Part c: Extension#
5. implement Jacobi iteration for sparse matrix#
using the following stopping criterion
\[
\frac{\|\vec{b} - A\vec{x} \|}{\|\vec{b}\|} < 10^{-6}.
\]
def norm(v):
'''
Insert code to compute Euclidean norm.
'''
def Jacobi_iteration_sparse(A_real, I_row, I_col, b, tol):
# To ensure that arrays are stored in double precision.
A_real = A_real.astype(np.float64)
b = b.astype(np.float64)
nonzero = len(A_real) # number of nonzero elements
n=len(b) # dimensions of the linear system of equations
D_diag=np.zeros(n) # diagonal of matrix A
n_diag = 0 # count the dimensions
for k in range(nonzero):
if (I_row[k] == I_col[k]):
'''
Insert code to find the diagonal array D_diag of matrix A_real
'''
if n_diag < n :
print('Zero diagonal element!')
return
P_real=np.zeros(nonzero) # matrix P = D^{-1}(L+U)
i = 0
for k in range(nonzero):
'''
Insert code to compute matrix P
'''
res = 2.*tol # initialise a residual
max_it = 1000 # set a maximal number of iterations in case the programme comes to a dead loop
it = 0 # record the total number of iterations
while (res > tol and it < max_it):
it = it + 1
xnew = p.copy()
for k in range(nonzero):
xnew[I_row[k]] = xnew[I_row[k]] - P_real[k] * x[I_col[k]]
x = xnew.copy()
Ax = np.zeros(n)
for k in range(nonzero):
'''
Insert code to compute Ax = A*x
'''
res = norm(Ax-b) / norm(b)
print ("res ===", res)
return x, it
6. Test Jacobi_iteration_sparse#
# Test different linear solvers starting from the above two-dimensional linear system
A = np.array(
[[3, 0, 1.5, 0], [2, 4, 2.2, 0], [0, 2, 6, 0], [0, 0, 4, -9]], dtype=np.float64
)
b = np.array([3, -1, 9, 0], dtype=np.float64)
# numpy linear solver
x0 = np.linalg.solve(A, b)
print("Solution by numpy solver:", x0)
A_real = np.array([3, 1.5, 2, 4, 2.2, 2, 6, 4, -9], dtype=np.float64)
I_row = np.array([0, 0, 1, 1, 1, 2, 2, 3, 3], dtype=np.int32)
I_col = np.array([0, 2, 0, 1, 2, 1, 2, 2, 3], dtype=np.int32)
tol = 1.0e-6
x, it = Jacobi_iteration_sparse(A_real, I_row, I_col, b, tol)
print(f"Solution by Jacobi iteration: {x} after {it} iterations")
7. Compare Jacobi_iteration_sparse() with Gaussian_elimination_pivoting()#
Notice that there is a condition for Jacobi iteration: non zero element on the diagonal of the system matrix.
n = 20
B = np.random.rand(n, n)
eps = 10
A = eps * np.eye(n) + B * B.T
b = np.random.rand(n)
# create a spare matrix
n = len(b)
A_real = np.array([], dtype=np.float64)
I_row = np.array([], dtype=np.int32)
I_col = np.array([], dtype=np.int32)
for i in range(n):
for j in range(n):
if A[i, j] < 0.5 and i != j:
A[i, j] = 0.0
else:
if i == j:
A[i, j] = A[i, j] + 2.0
A_real = np.append(A_real, A[i, j])
I_row = np.append(I_row, i)
I_col = np.append(I_col, j)
x, it = Jacobi_iteration_sparse(A_real, I_row, I_col, b, 1.0e-6)
print(f"Solution by Jacobi iteration: {x} after {it} iterations")
x = Gaussian_elimination_pivoting(A, b)
print("Gaussian elimination: ", x)
print("Residual: ", np.matmul(A, x) - b)