最速下降法高斯牛顿法LM共轭梯度法预条件共轭梯度法

最速下降法

原理

• 目标函数：
  最速下降法基于多元函数的梯度概念。对于一个目标函数 $f(x)$。其中 $x$ 是一个向量，例如 $x=[x_1,x_2,\cdots ,x_n{]}^T$，
  函数在某一点 $x_k$ 的梯度 $\nabla f(x_k)$ 表示函数在该点上升最快的方向。而最速下降法是沿着负梯度方向（即 $-\nabla f(x_k)$）去寻找函数的极小值点。

• 证明：

1. 方向导数
   函数$f(x)$在点$x$处沿方向$d$的变化率可以用方向导数来表达。对于可微函数，方向导数等于梯度与方向的内积。（方向导数是数值，而不是矢量。）

$$\begin{array}{rl}Df(x;d)& =\underset{\lambda \to 0^{+}}{\lim }\frac{f(x+\lambda d)-f(x)}{\lambda }\\ & =\underset{i\to 0^{+}}{\lim }\frac{f(x)+{\lambda }_{0}f(x)d+o\left(\lambda ||d||\right)-f(x)}{\lambda }\\ & =\nabla f(x{)}^{T}d+\underset{\lambda \to 0^{+}}{\lim }\frac{o(\lambda \Vert d\Vert )}{\lambda }=\nabla f(x{)}^{T}d.\end{array}$$

2. $f(x)在点$$x$处下降速度最快的方向是负梯度方向。
   $$||{\nabla }^{\mathrm{T}}f(x)d||\le ||\nabla f(x)||||d||\le ||\nabla f(x)||$$
   可知：
   $$-||\nabla f(x)||\le {\nabla }^{\mathrm{T}}f(x)d\le ||\nabla f(x)||$$
   从而可知，梯度方向是变化最快的方向，正梯度上升最快，负梯度下降最快。

3. 使用最速下降法，相邻的两个搜索方向是正交的。
   通常最速下降法，当初始值选取比较可靠，且函数凸特性比较明显时，最速下降法，高斯牛顿法均具有一次收敛性，即通过一次迭代便可达到局部极小值处。而当函数较为复杂等原因，不能通过最速下降法达到极小值时，相邻两次的迭代方向正交。
   设函数$f(x)$第$k$次迭代的参数更新量为$\Delta x_k={\lambda }_k{d}_k$，其中$d_k$为$f(x)$在k时刻变量$x_k$的负梯度方向，即，$d_k=-\nabla f(x_k)$。
   其下降步长
   $${\lambda }_k=\min \psi (\lambda )=f(x_k+\lambda d_k)$$
   对上式求导得：
   $${\psi }^{\prime }(\lambda )={\nabla }^{\mathrm{T}}f(x_k+\lambda d_k)d_k=0$$
   可知第$k+1$次迭代，迭代方向为：$\nabla f(x_k+{\lambda }_k{d}_k)$.
   从而可知：
   $${\nabla }^{\mathrm{T}}f(x_k+\lambda d_k)d_k=0\to d_{k+1}^{\mathrm{T}}{d}_k=0$$
   可知第$k$次迭代方向与第$k+1$次迭代方向相互正交。
   如下图所示：
   

• 更新规则：
  $$x_{k+1}=x_k-{\alpha }_k\nabla f(x_k)$$
  其中，${\alpha }_k$ 是步长，可以通过线搜索（如精确线搜索、非精确线搜索等方法）来确定，使得 $f(x_{k+1})<f(x_k)$。

伪代码

最速下降法通过梯度的方向来确定最陡的下降路径，即当前点的梯度方向。其伪代码如下：

输入: 目标函数f(x)，初始点x0，迭代终止条件（例如最大迭代次数max_iter或者梯度阈值epsilon）
输出: 近似极小值点x

1. 设置当前迭代点 x = x0
2. 设置迭代次数 k = 0
3. while k < max_iter 并且 ||梯度f(x)|| > epsilon  do
4.     计算目标函数f(x)在当前点x的梯度gradient = ∇f(x)
5.     通过线搜索确定步长alpha（例如可以使用精确线搜索或者非精确线搜索方法）
6.     更新当前点 x = x - alpha * gradient
7.     k = k + 1
8. end while
9. 返回x作为近似极小值点

