For multiplying large binary matrices (10Kx20K), what I usually to do is to convert the matrices to float ones and perform float matrix multiplication as integer matrix multiplication is pretty slow (have a look at here).
This time though, I'd need to perform over hundred thousands of these multiplications and even a millisecond performance improvement on average matters to me.
I want an int or float matrix as a result, because the product may have elements that aren't 0 or 1. The input matrix elements are all 0 or 1, so they can be stored as single bits.
In the inner-product between a row vector and a column vector (to produce one element of the output matrix), multiplication simplifies to bitwise AND. Addition is still addition, but we can add bits with a population-count function instead of looping over them individually.
Some other boolean/binary-matrix functions OR the bits instead of counting them, producing a bit-matrix result, but that's not what I need.
Here is a sample code showing that forming the problem as std::bitset, AND and count operations is faster than matrix multiplication.
#include <iostream>
using std::cout; using std::endl;
#include <vector>
using std::vector;
#include <chrono>
#include <Eigen/Dense>
using Eigen::Map; using Eigen::Matrix; using Eigen::MatrixXf;
#include <random>
using std::random_device; using std::mt19937; using std::uniform_int_distribution;
#include <bitset>
using std::bitset;
using std::floor;
const int NROW = 1000;
const int NCOL = 20000;
const float DENSITY = 0.4;
const float DENOMINATOR = 10.0 - (10*DENSITY);
void fill_random(vector<float>& vec) {
random_device rd;
mt19937 eng(rd());
uniform_int_distribution<> distr(0, 10);
int nnz = 0;
for (int i = 0; i < NROW*NCOL; ++i)
vec.push_back(floor(distr(eng)/DENOMINATOR));
}
void matmul(vector<float>& vec){
float *p = vec.data();
MatrixXf A = Eigen::Map<Eigen::Matrix<float, NROW, NCOL, Eigen::RowMajor>>(p);
cout << "Eigen matrix has " << A.rows() << " rows and " << A.cols() << " columns." << endl;
cout << "Total non-zero values : " << A.sum() << endl;
cout << "The density of non-zero values is " << A.sum() * 1.0 / (A.cols()*A.rows()) << endl;
auto start = std::chrono::steady_clock::now();
MatrixXf B = A.transpose() * A;
auto end = (std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::steady_clock::now() - start)).count();
cout << "Mat mul took " << end << " ms"<< endl;
// Just to make sure the operation is not skipped by compiler
cout << "Eigen coo ";
for (int i=0; i<10; ++i)
cout << B(0,i) << " ";
cout << endl;
}
void bitset_op(vector<float>& vec) {
// yeah it's not a great idea to set size at compile time but have to
vector<bitset<NROW>> col_major(NCOL);
// right, multiple par for isn't a good idea, maybe in a parallel block
// Doing this for simplicity to profile second loop timing
// converting row major float vec to col major bool vec
#pragma omp parallel for
for (int j=0; j < NCOL; ++j) {
for (int i=0; i < NROW; ++i) {
col_major[j].set(i, vec[i*NCOL + j] && 1);
}
}
auto start = std::chrono::steady_clock::now();
vector<int> coo;
coo.assign(NCOL*NCOL, 0);
#pragma omp parallel for
for (int j=0; j < NCOL; ++j) {
for (int k=0; k<NCOL; ++k) {
coo[j*NCOL + k] = (col_major[j]&col_major[k]).count();
}
}
auto end = (std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::steady_clock::now() - start)).count();
cout << "bitset intersection took " << end << " ms"<< endl;
// Just to make sure the operation is not skipped by compiler
cout << "biset coo ";
for (int i=0; i<10; ++i)
cout << coo[i] << " ";
cout << endl;
}
int main() {
// Saving to float instead of int to speed up matmul
vector<float> vec;
fill_random(vec);
matmul(vec);
bitset_op(vec);
}
Running this with:
g++ -O3 -fopenmp -march=native -I. -std=c++11 code.cpp -o code
I get:
Eigen matrix has 1000 rows and 20000 columns.
Total non-zero values : 9.08978e+06
The density of non-zero values is 0.454489
Mat mul took 1849 ms
Eigen coo 458 206 208 201 224 205 204 199 217 210
bitset intersection took 602 ms
biset coo 458 206 208 201 224 205 204 199 217 210
As you can see, matmul as set of bitset operations is about 3x faster than Eigen's float matmul, which makes sense.
I want to emphasize that I need to perform this operation over 100K (in the HPC or cloud) and a millisecond performance improvement on average would make a difference.
I'm not bound to any specific library, C++ standard, etc. So please feel free to answer with any solution that you think is faster other than those using GPU, as I can't use it for a number of reasons.