Sparse matrices

The numerical solution of differential equations often results in what is known as a ‘sparse linear’ system. A matrix is ‘sparse’ if most of its entries are zero, and most of the MATLAB functions you have used so far have versions that are specially optimised for sparse matrices, which can speed up your code immensely. The details of this speeding-up will be discussed further in the later unit ‘Software engineering’, but we introduce sparse matrices here as they will be used in exercises later in this unit.

The following will give you an introduction to the uses of sparse matrices.

Walkthrough

Create the matrix A using the following commands:

A=2*eye(5000)-diag(ones(4999,1),1)-diag(ones(4999,1),-1);

Run the following commands to investigate the structure of A:

A(1:5,1:5)
nnz(A)
spy(A)
spy(A(1:10,1:10))
Question

What is the structure of A?

Expand for solution
Solution

A is ‘tri-diagonal’, with 2 on the diagonal and −1 on the adjacent diagonals. Therefore, A has $5000+(2\times 4999)=14998$ non-zero entries, as found using nnz(A).

Calculate the inverse of A.

Note

The inverse of A is not sparse. This is true in general.

Ai=inv(A);
spy(abs(Ai)>1)
spy(abs(Ai)>5)
spy(abs(Ai)>10)
spy(abs(Ai)>100)

Calculate how long it takes to calculate $A^2$.

Info

You can use the commands tic and toc to record elapsed time in MATLAB.

tic; A*A; toc

Create the matrix B using the following commands. Note that A and B have the same non-zero entries:

B=spdiags(ones(5000,1)*[-1,2,-1],[-1,0,1],5000,5000);

The command spdiags creates a sparse matrix with the entries from the vector (passed as the first variable) on the diagonal indicated by the second variable. These variables can be passed individually or in groups as in the above example. The final two variables represent the size of the matrix. Look at the MATLAB help files for more information on spdiags.

Look at the structure of B:

B(1:5,1:5)
full(B(1:5,1:5))
nnz(B)
Question

Using the command whos, look at how the matrices A and B are stored. Which uses less memory?

Expand for solution
Solution
>> whos A
  Name         Size                  Bytes  Class     Attributes
  A         5000x5000            200000000  double              

>> whos B
  Name         Size               Bytes  Class     Attributes
  B         5000x5000            279976  double    sparse

So A used 2000000 bytes, and B only used 27976 bytes, saving about 98.6% of the space.

Question

Calculate how long it takes to calculate $B^2$:

tic; B*B; toc

Which is quicker, $A^2$ or $B^2$?

Expand for solution
Solution
>> tic; A*A; toc
Elapsed time is 0.768850 seconds.       

>> tic; B*B; toc
Elapsed time is 0.013921 seconds.

Calculating $B^2$ should be faster, and the proportion that it is faster should increase with the size of the matrix.

Question

Now compare the \ operator

b=ones(5000,1);
tic; x1=A\b; toc
tic; x2=B\b; toc

Which is quicker, A\b or B\b?

Expand for solution
Solution
>> b=ones(5000,1);
tic; x1=A\b; toc
tic; x2=B\b; toc
Elapsed time is 0.268863 seconds.
Elapsed time is 0.000719 seconds.

Calculating B\b should be faster, and the proportion that it is faster should increase with the size of the matrix.

Finally, you can also convert matrices from full to sparse using the following command

A=sparse(A);

Using the command whos, check that A is now sparse and uses the same memory as B.

Now clear the workspace

clear