Eigen linear algebra solvers seem slow

Question

I want to solve a linear algebraic equation Ax = b using Eigen solvers. In my case, A is a complex sparse matrix(26410*26410), b is a real vector (26410*1).

I use mex file in MATLAB to map the sparse matrix A and vector b to Eigen accepted format. The reason why I use Eigen solver is to hope it would be faster than solving directly in MATLAB using x = A\b.

However, after tried LDLT, SparseLU, CG and BiCGSTAB, I found the results are not very satisfying:

LDLT takes 1.462s with norm(A*x - b)/norm(b) = 331; SparseLU takes 37.994s with 1.5193e-4; BiCGSTAB takes 95.217s with 4.5977e-4; On the contrast, directly use x = A\b in MATLAB consumes 13.992s with norm of the error 2.606e-5.

I know it is a little stupid and also time consuming to map the sparse matrix A and vector b in MATLAB workspace to Eigen. But I am wondering whether the results I got are the best results which Eigen can give? Anyone can give me some pointers? Should I try some other linear equation solvers? Thanks a lot in advance! The following is the main part of codes.

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {

    //input vars

    //temp var
    size_t nrows;

    //output vars
    //double *x;

    //GetData
    /* check inputs 
    ...*/

    //"mxArray2SCM" is a sub-function for map the complex sparse matrix in Eigen
    SpCMat A = mxArray2SCM(prhs[0]);
    //SpMat A = mxArray2SM(prhs[0]);

    //"mxArray2ECV" is a sub-function for map the real vector in Eigen
    Eigen::VectorXcd b = mxArray2ECV(prhs[1]);
    //Eigen::VectorXd b = mxArray2EV(prhs[1]);   
    nrows = b.size();
    //Computation
    Eigen::VectorXcd x(nrows);
    //SparseLU<SparseMatrix<CD> > solver;
    BiCGSTAB<SparseMatrix<CD>,IncompleteLUT<CD> > BiCG;
    //BiCG.preconditioner().setDroptol(0.001);
    BiCG.compute(A);
    if(BiCG.info()!=Success){
        //decomposition failed
        return;
    }
    x = BiCG.solve(b);

    //Output results
    plhs[0] = ECV2mxArray(x);   
}

Kevyn

Do you have compiler optimizations enabled? – Gastrectomy

Yes. I use "mex -O -largeArrayDims xxx.cpp" command to compile. By the way, version of MatLAB is R2016b and compiler is VS2015 Professional.@Baum mit Augen – Kevyn

Why would you guess that Eigen would be faster than Matlab? Matlab relies on highly optimized and multi-threaded external libraries such as UmfPack and MKL. Eigen can also make use of them. For instance, if you use Eigen::UmfpackLU and link to UmfPack, you will probably get the same performance at MatLab. – Lehr

Since I have to do many iterations for my optimization, and in each iteration a linear algebra equation is solved, then I expect to find a way to reduce the computation time for each iteration. Thus, I tried Eigen to expect it could be faster. It seems that I was wrong. As you said, UmfPack maybe can provide the same performance as MatLab. Could you give me some tips for using that in Eigen? Thanks a lot! @Lehr – Kevyn

Based on what the matrix you describe, Matlab’s A\b will use a sparse LU solver. It’s not surprising LDLT fails (provides non-epsilon error) because Eigen’s LDLT is for square-Hermitian matrixes only. – Stradivari

As a sanity check, can you time how long it takes to convert Matlab format to Eigen format and back? I’d assume it’s 1 to 1.5 seconds (your LDLT example is around 1.5 seconds and I’m assuming that’s just a round-trip since the LDLT solver sees your matrix is non-square)? – Stradivari

@AhmedFasih For convert MatLAB sparse matrix to eigen and transfer back, it only takes about 0.07s. And my matrix A is square. I know LDLT is not quite appropriate to use, but I think it indeed solved the equation, though the solution is not very accurate. – Kevyn

Sorry, I totally misread the first dimension of A 😣! But like @Lehr says, Matlab’s linear algebra routines are optimized with Intel MKL and while it is disappointing Eigen is not faster, it’s not surprising either. It might be worth posting your findings on the Eigen mailing list—that project likes to have good benchmarks so a nearly 3× penalty might get their attention. – Stradivari

Do you know if your matrix has any properties that might make it amenable to preconditioners? Matlab’s pcg conjugate gradient solver can accept preconditioners to speed up convergence. – Stradivari

An error norm of 300 is so far beyond numerical precision that I wouldn’t say LDLT “solved” your equation 😋 it just returned something, which happened to be wrong! – Stradivari

For LDLT, your matrix needs to be at least hermitian, i.e., norm(A'-A)=0. – Lehr

@AhmedFasih Thanks! I will go have a try at pcg with preconditioners in Matlab, but not give it too much hope, after all it's a iteration solver. – Kevyn

@Lehr You are right. My matrix is not hermitian, norm(A'-A) is far from zero. Thanks! – Kevyn

Iterative solvers like pcg work really well with sparse arrays. Fast convergence with tunable errors, even if you don't have a preconditioner. Try it. – Stradivari

@AhmedFasih I have tried pcg in Matlab. Although it don't give any error, the solution x is all zeros. After checking the document of pcg, it is said that matrix A (A*x = b) must be symmetric and positive definite, however, in my case, A is neither symmetric nor definite. So I guess pcg does not work for me. – Kevyn

Omg, I am so sorry—my brain is obviously not working today, I’m going home and sleeping. – Stradivari

@AhmedFasih haha, never mind, thanks anyway! – Kevyn

Matlab is a matrix manipulation program. If you can beat the internal eigen solvers by knocking up your own little max function then Matlab's people can't be very good at their jobs. – Barbate

Answer 1

Have you considered using PetSc for Krylov solvers or SLEPc to compute eigenvalues?

Make sure you analyze the eigenspectrum before using a specific Krylov solver (CG works only for symmetric positive definite matrices).

PETSc has quite a few solvers that you can try out based on your eigenspectrum.

You can check Y. Saad's book on how these solvers work.

If your matrix is not symmetric positive definite GMRES is a good option.

Candlestand