Código Eigen para multiplicação de matrizes rodando mais devagar que a multiplicação em loop usando std::vector

Question

Estou aprendendo C++ e aprendizado de máquina, então decidi usar a biblioteca Eigen para multiplicação de matrizes. Eu estava treinando um perceptron para reconhecer um dígito do banco de dados MNIST. Para a fase de treinamento, defino o número de ciclos de treinamento (ou épocas) para T = 100.

A 'matriz de treinamento' é uma matriz de 10000 x 785. O elemento zero de cada linha contém o 'rótulo' que identifica o dígito para o qual os dados de entrada (os 784 elementos restantes da linha) são mapeados.

Há também um vetor de 'pesos' de 784 x 1 que contém os pesos para cada um dos 784 recursos. O vetor de pesos seria multiplicado com cada vetor de entrada (uma linha da matriz de treinamento excluindo o elemento zero) e seria atualizado a cada iteração, e isso aconteceria T vezes para cada uma das 10.000 entradas.

Escrevi o seguinte programa (que captura a essência do que estou fazendo), onde comparei a abordagem "baunilha" de multiplicar as linhas de uma matriz com o vetor de peso (usando std::vector e loops) com o que senti ser o melhor que pude fazer com uma abordagem Eigen. Não é realmente uma multiplicação de uma matriz com um vetor, na verdade estou fatiando a linha da matriz de treinamento e multiplicando isso com o vetor de peso.

A duração do período de treinamento para a abordagem std::vector foi de 160,662 ms e para o método Eigen foi geralmente superior a 10.000 ms.

Eu compilo o programa usando o seguinte comando:

clang++ -Wall -Wextra -pedantic -O3 -march=native -Xpreprocessor -fopenmp permute.cc -o perm -std=c++17

Estou usando um MacBook Pro "meados" de 2012 executando o macOS Catalina e tendo um i5 dual core de 2,5 GHz.

#include <iostream>
#include <algorithm>
#include <random>
#include <Eigen/Dense>
#include <ctime>
#include <chrono>
using namespace Eigen;

int main() {
    Matrix<uint8_t, Dynamic, Dynamic> m = Matrix<uint8_t, Dynamic, Dynamic>::Random(10000, 785);
    Matrix<double, 784, 1> weights_m = Matrix<double, 784, 1>::Random(784, 1);
    Matrix<uint8_t, 10000, 1> y_m, t_m;

    std::minstd_rand rng;
    rng.seed(time(NULL));
    std::uniform_int_distribution<> dist(0,1); //random integers between 0 and 1
    for (int i = 0; i < y_m.rows(); i++) {
        y_m(i) = dist(rng);
        t_m(i) = dist(rng);
    }

    int T = 100;
    int err;
    double eta;
    eta = 0.25; //learning rate
    Matrix<double, 1, 1> sum_wx_m;

    auto start1 = std::chrono::steady_clock::now(); //start of Eigen Matrix loop

    for (int iter = 0; iter < T; iter++) {
        for (int i = 0; i < m.rows(); i++) {
            sum_wx_m = m.block(i, 1, 1, 784).cast<double>() * weights_m;
        
            //some code to update y_m(i) based on the value of sum_wx_m which I left out
        
            err = y_m(i) - t_m(i);
            if (fabs(err) > 0) { //update the weights_m matrix if there's a difference between target and predicted
                weights_m = weights_m - eta * err * m.block(i, 1, 1, 784).transpose().cast<double>();
            } 
        }
    }

    auto end1 = std::chrono::steady_clock::now();
    auto diff1 = end1 - start1;
    std::cout << "Eigen matrix time is "<<std::chrono::duration <double, std::milli> (diff1).count() << " ms" << std::endl;

    //checking how std::vector form performs;

    std::vector<std::vector<uint8_t>> v(10000);
    std::vector<double> weights_v(784);
    std::vector<uint8_t> y_v(10000), t_v(10000);

    for (unsigned long i = 0; i < v.size(); i++) {
        for (int j = 0; j < m.cols(); j++) {
            v[i].push_back(m(i, j));
        }
    }

    for (unsigned long i = 0; i < weights_v.size(); i++) {
        weights_v[i] = weights_m(i);
    }

    for (unsigned long i = 0; i < y_v.size(); i++) {
        y_v[i] = dist(rng);
        t_v[i] = dist(rng);
    }

    double sum_wx_v;

    auto start2 = std::chrono::steady_clock::now(); //start of vector loop

    for (int iter = 0; iter < T; iter++) {
        for(unsigned long j = 0; j < v.size(); j++) {
            sum_wx_v = 0.0;
            for (unsigned long k = 1; k < v[0].size() ; k++) {
                sum_wx_v += weights_v[k - 1] * v[j][k];
            }
        
            //some code to update y_v[i] based on the value of sum_wx_v which I left out
        
            err = y_v[j] - t_v[j];
            if (fabs(err) > 0) {//update the weights_v matrix if there's a difference between target and predicted
                for (unsigned long k = 1; k < v[0].size(); k++) {
                    weights_v[k - 1] -= eta * err * v[j][k];
                }
            }
        }
    }

    auto end2 = std::chrono::steady_clock::now();
    auto diff2 = end2 - start2;
    std::cout << "std::vector time is "<<std::chrono::duration <double, std::milli> (diff2).count() << " ms" << std::endl;
}

