C++ Eigen LU schneller machen (meine Tests zeigen 2x langsamer als GSL)

Question

Ich vergleiche LU Dekomposition/Lösung von Eigen zu GSL, und finden Sie Eigen in der Größenordnung von 2x langsamer mit -O3 Optimierungen auf g ++/OSX. Ich isolierte das Timing des Zerlegens und des Lösens, aber finde, dass beide wesentlich langsamer sind als ihre GSL-Gegenstücke. Mache ich etwas Dummes, oder funktioniert Eigen für diesen Anwendungsfall nicht gut (z. B. sehr kleine Systeme?) Errichtet mit Eigen 3.2.8 und einer älteren GSL 1.15. Der Testfall ist sehr künstlich, spiegelt aber die Ergebnisse in einer nichtlinearen Anpassungssoftware wider, die ich schreibe - Eigen ist irgendwo von 1,5x - 2x + langsamer für die gesamte LU/Lösungsoperation.

#define NDEBUG 

#include "sys/time.h" 
#include "gsl/gsl_linalg.h" 
#include <Eigen/LU> 

// Ax=b is a 3x3 system for which soln is x=[8,2,3] 
// 
double avals_col[9] = { 4, 2, 2, 7, 5, 5, 7, 5, 9 }; 
    // col major 
double avals_row[9] = { 4, 7, 7, 2, 5, 5, 2, 5, 9 }; 
    // row major 
double bvals[9] = { 67, 41, 53 }; 

//----------- helpers 

void print_solution(double *x, int dim, char *which) { 
    printf("%s solve for x:\n", which); 
    for(int j=0; j<3; j++) { 
     printf("%g ", x[j]); 
    } 
    printf("\n"); 
} 

struct timeval tv; 
struct timezone tz; 
double timeNow() { 
    gettimeofday(&tv, &tz); 
    int _mils = tv.tv_usec/1000; 
    int _secs = tv.tv_sec; 
    return (double)_secs + ((double)_mils/1000.0); 
} 

//----------- 

void run_gsl(double *A, double *b, double *x, int dim, int reps) { 
    gsl_matrix_view gslA; 
    gsl_vector_view gslB; 
    gsl_vector_view gslX; 
    gsl_permutation *gslP; 
    int sign; 

    gslA = gsl_matrix_view_array(A, dim, dim); 
    gslP = gsl_permutation_alloc(dim); 
    gslB = gsl_vector_view_array(b, dim); 
    gslX = gsl_vector_view_array(x, dim); 

    int err; 
    double t, elapsed; 
    t = timeNow(); 
    for(int i=0; i<reps; i++) { 
     // gsl overwrites A during decompose, so we must copy the initial A each time. 
     memcpy(A, avals_row, sizeof(avals_row)); 
     err = gsl_linalg_LU_decomp(&gslA.matrix, gslP, &sign); 
    } 
    elapsed = timeNow() - t; 
    printf("GSL decompose (%dx) time = %g\n", reps, elapsed); 

    t = timeNow(); 
    for(int i=0; i<reps; i++) { 
     err = gsl_linalg_LU_solve(&gslA.matrix, gslP, &gslB.vector, &gslX.vector); 
    } 
    elapsed = timeNow() - t; 
    printf("GSL solve (%dx) time = %g\n", reps, elapsed); 

    gsl_permutation_free(gslP); 
} 

void run_eigen(double *A, double *b, double *x, int dim, int reps) { 
    Eigen::PartialPivLU<Eigen::MatrixXd> eigenA_lu; 

    Eigen::Map< Eigen::Matrix < double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > > ma(A, dim, dim); 
    Eigen::Map<Eigen::MatrixXd> mb(b, dim, 1); 

    int err; 
    double t, elapsed; 
    t = timeNow(); 
    for(int i=0; i<reps; i++) { 
     // This memcpy is not necessary for Eigen, which does not overwrite A in the 
     // decompose, but do it so that the time is directly comparable to GSL. 
     memcpy(A, avals_col, sizeof(avals_col)); 
     eigenA_lu.compute(ma); 
    } 
    elapsed = timeNow() - t; 
    printf("Eigen decompose (%dx) time = %g\n", reps, elapsed); 

    t = timeNow(); 
    Eigen::VectorXd _x; 
    for(int i=0; i<reps; i++) { 
     _x = eigenA_lu.solve(mb); 
    } 
    elapsed = timeNow() - t; 
    printf("Eigen solve (%dx) time = %g\n", reps, elapsed); 

    // copy soln to passed x 
    for(int i=0; i<dim; i++) { 
     x[i] = _x(i); 
    } 
} 

