Увеличение точности в SelfAdjointEigenSolver in Eigen

Question

Я пытаюсь определить собственные значения и собственные векторы разреженного массива в Eigen. Поскольку мне нужно вычислить все собственные векторы и собственные значения, и я не мог этого сделать, используя неподдерживаемый модуль ArpackSupport, я решил преобразовать систему в плотную матрицу и вычислить систему eigens с использованием SelfAdjointEigenSolver (я знаю, что моя матрица реальна и имеет вещественные собственные значения). Это хорошо работает, пока у меня не будет матриц размером 1024 * 1024, но затем я начну получать отклонения от ожидаемых результатов.

В документации этого модуля (https://eigen.tuxfamily.org/dox/classEigen_1_1SelfAdjointEigenSolver.html) от того, что я понял, что можно изменить число итераций макс:

Const ИНТ 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 

Как бы вы изменить его, чтобы изменить максимальное число итераций?

Кроме того, это решит мою проблему или я должен попытаться найти альтернативную функцию или алгоритм для решения системы eigensystem?

Мое спасибо заранее.

Answer 1

Я решил проблему, написав алгоритм Якоби адаптированный из книги Numerical Recipes:

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


       } 
      } 
     } 
    } 


} 

}

функция Джейкоби принимает размер квадратной матрицы п, матрицы мы хотим вычислить мы хотим решить (а) и матрицу, которая получит собственные векторы в каждом столбце и вектор, который будет получать собственные значения.Это немного медленнее, поэтому я попытался распараллелить его с помощью 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 поможет много!

Если это не поможет, пожалуйста, уточните «отклонения от ожидаемых результатов».