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

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,

Answer 1

[John D'Errico]

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.

[O45]

Dear John D'Errico, I appreciate your answer but in my case only one of the matrix is sparse and the other one is full. What I wanted to explain in my original post was that since matrix A is sparse so only a couple of rows are needed to generate each row of matrix C (I know that you know about matrix multiplication :) ).

[Guillaume]

And if you test John's example keeping one matrix always full, you'll see that matlab uses the knowledge that one is sparse to accelerate the muliplication. On my machine:

timeit(@() Af*Bf)    %both matrices full
ans =
   0.042679430433814

timeit(@() A*Bf)     %one sparse, one full
ans =
   0.008262039958923

In other words, matlab already does the right thing, as long as your sparse matrix is actually stored as sparse.

[John D'Errico]

So what? MATLAB understands how to take the product of a sparse and full matrix, and does so as efficiently as possible already!

Why not test these things yourself? It is trivial to do, and you get an answer immediately. I showed you how to do exactly that test.

A = sprand(1000,1000,0.005);
B = rand(1000,1000);
Af = full(A);

timeit(@() Af*B)
ans =
            0.051683486737

timeit(@() A*B)
ans =
            0.007581789737

So the product of a sparse and a full matrix is also faster than the product of two full matrices. MATLAB IS using an efficient algorithm, as long as your matrix is stored in sparse form. If you did not store the matrix as sparse, then no, MATLAB cannot possibly bother to look for zero elements inside every matrix multiply. That would make multiplies so much slower in general. It only works if you explicitly store your matrix as sparse.

Do you get as much of again from two fully, truly sparse matrices in the product. No. But then what can you expect?

Of course, if your matrix is not very sparse, then there is little gain. In fact, sparse methods can actually be slower when applied to nearly full matrices. But you have claimed that your matrices are indeed sparse. I don't know if that is true, since while you state that only up to 4 elements of each row are non-zero, I don't know how long the rows will be. Therefore your matrix may or may not fall into the realm of what I would consider sparse. As an example of this, consider two nearly full matrices.

A = sprand(1000,1000,0.5);
B = sprand(1000,1000,0.5);
timeit(@() A*B)
ans =
            0.345950141737

Af = full(A);
Bf = full(B);
timeit(@() Af*Bf)
ans =
            0.051654236737

So sparse matrices do indeed need to be truly sparse to benefit. Otherwise there is additional bookwork, as the sparse mtimes needs to "look for" the necessary elements in a multiply. That bookwork can overwhelm the actual matrix multiply.

One other thing. you need to consider that MATLAB is able to multi-thread matrix multiplies, but only when the matrices get sufficiently large that it sees a benefit. On small problems, the overhead of multithreading is simply not worth the effort.

For example, I executed the following test case on my CPU, while watching an activity monitor to tell me how many of my CPUS were being used:

A = rand(5000,5000);
B = rand(5000,5000);
C = A*B;

That multiply briefly woke up 4 of the CPUS on my machine.

Next, I tried it on a pair of big matrices, stored in sparse form. I explicitly made them not very sparse, to ensure that I would see if multiple cores were being used.

A = sprand(5000,5000,.5);
B = sprand(5000,5000,.5);
C = A*B;

That multiply never accessed more than one core. In fact, I believe that sparse multiplications NEVER go multi-threaded (at least in the current release of MATLAB. That may change in the future.)

The point is that a direct matrix multiply between full matrices can be very efficient on a machine with multiple CPUs, because on large problems it can use the full power of your machine. On small problems, there simply is not sufficient gain to try. (I have gotten as many as 8 cores, all running flat out, when computing millions of digits of pi.)

Remember that matrix multiples between full matrices also have the ability to fully use any pipelines on your processors. The low level routines can fill the pipelines very well then. But on sparse multiplies, your CPUs spend part of their time just lining up elements to then multiply and add.

So there are many issues to consider. But the fact is, A*B is as fast as you will get it to be, as long as you use sparsity properly.

[O45]

Thank you for your detailed response.

[D_coder]

what if I want to do element by element multiplication rather than product of two matrices , sparse performs slower than actual multiplication of full matrices

A = sprand(1000,1000,0.005);
B = rand(1000,1000);
Af = full(A);
timeit(@() Af.*B)

ans =

1.8185e-04

   timeit(@() A.*B)

   ans =

         0.0012