Que alterações devo fazer para obter melhores tempos de funcionamento?

Answer 1

Pode não ser a melhor solução, mas você pode tentar:

• Como a ordem de dados padrão de Eigen é Column-Major, você pode deixar sua matriz de treinamento ser 785x10000, de modo que cada rótulo/par de dados de treinamento seja contíguo na memória (altere também a linha em que sum_wx_m é calculado).
• Use a versão de tamanho fixo das operações de bloco, ou seja, você pode substituir m.block(i, 1, 1, 784) por m.block<1,784>(i, 1) (na ordem inversa se você já trocou sua matriz de treinamento layout ou você pode simplesmente mapear a parte de dados de sua matriz de treinamento e usar a referência .col() [veja o exemplo abaixo])

Aqui está o seu código modificado com base nessas ideias:

#include <iostream>
#include <algorithm>
#include <random>
#include <Eigen/Dense>
#include <ctime>
#include <chrono>
using namespace Eigen;

int main() {
    Matrix<uint8_t, Dynamic, Dynamic> m = Matrix<uint8_t, Dynamic, Dynamic>::Random(785, 10000);
    Map<Matrix<uint8_t, Dynamic, Dynamic>> m_data(m.data() + 785, 784, 10000);

    Matrix<double, 784, 1> weights_m = Matrix<double, 784, 1>::Random(784, 1);
    Matrix<uint8_t, 10000, 1> y_m, t_m;

    std::minstd_rand rng;
    rng.seed(time(NULL));
    std::uniform_int_distribution<> dist(0,1); //random integers between 0 and 1
    for (int i = 0; i < y_m.rows(); i++) {
        y_m(i) = dist(rng);
        t_m(i) = dist(rng);
    }

    int T = 100;
    int err;
    double eta;
    eta = 0.25; //learning rate
     Matrix<double, 1, 1> sum_wx_m;

    auto start1 = std::chrono::steady_clock::now(); //start of Eigen Matrix loop

    for (int iter = 0; iter < T; iter++) {
        for (int i = 0; i < m.cols(); i++) {
            sum_wx_m = weights_m.transpose() * m_data.col(i).cast<double>();
        
            //some code to update y_m(i) based on the value of sum_wx_m which I left out
        
            err = y_m(i) - t_m(i);
            if (fabs(err) > 0) { //update the weights_m matrix if there's a difference between target and predicted
                weights_m = weights_m - eta * err * m_data.col(i).cast<double>();
            } 
        }
    }

    auto end1 = std::chrono::steady_clock::now();
    auto diff1 = end1 - start1;
    std::cout << "Eigen matrix time is "<<std::chrono::duration <double, std::milli> (diff1).count() << " ms" << std::endl;

    //checking how std::vector form performs;

    std::vector<std::vector<uint8_t>> v(10000);
    std::vector<double> weights_v(784);
    std::vector<uint8_t> y_v(10000), t_v(10000);

    for (unsigned long i = 0; i < v.size(); i++) {
        for (int j = 0; j < m.rows(); j++) {
            v[i].push_back(m(j, i));
        }
    }

    for (unsigned long i = 0; i < weights_v.size(); i++) {
        weights_v[i] = weights_m(i);
    }

    for (unsigned long i = 0; i < y_v.size(); i++) {
        y_v[i] = dist(rng);
        t_v[i] = dist(rng);
    }

    double sum_wx_v;

    auto start2 = std::chrono::steady_clock::now(); //start of vector loop

    for (int iter = 0; iter < T; iter++) {
        for(unsigned long j = 0; j < v.size(); j++) {
            sum_wx_v = 0.0;
            for (unsigned long k = 1; k < v[0].size() ; k++) {
                sum_wx_v += weights_v[k - 1] * v[j][k];
            }
        
            //some code to update y_v[i] based on the value of sum_wx_v which I left out
        
            err = y_v[j] - t_v[j];
            if (fabs(err) > 0) {//update the weights_v matrix if there's a difference between target and predicted
                for (unsigned long k = 1; k < v[0].size(); k++) {
                    weights_v[k - 1] -= eta * err * v[j][k];
                }
            }
        }
    }

    auto end2 = std::chrono::steady_clock::now();
    auto diff2 = end2 - start2;
    std::cout << "std::vector time is "<<std::chrono::duration <double, std::milli> (diff2).count() << " ms" << std::endl;
}

Eu compilei este código no meu Ubuntu Desktop com i7-9700K:

