WS05: Sparse systems and pivoting

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.

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.0
z = np.zeros(dim)

for k in range(nonzero):
    # insert code to compute z=A*y
    continue

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):
        """
        Insert code to find the diagonal array D_diag of matrix A_real
        """
        if I_row[k] == I_col[k]:
            True
    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.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
    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)
Contents