优点与局限性

• 优点：

  • 计算简单，适用于大多数优化问题。
  • 梯度下降法在很多实际问题中有很好的表现。

• 局限性：

  • 如果目标函数的条件数较大（即矩阵的特征值差异大），算法收敛较慢。
  • 步长的选择对收敛速度有很大的影响。

C++ 示例代码

以下是一个简单的使用最速下降法求解二元函数最小值的示例代码（这里以简单的二次函数为例），假设目标函数为 $f(x,y)=x^2+y^2$：

#include <Eigen/Dense>
#include <iostream>

Eigen::VectorXd steepestDescent(const Eigen::MatrixXd& A, const Eigen::VectorXd& b, Eigen::VectorXd& x, int maxIter = 50, double tol = 1e-6) {
    Eigen::VectorXd r = b - A * x;  // Residual
    Eigen::VectorXd p = r;  // Search direction
    double rsold = r.dot(r);

    for (int i = 0; i < maxIter; ++i) {
        Eigen::VectorXd Ap = A * p;
        double alpha = rsold / p.dot(Ap);
        x += alpha * p;  // Update solution
        r -= alpha * Ap;  // Update residual

        double rsnew = r.dot(r);
        if (rsnew < tol) break;  // Convergence check
        p = r;  // Update search direction
        rsold = rsnew;
    }

    return x;
}

int main() {
    Eigen::MatrixXd A(3, 3);
    A << 4, 1, 2,
         1, 3, 0,
         2, 0, 3;
    Eigen::VectorXd b(3);
    b << 1, 2, 3;
    
    Eigen::VectorXd x = Eigen::VectorXd::Zero(3);  // Initial guess
    
    Eigen::VectorXd result = steepestDescent(A, b, x);
    std::cout << "Solution: " << result.transpose() << std::endl;
    return 0;
}

高斯牛顿法

原理

• 目标函数：
  考虑非线性最小二乘问题，目标是最小化函数：
  $$F(x)=\frac{1}{2}\sum _{i=1}^m{r}_i^2(x)$$
  其中，$r_i(x)$ 是残差函数（例如在拟合数据时，观测值与模型预测值的差值）。
  对 $F(x)$ 求一阶导数（梯度）并令其为 0，利用泰勒展开式将残差函数 $r_i(x)$ 在当前迭代点 $x_k$ 处展开：
  $$r_i(x_{k+1})\approx r_i(x_k)+J_i(x_k)\Delta x$$
  其中，$J_i(x_k)$ 是 $r_i(x)$ 在 $x_k$ 处的雅可比矩阵。将其代入 $F(x)$ 的梯度表达式，并令梯度为 0，经过推导可得高斯牛顿法的迭代公式：
  $$(J^T(x_k)J(x_k))\Delta x=-J^T(x_k)r(x_k)$$
  其中，$J(x_k)$ 是所有残差函数的雅可比矩阵组成的矩阵，$r(x_k)$ 是残差向量，求解 $\Delta x$ 后更新 $x_{k+1}=x_k+\Delta x$。
• 更新规则：
  $$x_{k+1}=x_k-{\left(J^{T}J\right)}^{-1}{J}^{T}r(x_k)$$

在处理非线性最小二乘拟合等问题时，普通的梯度下降法收敛速度可能较慢。同最速下降法，当初始值以及约束方程足够理想时，牛顿欧拉法也能够经历一次迭代到达极小值。高斯牛顿法利用了目标函数是残差平方和的特殊结构，通过近似二阶导数信息（用雅可比矩阵相关运算来近似海森矩阵），能在一定程度上加快收敛速度，更高效地找到使残差平方和最小的参数值。高斯-牛顿法通过泰勒展开对目标函数进行二次近似，利用雅可比矩阵 $J$ 来线性化非线性函数。

伪代码

输入: 残差函数向量r(x)，初始参数向量x0，迭代终止条件（例如最大迭代次数max_iter或者参数更新量阈值epsilon）
输出: 使残差平方和最小的参数向量x