g++ -Wall -Wextra -O3 -std=c++17
====================================
Eigen matrix time is 110.523 ms
std::vector time is 117.826 ms


g++ -Wall -Wextra -O3 -march=native -std=c++17
=============================================
Eigen matrix time is 66.3044 ms
std::vector time is 71.2296 ms

Answer 2

Após discussões com os usuários J. Schultke e puhu, fiz as seguintes alterações em meu código:

1. Eu mudei todas as chamadas m.block(i, 1, 1, 784) para m.block<1, 784>(i, 1) , isso reduz o tempo necessário para o loop da matriz Eigen em um terço . (sugerido pela primeira vez por J. Schultke)
2. Eu declarei minha matriz m como sendo armazenada na ordem RowMajor . Isso ocorre porque, por padrão, as matrizes Eigen são armazenadas na ordem ColMajor (coluna principal). Isso faria com que cada entrada em uma linha fosse armazenada de forma contígua. Então agora as chamadas m.block() , que eu uso para se referir a uma fatia de uma linha na matriz m, simplesmente buscariam todo o pedaço de memória de uma vez, reduzindo o tempo da "matriz própria" para abaixo do "std: :vetor" tempo. (sugerido por puhu)

Os tempos de execução médios agora são

cpp:Pro$ ./perm
Eigen matrix time is 134.76 ms
std::vector time is 155.574 ms

e o código modificado é:

#include <iostream>
#include <algorithm>
#include <random>
#include <Eigen/Dense>
#include <chrono>
#include <ctime>
using namespace Eigen;
int main() {
    Matrix<uint8_t, Dynamic, Dynamic, RowMajor> m = Matrix<uint8_t, Dynamic, Dynamic, RowMajor>::Random(10000, 785);
    Matrix<double, 784, 1> weights_m = Matrix<double, 784, 1>::Random(784, 1);
    Matrix<uint8_t, 10000, 1> y_m, t_m;
    std::minstd_rand rng;
    rng.seed(time(NULL));
    std::uniform_int_distribution<> dist(0,1); //random integers between 0 and 1
    for (int i = 0; i < y_m.rows(); i++) {
        y_m(i) = dist(rng);
        t_m(i) = dist(rng);
    }

    int T = 100;
    int err;
    double eta;
    eta = 0.25; //learning rate
    Matrix<double, 1, 1> sum_wx_m;

    auto start1 = std::chrono::steady_clock::now(); //start of Eigen Matrix loop

    for (int iter = 0; iter < T; iter++) {
        for (int i = 0; i < m.rows(); i++) {
            auto b = m.block<1, 784>(i, 1).cast<double>();
            sum_wx_m = b * weights_m;
    
            //some code to update y_m(i) based on the value of sum_wx_m which I left out
    
            err = y_m(i) - t_m(i);
            if (fabs(err) > 0) { //update the weights_m matrix if there's a difference between target and predicted
                weights_m = weights_m - eta * err * b.transpose();
            } 
        }
    }

    auto end1 = std::chrono::steady_clock::now();
    auto diff1 = end1 - start1;
    std::cout << "Eigen matrix time is "<<std::chrono::duration <double, std::milli> (diff1).count() << " ms" << std::endl;

    //checking how std::vector form performs;

    std::vector<std::vector<uint8_t>> v(10000);
    std::vector<double> weights_v(784);
    std::vector<uint8_t> y_v(10000), t_v(10000);

    for (unsigned long i = 0; i < v.size(); i++) {
        for (int j = 0; j < m.cols(); j++) {
            v[i].push_back(m(i, j));
        }
    }

    for (unsigned long i = 0; i < weights_v.size(); i++) {
        weights_v[i] = weights_m(i);
    }

    for (unsigned long i = 0; i < y_v.size(); i++) {
        y_v[i] = dist(rng);
        t_v[i] = dist(rng);
    } 

    double sum_wx_v;

    auto start2 = std::chrono::steady_clock::now(); //start of vector loop

    for (int iter = 0; iter < T; iter++) {
        for(unsigned long j = 0; j < v.size(); j++) {
            sum_wx_v = 0.0;
            for (unsigned long k = 1; k < v[0].size() ; k++) {
                sum_wx_v += weights_v[k - 1] * v[j][k];
            }
    
            //some code to update y_v[i] based on the value of sum_wx_v which I left out
    
            err = y_v[j] - t_v[j];
            if (fabs(err) > 0) {//update the weights_v matrix if there's a difference between target and predicted
                for (unsigned long k = 1; k < v[0].size(); k++) {
                    weights_v[k - 1] -= eta * err * v[j][k];
                }
            }
        }
    }

    auto end2 = std::chrono::steady_clock::now();
    auto diff2 = end2 - start2;
    std::cout << "std::vector time is "<<std::chrono::duration <double, std::milli> (diff2).count() << " ms" << std::endl;
}