用特征和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);}

然而，在许多方面，特征库似乎很好；使用起来更容易，并且使用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()));

然而，当我做一些测试时，艾根的表现要慢得多。

所以，我想知道是否还有其他方法可以更快地解决艾根的问题？也许是另一种得到行列式日志的方法？

编辑:注释后:我就是这样测试代码的：

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

在算法库中实现的LU算法只执行部分旋转，因此，在特征中，您应该使用等效的PartialPivLU类，它确实更快。此外，确保与编译器优化的台架。

Answer 2

所有以特征表示的计算时间都用于LU分解(称为Eigen::FullPivLU<Eigen::MatrixXcd> LU(m4))。其余的都可以忽略不计。你没什么可以改进的，AFAIK。