Увеличьте точность в SelfAdjointEigenSolver в 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
Как бы вы изменили его, чтобы изменить максимальное количество итераций?
Кроме того, решит ли это мою проблему или я должен попытаться найти альтернативную функцию или алгоритм для решения собственной системы?
Заранее спасибо
Решение
Увеличение количества итераций вряд ли поможет. С другой стороны, переходя от float в double очень поможет!
Если это не поможет, пожалуйста, уточните «отклонения от ожидаемых результатов».
Другие решения
m_maxIterations это static const int переменная, и как таковая она может рассматриваться как внутреннее свойство типа. Изменение такого свойства типа обычно выполняется через определенный параметр шаблона. В этом случае, однако, он установлен на постоянное число 30, так что это невозможно.
Поэтому вам остается только изменить значение в заголовочном файле и перекомпилировать вашу программу.
Однако, прежде чем сделать это, я бы попробовал Разложение по единственному значению. Согласно домашней странице, его точность «отлично-проверено». Более того, он может преодолеть проблемы из-за численно не полностью симметричных матриц.
Я решил проблему, написав алгоритм Якоби, адаптированный из книги «Численные рецепты»:
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, матрицу, которую мы хотим вычислить, которую мы хотим решить (а), и матрицу, которая будет получать собственные векторы в каждом столбце, и вектор, который будет получать собственные значения. Это немного медленнее, поэтому я попытался распараллелить его с OpenMp (см .: Распараллеливание алгоритма Якоби с использованием собственного c ++ с использованием openmp), но для матриц размером 4096×4096, что я не имел в виду улучшение времени вычислений, к сожалению.