与Alglib一样,使用Eigen快速获取行列式的对数

我需要一种快速的方法来计算复杂的行列式的对数,最好不要先求行列式,然后再求对数,因为这些数字可能真的很大,也可能真的很小(我后来使用了这些数字的比率,但只有当它们相似时才使用;所以它们的差异指数表现良好)。

到目前为止,我一直在使用alglib库;进行LU分解,然后沿对角线添加日志,然后添加i*pi乘以枢轴数。假设我有一个 alglib::complex_2d_array m 大小的 n ,我知道

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);

我在使用函数

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

然而,在许多方面,Eigen库看起来不错;更易于使用和使用 complex<double> 而不是它自己的复杂类。此外,我已经将其用于其他目的,这样可以简化工作。

我尝试以类似的方式使用它,假设 Eigen::MatrixXcd m 大小的 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()));

然而,当我做一些测试时,Eigen的表现要慢得多。

所以我想知道是否有另一种方法可以用Eigen做得更快?也许还有另一种方法可以完全得到行列式的对数?

编辑:评论后:这是我测试代码的方式:

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;

它是用 g++ -c -std=c++11 -O2 。典型运行给出:

Eigen time: 2.10524
Alglib time: 0.664027