1. 设置当前迭代点 x = x0
2. 设置迭代次数 k = 0
3. while k < max_iter 并且 ||参数更新量Δx|| > epsilon do
4.     计算残差向量 r = r(x)
5.     计算雅可比矩阵 J = J(x) （r(x)中每个残差函数关于参数向量x的偏导数组成的矩阵）
6.     计算近似海森矩阵 H = J^T * J
7.     计算梯度向量 g = J^T * r
8.     求解线性方程组 H * Δx = -g 得到参数更新量Δx（例如使用矩阵分解等方法求解）
9.     更新当前点 x = x + Δx
10.    k = k + 1
11. end while
12. 返回x作为最终参数向量

优点与局限性

• 优点：

  • 对于许多非线性最小二乘问题收敛较快。
  • 算法在处理数据拟合和非线性回归等问题时效果较好。

• 局限性：

  • 需要计算雅可比矩阵和二阶导数，可能导致计算量较大。
  • 当雅可比矩阵条件数较差时，算法可能收敛缓慢或发散。

C++ 示例代码

以下是一个简单的高斯牛顿法求解非线性最小二乘问题的示例代码框架（这里以简单的拟合曲线问题为例，假设残差函数形式简单）：

#include <Eigen/Dense>
#include <iostream>

Eigen::VectorXd gaussNewton(const Eigen::MatrixXd& A, const Eigen::VectorXd& b, Eigen::VectorXd& x, int maxIter = 50, double tol = 1e-6) { 

   // H: J^T*T r:-J^T*b
    Eigen::MatrixXd H = A.transpose() * A;
    Eigen::VectorXd r = b - A*x;
   //b - Ax
    for (int i = 0; i < maxIter; ++i) {

        Eigen::VectorXd Jb = A.transpose() * r;

        Eigen::VectorXd dx = H.ldlt().solve(Jb);
        x += dx;  // Update solution
        r -= A*dx;  // Update residual

        double rsnew = r.dot(r);
        if (rsnew < tol) break;  // Convergence check
    }   

    return x;
}

int main() {
    Eigen::MatrixXd A(4, 3); 
    A << 4, 1, 2,
         1, 3, 0,
         1.1, 2.9, 0.1,
         2, 0, 3;
    Eigen::VectorXd b(4);
    b << 1, 2, 2.1, 3;

    Eigen::VectorXd x = Eigen::VectorXd::Zero(3);  // Initial guess
    
    Eigen::VectorXd result = gaussNewton(A, b, x); 
    std::cout << "Solution: " << result.transpose() << std::endl;
    return 0;
}

//g++ gauss1.cpp `pkg-config eigen3 --libs --cflags`

LM 法（Levenberg-Marquardt 法）

原理

LM 法是高斯牛顿法的一种改进。在高斯牛顿法中，有时候近似的海森矩阵 $(J^{T}J)$ 可能是奇异的或者条件数很差，导致算法不稳定或者收敛缓慢。LM 法引入了一个阻尼因子 $\lambda$，对迭代公式进行修正，其迭代公式如下：

$$(J^T(x_k)J(x_k)+\lambda I)\Delta x=-J^T(x_k)r(x_k)$$

• 更新规则：
  $$x_{k+1}=x_k-{\left(J^{T}J+\lambda I\right)}^{-1}{J}^{T}r(x_k)$$

其中，$I$ 是单位矩阵。通过调整 $\lambda$ 的大小来平衡接近高斯牛顿法（$\lambda$ 较小时）的快速收敛特性和类似最速下降法（$\lambda$ 较大时）的全局收敛性，保证算法在不同情况下都能较好地收敛。

高斯牛顿法虽然在很多情况下对非线性最小二乘问题有较好的表现，但存在对初始值敏感、在某些情况下可能不收敛等问题。LM 法通过引入阻尼因子来改善算法的稳定性和收敛性，使其能够处理更广泛的非线性最小二乘问题，对初始值的要求相对没那么苛刻，在实际的曲线拟合、参数估计等应用中更加实用可靠。
莱文贝格-马夸尔特算法通过增加阻尼因子 ( \lambda ) 来控制步长，避免在条件数较差的情况下发散。

优点与局限性

• 优点：

  • 适用于非线性最小二乘问题。
  • 在非线性问题中比高斯-牛顿法更稳定，尤其是在接近最优解时。

