I love you guys for making such a great library!

Usually when I go from Matlab to Eigen, I get a speed up, but this time it was more than 10 times slower. Below are two programs that should to the same thing, Matlab takes less than a second but Eigen takes more than 12 seconds. I tested Eigen 3.2.4 and Dev-version with gcc5.2 -O3 and Matlab R2015b.

Am I doing something wrong?

BR,
Henrik

Eigen:

Code: Select all

#include <Eigen/Dense>
#include <Eigen/Eigenvalues>
#include <iostream>
#include <string>

int main(int /*argc*/, char *argv[]) {
  const unsigned int Vsize = std::stoul(std::string(argv[1]));
  const unsigned int N = std::stoul(std::string(argv[2]));

  Eigen::VectorXd x(Vsize);
  Eigen::MatrixXd C(Vsize, Vsize);

  for (int idx=0;idx!=10;idx++) {
    x.setRandom();
    C = x * x.transpose();

    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> E(C);

    double lambda = E.eigenvalues()[N];
    std::cout << "Eigenvalue " << N << " is " << lambda << '\n';
  }
}


Matlab:

Code: Select all

tic;
Vsize = 1000;
N = 1000;

for idx = 1:10
    x = rand(Vsize, 1);
    C = x * x';
    [V, D] = eig(C);
    fprintf('Eigenvalue %d is %f\n', N, D(N, N));
end

toc


An out-of-the-box version of the Eigen code would be

Code: Select all

#include <Eigen/Dense>
#include <Eigen/Eigenvalues>

#include <iostream>

int main() {
  const unsigned int Vsize = 1000;
  const unsigned int N = 999;

  Eigen::VectorXd x(Vsize);
  Eigen::MatrixXd C(Vsize, Vsize);

  for (int idx=0;idx!=10;idx++) {
    x.setRandom();
    C = x * x.transpose();

    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> E(C);

    double lambda = E.eigenvalues()[N];
    std::cout << "Eigenvalue " << N << " is " << lambda << '\n';
  }
}

It could be that Matlab sees that C is a rank-1 matrix, in which case lambda=|x|^2, all other eigenvalues are 0, and x/|n| is the last eigenvector. So the only work to do is to complete the remaining eigenvectors, e.g.,:

Eigen::MatrixXd V = x.householderQr().householderQ();
V.col(0).swap(V.col(N));

But that's probably not the explanation because in this case, Eigen takes only a few ms, so Matlab should also takes few ms, not an entire second.

Other possible explanations are a combination of
- in SelfAdjointEigenSolver, the computation of the eigenvectors does not scale well wrt matrix size and cache issues (this likely explain a factor 5 or more)
- matlab probably uses the D&C algorithm in MKL which is nicely multithreaded, whereas Eigen uses a slightly slower and non multithread algorithm.