数值计算：线性方程组的迭代解法 02 非静态迭代法

上一节记录了线性方程组的静态迭代方法
本节主要介绍非静态迭代方法中的最速下降法与共轭梯度法。

非静态迭代法的基本原理

基于变分原理，考虑以下函数

\[\begin{align} \Phi(\vec x)=\frac{1}{2}\vec{x}\cdot\boldsymbol A\vec{x}-\vec{x}\cdot \vec{b} \end{align} \]

若\(\boldsymbol A\)对称，则函数的梯度为

\[\begin{align} \nabla \Phi(\vec{x})=\boldsymbol A\vec{x}-\vec{b} \end{align} \]

若\(\boldsymbol A\)对称正定的，则解方程\(\boldsymbol A\vec{x}=\vec{b}\)等价于求函数\(\Phi(\vec{x})\)的最小值，令\(\vec{r}=\vec{b}-\boldsymbol A\vec{x}\)，则\(\vec{r}\)表示函数\(\Phi(\vec{x})\)的负梯度方向。

假设现已经得到\(k-1\)步迭代的解\(\vec{x}_{(k-1)}\)，在第\(k\)步中，沿着\(\vec{p_k}\)方向进行搜素，即

\[\begin{align} \vec{x}_{k}=\vec{x}_{k-1}+\alpha_{k}\cdot \vec{p_{k}} \end{align} \]

带入\(\Phi(\vec{x})\)函数表达式可以得到

\[\begin{align} \Phi(\vec{x}_{k})=\frac{1}{2}\vec{p_k}\cdot \boldsymbol A\vec{p_k}\alpha_k^2-\left(\vec{b}-\boldsymbol A\vec{x}_{k-1}\right)\cdot\vec{p_k}\alpha_k+\frac{1}{2}\vec{x}_{k-1}\cdot\left(\boldsymbol A\vec{x}_{k-1}-\vec{b}\right) \end{align} \]

可以视作\(\Phi(\alpha_k)\)函数取最小值时，\(\frac{\partial\Phi(\alpha_k)}{\partial \alpha_k}=0\)，则可以得到

\[\begin{equation} \begin{aligned} \alpha_k=\frac{\left(\vec{b}-\boldsymbol A\vec{x}_{k-1}\right)\cdot \vec{p_k}}{\vec{p_k}\cdot\boldsymbol A\vec{p_k}}&=\frac{\vec{r}_{k-1}\cdot \vec{p_k}}{\vec{p_k}\cdot\boldsymbol A\vec{p_k}}\\ \vec{r}_{k}&=\vec{b}-\boldsymbol A\vec{x}_{k}\\ &=\vec{b}-\boldsymbol A\cdot(\vec{x}_{k-1}+\alpha_k\vec{p}_k)\\ &=\vec{r}_{k-1}-\alpha_k\boldsymbol A\vec{p_k}\\ \end{aligned} \end{equation} \]

由以上推导，可以发现，迭代过程中最核心的两个部分分别是1、迭代搜索方向\(\vec{p}_k\)，2、迭代步长\(\alpha_k\)，其中迭代步长\(\alpha_k\)的最优解已经确定，它与负梯度方向\(\vec{r}_k\)与迭代搜索方向\(\vec{p}_k\)相关。

最速下降法

（steepest descent method）

对于最速下降法，迭代搜索方向简单得选择为负梯度方向，即：\(\vec{p_k}=\vec{r}_{k}\)，按照下降最快的方向去搜索。

A=[4 3 0;3 4 -1;0 -1 4];
b=[3;5;-5];
x0=[0;0;0];

[x,iter,Residual]=Steepest_Descent(A,b,x0,1e-6,100);

function [x,iter,Residual]=Steepest_Descent(A,b,x0,epsilon,iter_max)
%系数矩阵A,b,初始值x0,误差限制epsilon，最大迭代步数iter_max
iter=0;
x=x0;
r=b-A*x0;
while sqrt(dot(r,r))>=epsilon&&iter<iter_max
    iter=iter+1;
    Ar=A*r;
    alpha=dot(r,r)/dot(r, Ar);
    x_new=x+alpha*r;
    Residual=norm(x-x_new,2);
    x=x_new;
    r=r-alpha*Ar;
end
end

共轭梯度法

（conjugate gradient method）

最速下降法简单地选择负梯度方向为搜索方向，搜索效率较低，因为\(\{\vec{p}_k\}\)之间并没有相对限制关系，如果能够对迭代搜索方向加以限制，回大大提高收敛效率。

共轭梯度法是一种求解大型稀疏对称正定方程组的有效方法。搜索方向不再是简单的负梯度方向，共轭梯度法要求\(\{\vec{p}_k\}_{i=0}^{n}\)之间相互独立，理论上，对于\(n\)维空间，可以在\(n\)步之内得到收敛结果。

现在引入共轭的定义

1、若向量\(\vec{x},\vec{y}\ne0\)，且\(\vec{x},\vec{y}\)关于\(\boldsymbol A\)共轭，则有\(\vec{x}^T\boldsymbol A\vec{y}=0\)；

2、若\(\boldsymbol A\)对称正定，\(\{d_i\}\)关于\(\boldsymbol A\)共轭，则\(\{d_i\}\)相互独立；

---

同样的迭代格式

