c++ Large eigendecomposition speed

Question

As part of my pipeline I need to perform eigendecomposition of a big matrix in the order of 6000x6000. The matrix is dense, so except if I simplify the problem (sot sure if possible) no sparse method can be used.

At the moment I play with toy data. Using the Eigen library for a 513x513 matrix I need ~6.5 seconds, while for a 2049x2049 matrix I need ~130 seconds, which sounds prohibitive since the increase is not linear. This was achieved with Eigen::SelfAdjointEigenSolver, while with other methods like Eigen::EigenSolver or Eigen::ComplexEigenSolver I didn't get notable improvement. The same happened when I tried Armadillo with arma::eig_sym even with the option "dc"` that is supposed to give a faster but approximate result. Armadillo has some methods thatreturn only the first X eigenvalues for speedup, but this is only for sparse methods. At the moment I can probably get away with the first 10-20 eigenvalues.

Is there a way or a library/method that can give me a notable speedup?

Untutored

If you need only the few highest or smallest eigenvectors, then there are more efficient methods. – Ringler

This is exactly what I mention, that this might be a good enough workaround. Which are these methods? Any pointer please? – Untutored

Lapack provides such routines. Regarding your numbers, I obtained 7.5s only for a 2049x2049 matrix using Eigen::SelfAdjointEigenSolver, and 280s for a 6000x6000 matrix. Make sure you compiled with compiler optimizations ON. Of course, this is still prohibitive and better use a dedicated algorithm extracting only the first eigenvectors. – Foliar

Answer 1

Spectra is used to retrieve a few number of eigenvalues of a large matrix.

Sample code to calculate the largest and smallest 10 eigenvalues may look like:

#include <Eigen/Core>
#include <Eigen/Eigenvalues>
#include <MatOp/DenseGenMatProd.h>
#include <MatOp/DenseSymShiftSolve.h>
#include <SymEigsSolver.h>
#include <iostream>

using namespace Spectra;

int main()
{
    srand(0);
    // We are going to calculate the eigenvalues of M
    Eigen::MatrixXd A = Eigen::MatrixXd::Random(1000, 1000);
    Eigen::MatrixXd M = A.transpose() * A;

    // Matrix operation objects
    DenseGenMatProd<double> op_largest(M);
    DenseSymShiftSolve<double> op_smallest(M);

    // Construct solver object, requesting the largest 10 eigenvalues
    SymEigsSolver< double, LARGEST_MAGN, DenseGenMatProd<double> >
        eigs_largest(&op_largest, 10, 30);

    // Initialize and compute
    eigs_largest.init();
    eigs_largest.compute();

    std::cout << "Largest 10 Eigenvalues :\n" <<
        eigs_largest.eigenvalues() << std::endl;

    // Construct solver object, requesting the smallest 10 eigenvalues
    SymEigsShiftSolver< double, LARGEST_MAGN, DenseSymShiftSolve<double> >
        eigs_smallest(&op_smallest, 10, 30, 0.0);

    eigs_smallest.init();
    eigs_smallest.compute();
    std::cout << "Smallest 10 Eigenvalues :\n" <<
        eigs_smallest.eigenvalues() << std::endl;

    return 0;
}

Insinuate

Answer 2

I would recommend to try out Arpack-Eigen. I know from Octave/Matlab that it can compute the largest eigenvalue of a random 2049x2049 within a second, and the largest 10 within 5-20 seconds, eigs(rand(2049), 10). Now, its documentation help eigs points to ARPACK. Arpack-Eigen https://github.com/yixuan/arpack-eigen lets you request 10 eigenvalues from larger matrix like this: SymEigsSolver< double, LARGEST_ALGE, DenseGenMatProd<double> > eigs(&op, 10, 30);.

Ringler

ARPACK-Eigen is now superseded by Spectra. And the ncv parameter should be roughly two or three times the number of eigenvalues requested, so eigs(&op, 10, 30) is probably more appropriate. – Insinuate

@Insinuate Thanks for the correction, I edited my answer. – Ringler