MATLAB: How to reduce computation time of (sparse) matrix multiplication in Matlab
computation timematrix multiplication
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
Best Answer
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.
We understand what a matrix multiply means. :) But it sounds like you don't appreciate the use of sparse matrices in MATLAB. Just because your matrix has zero elements in it, does not make it a matrix stored in sparse form.
If your sparse matrix is indeed stored in sparse format, then MATLAB will AUTOMATICALLY use highly efficient multiplication.
A = sprand(1000,1000,0.005);
B = sprand(1000,1000,0.005);
Af = full(A);
Bf = full(B);
whos A B Af Bf
Name SizeBytesClassAttributes
A 1000x100087656doublesparse
Af 1000x10008000000double
B 1000x100087832doublesparse
Bf 1000x10008000000double
But sparse matrices are not only there to save space. MATLAB does use the known sparsity in a multiplication.
timeit(@() Af*Bf)
ans =
0.078021563737
timeit(@() A*B)
ans =
0.000990737737
If you want something faster than the already fast multiplication built into MATLAB that works with sparse matrices, then, no. Not gonna happen.
You can define anything you want to be sparse if you so desire. So sparse(ones(1000)) will produce a sparse matrix. ;-)
But seriously, this matrix really is reasonably sparse. Just read what is shown on the page you direct to!
We see that it is shown to be an 1899x1899 square matrix. There are 20296 non-zeros in the matrix.
20296/1899
ans =
10.688
So on average, roughly 10.7 non-zeros per row of a matrix, where each row has 1899 elements.
20296/1899/1899
ans =
0.0056281
Total density of around 0.6% non-zeros. Is that matrix sparse? Well, yes, it is. Not massively so, but it is sparse. Since I don't have that matrix, I'll do a little memory comparison on a random one with the same density.
As = sprand(1899,1899,0.0056);
Af = full(As);
whos As Af
Name SizeBytesClassAttributes
Af 1899x189928849608double
As 1899x1899337568doublesparse
So we see that storing the matrix as sparse gives me a decrese in memory of almost a factor of 100-1. (Sparse matrix storage also requires we store the location of those elements, so it is not quite as efficient as we might like.)
b = rand(1899,1);
timeit(@() Af*b)
ans =
0.0023482
timeit(@() As*b)
ans =
6.1858e-05
And working with the matrix with a matrix multiply does give a significant savings.
Much of the time, when you work with a sparse matrix, you will perform a decomposition of some sort. Odds are that will not produce a speed increase here, because the matrix itself is not a nicely tightly banded matrix. I think much of the time, people think of a sparse matrix in the form of something like a tridiagonal matrix. A matrix factorization of a tridiagonal matrix will produce no fill-in at all. On the matrix you show, if we tried to do a Cholesky or LU factoriation, for example, the result will be a virtually completely full triangular matrix.
tic,[L,U] = lu(Af);toc
Elapsed time is 0.205769 seconds.
tic,[L,U] = lu(As);toc
Elapsed time is 3.032499 seconds.
As you see, not all computations on such a matrix will see a gain in speed. Here, we see the result ends up as a nearly full triangular matrix for the factors. So treating that matrix as sparse, IF you intend to do a decomposition of the matrix will be a time loss, because the overhead from the sparse algorithms is too costly.
nnz(L)
ans =
1177842
nnz(U)
ans =
1208969
numel(L)
ans =
3606201
That is to be expected on this matrix. It does not mean the original was not sparse. It is indeed sparse. Just perhaps not as massively sparse as you may want a matrix to be. Sometimes sparse must be measured by the eye of the beholder, in proper context.
Best Answer