EigenのSelfAdjointEigenSolverの精度を高める

Question

私はEigenの疎な配列の固有値と固有ベクトルを決定しようとしています。私はすべての固有ベクトルと固有値を計算する必要があるので、サポートされていない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 

どのように反復の最大数を変更するためにそれを修正するのでしょうか？

また、これは私の問題を解決するか、または固有のシステムを解決するための代替機能やアルゴリズムを見つけようとすべきですか？

事前に感謝します。

Answer 1

私はブック数値レシピから適応ヤコビアルゴリズムを記述することによって、問題を解決：

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); 


       } 
      } 
     } 
    } 


} 

}

を機能ジャコビーは、正方行列nのサイズ、我々がしたい行列を受け取ります（a）と、各列の固有ベクトルを受け取る行列と、固有値を受け取るベクトルとを求める。私はOpenMp（Parallelization of Jacobi algorithm using eigen c++ using openmp参照）と並列化しようとしましたが、4096x4096のサイズの行列については、計算時間の改善を意味しませんでした。残念なことに、

Answer 2

m_maxIterationsは、static const intであり、このように、この型の固有の性質と考えることができる。このような型プロパティを変更するには、通常、特定のテンプレートパラメータを使用します。ただし、この場合は定数番号30に設定されているため、不可能です。

したがって、ヘッダーファイルの値を変更してプログラムを再コンパイルするだけです。

しかし、それを行う前に、私はSingular Value Decompositionを試してみましょう。ホームページによると、その精度は「優れている」と評価されています。さらに、数値的に完全に対称な行列ではない問題を克服することができます。

Answer 3

繰り返しの回数を増やすことは役立ちません。一方、floatからdoubleに移動すると、多くの助けになります！

それが役に立たない場合は、「期待された結果からの逸脱」を具体的にご指定ください。