Back to tutorial index

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

  1. 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.
  2. 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] $$
  3. 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

Download this file

You can download this code here lab_27.m

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)

Powered by MathJax

Published with MATLAB® 8.4