\[\begin{align} \vec{x}_{k}&=\vec{x}_{k-1}+\alpha_{k}\vec{p}_{k}\\ \end{align} \]

假设现在已经存在线性无关的向量\(\vec{p}_1,\vec{p}_2,\cdots,\vec{p}_{k}\)关于\(\boldsymbol A\)共轭，现在要找到\(\vec{p}_{k+1}\)方向，可以认为

\[\begin{align} \vec{p}_{k+1}=\beta_k \vec{p}_k+\vec{r}_k \end{align} \]

迭代公式含义可以理解为，\(\vec{p}_{k+1}\)等于\(\vec{r}_{k}\)减去\(\vec{p}_{k}\)的分量，\(\beta\)好比是\(\vec{r}_{k}\)在\(\vec{p}_{k}\)上的负投影.

\(\vec{p}_{k+1}\)与\(\vec{p}_1,\vec{p}_2,\cdots,\vec{p}_{k}\)关于\(\boldsymbol A\)共轭，等效于\(\vec{p}_{k+1}\)与\(\vec{p}_{k}\)关于\(\boldsymbol A\)共轭。则存在

\[\begin{equation} \begin{aligned} \vec{p}_{k+1}^T\boldsymbol A\vec{p}_k&=0\\ (\beta_k \vec{p}_k+\vec{r}_k)^T\boldsymbol A\vec{p}_k&=0\\ \beta_k&=-\frac{\vec{r}_k^T\boldsymbol A\vec{p}_k}{\vec{p}_k^T\boldsymbol A\vec{p}_k} \end{aligned} \end{equation} \]

之前推导的\(\alpha_k\)可以作进一步简化

\[\begin{equation} \begin{aligned} \alpha_k&=\frac{\vec{r}_{k-1}\cdot \vec{p_k}}{\vec{p_k}\cdot\boldsymbol A\vec{p_k}}\\ &=\frac{\vec{r}_{k-1}\cdot ({\beta_k \vec p_{k-1}+\vec{r}_{k-1}})}{\vec{p_k}\cdot\boldsymbol A\vec{p_k}}\\ &=\frac{\vec{r}_{k-1}\cdot\vec{r}_{k-1}}{\vec{p_k}\cdot\boldsymbol A\vec{p_k}}\\ \end{aligned} \end{equation} \]

到此，便是完整的梯度下降法，详细的数学推导公式见参考文章

---

伪代码:

\[\begin{equation}\begin{aligned} &\boldsymbol x_0=0\\ &\boldsymbol r_0=b-\boldsymbol Ax_0\\ &\boldsymbol p_1=r_0\\ &\text{Computing} \quad \boldsymbol A\boldsymbol p_1\\ &\text{for} \quad k=1,2,\cdots\\ &\quad \text{precise}=||r_{k-1}||_2/||b||_2\\ &\quad \text{if} (\text{precise<}\varepsilon) \text {break;}\\ &\quad a_k=\frac{(\boldsymbol r_{k-1},\boldsymbol r_{k-1})}{(\boldsymbol p_k,\boldsymbol A\boldsymbol p_k)}\\ &\quad \boldsymbol x_k=\boldsymbol x_{k-1}+\alpha_k\boldsymbol p_k\\ &\quad \boldsymbol r_k=\boldsymbol r_{k-1}-\alpha_k\boldsymbol A p_k\\ &\quad\text{Computing} \quad \boldsymbol A\boldsymbol r_k\\ &\quad \beta_k=-\frac{(\boldsymbol p_k,\boldsymbol {A r}_k)}{(\boldsymbol p_k,\boldsymbol {A p}_k)} \\ &\quad \boldsymbol p_{k+1}=\boldsymbol r_k+\beta_k \boldsymbol p_k\\ &\quad \boldsymbol A \boldsymbol p_{k+1}=\boldsymbol A \boldsymbol r_{k}+\beta_k\boldsymbol A \boldsymbol p_{k} \\&\text{end} \end{aligned}\end{equation} \]

---

clc;clear;
A=[4 3 0;3 4 -1;0 -1 4];
b=[3;5;-5];
x0=[0;0;0];

[x,iter,Residual]=Conjugate_Gradient(A,b,x0,1e-6,100);

function [x,iter,Residual]=Conjugate_Gradient(A,b,x0,epsilon,iter_max)
%系数矩阵A,b,初始值x0,误差限制epsilon，最大迭代步数iter_max
iter=0;
x=x0;
r=b-A*x0;
p=r;
Ap=A*p;
while sqrt(dot(r,r))>=epsilon&&iter<iter_max
    iter=iter+1;
    alpha=dot(r,r)/dot(p,Ap);
    x_new=x+alpha*p;
    Residual=norm(x-x_new,2);
    x=x_new;
    r=r-alpha*Ap;
    Ar=A*r;
    beta=-dot(p,Ar)/dot(p,Ap);
    p=r+beta*p;
    Ap=Ar+beta*Ap;
end
end

除此之外，还有基于共轭梯度法的改进方法，等用到了再来续写。

双共轭梯度法

（Bi-conjugate gradient method）

显然，共轭梯度法针对的是对称正定矩阵，而双共轭梯度法是针对于非对称矩阵的一类改进方法

稳定双共轭梯度法

（Bi-conjugate gradient stabilized method）