int main() { 
    // solve a 3x3 system many times 

    double A[9], b[3], x[3]; 
    int dim = 3; 
    int reps = 1000000; 

    memcpy(b, bvals, sizeof(bvals)); 
     // init b vector, A is copied multiple times in run_gsl/run_eigen 

    run_eigen(A, b, x, dim, reps); 
    print_solution(x, dim, "Eigen"); 

    run_gsl(A, b, x, dim, reps); 
    print_solution(x, dim, "GSL"); 

    return 0; 
} 

Dies erzeugt, zum Beispiel:

Eigen decompose (1000000x) time = 0.242 
Eigen solve (1000000x) time = 0.108 
Eigen solve for x: 
8 2 3 
GSL decompose (1000000x) time = 0.049 
GSL solve (1000000x) time = 0.075 
GSL solve for x: 
8 2 3 

Answer 1

Ihre Benchmark nicht wirklich fair ist, wie Sie die Kopie der Eingangsmatrix zweimal in der Eigen-Version tun: eine manuell durch memcpy, und ein innerhalb PartialPivLU. Sie lassen auch Eigen wissen, dass mb ein Vektor ist, indem Sie es als Map<Eigen::Vectord> deklarieren. Dann bekomme ich (GCC5, O3, Eigen3.3):

Eigen decompose (1000000x) time = 0.087 
Eigen solve (1000000x) time = 0.036 
Eigen solve for x: 
8 2 3 
GSL decompose (1000000x) time = 0.032 
GSL solve (1000000x) time = 0.062 
GSL solve for x: 
8 2 3 

Außerdem Eigens PartialPivLU ist nicht wirklich entworfen für eine solche extrem kleinen Matrizen (siehe unten). Für 3x3-Matrizen sollte die Inverse besser explizit berechnet werden (für Matrizen bis 4x4 ist es normalerweise ok, aber nicht für größere!). In diesem Fall müssen Sie die Größen zur Compile-Zeit beheben:

Eigen::PartialPivLU<Eigen::Matrix3d> eigenA_lu; 
Eigen::Map<Eigen::Matrix3d> ma(avals_col); 
Eigen::Map<Eigen::Vector3d> mb(b); 
Eigen::Matrix3d inv; 
Eigen::Vector3d _x; 
double t, elapsed; 
t = timeNow(); 
for(int i=0; i<reps; i++) { 
    inv = ma.inverse(); 
} 
elapsed = timeNow() - t; 
printf("Eigen decompose (%dx) time = %g\n", reps, elapsed); 

t = timeNow(); 
for(int i=0; i<reps; i++) { 
    _x.noalias() = inv * mb; 
} 
elapsed = timeNow() - t; 
printf("Eigen solve (%dx) time = %g\n", reps, elapsed); 

das gibt mir:

Eigen inverse and solve (1000000x) time = 0.0209999 
Eigen solve (1000000x) time = 0.000999928 

so viel schneller.

Nun, wenn wir ein viel größeres Problem versuchen, wie 3000 x 3000, erhalten wir mehr als eine Größenordnung Unterschied zugunsten von Eigen:

Eigen decompose (1x) time = 0.411 
GSL decompose (1x) time = 6.073 

Dies ist in der Regel die Optimierungen, die für eine solche Leistung erlaubt große Probleme, die auch einen Overhead für sehr kleine Matrizen verursachen.