• 局限性：

  • 需要调整阻尼因子，选择不当可能影响收敛速度。

伪代码

输入: 残差函数向量r(x)，初始参数向量x0，初始阻尼因子lambda，迭代终止条件（例如最大迭代次数max_iter或者参数更新量阈值epsilon），阻尼因子调整倍数lambda_factor
输出: 使残差平方和最小的参数向量x

1. 设置当前迭代点 x = x0
2. 设置迭代次数 k = 0
3. while k < max_iter 并且 ||参数更新量Δx|| > epsilon do
4.     计算残差向量 r = r(x)
5.     计算雅可比矩阵 J = J(x)
6.     计算近似海森矩阵 H = J^T * J + lambda * I （I为单位矩阵）
7.     计算梯度向量 g = J^T * r
8.     求解线性方程组 H * Δx = -g 得到参数更新量Δx（例如使用矩阵分解等方法求解）
9.     计算新的目标函数值 new_f = 0.5 * ||r(x + Δx)||^2
10.    计算当前目标函数值 current_f = 0.5 * ||r(x)||^2
11.    if new_f < current_f then
12.        lambda = lambda / lambda_factor  // 减小阻尼因子，更接近高斯牛顿法
13.        更新当前点 x = x + Δx
14.    else
15.        lambda = lambda * lambda_factor  // 增大阻尼因子，增强稳定性
16.    end if
17.    k = k + 1
18. end while
19. 返回x作为最终参数向量

C++ 示例代码

以下是一个简单的 LM 法求解非线性最小二乘问题的示例代码（基于前面高斯牛顿法代码进行修改）：

#include <Eigen/Dense>
#include <iostream>

Eigen::VectorXd gaussNewton(const Eigen::MatrixXd& A, const Eigen::VectorXd& b, Eigen::VectorXd& x, int maxIter = 50, double tol = 1e-6) { 

    double lambda = 0.1;
   // H: J^T*T r:-J^T*b
    Eigen::MatrixXd H = A.transpose() * A + lambda * Eigen::MatrixXd::Identity(A.cols(), A.cols());
    Eigen::VectorXd r = b - A*x;
   //b - Ax
    for (int i = 0; i < maxIter; ++i) {

        Eigen::VectorXd Jb = A.transpose() * r;

        Eigen::VectorXd dx = H.ldlt().solve(Jb);
        x += dx;  // Update solution
        r -= A*dx;  // Update residual

        double rsnew = r.dot(r);
        if (rsnew < tol) break;  // Convergence check
    }

    return x;
}

int main() {
    Eigen::MatrixXd A(4, 3);
    A << 4, 1, 2,
         1, 3, 0,
         1.1, 2.9, 0.1,
         2, 0, 3; 
    Eigen::VectorXd b(4);
    b << 1, 2, 2.1, 3;

    Eigen::VectorXd x = Eigen::VectorXd::Zero(3);  // Initial guess
        
    Eigen::VectorXd result = gaussNewton(A, b, x);
    std::cout << "Solution: " << result.transpose() << std::endl;
    return 0;
}

//g++ lm.cpp `pkg-config eigen3 --libs --cflags`

共轭梯度法

提出

共轭梯度法（Conjugate Gradient Method，CG法）最早由 Magnus Hestenes 和 Eduard Stiefel 在1952年提出，旨在求解大型线性系统 $A\mathbf{x}=\mathbf{b}$，其中 $A$ 是对称正定矩阵。这类问题在科学计算和工程优化中经常出现，尤其是求解大规模稀疏系统时，直接解法（如高斯消元法或LU分解）由于计算和存储需求过高而变得不可行。共轭梯度法是针对这一问题的高效迭代方法。
与梯度下降法不同，梯度下降法是沿着目标函数的梯度方向进行更新，但这并不总是最有效的。梯度下降法在某些情况下，尤其是在具有“长尾（梯度过小）”或“陡峭（梯度过大）”梯度的情况下，收敛速度较慢。而共轭梯度法通过利用矩阵的特征值分解，选择更加合适的方向，能显著提高收敛效率。

原理

