# Fastest way to perform element-wise multiplication on a sparse matrix

I have two large-ish matrices (~100K cols x ~100K rows).

They are sparse and symmetrical (about 0.1% of them values are non-zero).

I want to do element-wise multiplication between them.

Also, I plan to perform this operation about 1,000,000 times, so speed is is definitely going to be an issue.

Does anyone know of a good library/software that could help me accomplish this?

Also, I've stumbled unto SuiteSparse. If anyone knows more about it, do you think I could use it for this?

Best.

• If all your matrices have the same sparsity pattern, then you can simply multiply the VAL array element-wise. There is probly no API function to do that, but it is very easy (just a for loop). In addition, you can "#pragma omp parallel for" the loop to gain a bit more. Oct 9 '16 at 11:26
• Do your matrices all have the same sparsity pattern? What form are they stored in (e.g. CRS or CCS)? Oct 9 '16 at 15:06

I suggest you take a look at the Eigen C++ matrix class library, http://eigen.tuxfamily.org/index.php?title=Main_Page.

I did a small test with sparse matrices of the size and sparsity you state and it takes about 1 ms per matrix multiplication on my moderate-power windows machine.

The code for my experiment is below. As you can see, most of the code is for setting up the test matrices. The actual matrix multiply is a simple one-liner. (I multiplied the matrices 500 times to make the timing more reliable.)

#include <iostream>
using std::cout;
using std::endl;
#include <random>

#include <Eigen/Core>
#include <Eigen/Sparse>

#include <boost/timer.hpp>

int main()
{
const int n = 100000;
const int nnz = static_cast<int>(.001*n);
typedef Eigen::Triplet<double> Triplet;
std::vector<Triplet> at(nnz), bt(nnz);
std::random_device rd; // obtain a random number from hardware
std::mt19937 eng(rd()); // seed the generator
std::uniform_int_distribution<> distr(0, nnz-1); // define the range
for (int i = 0; i < nnz; i++) {
int ii = distr(eng), jj = distr(eng);
at[i] = Triplet(ii, jj, 1);
bt[i] = Triplet(jj, ii, 3);
}
Eigen::SparseMatrix<double> as(n, n), bs(n, n), cs(n, n);
as.setFromTriplets(at.begin(), at.end());
bs.setFromTriplets(bt.begin(), bt.end());
boost::timer timer;
const int nt = 500;
for (int i = 0; i < nt; i++)
cs = as.cwiseProduct(bs);
cout << "elapsed time=" << timer.elapsed()/(double) nt  << endl;
}

• Thank you, I ran the code as well and it's fast enough for my purposes. Just a small remark, were you compiling with optimization enabled? My code ran about 10 times slower than yours without optimization -O0, but when using -O3 it's actually down to 0.5ms / step. I'm asking to compare apples w/ apples. Oct 9 '16 at 21:13
• Yes, I was running with optimization turned on; my computer does have SSE2 but is actually rather slow compared with current machines. In general, Eigen is much slower in debug mode because the steps they take to optimize the code, including supporting vector instructions sets like SSE, is very good. Oct 9 '16 at 21:37
• By the way, you may have to explicitly enable vector instructions with the appropriate compiler flag if you want Eigen to exploit those. Oct 9 '16 at 21:44
• I explicitly added -msse3 and it became a wee faster hehe, thanks. Oct 9 '16 at 23:59

If there is a choice in programming language, one option can be to use Julia, which has built in support for sparse matrices (via Suitsparse). The timing come out to about one and a half milliseconds on my laptop, and you get to use an interactive, dynamic language, which may be useful in certain situations.

julia> a=sprand(1000, 1000, 0.1)
1000x1000 sparse matrix with 99749 Float64 entries:
[5   ,    1]  =  0.725824
[17  ,    1]  =  0.420022
[19  ,    1]  =  0.404282
[21  ,    1]  =  0.0307138
[52  ,    1]  =  0.453376
[55  ,    1]  =  0.30054
[69  ,    1]  =  0.360203
[74  ,    1]  =  0.346881
[94  ,    1]  =  0.312849
⋮
[932 , 1000]  =  0.978966
[933 , 1000]  =  0.149551
[954 , 1000]  =  0.417852
[959 , 1000]  =  0.722707
[964 , 1000]  =  0.519931
[967 , 1000]  =  0.567152
[971 , 1000]  =  0.964192
[979 , 1000]  =  0.88494
[987 , 1000]  =  0.286723
[988 , 1000]  =  0.24282

julia> b=sprand(1000, 1000, 0.1)
1000x1000 sparse matrix with 99998 Float64 entries:
[1   ,    1]  =  0.920533
[3   ,    1]  =  0.879179
[7   ,    1]  =  0.267203
[25  ,    1]  =  0.522407
[34  ,    1]  =  0.656031
[41  ,    1]  =  0.280885
[44  ,    1]  =  0.735824
[68  ,    1]  =  0.433098
[69  ,    1]  =  0.124862
⋮
[932 , 1000]  =  0.505959
[939 , 1000]  =  0.983413
[947 , 1000]  =  0.418157
[949 , 1000]  =  0.884657
[963 , 1000]  =  0.412645
[964 , 1000]  =  0.544348
[966 , 1000]  =  0.709398
[983 , 1000]  =  0.260483
[989 , 1000]  =  0.1218
[1000, 1000]  =  0.468975

julia> a.*b
1000x1000 sparse matrix with 9876 Float64 entries:
[69  ,    1]  =  0.0449757
[102 ,    1]  =  0.0340867
[137 ,    1]  =  0.0794594
[247 ,    1]  =  0.108002
[376 ,    1]  =  0.248346
[609 ,    1]  =  0.241789
[633 ,    1]  =  0.224115
[658 ,    1]  =  0.379804
[754 ,    1]  =  0.272618
⋮
[224 , 1000]  =  0.0122434
[301 , 1000]  =  0.163899
[309 , 1000]  =  0.0972784
[403 , 1000]  =  0.0245688
[659 , 1000]  =  0.0801249
[700 , 1000]  =  0.158823
[760 , 1000]  =  0.388442
[926 , 1000]  =  0.193808
[932 , 1000]  =  0.495317
[964 , 1000]  =  0.283024

julia> @time a.*b
0.001649 seconds (25 allocations: 3.056 MB)