在 Eigen 中提高 SelfAdjointEigenSolver 的精度

Question

我正在尝试确定 Eigen 中稀疏数组的特征值和特征向量。由于我需要计算所有特征向量和特征值，并且我无法使用不受支持的 ArpackSupport 模块工作来完成此操作，因此我选择将系统转换为密集矩阵并使用 SelfAdjointEigenSolver 计算特征系统（我知道我的矩阵是真实的并且具有实特征值）。在我拥有大小为 1024*1024 的矩阵之前，此方法运行良好，但随后我开始与预期结果出现偏差。

根据我的理解，在该模块 (https://eigen.tuxfamily.org/dox/classEigen_1_1SelfAdjointEigenSolver.html) 的文档中，可以更改最大迭代次数：

const int m_maxIterations static Maximum number of iterations.

The algorithm terminates if it does not converge within m_maxIterations * n iterations, where n denotes the size of the matrix. This value is currently set to 30 (copied from LAPACK).

但是，我不明白你是如何实现的，使用他们的例子：

SelfAdjointEigenSolver<Matrix4f> es;
Matrix4f X = Matrix4f::Random(4,4);
Matrix4f A = X + X.transpose();
es.compute(A);
cout << "The eigenvalues of A are: " <<   es.eigenvalues().transpose() << endl;
es.compute(A + Matrix4f::Identity(4,4)); // re-use es to compute eigenvalues of A+I
cout << "The eigenvalues of A+I are: " << es.eigenvalues().transpose() << endl

您将如何修改它以更改最大迭代次数？

此外，这会解决我的问题还是我应该尝试找到替代函数或算法来求解特征系统？

提前致谢。

Answer 1

m_maxIterations 是一个 static const int 变量，因此它可以被认为是类型的内在属性。更改这样的类型属性通常会通过特定的模板参数来完成。然而，在本例中，它被设置为常量 30，因此这是不可能的。

因此，您唯一的选择就是更改头文件中的值并重新编译您的程序。

但是，在这样做之前，我会尝试 Singular Value Decomposition。根据主页，它的准确度是"Excellent-Proven"。此外，它可以克服由于数值上不完全对称矩阵而导致的问题。

Answer 2

增加迭代次数不太可能有帮助。另一方面，从 float 移动到 double 会有很大帮助！

如果这没有帮助，请更具体地说明 "deviations from the expected results"。

Answer 3

我通过编写改编自 Book Numerical Recipes 的 Jacobi 算法解决了这个问题：

void ROTATy(MatrixXd &a, int i, int j, int k, int l, double s, double tau)
{
double g,h;
g=a(i,j);
h=a(k,l);
a(i,j)=g-s*(h+g*tau);
a(k,l)=h+s*(g-h*tau);

}

void jacoby(int n, MatrixXd &a, MatrixXd &v, VectorXd &d )
{
int j,iq,ip,i;
double tresh,theta,tau,t,sm,s,h,g,c;


VectorXd b(n);
VectorXd z(n);

v.setIdentity();    
z.setZero();


for (ip=0;ip<n;ip++)
{   
    d(ip)=a(ip,ip);
    b(ip)=d(ip);
}


for (i=0;i<50;i++) 
{
    sm=0.0;
    for (ip=0;ip<n-1;ip++) 
    {

        for (iq=ip+1;iq<n;iq++)
            sm += fabs(a(ip,iq));
    }

    if (sm == 0.0) {
        break;
    }

    if (i < 3)
    tresh=0.2*sm/(n*n); 
    else
    tresh=0.0;  


    for (ip=0;ip<n-1;ip++)
    {
        for (iq=ip+1;iq<n;iq++)
        {
            g=100.0*fabs(a(ip,iq));
            if (i > 3 && (fabs(d(ip))+g) == fabs(d[ip]) && (fabs(d[iq])+g) == fabs(d[iq]))
            a(ip,iq)=0.0;
            else if (fabs(a(ip,iq)) > tresh)
            {
                h=d(iq)-d(ip);
                if ((fabs(h)+g) == fabs(h))
                {
                    t=(a(ip,iq))/h;
                }   
                else 
                {
                    theta=0.5*h/(a(ip,iq));
                    t=1.0/(fabs(theta)+sqrt(1.0+theta*theta));
                    if (theta < 0.0)
                    {
                        t = -t;
                    }
                    c=1.0/sqrt(1+t*t);
                    s=t*c;
                    tau=s/(1.0+c);
                    h=t*a(ip,iq);


                    z(ip)=z(ip)-h;
                    z(iq)=z(iq)+h;
                    d(ip)=d(ip)- h;
                    d(iq)=d(iq) + h;
                    a(ip,iq)=0.0;


                    for (j=0;j<ip;j++)
                        ROTATy(a,j,ip,j,iq,s,tau);
                    for (j=ip+1;j<iq;j++)
                        ROTATy(a,ip,j,j,iq,s,tau);
                    for (j=iq+1;j<n;j++)
                        ROTATy(a,ip,j,iq,j,s,tau);
                    for (j=0;j<n;j++)
                        ROTATy(v,j,ip,j,iq,s,tau);


                }
            } 
        }
    }


}

}

函数 jacoby 接收方阵 n 的大小、我们要计算的矩阵、我们要求解的矩阵 (a) 和一个将接收每一列中的特征向量的矩阵和一个将要计算的向量接收特征值。它有点慢，所以我尝试将它与 OpenMp 并行化（参见：），但不幸的是，对于 4096x4096 大小的矩阵，我并不是说计算时间有所改善。

在 Eigen 中提高 SelfAdjointEigenSolver 的精度

Increase precision in SelfAdjointEigenSolver in Eigen

c++

eigenvalue

eigenvector

eigen

eigen3