How to reduce computation time of (sparse) matrix multiplication in Matlab?

In my Matlab script, I have one line of code which uses the chain rule to calculate the derivative of a vector A (or an array) with respect to another vector B, it's something like:

dA/dB = dP/dB * dQ/dP * dA/dQ.

The 3 terms on the right hand side of the equation are all sparse (non-square non-diagonal) matrices, and in each matrix, the number of rows and of columns are both about 4000. This line of code does give correct derivatives, but the matrix multiplication is too slow! It's so slow that I can't really use this code..

In Matlab Workspace, the 3 matrices are shown to be "sparse double". Is it strange that sparse matrix multiplication takes so long in Matlab?

I would appreciate any suggestion, thank you!

-Paula

[Walter Roberson]

How sparse are they individually? nnz()/numel()

[Steven Lord]

Walter's question is a good one, since storing a matrix that's not sparsely populated in sparse storage format can actually take up more memory than storing it using full storage!

A = magic(100);
S = sparse(A);
fprintf("S has %g%% non-zero elements", 100*nnz(S)/numel(S))

S has 100% non-zero elements

whos

  Name        Size              Bytes  Class     Attributes

  A         100x100             80000  double              
  S         100x100            160808  double    sparse    

But if the matrix is truly sparsely populated you can save storage and potentially computation time:

S2 = bucky;
A2 = full(S2);
fprintf("S2 has %g%% non-zero elements", 100*nnz(S2)/numel(S2))

S2 has 5% non-zero elements

whos A2 S2

  Name       Size            Bytes  Class     Attributes

  A2        60x60            28800  double              
  S2        60x60             3368  double    sparse    

[Paula]

Hi Walter and Steven,

dA/dB = dP/dB * dQ/dP * dA/dQ.

nnz(dP/dB) = nnz(dQ/dP) = 14023800, nnz(dA/dQ) = 14700.

dP/dB: 3816 by 3675, (which has 14023800 components, which means this matrix are all non-zero entries!)

dQ/dP: 3675 by 3816, (same, all non-zero entries)

dA/dQ: 3816 by 3675.

By doing this calculation I just realized that matrices dP/dB and dQ/dP don't really have zero entries.. So I guess this is probably the reason why the multiplication is so slow even if they are defined to be "sparse"?

I'll define the 3 matrices to be full matrices and try again, as @John D'Errico explained.

There are a lot of other "intermediate" matrices involved in this code, it's taking me some time to identify which matrices should become "full".

Thanks!

Answer 1

[John D'Errico]

The answer is easy. It depends!

What does it depend on? How sparse are your matrices? A lot of people think if half the entries in a matrix are zero, then using sparse form is right. WRONG.You don't really gain much there. In fact, it probably runs MUCH more slowly. As a test on a relatvely small pair of matrices...

As = sprand(1000,1000,0.5); % 50% non-zero, so 500 non-zeros per row
Af = full(As);
Bs = sprand(1000,1000,0.5);
Bf = full(Bs);

timeit(@() As*Bs)

ans = 0.2835

timeit(@() Af*Bf)

ans = 0.0289

And that was only for a 1000x1000 pair of matrices.

A significant problem is that sparse matrix multiplication does not use multiple processors. On my own computer, I have 8 real cores. But that sparse multiply uses only ONE of them, whereas the full matrix multiply uses all 8 in parallel.

Now, let me try it again, but for larger matrices that are truly sparse.

As = sprand(5000,5000,0.001); % 0.1% non-zero, roughly 5 non-zeros per row. 
Af = full(As);
Bs = sprand(5000,5000,0.001);
Bf = full(Bs);

timeit(@() As*Bs)

ans = 0.0039

timeit(@() Af*Bf)

ans = 3.4139

As you see in the second test, the sparse multiply whizzes past the full multiply. We are talking tortoise and hare here.

So the real problem is you most likely do not have matrices that are even remotely sparse. Hey, I have a lot of zeros. So sparse must be good. Just having a lot of zeros is not enough. A matrix needs to be seriously sparse for you to see a gain. But when that is the case, sparse is a HUGE benefit.

You may want to do some reading about sparse matrices and the concept of fill-in to understand why all of this works as it does.

For example, consider the matrices:

A1 = sprand(1000,1000,0.5);
B1 = sprand(1000,1000,0.5);
[numel(A1),nnz(A1),nnz(B1),nnz(A1*B1)]

ans = 1×4

1000000 393668 393681 1000000

A2 = sprand(5000,5000,0.001);
B2 = sprand(5000,5000,0.001);
[numel(A2),nnz(A2),nnz(B2),nnz(A2*B2)]

ans = 1×4

25000000 24985 24989 124732

Do you see that the product A1*B1 is a completely full matrix? So even though A1 and B1 were sparse in theory, the product A1*B1 is completely full. It is still stored in sparse form of course. And if you then multiply that matrix by ANOTHER sparse matrix, the product will be even slower to compute.

[Paula]

Hi John,

Thanks a lot for your explanation! I just realized by doing nnz() that I don't really have zeros in two of the matrices.. Clearly I didn't understand how to use "sparse matrix".

I'll read on sparse matrix, and I'll define the 3 matrices on the right hand side of the equation to be full matrices, and let you know.

[Paula]

I modified the code, and it worked! So for my matrix multiplication: P = A * B * C, matrices A and B have only non-zero entries, matrix C is sparse enough: nnz(C)/numel(C) = 0.001.

So I changed this line to: P = full(A) * full(B) * C, and made corresponding changes to the other matrices, and the elapsed time of this line of code reduced remarkably (1/60 of the previous time).

I learned a lot from this question, thank you so much for your help, @Walter Roberson @Steven Lord @John D'Errico !

-Paula

[John D'Errico]

Super. A question where you learn something is a good question!

Used properly, sparse matrices are an IMMENSELY valuable tool. I recall figuring out on one problem, before we had sparse matrices, that it would take days to solve a problem we needed to solve on the biggest, fastest computer we had available. With sparse matrices, it just happens in seconds. But your matrices need to be truly sparse to show a real gain.