Eigen库的共轭梯度法(Conjugate Gradient)代码分析

ConjugateGradient类

Eigen使用ConjugateGradient类进行共轭梯度法计算。一个代码示例为(https://eigen.tuxfamily.org/dox/classEigen_1_1ConjugateGradient.html ):

“`c++
int n = 10000;
VectorXd x(n), b(n);
SparseMatrix A(n,n);
// fill A and b
ConjugateGradient, Lower|Upper> cg;
cg.compute(A);
x = cg.solve(b);
std::cout << "#iterations: " << cg.iterations() << std::endl; std::cout << "estimated error: " << cg.error() << std::endl; // update b, and solve again x = cg.solve(b);

这里提到,`solve()`函数会默认使用零向量作为初始猜测。

# IterativeSolverBase类
ConjugateGradient类是IterativeSolverBase的一个派生类。它的文档在:https://eigen.tuxfamily.org/dox/classEigen_1_1IterativeSolverBase.html

下面我们来看这个类。它给了两个控制收敛的判据,一个是相对精度,一个是迭代次数。相对精度的意思是,对于系统Ax=b,只要|Ax-b|/|b|<\epsilon就判定收敛,这里的\epsilon是设置的相对精度,其默认值为`NumTraits<Scalar>::epsilon()`,其中`Scalar`是矩阵里面的数的类型,比如`double`或者`float`。这个值,根据测试,其实就是机器精度`std::numeric_limits<Scalar>::epsilon()`。

而迭代次数,如果未设置的话,初始值为矩阵列数的两倍。

# ConjugateGradient::comput()函数

根据之前的示例代码,ConjugateGradient类先后调用两个函数,首先是`compute()`,然后是`solve()`。我们先来看`compute()`(代码见:https://eigen.tuxfamily.org/dox/IterativeSolverBase_8h_source.html ):

```c++
   template<typename MatrixDerived>
   Derived& compute(const EigenBase<MatrixDerived>& A)
   {
     grab(A.derived());
     m_preconditioner.compute(matrix());
     m_isInitialized = true;
     m_analysisIsOk = true;
     m_factorizationIsOk = true;
     m_info = m_preconditioner.info();
     return derived();
   }
</code></pre>

第4行<code>grab()</code>好像是某种复制,先不管它。这个函数实际上就是调用了preconditioner的<code>compute()</code>函数,接下来我们来看这个。

<h1>DiagonalPreconditioner::compute()函数</h1>

ConjugateGradient的默认preconditioner是这个DiagonalPreconditioner:https://eigen.tuxfamily.org/dox/classEigen_1_1DiagonalPreconditioner.html

下面我们来看这个类:https://eigen.tuxfamily.org/dox/BasicPreconditioners_8h_source.html

它的<code>compute()</code>函数的唯一作用是调用<code>factorize()</code>函数:


```c++
     template
     DiagonalPreconditioner& factorize(const MatType& mat)
     {
       m_invdiag.resize(mat.cols());
       for(int j=0; j其中关键是`m_invdiag`向量,它的长度是矩阵的列数,也就是`mat.cols()`。在这个`factorize()`完成之后,这个向量就储存了矩阵A对角线的倒数。为了看懂这个,我们需要搞清楚Eigen是怎么存的稀疏矩阵。Eigen的文档称(https://eigen.tuxfamily.org/dox/classEigen_1_1SparseMatrix.html#a5ff54ffc10296f9466dc81fa888733fd )它的默认存储是列优先,也就是CSC格式,不过可以改成CSR(行优先)格式,这个比较符合直觉。

为了方便起见,我们先来研究行优先格式。可以参照:https://zhuanlan.zhihu.com/p/37525925

假设一共有nnz个非零元素(顺带一提,nnz在Eigen的稀疏矩阵里面可以用`nonZeros()`函数获取),那么CSR格式有两个长度为nnz的数组,一个是inner index(可以用innerIndexPtr()获取),另外一个是values(可以用valuePtr()获取)。这个inner index数组里面存的,就是每个非零元素的列号(注意我们现在谈论的是CSR,Eigen库可以同时支持CSR和CSC,所以它在取名上用了inner/outer作为区分),而values里面存的就是每个非零元素的值。最后是outer index,它的长度是行数+1,代表着每一行的元素起始点。比如说第k行的第一个非零元在inner index里面的下标就是outer_index[k].

所以,这样看中间那个循环就很明显了。首先是j遍历每一行(我们假设是CSR),然后在其中用it找到相应的j行j列的元素。如果没有或者是0,那么就把`m_invdiag(j)`设成1,否则就设成这个元素的倒数。

# conjugate_gradient()函数

对于求解,我们直接看`ConjugateGradient.h`里面的`conjugate_gradient()`函数:

```c++
 template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
 EIGEN_DONT_INLINE
 void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x,
                         const Preconditioner& precond, Index& iters,
                         typename Dest::RealScalar& tol_error)
 {
   using std::sqrt;
   using std::abs;
   typedef typename Dest::RealScalar RealScalar;
   typedef typename Dest::Scalar Scalar;
   typedef Matrix<Scalar,Dynamic,1> VectorType;

   RealScalar tol = tol_error;
   Index maxIters = iters;

   Index n = mat.cols();

   VectorType residual = rhs - mat * x; //initial residual
   //注意,如果初始猜测为0,那它其实就是rhs

   RealScalar rhsNorm2 = rhs.squaredNorm();
   if(rhsNorm2 == 0) 
   {
     x.setZero();
     iters = 0;
     tol_error = 0;
     return;
   }
   const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
   RealScalar threshold = numext::maxi(RealScalar(tol*tol*rhsNorm2),considerAsZero);
   RealScalar residualNorm2 = residual.squaredNorm();
   if (residualNorm2 < threshold)
   {
     iters = 0;
     tol_error = sqrt(residualNorm2 / rhsNorm2);
     return;
   }

   VectorType p(n);
   p = precond.solve(residual);      // initial search direction

   VectorType z(n), tmp(n);
   RealScalar absNew = numext::real(residual.dot(p));  // the square of the absolute value of r scaled by invM
   Index i = 0;
   while(i < maxIters)
   {
     tmp.noalias() = mat * p;                    // the bottleneck of the algorithm

     Scalar alpha = absNew / p.dot(tmp);         // the amount we travel on dir
     x += alpha * p;                             // update solution
     residual -= alpha * tmp;                    // update residual

     residualNorm2 = residual.squaredNorm();
     if(residualNorm2 < threshold)
       break;

     z = precond.solve(residual);                // approximately solve for "A z = residual"

     RealScalar absOld = absNew;
     absNew = numext::real(residual.dot(z));     // update the absolute value of r
     RealScalar beta = absNew / absOld;          // calculate the Gram-Schmidt value used to create the new search direction
     p = z + beta * p;                           // update search direction
     i++;
   }
   tol_error = sqrt(residualNorm2 / rhsNorm2);
   iters = i;
 }

发表评论

您的电子邮箱地址不会被公开。 必填项已用*标注