在 description of the CG method在 Eigen 库中,可以找到以下语句:
This class allows to solve for A.x = b linear problems using an iterative conjugate gradient algorithm. The matrix A must be selfadjoint.
然而,在文献中,共轭梯度法通常用于实对称正定矩阵。
示例表明 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 失败。
implementation看起来类似于经典的 CG 实现。
问题:
Eigen 中是否有任何技巧可以处理非正定矩阵?
示例 2 是否不遵守某些要求以便由 Eigen 用 CG 解决?
最佳答案
CG 可用于求解具有非正定和对称矩阵的系统,方法如下: CG 算法必须应用于系统 [A]T[A]x=[A]Tb 其中 [A]T 代表转置矩阵。在那种情况下 [A]T[A] 是对称且正定的,除非 [答] 是单数。一个缺点是 [A]T[A] 具有原始矩阵条件比率的平方,因此如果 cond( [A] ) 超过约。 10e7,CG 迭代根本不可能收敛,和/或结果向量 x 可能没有任何有效数字。如果您的矩阵在数值意义上相当“好”,例如 cond( [A] ) 不超过大约 10e3 或 10e4,您可能期望迭代收敛,并且解将有几个有效数字.以下出版物包含实现此类算法的源代码:https://www.amazon.com/Solution-Systems-Algebraic-Equations-Matrices/dp/0646990454
关于matrix - 非正定矩阵的 Eigen 共轭梯度,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/53010866/