Search code examples
rrcppeigeneigenvector

R eigen solver is faster than Eigen's one?

I tried using the eigen solver of the Eigen library in R to improve performance:

// [[Rcpp::export]]
MatrixXd Eigen4(const Map<MatrixXd> bM) {
  SelfAdjointEigenSolver<MatrixXd> es(bM);

  return(es.eigenvectors());
}

Yet, when comparing on a 2000x2000 matrix:

n <- 5e3
m <- 2e3
b <- crossprod(matrix(rnorm(n*m), n))

print(system.time(test <- Eigen4(b))) # 18 sec
print(system.time(test2 <- eigen(b, symmetric = TRUE))) # 8.5 sec

For the result of microbenchmark:

Unit: seconds
    expr                         min        lq      mean    median        uq       max neval
Eigen4(b)                  18.614694 18.687407 19.136380 18.952063 19.292021 20.812116    10
eigen(b, symmetric = TRUE)  8.652628  8.663302  8.696543  8.676914  8.718517  8.831664    10

R is twice as fast as Eigen ? I'm using latest versions of R and RcppEigen.

Am I doing something wrong ?

Solution

  • R's eigen is an interface to Fortran functions from LAPACK. Eigen uses its generic C++ code by default, although it can be configured to use external BLAS/LAPACK backends for certain dense matrix operations, including eigendecomposition. Depending on your architecture and compilers, R's default LAPACK may well be faster. If you configure both R and Eigen to use the same highly optimized platform-specific BLAS/LAPACK (e.g. MKL on Intel) you should get virtually identical (and better) results.