El código propio para la multiplicación de matrices se ejecuta más lentamente que la multiplicación en bucle usando std::vector

Question

Estoy aprendiendo C++ además de aprendizaje automático, así que decidí usar la biblioteca Eigen para la multiplicación de matrices. Estaba entrenando un perceptrón para reconocer un dígito de la base de datos MNIST. Para la fase de entrenamiento establecí el número de ciclos de entrenamiento (o épocas) en T = 100.

La 'matriz de entrenamiento' es una matriz de 10000 x 785. El elemento cero de cada fila contiene la 'etiqueta' que identifica el dígito al que se asignan los datos de entrada (los 784 elementos restantes de la fila).

También hay un vector de 'pesos' de 784 x 1 que contiene los pesos para cada una de las 784 características. El vector de pesos se multiplicaría con cada vector de entrada (una fila de la matriz de entrenamiento que excluye el elemento cero) y se actualizaría en cada iteración, y esto sucedería T veces para cada una de las 10000 entradas.

Escribí el siguiente programa (que captura la esencia de lo que estoy haciendo), donde comparé el enfoque "vainilla" de multiplicar las filas de una matriz con el vector de peso (usando std::vector y loops) con lo que sentí que era lo mejor que pude hacer con un enfoque Eigen. No es realmente una multiplicación de una matriz con un vector, en realidad estoy cortando la fila de la matriz de entrenamiento y multiplicándola con el vector de peso.

La duración del período de entrenamiento para el enfoque std::vector fue de 160,662 ms y para el método Eigen fue generalmente de más de 10,000 ms.

Compilo el programa usando el siguiente comando:

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

Estoy usando un MacBook Pro de "mediados" de 2012 con macOS Catalina y un i5 de doble núcleo a 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;
}

¿Qué cambios debo hacer para obtener mejores tiempos de ejecución?

Answer 1

Puede que no sea la mejor solución, pero puedes probar:

• Dado que el orden de datos predeterminado de Eigen es Column-Major, puede dejar que su matriz de entrenamiento sea 785x10000, de modo que cada etiqueta de entrenamiento/par de datos sea contiguo en la memoria (también cambie la línea donde se calcula sum_wx_m).
• Use la versión de tamaño fijo de las operaciones de bloque, es decir, puede reemplazar m.block(i, 1, 1, 784) con m.block<1,784>(i, 1) (en orden inverso si ya ha cambiado su matriz de entrenamiento diseño o simplemente puede mapear la parte de datos de su matriz de entrenamiento y usar la referencia .col() [vea el ejemplo a continuación])

Aquí está su código modificado basado en estas ideas:

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

Compilé este código en mi Ubuntu Desktop con 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

Después de conversaciones con los usuarios J. Schultke y puhu, realicé los siguientes cambios en mi código:

1. Cambié todas las llamadas m.block(i, 1, 1, 784) a m.block<1, 784>(i, 1) , esto reduce el tiempo requerido para el bucle de matriz Eigen en un tercio . (sugerido por primera vez por J. Schultke)
2. He declarado que mi matriz m está almacenada en el orden RowMajor . Esto se debe a que, de forma predeterminada, las matrices Eigen se almacenan en el orden ColMajor (columna principal). Esto haría que cada entrada en una fila se almacenara de forma contigua. Así que ahora las llamadas m.block() , que utilizo para referirme a una porción de una fila en la matriz m, simplemente recuperarían toda la memoria a la vez, reduciendo el tiempo de la "matriz Eigen" por debajo del "std: :vector" tiempo. (sugerido por puhu)

Los tiempos de ejecución promedio ahora son

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

y el código modificado es:

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