ConjugateGradient类

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

int n = 10000;
VectorXd x(n), b(n);
SparseMatrix<double> A(n,n);
// fill A and b
ConjugateGradient<SparseMatrix<double>, 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 ）：

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

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

DiagonalPreconditioner::compute()函数

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

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

它的compute()函数的唯一作用是调用factorize()函数：

     template<typename MatType>
     DiagonalPreconditioner& factorize(const MatType& mat)
     {
       m_invdiag.resize(mat.cols());
       for(int j=0; j<mat.outerSize(); ++j)
       {
         typename MatType::InnerIterator it(mat,j);
         while(it && it.index()!=j) ++it;
         if(it && it.index()==j && it.value()!=Scalar(0))
           m_invdiag(j) = Scalar(1)/it.value();
         else
           m_invdiag(j) = Scalar(1);
       }
       m_isInitialized = true;
       return *this;
     }

其中关键是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()函数：

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