用于求解大型稀疏对称正定系统 $Ax=b$（其中 $A$ 是对称正定矩阵），共轭梯度法的核心思想是构造一组关于矩阵A的共轭向量 $p_0,p_1,\cdots ,p_{n-1}$（满足 $p_i^{T}A{p}_j=0$，$i\ne j$），并沿着这些共轭方向逐步搜索来逼近方程组的解。

共轭梯度方向是上一次迭代方向与当前梯度方向的一个融合，是有限次的迭代方法。适合用于大型稀疏系统，因为其不需要计算hessian矩阵的逆，而计算的是$p_i^{T}A{p}_j$的逆。

对于n维变量求解，共轭梯度法可以在 n 次迭代内，解会完全收敛，也就是说，经过 n 次迭代后，解的误差已经足够小，达到精度要求。这是因为：
每次迭代时，共轭梯度法都会在一个新的共轭方向上进行更新，随着迭代次数增加，解逐渐逼近真实解。
共轭方向在几何上形成了一个与目标函数的特征值相关的正交基，这使得每个方向的贡献是独立的，因此最坏情况下，最多n 次迭代就能得到精确解。可以理解为分别在每一个维度上逐个优化，从而得到最终的极小值。当引入误差等特殊情况，则需要更多次的迭代。

• 目标问题：
  $$\underset{x}{\min }f(x)=\underset{x}{\min }\left(\frac{1}{2}||r(x)||\right)=\underset{x}{\min }\left(\frac{1}{2}(b-Ax{)}^{\mathrm{T}}(b-Ax)\right)=\underset{x}{\min }\left(\frac{1}{2}{x}^{T}Ax-b^{T}x\right)$$
  从而可知，其梯度：
  $$\nabla f(x)=g(x)=Ax-b=-r(x)\to r(x)=-g(x)$$

• 关于矩阵共轭的几何意义：
  满足 $p_i^{T}A{p}_j=0$的向量，其中A为正定对称矩阵，$p_i$和$p_j$成为关于矩阵A的共轭向量，此时，$p_i$与经过A变换后的$p_j$相互正交。

• 更新规则：
  变量更新：
  $$x_{k+1}=x_k+{\alpha }_k{p}_k$$ 残差更新：（残差更新即是梯度更新$r=g$）$$r_{k+1}=r_k-{\alpha }_{k}A{p}_k$$方向更新:$${\beta }_k=\frac{r_{k+1}^T{r}_{k+1}}{r_k^T{r}_k}$$$$p_{k+1}=r_{k+1}+{\beta }_k{p}_k$$

• 性质：
  对于该正定二次函数，使用一维搜索的共轭下降法，下列关系成立
  $$p_i^{\mathrm{T}}A{p}_j=0,(j=1,2,...,i-1)$$
  $$g_i^{\mathrm{T}}{g}_j=0,(j=1,2,...,i-1)$$即$$r_i^{\mathrm{T}}{r}_j=0,(j=1,2,...,i-1)$$
  $$g_i^{\mathrm{T}}{p}_i=-g_i^{\mathrm{T}}{g}_i\to r_i^{\mathrm{T}}{p}_i=r_i^{\mathrm{T}}{r}_i$$

• 推导：

• ${\alpha }_k$的计算
  $$-\nabla f(x)=b-Ax=r(x)$$$$\underset{{\alpha }_k\ge 0}{\min }\psi (x)=f(x_k+{\alpha }_k{p}_k)$$
  $$\begin{gathered} {\psi }^{\prime }({\alpha }_k)={\nabla }^{\mathrm{T}}f(x_k+{\alpha }_k{p}_k)p_k=0\Rightarrow g_{k+1}^{\mathrm{T}}{p}_k=0\Rightarrow r_{k+1}^{\mathrm{T}}{p}_k=0 \\ {\left(b-A(x_k+{\alpha }_k{p}_k)\right)}^{\mathrm{T}}{p}_k=0 \\ (b-A{x}_k-{\alpha }_{k}A{p}_k{)}^{\mathrm{T}}{p}_k=(r_k-{\alpha }_{k}A{p}_k{)}^{\mathrm{T}}{p}_k=0 \end{gathered}$$

