Getting log of determinant as fast with Eigen as with Alglib

Question

I need a quick way to get logarithms of complex determinants, preferably not by getting the determinants first and then taking the log, as the numbers can get really big or really small (and I later use ratios of such numbers but only when they are similar; so the exponential of their difference is well-behaved).

Up until now I have been using the alglib library; doing an LU decomposition and then adding the logs along the diagonal and then adding i*pi times the number of pivots. Assuming I have an alglib::complex_2d_array m of size n, I do

alglib::integer_1d_array pivots;
cmatrixlu(m, n, n, pivots);
int nopivs=0;
for(int j=0;j<n;j++) nopivs+=(pivots(j)!=j);
complex<double> aldet=0;
for(int i=0;i<n;i++) aldet+=log(m[i][i]);
aldet+=complex<double>(0, nopivs*pi);

where I'm using a function

complex<double> log(alglib::complex a) {return log(complex<double>(a.x,a.y);}

However in many ways the Eigen library seems nice; easier to use and uses complex<double> instead of its own complex class. Also I'm already using it for some other purposes so this would streamline things.

I try to use it in a similar way, assuming a Eigen::MatrixXcd m of size n:

Eigen::FullPivLU<Eigen::MatrixXcd> LU(m);
Eigen::MatrixXcd U=LU.matrixLU().triangularView<Eigen::Upper>();
complex<double> Edet=0;
for(int i=0;i<n;i++) Edet+=log(U(i,i));
Edet+=log(CD(LU.permutationP().determinant()*LU.permutationQ().determinant()));

However when I do some testing, Eigen performs much slower.

So I'm wondering if there is another way of doing this with Eigen that would be faster? Maybe another way of getting the log of the determinant altogether?

EDIT: After comment: This is how I test the code:

int n=20, k=5000;
Eigen::MatrixXcd m(n, n);
srand((unsigned int) time(0));
m.setRandom();
alglib::complex_2d_array m2=Eigen2AL_2d(m);

Eigen::FullPivLU<Eigen::MatrixXcd> LU(m);
CD Edet=0.0, aldet=0.0, test=LU.determinant();

clock_t starttime=clock();
for(int i=0;i<k;i++) {
    Eigen::MatrixXcd m4=m;
    Eigen::FullPivLU<Eigen::MatrixXcd> LU(m4);
    Eigen::MatrixXcd U=LU.matrixLU().triangularView<Eigen::Upper>();
    Edet=0;
    for(int i=0;i<n;i++) Edet+=log(U(i,i));
        Edet+=log(CD(LU.permutationP().determinant()*LU.permutationQ().determinant()));
}
cout << "Eigen time: " << (clock()-starttime)/(double)CLOCKS_PER_SEC << endl;

starttime=clock();
for(int i=0;i<k;i++) {
     alglib::integer_1d_array pivots;
    alglib::complex_2d_array m3=m2;
    cmatrixlu(m3, n, n, pivots);
    int nopivs=0;
    for(int j=0;j<n;j++) nopivs+=(pivots(j)!=j);
    aldet=0;
    for(int i=0;i<n;i++) aldet+=log(m3[i][i]);
    aldet+=CD(0, nopivs*pi);
}
cout << "Alglib time: " << (clock()-starttime)/(double)CLOCKS_PER_SEC << endl;

cout << "det = " << test << " " << exp(aldet) << " " << exp(Edet) << endl;

It is compiled with g++ -c -std=c++11 -O2. A typical run gives:

Eigen time: 2.10524
Alglib time: 0.664027

Answer 1

The LU algorithm implemented in alglib performs partial pivoting only, therefore, in Eigen you should use the equivalent PartialPivLU class which is indeed order of magnitude faster. Also, make sure to bench with compiler optimization ON.