（五）Richardson 迭代法

1. Richardson 迭代法

求解线性方程组：
$$\mathbf{A}\vec{x}=\vec{b}(\ast )$$
其中，$\mathbf{A}$ 若为 $n$ 阶对称矩阵则要求正定，采用迭代格式：
$${\vec{x}}_{i+1}={\vec{x}}_i+\alpha (\vec{b}-\mathbf{A}{\vec{x}}_i)(i=0,1,\dots ;\alpha =constant>0)$$
进行求解的方法称作Richardson 迭代法。

若 $\mathbf{A}$ 为对称正定矩阵 $(SPD)$，此时求解线性方程组 $(\ast )$ 与寻找如下多元二次函数 $g(\vec{x})$ 的最小值点的问题等价：
$$g(\vec{x})=\frac{1}{2}{\vec{x}}^T\mathbf{A}\vec{x}-{\vec{b}}^T\vec{x}$$

说明如下：

一方面，由 $\vec{u}$ 处取得二次函数最小值的必要条件：
$$▽g(\vec{u})=\mathbf{A}\vec{u}-\vec{b}=0⟹\mathbf{A}\vec{u}=\vec{b}$$
另一方面，若
$$\mathbf{A}\vec{u}=\vec{b}⟹▽g(\vec{u})=\mathbf{A}\vec{u}-\vec{b}=0$$
且
$$H(g)=\mathbf{A}(Hessian矩阵)$$
在 $\vec{x}=\vec{u}$ 处正定。故 $g(\vec{x})$ 在 $\vec{u}$ 处取得极小值，又由于系数矩阵正定，故使得梯度为零的点存在且唯一，故 $\vec{u}$ 为最小值点。

采用最速下降法求解该最小值问题：
$${\vec{x}}_{i+1}={\vec{x}}_i-\alpha ▽g({\vec{x}}_i)={\vec{x}}_i+\alpha (\vec{b}-\mathbf{A}{\vec{x}}_i)$$
其中，$\alpha >0$反映沿负梯度方向前进的长度。上式实际上就是Richardson 迭代法所采用的迭代格式，说明了Richardson 迭代法的合理性，解释了为何要求 $\alpha >0$。

2. Richardson 迭代法的收敛条件

