提高SelfAdjointEigenSolver在本征中的精度

提问者：小点点

提高SelfAdjointEigenSolver在本征中的精度

我试图确定一个稀疏阵列的本征值和本征向量。由于我需要计算所有的特征向量和特征值，而我无法使用不支持的ArpackSupport模块来完成这一工作，所以我选择将系统转换为密集矩阵，并使用SelfAdjointEigenSolver来计算本征系统（我知道我的矩阵是实的，并且具有实的特征值）。在我的矩阵大小为1024*1024之前，这种方法工作得很好，但随后我开始得到与预期结果有偏差的结果。

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

const int m_maxIterations静态最大迭代数。

如果算法在m_maxIterations*n次迭代内不收敛，则该算法终止，其中n表示矩阵的大小。此值当前设置为30（从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

为了改变最大迭代次数，您将如何修改它？

另外，这会解决我的问题吗，还是我应该尝试找到一个替代函数或算法来解决本征系统？

我提前表示感谢。

共3个答案

匿名用户

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

如果这无济于事，请更具体地说明“与预期结果的偏差”。

匿名用户

是一个变量，因此可以将其视为该类型的固有属性。更改这样的类型属性通常可以通过特定的模板参数来完成。但是，在本例中，它被设置为常量，因此它是不可能的。

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

但是，在这样做之前，我要尝试一下奇异值分解。据主页介绍，它的精确度是“卓越--出处”（excellent-proven)。此外，它还能克服由于矩阵在数值上不完全对称而引起的问题。

匿名用户

我通过编写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的大小，方阵n是我们要计算的矩阵，我们要求解的矩阵(a)，一个矩阵将接收每一列的特征向量，一个向量将接收特征值。它有点慢，所以我尝试用OpenMp并行化它（请参阅:Jacobi算法的并行化使用eigen C++使用OpenMp），但不幸的是，对于4096x4096大小的矩阵，我并不意味着计算时间的改善。