Lecture 08: Sparse matrices and stopping criteria

Lecture 08: Sparse matrices and stopping criteria#

Sparse Matrices#

There are two main ways in which sparse matrices can be exploited in order to obtain benefits within iterative methods.

The storage must be reduced for \(O(n^2)\).
The cost per iteration must be reduced from \(O(n^2)\).

Recall that a sparse matrix is defined to be such that it has at most \(\alpha n\) non-zero entries (where \(\alpha\) is independent of \(n\)).

Typically this happens when we know there are at most \(\alpha\) non-zero entries in any row.

Sparse matrix storage#

The simplest way in which a sparse matrix is stored is using three arrays:

an array of floating point numbers (A_real say) that stores the non-zero entries;
an array of integers (I_row say) that stores the row number of the corresponding entry in the real array;
an array of integers (I_col say) that stores the column numbers of the corresponding entry in the real array.

This requires just \(3 \alpha n\) units of storage - i.e. \(O(n)\).

Sparse matrix multiplication#

Given the above storage pattern, the following algorithm will execute a sparse matrix-vector multiplication (\(\vec{z} = A \vec{y}\)) in \(O(n)\) operations:

z = np.zeros((n, 1))
for k in range(nonzero):
    z[I_row[k]] = z[I_row[k]] + A_real[k] * y[I_col[k]]

Here nonzero is the number of non-zero entries in the matrix.
Note that the cost of this operation is \(O(n)\) as required.

Convergence of an iterative method#

We have discussed the construction of iterations which aim to find the solution of the equations \(A \vec{x} = \vec{b}\) through a sequence of better and better approximations \(\vec{x}^{(k)}\).

In general the iteration takes the form

\[ \vec{x}^{(k+1)} = \vec{F}(\vec{x}^{(k)}) \]

where \(\vec{x}^{(k)}\) is a vector of values and \(\vec{F}\) is some vector-valued function which we have defined.
How can we decide if this iteration has converged?
- We need \(\vec{x} - \vec{x}^{(k)}\) to be small.

Magnitude of a vector#

How do we decide that a vector/array is small?

We need a way of defining the size/magnitude of an array.
The most common measure is to use the “Euclidean norm” of an array.
This is defined to be the square root of the sum of squares of the entries of the array:

\[ \| \vec{r} \| = \sqrt{ \sum_{i=1}^n r_i^2 } \]

where \(\vec{r}\) is a vector with \(n\) entries.

Examples#

Consider the following sequence \(\vec{x}^{(k)}\):

\[\begin{split} \begin{pmatrix} 1 \\ -1 \end{pmatrix}, \begin{pmatrix} 1.5 \\ 0.5 \end{pmatrix}, \begin{pmatrix} 1.75 \\ 0.25 \end{pmatrix}, \begin{pmatrix} 1.875 \\ 0.125 \end{pmatrix}, \begin{pmatrix} 1.9375 \\ -0.0625 \end{pmatrix}, \begin{pmatrix} 1.96875 \\ -0.03125 \end{pmatrix}, \ldots \end{split}\]

What is \(\|\vec{x}^{(1)} - \vec{x}^{(0)}\|\)?
What is \(\|\vec{x}^{(5)} - \vec{x}^{(4)}\|\)?

Let \(\vec{x} = \begin{pmatrix} 2 \\ 0 \end{pmatrix}\).

What is \(\|\vec{x} - \vec{x}^{(3)}\|\)?
What is \(\|\vec{x} - \vec{x}^{(4)}\|\)?
What is \(\|\vec{x} - \vec{x}^{(5)}\|\)?

Convergence detection#

Rather than decide in advance how many iterations (of the Jacobi or Gauss-Seidel methods) to use stopping criteria:

This could be a maximum number of iterations.
This could be the change in values is small enough:

\[ \|x^{(k+1)} - \vec{x}^{(k)}\| < tol, \]
This could be the norm of the residual is small enough:

\[ \| \vec{r} \| = \| \vec{b} - A \vec{x}^{(k)} \| < tol \]
In both cases, we call \(tol\) the convergence tolerance.
The choice of \(tol\) will control the accuracy of the solution.

Discussion#

What is a good convergence tolerance?

Failure to converge#

In general there are two possible reasons that an iteration may fail to converge.

It may diverge - this means that \(\|\vec{x}^{(k)}\| \to \infty\) as \(k\) (the number of iterations) increases, e.g.:

\[\begin{split} \begin{pmatrix} 1 \\ 1 \end{pmatrix}, \begin{pmatrix} 4 \\ 2 \end{pmatrix}, \begin{pmatrix} 16 \\ 4 \end{pmatrix}, \begin{pmatrix} 64 \\ 8 \end{pmatrix}, \begin{pmatrix} 256 \\ 16 \end{pmatrix}, \begin{pmatrix} 1024 \\ 32 \end{pmatrix}, \ldots \end{split}\]
It may neither converge nor diverge, e.g.:

\[\begin{split} \begin{pmatrix} 1 \\ 1 \end{pmatrix}, \begin{pmatrix} 2 \\ 0 \end{pmatrix}, \begin{pmatrix} 3 \\ 1 \end{pmatrix}, \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \begin{pmatrix} 2 \\ 1 \end{pmatrix}, \begin{pmatrix} 3 \\ 0 \end{pmatrix}, \ldots \end{split}\]

Implementation#

In addition to testing for convergence it is also necessary to include tests for failure to converge.

Divergence may be detected by monitoring \(\|\vec{x}^{(k)}\|\).
Impose a maximum number of iterations to ensure that the loop is not repeated forever!