# Sparse and Banded Matrices

## Topics in this lab

## Introduction

Thus far we have been thinking of matrices as tables (or arrays) of
numbers. However, matrices that arise from many real-world applications
are typically extremely large and very **sparse**, which means that
most entries in the matrix are zero.

`clear all`

Back to the top

## Example : Adjacency matrices

Consider the connectivity of the vertices in the graph below

nv = 12; P = haltonset(2); V = net(P,nv); DT = DelaunayTri(V); h = triplot(DT,'k-','Marker','.','MarkerSize',6); axis square; axis off text(V(:,1),V(:,2),cellstr(num2str((1:nv)')),'FontSize',14,'Color','r');

We can represent the connectivity of the graph with an adjacency matrix $A$. Here entry $a_{ij}$ is 1 if there is a line (edge) from vertex $i$ to vertex $j$, otherwise the entry is zero. For the graph above the adjacency matrix is given by

A = zeros(nv); for j=1:size(DT,1) A(DT(j,1),DT(j,:)) = 1; A(DT(j,2),DT(j,:)) = 1; A(DT(j,3),DT(j,:)) = 1; end A

A = 1 0 0 0 1 0 1 0 1 1 0 0 0 1 1 1 0 1 1 1 0 1 1 0 0 1 1 0 1 1 0 0 1 0 1 0 0 1 0 1 0 0 0 1 0 1 0 0 1 0 1 0 1 0 1 0 1 0 1 0 0 1 1 0 0 1 0 1 1 0 0 1 1 1 0 0 1 0 1 0 0 1 1 0 0 1 0 1 0 1 0 1 0 0 0 1 1 0 1 0 1 1 0 0 1 0 0 0 1 1 0 1 0 0 1 0 0 1 0 0 0 1 1 0 1 0 1 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 1

For a sparse matrix it is extremely wasteful (or impossible if the matrix is too large) to store every entry of the matrix. Indeed all we really care about is the nonzero entries and their locations in the matrix. For this purpose, MATLAB has a special data type called sparse. The sparse representation for the above adjacency matrix is

sparse(A)

ans = (1,1) 1 (5,1) 1 (7,1) 1 (9,1) 1 (10,1) 1 (2,2) 1 (3,2) 1 (4,2) 1 (6,2) 1 (7,2) 1 (8,2) 1 (10,2) 1 (11,2) 1 (2,3) 1 (3,3) 1 (5,3) 1 (6,3) 1 (9,3) 1 (11,3) 1 (2,4) 1 (4,4) 1 (8,4) 1 (10,4) 1 (1,5) 1 (3,5) 1 (5,5) 1 (7,5) 1 (9,5) 1 (11,5) 1 (2,6) 1 (3,6) 1 (6,6) 1 (8,6) 1 (9,6) 1 (12,6) 1 (1,7) 1 (2,7) 1 (5,7) 1 (7,7) 1 (10,7) 1 (11,7) 1 (2,8) 1 (4,8) 1 (6,8) 1 (8,8) 1 (12,8) 1 (1,9) 1 (3,9) 1 (5,9) 1 (6,9) 1 (9,9) 1 (1,10) 1 (2,10) 1 (4,10) 1 (7,10) 1 (10,10) 1 (2,11) 1 (3,11) 1 (5,11) 1 (7,11) 1 (11,11) 1 (6,12) 1 (8,12) 1 (12,12) 1

The full
function is used to convert a matrix in sparse format back to
the regular table of numbers format.

While one can use the sparse
function to *sparsify* a full matrix, this is
usually not the right way to create sparse matrices. Instead MATLAB
provides many alternatives.

Back to the top

## Creating sparse matrices

The most basic way to create a sparse matrix is to simply give the sparse function the location of the non-zero entries and the values of these entries. Consider the following example, which creates non-zero entries on the diagonal and first superdiagonal:

A = sparse([1 1 2 2 3 3 4 4 5],[1 2 2 3 3 4 4 5 5],[1 2 3 4 5 6 7 8 9])

A = (1,1) 1 (1,2) 2 (2,2) 3 (2,3) 4 (3,3) 5 (3,4) 6 (4,4) 7 (4,5) 8 (5,5) 9

In **full** form this matrix is given by

full(A)

ans = 1 2 0 0 0 0 3 4 0 0 0 0 5 6 0 0 0 0 7 8 0 0 0 0 9

For constructing sparse, **banded** matrices MATLAB provides the useful
command spdiags.
This function is similar to the diag function in
that it lets you specify the entries that are to go on one of the
diagonals of the sparse matrix.

The example below, shows how to construct a 10-by-10 sparse matrix that has the
constant 1 on the diagonal, -2 on the second subdiagonal, and 3 on the
third superdiagonal.

e = ones(10,1)*[-2 1 3]; A = spdiags(e,[-2 0 3],10,10);

Here is the matrix is **full** form:

full(A)

ans = 1 0 0 3 0 0 0 0 0 0 0 1 0 0 3 0 0 0 0 0 -2 0 1 0 0 3 0 0 0 0 0 -2 0 1 0 0 3 0 0 0 0 0 -2 0 1 0 0 3 0 0 0 0 0 -2 0 1 0 0 3 0 0 0 0 0 -2 0 1 0 0 3 0 0 0 0 0 -2 0 1 0 0 0 0 0 0 0 0 -2 0 1 0 0 0 0 0 0 0 0 -2 0 1

Back to the top

## Properties of sparse matrices.

To determine how sparse a matrix is we simply compute the fraction of
entries of the matrix that are zeros, which we call the
**sparsity ratio**. The nnz
function of MATLAB makes this calculation easy. It returns the
number of non-zero entries in a sparse matrix.

For the sparse matrix from the previous example the sparsity ratio is

sparsity_ratio = 1 - nnz(A)/numel(A)

sparsity_ratio = 0.750000000000000

We consider a matrix to be sparse if its sparsity ratio is close to one.
A dense matrix has a sparsity ratio near zero.

The bandwidth of a matrix is defined to be the
maximum distance of the non-zero elements from the main diagonal. We can
compute the bandwidth using the find and max functions.

For the sparse matrix from the previous example the bandwidth is

[i,j] = find(A); bandwidth = max(abs(i-j))

bandwidth = 3

which matches what we expect since we put non-zero entries on the third
subdiagonal.

A matrix is said to *banded* if its bandwidth is small. For matrices
with small bandwidth and high sparsity, Gaussian elimination can be performed
in far fewer operations than the standard
$N^3$.

Back to the top

## Visualizing sparse matrices

Many times we want to see the structure of a sparse matrix with out
looking at the individual entries. We can use the
spy command in MATLAB to do this.

For example, we can see the structure of the sparse matrix from the
previous examples using

spy(A)

The code below generates a sparse matrix by randomly filling in non-zero entries

ii = randi(100,[500 1]); jj = randi(100,[500 1]); A = sparse(ii,jj,rand(500,1),100,100);

We can view the structure of this matrix using

spy(A)

While this matrix is sparse, it is certainly not banded.

Here is one more neat example. It first plots a Buckminster Fuller
geodesic dome (Bucky Ball) and then shows the structure of the adjacency
matrix.

[B,V] = bucky; clf; [i,j] = find(B); [~, p] = sort(max(i,j)); i = i(p); j = j(p); X = [V(i,1) V(j,1)]'; Y = [V(i,2) V(j,2)]'; Z = [V(i,3) V(j,3)]'; X = [X; NaN(size(i))']; Y = [Y; NaN(size(i))']; Z = [Z; NaN(size(i))']; % Plot the vertices plot3(X,Y,Z,'k-','LineWidth',2); hold on; plot3(V(:,1),V(:,2),V(:,3),'b.','MarkerSize',40); % Add a sphere sphere(100); % Lighting and shading effects on the sphere colormap([0.4 0.4 0.4]) alpha(0.5); % transparency shading interp; lighting phong; camlight left; % Lighting % Adjust aspect ratio, view angle, add title. daspect([1 1 1]) view(42,10); axis off; hold off; title('Bucky Ball','fontsize',18); snapnow;

We can visualize the sparsity structure using the spy command.

spy(B); title('Adjacency matrix of the Bucky Ball','fontsize',18); snapnow;

Back to the top

## Operations on sparse matrices

Most operations that can be done on **full** matrices can also be done on
**sparse** matrices. For operations on sparse matrices like matrix-vector
multiply, amd matrix-matrix multiply, special algorithms are used that
exploit the non-zeros of the matrices. Most of the time, the sparse
versions of these operations will be many, many times faster than the
full-versions.

As mentioned above, solving linear systems involving sparse and banded
matrices will also be much faster than full matrices. There are no
special functions you have to call to do this, you only have to use the
backslash operator on a matrix in sparse format.
MATLAB will
automatically detect this and use the approrpriate algorithm.

The following example illustrates the difference in timing for sparse
matrix solve and a full matrix solve. We use the 2D Poisson matrix,
which arises when solving the Poisson equation in 2D dimensions with
finite-differences.

A = gallery('poisson',50,50); % Load in the Poisson matrix for 50x50 grid b = ones(size(A,1),1); % Create an artificial right hand side.

Here is the structure of the Poisson matrix:

spy(A)

Now lets see how long it takes to solve a linear system with the sparse version of this matrix and the full version.

tic x = A\b; % Sparse version tms = toc; A = full(A); % Convert the matrix to full tic x = A\b; tmf = toc; fprintf('Time for sparse matrix solve is %1.4f\n',tms); fprintf('Time for full matrix solve is %1.4f\n',tmf);

Time for sparse matrix solve is 0.0025 Time for full matrix solve is 0.1520

Clearly the sparse version is significantly faster.

Back to the top

## Exercises

- Load the wathen
matrix into MATLAB using the following command
A = gallery('wathen',40,40);

Determine the sparsity ratio and bandwidth of the matrix $A$ and produce a plot showing the matrix structure. - Use spdiags
to produce the following two 10-by-10 matrices
$$
A =
\left[
\begin{array}{rrrrrrrrrr}
-2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
1 & -2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & -2 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2
\end{array}
\right]
$$

$$ B = \left[ \begin{array}{rrrrrrrrrr} 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 \\ -1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & -1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -1 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & -1 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & -1 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 1 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 0 \end{array} \right] $$ - Load in the Poisson matrix for a 20-by-20 grid
using the command
A = gallery('poisson',20,20);

Make a series of spy plots of $A$, $A^2$, ..., $A^8$ to see how the number of nonzero entries increases with each new matrix multiplication. This phenomenon is called*fill-in*and occurs for many matrix operations. Also make a spy plot of the inv(A). What is the sparsity ratio and bandwidth of inv(A)?

Compare your answers with the solutions.

Back to the top

## Get the code

Do you want to try the above code fragments on your own? Download the Matlab script that produces this page here. (lab_27.m)

Published with MATLAB® 8.4