Richardson 迭代格式为：
$${\vec{x}}_{i+1}={\vec{x}}_i+\alpha (\vec{b}-\mathbf{A}{\vec{x}}_i)=(\mathbf{E}-\alpha \mathbf{A}){\vec{x}}_i+\alpha \vec{b}(i=0,1,\dots )(1)$$
考虑该迭代格式的收敛性，等价于考虑迭代过程中误差的收敛性。设方程组有精确解 $\vec{x}$，则
$$\vec{x}=(\mathbf{E}-\alpha \mathbf{A})\vec{x}+\alpha \vec{b}(2)$$
由 $(1)-(2)$ 得：
$${\vec{e}}_{i+1}=(\mathbf{E}-\alpha \mathbf{A}){\vec{e}}_i(i=0,1,\dots )$$
则有：
$${\vec{e}}_i=(\mathbf{E}-\alpha \mathbf{A}){\vec{e}}_{i-1}=(\mathbf{E}-\alpha \mathbf{A}{)}^i{\vec{e}}_0$$
假设 $\mathbf{A}$ 为实对称矩阵，则 $\mathbf{B}\triangleq \mathbf{E}-\alpha \mathbf{A}$ 为实对称矩阵，其特征值 ${\lambda }_i^B(i=1,2,\dots ,n)$ 为实数，特征向量 ${\vec{u}}_i(i=1,2,\dots ,n)$ 两两正交，将误差在 $\mathbf{B}$ 特征向量构成的基上展开，有：
$$\begin{gathered} {\vec{e}}_i=\sum _{j=1}^n{c}_j^i{\vec{u}}_j=(\mathbf{E}-\alpha \mathbf{A}{)}^i{\vec{e}}_0=(\mathbf{E}-\alpha \mathbf{A}{)}^i\sum _{j=i}^n{c}_j^0{\vec{u}}_j=\sum _{j=1}^n{c}_j^0({\lambda }_j^B{)}^i{\vec{u}}_j \\ ⟹c_j^i=c_j^0({\lambda }_j^B{)}^i(j=1,2,\dots ,n) \end{gathered}$$
显然，想要随着迭代过程的进行，误差逐渐收敛，应有
$$-1<{\lambda }_j^B<1(i=1,2,\dots ,n)$$
即，
$$\rho (\mathbf{E}-\alpha \mathbf{A})<1⟹\{\begin{array}{l}1-\alpha {\lambda }_{min}^A<1⟹{\lambda }_{min}^A>0\\ \\ 1-\alpha {\lambda }_{max}^A>-1⟹0<\alpha <\frac{2}{{\lambda }_{max}^A}\end{array}$$
上式给出了Richardson 迭代的收敛条件：$\rho (\mathbf{E}-\alpha \mathbf{A})<1$，或者具体来说

(1) 系数矩阵 $\mathbf{A}$ 对称正定；

(2) $\alpha \in (0,\frac{2}{{\lambda }_{max}^A})$；

此外，上述讨论还说明 $\rho (\mathbf{E}-\alpha \mathbf{A})$ 越小，误差收敛越快，而 $\mathbf{B}$ 的收敛半径与参数 $\alpha$ 的选取有关，这意味着存在最佳的 $\alpha$ 使得Richardson 迭代的收敛速度最快，又：
ρ ( E − α A ) = m a x { ∣ 1 − α m i n A ∣ , ∣ 1 − α m a x A ∣ } \rho(\bold E-\alpha\bold A)=max\{|1-\alpha_{min}^A|,|1-\alpha_{max}^A|\} ρ(E−αA)=max{∣1−αminA​∣,∣1−αmaxA​∣}

由图可知：
$$min(\rho (\mathbf{E}-\alpha \mathbf{A}))=\frac{{\lambda }_{max}^A-{\lambda }_{min}^A}{{\lambda }_{max}^A+{\lambda }_{min}^A}$$
且此时，
$${\alpha }_{opt}=\frac{2}{{\lambda }_{max}^A+{\lambda }_{min}^A}$$

3. 算法实现与算例验证

function [x] = Richardson(A,b,x0,alpha,N)
% 函数功能：利用 Richardson 迭代法求解 n 元线性方程组 Ax=b
% 变量说明：
%         (1) A      = [n,n]  -   线性方程组的系数矩阵；              
%         (2) b      = [n,1]  -   线性方程组的常数项；                
%         (3) x      = [n,1]  -   线性方程组的解；                    
%         (4) alpha           -   Richardson 迭代法中的参数（>0）;
%                    = "数值" -   按设定的值进行计算（不满足方法条件时，警告）；
%                    = "best" -   按迭代速度最佳的alpha值进行计算；
%         (5) n               -   线性方程组未知量的个数；
%         (6) x0              -   解的初始猜测值；
%         (7) N               -   迭代次数；

% 检验输入参数是否符合 Richardson 迭代法求解的要求：
if isnumeric(N)
    if N <= 0
       error("错误：迭代步数[N]应为正整数")
    end
    if rem(N,1) ~= 0
       error("错误：迭代步数[N]应为正整数")
    end
else
    error("错误：迭代步数[N]应为正整数")
end
n = size(A,1);
if n ~= size(A,2)
   error("错误：系数矩阵[A]应为方阵")
end
if size(b,2) ~= 1
   error("错误：常数项[b]应为列向量")
end
if n ~= size(b,1)
   error("错误：系数矩阵[A]与常数项[b]维度不匹配")
end
if size(x0,2) ~= 1
   error("错误：初值[x0]应为列向量")
end
if n ~= size(x0,1)
   error("错误：初值[x0]维度不匹配")
end

% 检验线性方程组的系数矩阵是否为对称矩阵，若对称则必须正定:
flag = 0;
for i = 2:n
    for j = 1:i-1
        if A(i,j) ~= A(j,i)
           flag = 1;
        end
    end
end
if flag == 1
   for i = 1:n
       if det(A(1:i,1:i)) <= 0
          error("错误：系数矩阵[A]对称时应正定")
       end
   end
end

% 检验输入指定数值参数 alpha 时是否满足要求：
if isnumeric(alpha)
   if alpha > 2/vrho(A)
      warning("警告：参数[alpha]设置不满足Richardson 迭代法的收敛条件")
   end
elseif isstring(alpha)
   if strcmp(alpha,"best")
      alpha = 2/(max(eig(A))+min(eig(A)));
   else
      error("错误：参数[alpha]只能为数值或关键字[best]")
   end
else
   error("错误：参数[alpha]只能为数值或关键字[best]")
end

% Richardson 迭代:
i = 0;
while i < N
   x  = x0+alpha*(b-A*x0);
   i  = i+1;
   x0 = x;
end
end

利用编写完成的 Richardson 函数求解如下线性方程组 $\mathbf{A}\vec{x}=\vec{b}$ ：
$$\left[\begin{array}{cc}6& 3\\ 3& 4\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]=\left[\begin{array}{c}-3\\ -9\end{array}\right]$$
该线性方程组的精确解为：
$${\vec{x}}_{ex}=[1-3{]}^T$$
系数矩阵的特征值为：
$${\vec{\lambda }}^A={\left[\begin{array}{cc}1.8377& 8.1623\end{array}\right]}^T$$
由前面的讨论可知：

使得Richardson迭代法收敛的参数 $\alpha$ 的取值范围为：
$$0<\alpha <\frac{2}{{\lambda }_{max}^A}=0.245$$
且有：
$$\begin{gathered} {\alpha }_{opt}=\frac{2}{{\lambda }_{max}^A+{\lambda }_{min}^A}=0.2 \\ {\rho }_{opt}=\frac{{\lambda }_{max}^A-{\lambda }_{min}^A}{{\lambda }_{max}^A+{\lambda }_{min}^A}=0.6325 \end{gathered}$$
下面实际操作来进行下测试：

• 初值：${\vec{x}}_0=[0.00.0{]}^T$；
• $\alpha ={\left[\begin{array}{cccccc}0.06& 0.1& 0.2& 0.22& 0.24& 0.4\end{array}\right]}^T$;

% 程序功能：
%         Richardson 函数的测试算例

clear;clc;close all;
% 方程的系数矩阵与常数项：
A = [6,3;3,4];
b = [-3;-9];
% 该方程的精确解：
x_exc = [1;-3];
% 前50步的收敛曲线：
figure(1)
step = 50;
e = zeros(length(alpha),step);
for i = 1:length(alpha)
    for j = 0:step
        if j == 0
            e(i,j+1) = norm(x0-x_exc);
        else
            [x] = Richardson(A,b,x0,alpha(i),j);
            e(i,j+1) = norm(x-x_exc);
        end
    end
end

semilogy( 1:51,e(1,:), ...
         1:51,e(2,:), ...
         1:51,e(3,:), ...
         1:51,e(4,:), ...
         1:51,e(5,:), ...
         1:51,e(6,:),'LineWidth',1.8);
legend('\alpha = 0.06','\alpha = 0.10', ...
       '\alpha = 0.20','\alpha = 0.22', ...
       '\alpha = 0.24','\alpha = 0.40', ...
       'Location','northwest',          ...
       'FontName','Times New Roman',    ...
       'FontSize',12,                   ...
       'NumColumns',2)
legend('boxoff')
xlabel('Steps  \it{N}','FontSize',12,'FontName','Times New Roman')
ylabel('$||\vec{x}-\vec{x}_{ex}||_2$','FontSize',12, ...
       'FontName','Times New Roman','Interpreter','latex')
grid on

计算结果，如图所示：

结果与预期一致：

• 当 $\alpha =0.4>0.245$ 时，Richardson迭代法发散，其余算例均收敛；
• 当 $\alpha =0.2$ 时，Richardson迭代法收敛速度较其它算例更快;

观察收敛曲线可知，随着迭代的进行误差的对数与步数间逐渐线性化，这该如何解释？

实事上，收敛曲线的斜率满足：

从收敛曲线后半段斜率的几何意义上，我们也可以更好的理解为什么要求：
$$\rho (\mathbf{B})=\rho (\mathbf{E}-\alpha \mathbf{A})<1$$
且越小越好。为说明Richardson迭代法也适用于非对称系数矩阵的线性方程组求解，对上述算例进行雅可比左预处理，即
$$\left[\begin{array}{cc}\frac{1}{6}& 0\\ & \\ 0& \frac{1}{4}\end{array}\right]\left[\begin{array}{cc}6& 3\\ & \\ 3& 4\end{array}\right]\left[\begin{array}{c}{x}_1\\ \\ x_2\end{array}\right]=\left[\begin{array}{cc}\frac{1}{6}& 0\\ & \\ 0& \frac{1}{4}\end{array}\right]\left[\begin{array}{c}-3\\ \\ -9\end{array}\right]⟺\left[\begin{array}{cc}1& \frac{1}{2}\\ & \\ \frac{3}{4}& 1\end{array}\right]\left[\begin{array}{c}{x}_1\\ \\ x_2\end{array}\right]=\left[\begin{array}{c}-\frac{1}{2}\\ \\ -\frac{9}{4}\end{array}\right]$$
显然，这个与原方程有相同解的新方程组的系数矩阵非对称。对于新方程而言，系数矩阵的特征值为：
$${\vec{\lambda }}^{A^{\prime }}={\left[\begin{array}{cc}1.6124& 0.3876\end{array}\right]}^T$$
使得Richardson迭代法收敛的参数 $\alpha$ 的取值范围为：
$$0<\alpha <\frac{2}{{\lambda }_{max}^{A^{\prime }}}=1.24$$
且有：
$$\begin{gathered} {\alpha }_{opt}=\frac{2}{{\lambda }_{max}^{A^{\prime }}+{\lambda }_{min}^{A^{\prime }}}=1 \\ {\rho }_{opt}=\frac{{\lambda }_{max}^{A^{\prime }}-{\lambda }_{min}^{A^{\prime }}}{{\lambda }_{max}^{A^{\prime }}+{\lambda }_{min}^{A^{\prime }}}=0.6124 \end{gathered}$$

实际计算结果如图：