MATLAB: Is there a faster way to multiply a sparse and full matrix than standard multiplication in Matlab

matrix manipulation

I have a sparse matrix (say A) which has at most 3 to 4 non-zero entries regardless of matrix size. I want to multiply it with a full matrix (Say B) of the same dimension as matrix A.

C = A*B

One way to interpret this product is that each row in C is actually a linear combination of rows of B and the coefficients of the linear combinations comes from rows of A. Since the rows of A has only few (3-4) non-zero entries so it means we only need 3 to 4 rows of B matrix to produce each row of Matrix C. If we can somehow avoid the operations performed by the zero entries and perform the multiplication for many rows in parallel then can we save a implementation time? If yes, I want to know is there any command or any way we can implement this in Matlab?

Thanks,

Best Answer

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         Size                Bytes  Class     Attributes
A         1000x1000              87656  double    sparse    
Af        1000x1000            8000000  double              
B         1000x1000              87832  double    sparse    
Bf        1000x1000            8000000  double

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.

Related Solutions

MATLAB: Why tis matrix is sparse

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         Size                 Bytes  Class     Attributes
  Af        1899x1899            28849608  double              
  As        1899x1899              337568  double    sparse

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.

spy(As)

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

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.

Best Answer

Related Solutions

MATLAB: Why tis matrix is sparse

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

Related Question