使用 Eigen 与使用 Alglib 一样快速地获取行列式的日志

Question

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

到目前为止，我一直在使用 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> 而不是其自身的复杂 class。此外，我已经将它用于其他一些目的，所以这会简化事情。

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

Answer 1

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

Answer 2

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

使用 Eigen 与使用 Alglib 一样快速地获取行列式的日志

Getting log of determinant as fast with Eigen as with Alglib

performance

logarithm

determinants

eigen

alglib