Hello,

I'm in the process of replacing our own linear algebra implementation with Eigen. An interpolation algorithm that uses radial basis functions performs considerably slower with Eigen. In both implementations most of the time is spent in the LU decomposition:

matrixCLU is an n x n matrix (MatrixXd in Eigen).
_lu is an Eigen::FullPivLU<Eigen::MatrixXd>()

Our old code:

Code: Select all

lu(_matrixCLU, _pivotsCLU);

Eigen:

Code: Select all

_lu = matrixCLU.fullPivLu();

matrixCLU is dense and rather ill conditioned.

The _lu.solve(RHS) resp. forward- and backwardSubsitutions with different right-hand-sides takes almost no time after that.

I have also tried partialPivLU but that only performed insignificantly better (not even in all cases). Furthermode I need to know the rank of the matrixCLU which we do like that:

Code: Select all

  int rankDeficiency = 0;
  for (int i=0; i < n; i++){
    if (equals(_matrixCLU(i,i), 0.0)){
      rankDeficiency++;
    }
  }

resp. in Eigen:

Code: Select all

 int rankDeficiency = _lu.rank() - n;

For n = 4000 Eigen took 134s, our implementation took 69s.

Our LU implementation is pretty straighforward. I have uploaded it at http://pastebin.com/rQwuXmeQ (less then 70 lines).

What I can I do to speed up Eigen for this problem?

Thanks a lot,
Florian

Make sure you compiled with compiler optimizations ON. For n=4k, PartialPivLU should be 1 or 2 order of magnitude faster than FullPivLU. You might also consider (ColPiv)HouseholderQR that should be faster than FullPivLU.

Hello,

I did compile with -O3.

I just redid my benchmark runs and I wasn't able to reproduce the results.

My old results are there, probably faulty: http://xgm.de/upload/EigenOld.png Sorry for that!

New ones are there: http://xgm.de/upload/Eigen.png I think this looks more like what one would expect.

Problem that is still there: Is there a way to determine the rank deficiency with partialPivLU, like that

int rankDeficiency = _lu.rank() - n; // with FullPivLU, n being the size of the matrix.

Thanks and sorry for my confusion!

Florian

Partial pivoting LU is not rank-revealing. I guess that the best approximation you could do would be to count the number of significant pivots afterwards, eg.:

Code: Select all

auto abs_diag = lu.matrixLU().diagonal().array().abs();
int rank_approx = (abs_diag>=abs_diag.maxCoeff()*NumTraits<double>::epsilon()).count();