множения больших двоичных матриц (10Kx20K) я обычно преобразую матрицы в числа с плавающей запятой и выполняю умножение матрицы с плавающей запятой, так как умножение целочисленной матрицы выполняется довольно медленно (посмотрите здесь).

На этот раз мне нужно выполнить более ста тысяч этих умножений идаже улучшение производительности в миллисекундах в среднем имеет для меня значение.

Я хочуint или жеfloat матрица в результате, потому что продукт может иметь элементы, которые не равны 0 или 1. Все входные матричные элементы равны 0 или 1, поэтому они могут храниться как отдельные биты.

Во внутреннем произведении между вектором строки и вектором столбца (для получения одного элемента выходной матрицы) умножение упрощается до побитового И. Добавление по-прежнему является дополнением, но мы можем добавить биты с помощью функции подсчета населения вместо того, чтобы циклически обрабатывать их по отдельности.

Некоторые другие функции булевой / двоичной матрицы ИЛИ биты вместо того, чтобы считать их, производя результат битовой матрицы, но это не то, что мне нужно.

Вот пример кода, показывающий, что формирование проблемы какstd::bitset, AND а такжеcount операции быстрее, чем умножение матриц.

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

Запуск этого с:

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

Я получил:

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

Как вы можете видеть, matmul как набор битовых операций примерно в 3 раза быстрее, чем Eigen float matmul, что имеет смысл.

Хочу подчеркнуть, что мне нужно выполнить эту операцию более 100К (в HPC или облаке) и повышение производительности в миллисекундах в среднем будет иметь значение.

Я не привязан к какой-либо конкретной библиотеке, стандарту C ++ и т. Д. Поэтому, пожалуйста, не стесняйтесь отвечать любым решением, которое, по вашему мнению, работает быстрее, чем те, которые используют GPU, поскольку я не могу использовать его по ряду причин.