得：
$${\alpha }_k=\frac{r_k^{\mathrm{T}}{p}_k}{p_k^{\mathrm{T}}A{p}_k}=\frac{r_k^{\mathrm{T}}{r}_k}{p_k^{\mathrm{T}}A{p}_k}$$

*${\beta }_k$的计算：新k+1的方向为：k+1的负梯度方向与上一次迭代方向的融合
$$p_{k+1}=-g_{k+1}+{\beta }_k{d}_k=r_{k+1}+{\beta }_k{d}_k$$
令：
$$\begin{array}{rl}{p}_k^{T}A{p}_{k+1}& =p_k^{T}A(r_{k+1}+{\beta }_k{p}_k)\\ & =p_k^{\mathrm{T}}A{r}_{k+1}+{\beta }_k{p}_k^{\mathrm{T}}A{p}_k=0\\ & \Rightarrow {\beta }_k=-\frac{p_k^{T}A{r}_{k+1}}{p_k^{T}A{p}_k}\end{array}$$
将$$r_{k+1}=r_k-{\alpha }_{k}A{p}_k$$带入上式，可得：
$${\beta }_k=\frac{r_{k+1}^T{r}_{k+1}}{r_k^T{r}_k}$$

如果令：
$$\begin{array}{rl}{p}_{k+1}^{T}A{p}_k& =(r_{k+1}+{\beta }_k{p}_k{)}^{\mathrm{T}}A{p}_k\\ & =r_{k+1}^{\mathrm{T}}A{p}_k+{\beta }_k{p}_k^{\mathrm{T}}A{p}_k=0\\ & \Rightarrow {\beta }_k=-\frac{r_{k+1}^{\mathrm{T}}A{p}_k}{p_k^{T}A{p}_k}\end{array}$$

以上三式等价，其中使用$r_k$及$r_{k+1}$来组合${\beta }_k$效率最高，因而是最常用的迭代方式。

共轭梯度法通过构造共轭方向，并逐步更新解，最终收敛于最优解。每次更新时，新的搜索方向和之前的方向共轭。其迭代过程如下：

初始化：选择初始解向量 $x_0$，计算初始残差向量 $r_0=b-A{x}_0$，令 $p_0=r_0$。
迭代步骤：对于 $k=0,1,\cdots ,n-1$，计算步长：
$${\alpha }_k=\frac{r_k^T{r}_k}{p_k^{T}A{p}_k}$$
更新解向量：
$$x_{k+1}=x_k+{\alpha }_k{p}_k$$
更新残差向量：
$$r_{k+1}=r_k-{\alpha }_{k}A{p}_k$$
计算新的共轭方向：
$${\beta }_k=\frac{r_{k+1}^T{r}_{k+1}}{r_k^T{r}_k}$$$$p_{k+1}=r_{k+1}+{\beta }_k{p}_k$$

在求解大型稀疏线性方程组时（例如在数值模拟、有限元分析等领域经常遇到），直接求解（如 LU 分解等方法）计算量过大，存储需求也很高。共轭梯度法利用了矩阵的对称正定性质，通过迭代的方式逐步逼近精确解，并且具有收敛速度较快、存储需求相对较低（只需存储几个向量即可进行迭代计算）的优点，特别适用于处理大规模问题。

伪代码

输入: 系数矩阵A（对称正定矩阵），右端向量b，初始解向量x0，迭代终止条件（例如最大迭代次数max_iter或者残差范数阈值epsilon）
输出: 线性方程组Ax = b的解向量x

1. 设置当前迭代点 x = x0
2. 计算初始残差向量 r = b - A * x
3. 设置搜索方向向量 p = r
4. 设置初始残差的平方范数 rsold = r^T * r
5. 设置迭代次数 k = 0
6. while k < max_iter 并且 √(rsnew) > epsilon do
7.     计算矩阵向量乘积 Ap = A * p
8.     计算步长 alpha = rsold / (p^T * Ap)
9.     更新当前点 x = x + alpha * p
10.    更新残差向量 r = r - alpha * Ap
11.    计算新的残差的平方范数 rsnew = r^T * r
12.    计算共轭系数 beta = rsnew / rsold
13.    更新搜索方向向量 p = r + beta * p
14.    更新 rsold = rsnew
15.    k = k + 1
16. end while
17. 返回x作为线性方程组的解

