4

I am comparing the performance of Eigen (v3.2.8) and GSL (v2.1) libraries in various linear algebra operations. While in most operations Eigen wins by a large margin (factor of few), in the computation of singular-value decomposition it lags behind. I use two implementations of SVD - JacobiSVD and BDCSVD (the latter is still in the 'unsupported' modules as of this version), and tested the following code on several platforms:

#include <gsl/gsl_rng.h>
#include <gsl/gsl_linalg.h>
#include <unsupported/Eigen/SVD>
#include <ctime>
#include <iostream>

void testSVD(const int size)
{
clock_t tbegin;
gsl_rng* rng = gsl_rng_alloc(gsl_rng_default);
gsl_rng_set(rng, 42);
Eigen::MatrixXd mat(size, size);
for(int i=0; i<size; i++)
    for(int j=0; j<size; j++)
        mat(i, j) = gsl_rng_uniform(rng);
gsl_matrix_view matv = gsl_matrix_view_array(mat.data(), size, size);
gsl_matrix* mat1 = &matv.matrix;
gsl_matrix* mat2 = gsl_matrix_alloc(size, size);
gsl_vector* vec1 = gsl_vector_alloc(size);
gsl_vector* vec2 = gsl_vector_alloc(size);
tbegin = std::clock();
Eigen::JacobiSVD<Eigen::MatrixXd> jac(mat, Eigen::ComputeThinU | Eigen::ComputeThinV);
double tjac = (std::clock() - tbegin)*1.0/CLOCKS_PER_SEC;
tbegin = std::clock();
Eigen::BDCSVD<Eigen::MatrixXd> bdc(mat, Eigen::ComputeThinU | Eigen::ComputeThinV);
double tbdc = (std::clock() - tbegin)*1.0/CLOCKS_PER_SEC;
tbegin = std::clock();
gsl_linalg_SV_decomp(mat1, mat2, vec1, vec2);
double tgsl = (std::clock() - tbegin)*1.0/CLOCKS_PER_SEC;
std::cout << "Matrix size: "<<size<<", GSL SVD: "<<tgsl<<" s, JacobiSVD: "<<tjac<<" s, BDCSVD: "<<tbdc<<" s\n";
gsl_vector_free(vec2);
gsl_vector_free(vec1);
gsl_matrix_free(mat2);
}

Machine 1 (old, gcc 4.2.1, compiler options -O3):

Matrix size: 125, GSL SVD: 0.054497 s, JacobiSVD: 0.143428 s, BDCSVD: 0.177061 s
Matrix size: 250, GSL SVD: 0.457591 s, JacobiSVD: 1.16189 s, BDCSVD: 1.291 s
Matrix size: 500, GSL SVD: 6.01799 s, JacobiSVD: 15.2217 s, BDCSVD: 12.1353 s
Matrix size: 1000, GSL SVD: 67.3727 s, JacobiSVD: 155.501 s, BDCSVD: 114.372 s

Machine 2 (gcc 4.8.3, -O3):

Matrix size: 125, GSL SVD: 0.008666 s, JacobiSVD: 0.029863 s, BDCSVD: 0.026779 s
Matrix size: 250, GSL SVD: 0.082332 s, JacobiSVD: 0.318304 s, BDCSVD: 0.235657 s
Matrix size: 500, GSL SVD: 0.807378 s, JacobiSVD: 2.93758 s, BDCSVD: 1.98003 s
Matrix size: 1000, GSL SVD: 23.1846 s, JacobiSVD: 59.1015 s, BDCSVD: 27.8348 s

Machine 3 (icc 13.1.3, -O3):

Matrix size: 125, GSL SVD: 0.02 s, JacobiSVD: 0.05 s, BDCSVD: 0.04 s
Matrix size: 250, GSL SVD: 0.15 s, JacobiSVD: 0.4 s, BDCSVD: 0.28 s
Matrix size: 500, GSL SVD: 2.18 s, JacobiSVD: 3.51 s, BDCSVD: 2.41 s
Matrix size: 1000, GSL SVD: 10 s, JacobiSVD: 42.87 s, BDCSVD: 23.4 s

The question is, why is Eigen not beating GSL?

Jongware
  • 22,200
  • 8
  • 54
  • 100
user3097263
  • 63
  • 1
  • 5

1 Answers1

5

JacobiSVD is expected to be slow for large matrices. Regarding BDCSVD, please try the one in the 3.3-beta1 release (or devel branch). It is now in the official Eigen/SVD module and it has been considerably improved.

For the record, on OSX, using MacPorts's GSL, -O3 -mavx:

Matrix size: 100, GSL SVD: 0.006717 s, JacobiSVD: 0.01651 s, BDCSVD: 0.007867 s
Matrix size: 200, GSL SVD: 0.060359 s, JacobiSVD: 0.179859 s, BDCSVD: 0.031222 s
Matrix size: 300, GSL SVD: 0.228808 s, JacobiSVD: 0.806909 s, BDCSVD: 0.073439 s
Matrix size: 400, GSL SVD: 0.644856 s, JacobiSVD: 1.98206 s, BDCSVD: 0.132627 s
Matrix size: 600, GSL SVD: 4.24252 s, JacobiSVD: 7.24865 s, BDCSVD: 0.345525 s
Matrix size: 1000, GSL SVD: 33.7155 s, JacobiSVD: 78.6189 s, BDCSVD: 1.31173 s
ggael
  • 28,425
  • 2
  • 65
  • 71
  • @user3097263, for small matrix, `GSL` may be after, but `BDCSVD` in `Eigen` is faster for larger matrix. My benchmark is: `Matrix size: 100, GSL SVD: 0.008633 s, JacobiSVD: 0.058554 s, BDCSVD: 0.0107 s` and `Matrix size: 1000, GSL SVD: 8.68049 s, JacobiSVD: 35.4024 s, BDCSVD: 1.32387 s`. I used `g++ -DNDEBUG -O2 -o tt -msse2 tt.cpp -I . -lgsl -lcblas -lgslcblas` to compile. – zhanxw Nov 21 '16 at 23:57
  • I found `JacobiSVD` to give a much more accurate solution than `BDCSVD`. – Felix Crazzolara May 24 '21 at 19:39