Сопряженный градиент в собственных для неположительно определенных матриц - PullRequest
Новая
       16

Сопряженный градиент в собственных для неположительно определенных матриц

0 голосов
/ 26 октября 2018

В описании метода CG в библиотеке Eigen можно найти утверждение:

Этот класс позволяет решать линейные задачи A.x = b с использованием алгоритма итеративного сопряженного градиента. Матрица A должна быть самосопряженной.

Однако в литературе метод сопряженных градиентов обычно представлен для реальных симметричных положительно определенных матриц. Примеры показывают, что Eigen CG действительно работает для неположительно определенных матриц, которые не может обработать matlab pcg .

Например, запустив код: 

#include <cstdio>
#include <iostream>
#include <vector>
#include "Eigen/Dense"
#include "Eigen/IterativeLinearSolvers" 
#include "Eigen/Eigenvalues"

int main()
{
    srand(static_cast<unsigned int>(time(0)));
    const int N = 10;
    Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor>  S(N,N);
    const Eigen::Matrix<double,Eigen::Dynamic,2> sources = Eigen::MatrixXd::Random(N,2);

    for(size_t iEx = 0; iEx < 4; iEx++ )
    {
        std::cout<<"EX "<<iEx<<":\n";
        if(iEx == 0)
            for(int i = 0; i < N; i++)
                for(int j = i; j < N; j++)
                    S(i,j) = S(j,i) = 1./std::sqrt((sources.row(i) - sources.row(j)).squaredNorm() +1.);
        if(iEx == 1)
            for(int i = 0; i < N; i++)
                for(int j = i; j < N; j++)
                    S(i,j) = S(j,i) = (sources.row(i) - sources.row(j)).norm();
        if(iEx == 2)
            for(int i = 0; i < N; i++)
                for(int j = i; j < N; j++)
                    S(i,j) = S(j,i) = sources.row(i).dot(sources.row(j));

        if(iEx == 3)
            S = Eigen::MatrixXd::Random(N,N).selfadjointView<Eigen::Lower>();

        Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> Sadj = S.selfadjointView<Eigen::Lower>();
        std::cout<<"\tIS SELFADJOINT: "<<((Sadj.array() == S.array()).all()?"YES\n":"NO\n");
        Eigen::EigenSolver< Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> > eigensolver(S);
        std::cout<<"\tNUMBER OF NEGATIVE EIGEN VALUES: "<<(eigensolver.eigenvalues().real().array() < 0.).count()<<" OUT OF "<<N<<"\n";

        const Eigen::Matrix<double,Eigen::Dynamic,1> xExact = Eigen::VectorXd::Ones(N);
        const Eigen::Matrix<double,Eigen::Dynamic,1> b = S * xExact;

        Eigen::ConjugateGradient< Eigen::MatrixXd, Eigen::Lower|Eigen::Upper, Eigen::IdentityPreconditioner> cg(S);
        cg.setMaxIterations(3000);
        cg.setTolerance(1e-10);

        const Eigen::Matrix<double,Eigen::Dynamic,1> xSol = cg.solve(b);
        std::cout<<"\tITERATIONS       : " << cg.iterations() << "\n";
        std::cout<<"\tESTIMATED ERROR  : " << cg.error()      << "\n";

        std::cout<<"\tNORM 2 ERROR     : "<<(xExact-xSol).norm()<<"\n";
        std::cout<<"\tNORM 2 AVG ERROR : "<<(xExact-xSol).norm()/static_cast<double>(N)<<"\n";
        std::cout<<"\tNORM INF ERROR   : "<<(xExact-xSol).lpNorm<Eigen::Infinity>()<<"\n";
        std::cout<<std::flush;
    }
}

Дает вывод: 

EX 0:
        IS SELFADJOINT: YES
        NUMBER OF NEGATIVE EIGEN VALUES: 0 OUT OF 10
        ITERATIONS       : 11
        ESTIMATED ERROR  : 1.01319e-11
        NORM 2 ERROR     : 2.49293e-10
        NORM 2 AVG ERROR : 2.49293e-11
        NORM INF ERROR   : 1.20759e-10
EX 1:
        IS SELFADJOINT: YES
        NUMBER OF NEGATIVE EIGEN VALUES: 9 OUT OF 10
        ITERATIONS       : 10
        ESTIMATED ERROR  : 2.43788e-12
        NORM 2 ERROR     : 1.77969e-11
        NORM 2 AVG ERROR : 1.77969e-12
        NORM INF ERROR   : 8.2061e-12
EX 2:
        IS SELFADJOINT: YES
        NUMBER OF NEGATIVE EIGEN VALUES: 4 OUT OF 10
        ITERATIONS       : 1
        ESTIMATED ERROR  : 1.72812e-16
        NORM 2 ERROR     : 2.97281
        NORM 2 AVG ERROR : 0.297281
        NORM INF ERROR   : 1.45547
EX 3:
        IS SELFADJOINT: YES
        NUMBER OF NEGATIVE EIGEN VALUES: 5 OUT OF 10
        ITERATIONS       : 9
        ESTIMATED ERROR  : 7.73713e-14
        NORM 2 ERROR     : 8.55003e-14
        NORM 2 AVG ERROR : 8.55003e-15
        NORM INF ERROR   : 5.29576e-14

Пример 0 является положительно определенной матрицей. Примеры 1 2 3 являются симметричными НЕ положительно определенными матрицами. Примеры 1 и 3 решены правильно, тогда как пример 2 терпит неудачу.

Реализация выглядит аналогично классическим реализациям CG.

ВОПРОСЫ: Есть ли в Eigen трюк, позволяющий работать с неположительно определенными матрицами? В примере 2 не соблюдаются некоторые требования, которые должны быть решены Eigen с помощью CG?

1 Ответ

0 голосов
/ 27 октября 2018

CG может использоваться для решения систем с матрицами, которые не являются положительно определенными и симметричными, следующим образом: алгоритм CG должен применяться к системе [A] T [A] x = [A] Tb , где [A] T обозначает транспонированную матрицу.В этом случае [A] T [A] является симметричным и положительно определенным, если [A] не является единственным.Недостатком является то, что [A] T [A] имеет квадрат исходного отношения условий матрицы, поэтому, если cond ( [A] ) превышает appox.10e7, итерации CG вряд ли будут сходиться вообще, и / или результирующий вектор x может не иметь каких-либо значимых цифр.Если ваша матрица достаточно «хороша» в числовом смысле, скажем, cond ( [A] ) не превышает примерно 10e3 или 10e4, вы можете ожидать, что итерации сходятся, и решение будет иметь несколько значащих цифр,Следующая публикация содержит исходные коды, реализующие такой алгоритм: https://www.amazon.com/Solution-Systems-Algebraic-Equations-Matrices/dp/0646990454

...