WS05: Sparse systems and pivoting - partial solutions#
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.
These are partial solutions. Please create Issues and Pull requests with your solutions.
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{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}. \)$
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.0
z = np.zeros(dim)
for k in range(nonzero):
z[I_row[k]] = z[I_row[k]] + A_real[k] * y[I_col[k]]
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
"""
# we should take care to ensure that arrays are stored with the correct type - float!
A = A.astype(np.float64)
b = b.astype(np.float64)
# check sizes of A and b match appropriately
nb = len(b)
n, m = A.shape
if n != m:
raise ValueError(f"A is not a square matrix! {A.shape=}")
if n != nb:
raise ValueError(f"shapes of A and b do not match! {A.shape=} {b.shape=}")
# check A is upper triangular
for i in range(n):
for j in range(0, i):
if not np.isclose(A[i, j], 0.0):
raise ValueError(f"A is not upper triangular! {A[i, j]=}")
# checks whether A has zero diagonal element
for i in range(n):
if np.isclose(A[i, i], 0.0):
raise ValueError(f"A[{i}, {i}] is zero")
# create a new array to store the results
x = np.empty_like(b)
# perform backwards subsitution
x[n - 1] = b[n - 1] / A[n - 1, n - 1]
for i in range(2, n + 1):
x[n - i] = b[n - i] / A[n - i, n - i]
for j in range(n - i + 1, n):
x[n - i] = x[n - i] - A[n - i, j] * x[j] / A[n - i, n - i]
return x
def Gaussian_elimination_pivoting(A, b, verbose=False):
# To ensure that arrays are stored in double precision.
A = A.astype(np.float64)
b = b.astype(np.float64)
# size of solution vector / the square matrix A
n = len(b) # or n, n = A.shape
# check sizes of A and b match appropriately
nb = len(b)
n, m = A.shape
if n != m:
raise ValueError(f"A is not a square matrix! {A.shape=}")
if n != nb:
raise ValueError(f"shapes of A and b do not match! {A.shape=} {b.shape=}")
if verbose:
print("starting system\n", A, b)
# perform forward elimination
for i in range(n):
# find the index of the maximal value in column i on or below
# the diagonal of A
maximum = abs(A[i, i])
max_index = i
for j in range(i + 1, n):
if abs(A[j, i]) > maximum:
maximum = abs(A[j, i])
max_index = j
# swap two max_indexs: i and max_index[i]
temp = b[i]
b[i] = b[max_index]
b[max_index] = temp
for j in range(n):
temp = A[i, j]
A[i, j] = A[max_index, j]
A[max_index, j] = temp
# check diagonal
if np.isclose(A[i, i], 0.0):
raise ValueError(f"A has zero on diagonal! A[{i}, {i}] = 0")
# row j <- row j - (a_{ji} / a_{ii}) row i
for j in range(i + 1, n):
if verbose:
print(f"row {j} <- row {j} - {A[j, i] / A[i, i]} row {i}")
factor = A[j, i] / A[i, i]
for k in range(0, n):
A[j, k] = A[j, k] - factor * A[i, k]
b[j] = b[j] - factor * b[i]
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)
# x0 = np.linalg.inv(A).dot(b)
print("Solution by numpy solver:", x0)
x = Gaussian_elimination_pivoting(A, b)
print("Gaussian elimination: ", x)
print("Residual: ", np.matmul(A, x) - b)
# extra test to ensure we have solved the problem
np.testing.assert_almost_equal(np.linalg.norm(b - A @ x), 0.0)
n = 20
A = np.random.rand(n, n)
b = np.random.rand(n)
# numpy linear solvers
x0 = np.linalg.solve(A, b)
# x0 = np.linalg.inv(A).dot(b)
print("Solution by numpy solver:", x0)
x = Gaussian_elimination_pivoting(A, b)
print("Gaussian elimination: ", x)
print("Residual: ", np.matmul(A, x) - b)
# extra test to ensure we have solved the problem
np.testing.assert_almost_equal(np.linalg.norm(b - A @ x), 0.0)
Part c: Extension#
5. implement Jacobi iteration for spare 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.
"""
n = len(v) # length of vector
sum2 = 0.0
for i in range(n):
sum2 += v[i] * v[i]
return np.sqrt(sum2)
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
"""
D_diag[I_row[k]] = A_real[k]
n_diag += 1
if n_diag < n:
print("Zero diagonal element!")
return
P_real = np.zeros(nonzero) # matrix P = D^{-1}(L+U)
p = np.zeros(n) # vector p = D^{-1} b
i = 0
for k in range(nonzero):
"""
Insert code to compute matrix P
"""
if I_col[k] != I_row[k]: # remove diagonal
P_real[k] = A_real[k] / D_diag[I_row[k]]
for i in range(n):
p[i] = b[i] / D_diag[i]
res = 2.0 * 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
x = np.zeros_like(b)
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
"""
Ax[I_row[k]] = Ax[I_row[k]] + A_real[k] * x[I_col[k]]
res = norm(Ax - b) / norm(b)
print(f"{res=:.8e}")
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")
np.testing.assert_almost_equal(np.linalg.norm(b - A @ x), 0.0, decimal=5)
np.testing.assert_allclose(x, x0, rtol=1.0e-4)
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)