0

我需要一种快速的方法来获取复杂行列式的对数，最好不要先获取行列式然后取对数，因为数字可以变得非常大或非常小（我后来使用这些数字的比率，但仅在它们相似时; 所以他们的差异的指数表现良好）。

到目前为止，我一直在使用 alglib 库；进行 LU 分解，然后沿对角线添加对数，然后添加 i*pi 乘以枢轴点数。假设我有一个alglib::complex_2d_array msize 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>而不是它自己的复杂类。此外，我已经将它用于其他目的，因此这将简化事情。

我尝试以类似的方式使用它，假设 aEigen::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
4

2 回答 2

2

在 alglib 中实现的 LU 算法仅执行部分旋转，因此，在 Eigen 中，您应该使用等效的 PartialPivLU 类，它确实要快一个数量级。此外，请确保在编译器优化开启的情况下进行测试。

于 2015-12-22T20:03:49.523 回答
1

Eigen 中的所有计算时间都花在了 LU 分解上（称为 when Eigen::FullPivLU<Eigen::MatrixXcd> LU(m4)）。其余的可以忽略不计。AFAIK，您无能为力来改善这一点。

于 2015-12-22T15:37:37.543 回答