C++ 示例代码

以下是一个简单的共轭梯度法求解线性方程组的示例代码（假设矩阵 $A$ 和向量 $b$ 已知）：

#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <iostream>

Eigen::VectorXd conjugateGradient(const Eigen::SparseMatrix<double>& A, const Eigen::VectorXd& b, Eigen::VectorXd& x, int maxIter = 1000, double tol = 1e-6) {
    Eigen::VectorXd r = b - A * x;  // Residual
    Eigen::VectorXd p = r;  // Search direction
    double rsold = r.dot(r);

    for (int i = 0; i < maxIter; ++i) {
        Eigen::VectorXd Ap = A * p;
        double alpha = rsold / p.dot(Ap);
        x += alpha * p;  // Update solution
        r -= alpha * Ap;  // Update residual

        double rsnew = r.dot(r);
        if (rsnew < tol) break;  // Convergence check
        p = r + (rsnew / rsold) * p;  // Update direction
        rsold = rsnew;
    }

    return x;
}

int main() {
    Eigen::SparseMatrix<double> A(3, 3);
    A.insert(0, 0) = 4;
    A.insert(1, 1) = 4;
    A.insert(2, 2) = 4;
    A.insert(0, 1) = 1;
    A.insert(1, 0) = 1;
    Eigen::VectorXd b(3);
    b << 1, 2, 3;

    Eigen::VectorXd x = Eigen::VectorXd::Zero(3);  // Initial guess
    
    Eigen::VectorXd result = conjugateGradient(A, b, x);
    std::cout << "Solution: " << result.transpose() << std::endl;
    return 0;
}
//g++ cg.cpp `pkg-config eigen3 --libs --cflags`

预条件共轭梯度法

原理

预条件共轭梯度法是在共轭梯度法的基础上改进而来。共轭梯度法是求解对称正定系数的线性方程组的极为有效的方法，并且指出：对于n阶线性方程组，通常最多n+1次迭代可以获得收敛。在实践中，需要求解的经常是非常大，极端大的线性方程组，此时，采用共轭梯度法，有时可能需要上万次的迭代才能收敛，虽然每次收敛的计算量并不大，但是整体求解也会花费较多时间。
实际上，共轭梯度法能否快速收敛与系数矩阵A密切相关。对于共轭梯度法存在以下定理：如果$A=I+B$为 $n\ast n$ 的对称正定矩阵，且$rank(B)=r$，则共轭梯度法至多 $r+1$ 步收敛。其中 $I$ 指单位矩阵，rank是矩阵的秩。

从该定理可知，共轭梯度法能否快速收敛，主要取决于“A”是否“接近于”单位矩阵。特别地，当A就是单位矩阵时，$rank(B)=0$，此时一步即可收敛，当然这是很显然的。基于上述定理，多种基于共轭梯度法法衍生的预处理共轭梯度法（PCG）得以广泛地应用。

其核心是引入一个预条件矩阵 $M$，对原线性方程组 $Ax=b$ 进行变换，使得新的方程组 $M^{-1}Ax=M^{-1}b$ 的条件数得到改善（通常希望变换后的矩阵更接近单位矩阵，这样迭代收敛会更快），从而构成新的方程形式$\hat{A}x=\hat{b},\hat{A}=M^{-1}A,\hat{b}=M^{-1}b$。

• 目标问题：
  $$A\mathbf{x}=\mathbf{b}$$

• 更新规则：
  $$\mathbf{z}=M^{-1}\mathbf{r}$$
  之后使用共轭梯度法求解预处理后的系统。

在迭代过程中，将共轭梯度法中的一些计算（如计算步长、更新残差和共轭方向等）涉及的矩阵运算都用预条件矩阵进行相应的修改。例如，计算步长时变为：

$${\alpha }_k=\frac{(M^{-1}{r}_k{)}^T{r}_k}{(M^{-1}{p}_k{)}^{T}A{p}_k}$$

其他步骤也做类似调整，通过合适的预条件矩阵选择，可以加快收敛速度。

