Para multiplicar matrizes binárias grandes (10Kx20K), o que costumo fazer é converter as matrizes em unidades flutuantes e executar a multiplicação da matriz flutuante, pois a multiplicação da matriz inteira é bem lenta (dê uma olhada aqui)

Dessa vez, porém, eu precisaria realizar mais de cem milhares dessas multiplicações emesmo uma melhoria de desempenho de milissegundos em média é importante para mim.

Eu quero umint oufloat como resultado, porque o produto pode ter elementos que não são 0 ou 1. Os elementos da matriz de entrada são todos 0 ou 1, para que possam ser armazenados como bits únicos.

No produto interno entre um vetor de linha e um vetor de coluna (para produzir um elemento da matriz de saída), a multiplicação simplifica para AND bit a bit. Adição ainda é adição, mas podemos adicionar bits com uma função de contagem de população em vez de fazer um loop sobre eles individualmente.

Algumas outras funções de matriz booleana / binária OU os bits, em vez de contá-los, produzindo um resultado de matriz de bits, mas não é disso que eu preciso.

Aqui está um código de exemplo mostrando que a formação do problema comostd::bitset, AND ecount operações é mais rápida que a multiplicação de matrizes.

#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);
}

Executando isso com:

g++ -O3 -fopenmp -march=native -I. -std=c++11 code.cpp -o code

Eu recebo:

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

Como você pode ver, matmul como conjunto de operações de bits é cerca de 3x mais rápido que o float matigen de Eigen, o que faz sentido.

Quero enfatizar que preciso executar esta operação acima de 100K (no HPC ou na nuvem) e uma melhoria de desempenho de milissegundos em média faria a diferença.

Não estou vinculado a nenhuma biblioteca específica, padrão C ++, etc. Portanto, sinta-se à vontade para responder com qualquer solução que considere mais rápida do que aquelas que usam GPU, pois não posso usá-lo por vários motivos.