Para multiplicar matrices binarias grandes (10Kx20K), lo que suelo hacer es convertir las matrices en flotantes y realizar una multiplicación de matriz flotante ya que la multiplicación de matriz entera es bastante lenta (mira aquí)

Esta vez, sin embargo, necesitaría realizar más de cien mil de estas multiplicaciones yincluso una mejora de rendimiento de milisegundos en promedio me importa.

Quiero unint ofloat como resultado, porque el producto puede tener elementos que no son 0 o 1. Los elementos de la matriz de entrada son todos 0 o 1, por lo que pueden almacenarse como bits individuales.

En el producto interno entre un vector de fila y un vector de columna (para producir un elemento de la matriz de salida), la multiplicación se simplifica a Y bit a bit. La suma sigue siendo la suma, pero podemos agregar bits con una función de conteo de población en lugar de recorrerlos individualmente.

Algunas otras funciones booleanas / de matriz binaria O los bits en lugar de contarlos, producen un resultado de matriz de bits, pero eso no es lo que necesito.

Aquí hay un código de muestra que muestra que formando el problema comostd::bitset, AND ycount Las operaciones son más rápidas que la multiplicación de matrices.

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

Ejecutando esto con:

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

Yo obtengo:

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 puede ver, el matmul como conjunto de operaciones de conjunto de bits es aproximadamente 3 veces más rápido que el matmul flotante de Eigen, lo cual tiene sentido.

Quiero enfatizar que necesito realizar esta operación a más de 100K (en HPC o en la nube) y una mejora de rendimiento de milisegundos en promedio marcaría la diferencia.

No estoy obligado a ninguna biblioteca específica, estándar C ++, etc. Por lo tanto, siéntase libre de responder con cualquier solución que considere más rápida que las que usan GPU, ya que no puedo usarla por varias razones.