I am trying to compute Eigen vectors for a matrix in C++ in a very efficient way. The problem is that the most representative C++ libraries OpenCV, Eigen and Armadillo are considerably slower than MATLAB equivalent eigen function. This is devastating and very concerning for me, I would never expect that MATLAB could beat C++ code performance, specially if it is a very well known/used library like the aforementioned ones. To give an idea of the performance difference, below find all the computation times for eigen vector computation in different libraries (units in seconds, some of them are in hours and days! yes days!) (all of them give me the expect results including MATLAB, it just the execution time is vastly difference)
The unfortunate part is that for my purposes I need to compute the Eigen values for a matrix whose dimensions are (8192, 8192). From the table you can see that MATLAB only takes 17 secs, and the second best (armadillo) takes 47 secs, might not sound like a big deal, but I need to repeat this operation thousands of times, so you can imagine how this will add up time quickly. I would tremendously appreciate if someone can tell me what I am doing wrong or this is just the sad reality I have to face and just have slow C++ code (at least slower than MATLAB) (by the way, props to the MATLAB coders that manage to have the eig function considerably faster than any other C++ library). For those who are interested in looking at the code and way how I compute the eigen vectors, I will leave it down here.
PD: All these methods were tested using release builds for the libraries (including libopenblas.dll for armadillo, excepting eigen that is a header library)
PD: All these computation times were obtained using same computer
// COV is a covariance matrix, meaning is a symmetric matrix
// has just a bunch of double numbers (CV_64FC1).
// For those interesting this covariance matrix COV is just
// the result of this operation COV = X' * X, where M is just literally
// any double (CV_64FC1) matrix dimensions (m, 8192), m can be any positive
// value to be honest.
cv::Mat X = read_opencv_matrix() // dim(X) = (8192, m)
cv::Mat COV, MEAN;
// dim(COV) = (8192, 8192)
cv::calcCovarMatrix(X.t(), COV, MEAN, cv::COVAR_NORMAL | cv::COVAR_ROWS);
int numRows = X.rows; // should be 8192
int numCols = X.cols; // can be anything to be honest
// computing eigen values using different libraries (opencv, armadillo, eigen)
// all of them give me the same results (including MATLAB function eig())
// the problem is that all of them are considerably slower than MATLAB
///////////////////////////////// OPENCV //////////////////////////
// opencv + eigen
if (do_opencv_eigen)
{
cv::Mat D, V;
std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();
cv::eigen(cov, D, V);
std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
auto count = std::chrono::duration_cast<std::chrono::milliseconds>(end - begin).count();
std::cout << "[TIME] OpenCV + cv::eigen = " << static_cast<double>(count) / 1000.0 << "[sec]" << std::endl;
}
// opencv + SVD
if (do_opencv_eigen)
{
cv::Mat u, w, vt;
std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();
cv::SVD::compute(cov, w, u, vt);
std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
auto count = std::chrono::duration_cast<std::chrono::milliseconds>(end - begin).count();
std::cout << "[TIME] OpenCV + cv::SVD = " << static_cast<double>(count) / 1000.0 << "[sec]" << std::endl;
}
// opencv + SVD + MODIFY_A flag
if (do_opencv_svd_mod_a)
{
cv::Mat u2, w2, vt2;
std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();
cv::SVD::compute(cov, w2, u2, vt2, cv::SVD::MODIFY_A);
std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
auto count = std::chrono::duration_cast<std::chrono::milliseconds>(end - begin).count();
std::cout << "[TIME] OpenCV + cv::SVD::MODIFY_A = " << static_cast<double>(count) / 1000.0 << "[sec]" << std::endl;
}
///////////////////////// ARMADILLO /////////////////////////////
arma::mat arma_cov = Utils::opencv2arma(cov); // helper function to convert cv::mat to arma::mat
arma::mat arma_X = Utils::opencv2arma(X);
arma::mat arma_col_mean_rep = Utils::opencv2arma(col_mean_rep);
// compute arma eigen gen vectors
if (do_arma_eigen_gen)
{
arma::cx_vec arma_Dc;
arma::cx_mat arma_Vc;
std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();
arma::eig_gen(arma_Dc, arma_Vc, arma_cov);
std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
auto count = std::chrono::duration_cast<std::chrono::milliseconds>(end - begin).count();
std::cout << "[TIME] Arma + arma::eig_gen = " << static_cast<double>(count) / 1000.0 << "[sec]" << std::endl;
}
// compute arma eigen gen vectors
if (do_arma_eigen_sym)
{
arma::vec arma_D;
arma::mat arma_V;
std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();
arma::eig_sym(arma_D, arma_V, arma_cov);
std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
auto count = std::chrono::duration_cast<std::chrono::milliseconds>(end - begin).count();
std::cout << "[TIME] Arma + arma::eig_sym = " << static_cast<double>(count) / 1000.0 << "[sec]" << std::endl;
}
// armadillo + svd
if (do_arma_svd)
{
arma::mat arma_U2;
arma::vec arma_s2;
arma::mat arma_V2;
std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();
arma::svd(arma_U2, arma_s2, arma_V2, arma_cov);
std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
auto count = std::chrono::duration_cast<std::chrono::milliseconds>(end - begin).count();
std::cout << "[TIME] Arma + arma::svd = " << static_cast<double>(count) / 1000.0 << "[sec]" << std::endl;
}
// armadillo + svd + econ
if (do_arma_svd_econ)
{
arma::mat arma_U2_econ;
arma::vec arma_s2_econ;
arma::mat arma_V2_econ;
std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();
arma::svd_econ(arma_U2_econ, arma_s2_econ, arma_V2_econ, arma_cov);
std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
auto count = std::chrono::duration_cast<std::chrono::milliseconds>(end - begin).count();
std::cout << "[TIME] Arma + arma::svd_econ = " << static_cast<double>(count) / 1000.0 << "[sec]" << std::endl;
}
/////////////////// EIGEN /////////////////////////////
Eigen::Matrix eig_cov = Utils::opencv2eigen(cov); // helper function to convert cv::mat to Eigen::Matrix
Eigen::Matrix eig_X = Utils::opencv2eigen(X);
Eigen::Matrix eige_col_mean_rep = Utils::opencv2eigen(col_mean_rep);
//Eigen general eigen function
if (do_eigen_eigen)
{
std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();
Eigen::EigenSolver<Eigen::MatrixXd> es(eig_cov);
Eigen::MatrixXcd eig_VC = es.eigenvectors();
Eigen::MatrixXd eig_V = eig_VC.real();
std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
auto count = std::chrono::duration_cast<std::chrono::milliseconds>(end - begin).count();
std::cout << "[TIME] Eigen + Eigen::EigenSolver = " << static_cast<double>(count) / 1000.0 << "[sec]" << std::endl;
}
// eigen library + SVD
if (do_eigen_SVD)
{
std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();
Eigen::BDCSVD<Eigen::MatrixXd> SVD(eig_cov, Eigen::ComputeThinU | Eigen::ComputeThinV);
Eigen::MatrixXd eig_V2 = SVD.matrixV();
std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
auto count = std::chrono::duration_cast<std::chrono::milliseconds>(end - begin).count();
std::cout << "[TIME] Eigen + Eigen::BDCSVD = " << static_cast<double>(count) / 1000.0 << "[sec]" << std::endl;
}
//Eigen library + SelfAdjointEigenSolver
if (do_eigen_SelfAdjoin)
{
std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> esa;
esa.compute(eig_cov);
std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
auto count = std::chrono::duration_cast<std::chrono::milliseconds>(end - begin).count();
std::cout << "[TIME] Eigen + Eigen::SelfAdjointEigenSolver = " << static_cast<double>(count) / 1000.0 << "[sec]" << std::endl;
}
UPDATE 1 OS : Windows 10 IDE : Visual Studio 2022 CPU : Intel Xeon CPU E5 2687W v4 OpenCV Version: 4.5.3 (compiled with the official intrusction from opencv https://docs.opencv.org/4.x/d3/d52/tutorial_windows_install.html) Armadillo: uses the default libopenblas.dll that comes with it Eigen: I downloaded the lastest and greatest version up today
