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.
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;
}
-O0, but when using -O3 it's actually down to 0.5ms / step. I'm asking to compare apples w/ apples.
$\endgroup$
Commented
Oct 9, 2016 at 21:13
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)