共轭梯度法虽然在很多情况下能有效求解线性方程组，但对于一些条件数较大的矩阵（即矩阵的病态程度较高），收敛速度可能会变得很慢。预条件共轭梯度法通过引入预条件矩阵来改善矩阵的条件数，从而进一步提高算法的收敛速度，使其在处理一些较难收敛的线性方程组问题时更具优势，能够更快地得到满足精度要求的解。

优点与局限性

• 优点：

  • 通过选择合适的预条件器加速共轭梯度法的收敛速度。

• 局限性：

  • 预条件器的选择直接影响算法性能，选择不当可能导致收敛速度降低。

伪代码

输入: 系数矩阵A，右端向量b，初始解向量x0，迭代终止条件（例如最大迭代次数max_iter或者残差范数阈值epsilon），预条件矩阵构造函数M(A)
输出: 线性方程组Ax = b的解向量x

1. 设置当前迭代点 x = x0
2. 计算初始残差向量 r = b - A * x
3. 计算预条件后的向量 z = M(A)^(-1) * r （通过求解线性方程组M(A) * z = r实现）
4. 设置搜索方向向量 p = z
5. 设置初始残差的平方范数 rsold = r^T * r
6. 设置迭代次数 k = 0
7. while k < max_iter 并且 √(rsnew) > epsilon do
8.     计算矩阵向量乘积 Ap = A * p
9.     计算步长 alpha = (z^T * r) / (p^T * Ap) （这里z是预条件后的向量）
10.    更新当前点 x = x + alpha * p
11.    更新残差向量 r = r - alpha * Ap
12.    计算预条件后的新向量 new_z = M(A)^(-1) * r
13.    计算新的残差的平方范数 rsnew = r^T * r
14.    计算共轭系数 beta = rsnew / rsold
15.    更新搜索方向向量 p = new_z + beta * p
16.    更新 rsold = rsnew
17.    k = k + 1
18. end while
19. 返回x作为线性方程组的解

C++ 示例代码

以下是一个简单的预条件共轭梯度法求解线性方程组的示例代码（这里简单地以对角矩阵作为预条件矩阵示例，实际中预条件矩阵选择更复杂）：

#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <iostream>

Eigen::VectorXd preconditionedConjugateGradient(
    const Eigen::SparseMatrix<double>& A,
    const Eigen::VectorXd& b,
    const Eigen::VectorXd& preconditioner,
    Eigen::VectorXd& x,
    int max_iterations = 1000,
    double tolerance = 1e-6) {

    Eigen::VectorXd r = b - A * x;
    Eigen::VectorXd z = r.cwiseQuotient(preconditioner);
    Eigen::VectorXd p = z;

    double rsold = r.dot(z);
    for (int i = 0; i < max_iterations; ++i) {
        Eigen::VectorXd Ap = A * p;
        double alpha = rsold / p.dot(Ap);
        x += alpha * p;
        r -= alpha * Ap;

        if (r.norm() < tolerance) break;

        z = r.cwiseQuotient(preconditioner);
        double rsnew = r.dot(z);
        p = z + (rsnew / rsold) * p;
        rsold = rsnew;
    }

    return x;
}

int main() {
    Eigen::SparseMatrix<double> A(3, 3);
    A.insert(0, 0) = 4;
    A.insert(1, 1) = 4;
    A.insert(2, 2) = 4;
    A.insert(0, 1) = 1;
    A.insert(1, 0) = 1;
    Eigen::VectorXd b(3);
    b << 1, 2, 3;

    Eigen::VectorXd preconditioner(3);
    preconditioner << 1, 1, 1;  // Diagonal preconditioner

    Eigen::VectorXd x = Eigen::VectorXd::Zero(3);  // Initial guess

    Eigen::VectorXd result = preconditionedConjugateGradient(A, b, preconditioner, x);
    std::cout << "Solution: " << result.transpose() << std::endl;
    return 0;
}

//g++ pcg.cpp `pkg-config eigen3 --libs --cflags`

参考链接

https://github.com/QiangLong2017/Optimization-Theory-and-Algorithm
https://www.bilibili.com/video/BV1e64y1Y7Sr/?spm_id_from=333.788
https://blog.csdn.net/weixin_41469272/article/details/121994485
https://mp.weixin.qq.com/s/sRDDSZrCmUWHvg7XxERJaw
https://blog.csdn.net/weixin_43940314/article/details/121125847