title: Eigen library: SVD is slow compared to GSL
body: <p>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:</p>
  
  <pre><code>#include &lt;gsl/gsl_rng.h&gt;
  #include &lt;gsl/gsl_linalg.h&gt;
  #include &lt;unsupported/Eigen/SVD&gt;
  #include &lt;ctime&gt;
  #include &lt;iostream&gt;
  
  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&lt;size; i++)
      for(int j=0; j&lt;size; j++)
          mat(i, j) = gsl_rng_uniform(rng);
  gsl_matrix_view matv = gsl_matrix_view_array(mat.data(), size, size);
  gsl_matrix* mat1 = &amp;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&lt;Eigen::MatrixXd&gt; jac(mat, Eigen::ComputeThinU | Eigen::ComputeThinV);
  double tjac = (std::clock() - tbegin)*1.0/CLOCKS_PER_SEC;
  tbegin = std::clock();
  Eigen::BDCSVD&lt;Eigen::MatrixXd&gt; 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 &lt;&lt; "Matrix size: "&lt;&lt;size&lt;&lt;", GSL SVD: "&lt;&lt;tgsl&lt;&lt;" s, JacobiSVD: "&lt;&lt;tjac&lt;&lt;" s, BDCSVD: "&lt;&lt;tbdc&lt;&lt;" s\n";
  gsl_vector_free(vec2);
  gsl_vector_free(vec1);
  gsl_matrix_free(mat2);
  }
  </code></pre>
  
  <p>Machine 1 (old, gcc 4.2.1, compiler options -O3):</p>
  
  <pre>
  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
  </pre>
  
  <p>Machine 2 (gcc 4.8.3, -O3):</p>
  
  <pre>
  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
  </pre>
  
  <p>Machine 3 (icc 13.1.3, -O3):</p>
  
  <pre>
  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
  </pre>
  
  <p>The question is, why is Eigen not beating GSL?</p>