迭代法（逐次逼近法）

迭代法（逐次逼近法）
在线性代数中我们常看到方程组被写为这样的形式：
$$Ax=b$$
其中A是非奇异矩阵（行列式不等于0）。本科阶段，我们求解的方程组阶数都不高，一般使用主元消去法求解。但对于A的阶数很大，而且零元素很多的大型稀疏矩阵方程组，例如，训练一个包含几十MB乃至几百MB的数据集时，主元消去法就显得力不从心了，而一般要选用逐次逼近法（或称为迭代法）求解。
为了便于说明，下面我们举一个求解线性方程组的迭代法例子。
$$\{\begin{array}{l}8x_1-3x_2+2x_3=20\\ 4x_1+11x_2-x_3=33\\ 6x_1+3x_2+12x_3=36\end{array}$$
如果记为Ax=b，其中：
$$A=\left[\begin{array}{ccc}8& -3& 2\\ 4& 11& -1\\ 6& 3& 12\end{array}\right]x=\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]b=\left[\begin{array}{c}20\\ 33\\ 36\end{array}\right]$$
方程组的精确解是：
$$x^{\ast }=(3,2,1{)}^{\mathrm{T}}$$
如果记为另一种形式：
$$\{\begin{array}{c}x=B_{0}x+f\\ x_1=\frac{1}{8}\left(3x_2-2x_3+20\right)\\ x_2=\frac{1}{11}\left(-4x_1+x_3+33\right)\\ x_3=\frac{1}{12}\left(-6x_1-3x_2+36\right)\end{array}$$
转换为矩阵的形式：
$$B_0=\left[\begin{array}{ccc}0& \frac{3}{8}& \frac{-2}{8}\\ \frac{-4}{11}& 0& \frac{1}{11}\\ \frac{-6}{12}& \frac{-3}{12}& 0\end{array}\right]f=\left[\begin{array}{c}\frac{20}{8}\\ \frac{33}{11}\\ \frac{36}{12}\end{array}\right]$$
任取初始值，例如取$x^{(0)}=(0,0,0{)}^{\mathrm{T}}$。将这些值代入公式(5)右边，即求得方程组的第一次迭代方程组的解，得到新的值。
$$x^{(1)}={\left(x_1^{(1)},x_2^{(1)},x_3^{(1)}\right)}^{\mathrm{T}}=(2.5,3,3{)}^{\mathrm{T}}$$
再将$x^{(1)}$的分量代入公式(5)右边得到$x^{(2)}$。反复利用这个计算程序，得到一个向量序列和一般的计算公式（迭代公式）简写为：
$${\mathbf{x}}^{(k+1)}={\mathbf{B}}_0{\mathbf{x}}^{(k)}+\mathbf{f}$$
其中k为迭代次数（k=0，1，2,…）。
迭代10次之后得到：
$$x^{(10)}=(3.000032,1.999838,0.9998813{)}^{\mathrm{T}}$$
误差向量范数：
$${\Vert {\epsilon }^{(10)}\Vert }_{\infty }=0.000187\left({\epsilon }^{(10)}=x^{(10)}-x^{\ast }\right)$$
代码实现：

import numpy as np
from numpy import linalg
import matplotlib.pyplot as plt

A = np.mat([[8,-3,2],[4,11,-1],[6,3,12]])
b = np.mat([20,33,36])
result = linalg.solve(A,b.T)
print(result)

# 迭代法
B0 = np.mat([[0,3/8,-2/8],[-4/11,0,1/11],[-6/12,-3/12,0]])
f = np.mat([20/8,3,3])
# x = B0x+f.T
error = 1.0e-10 # 误差阈值
steps = 100 # 迭代次数
xk = np.zeros((3,1)) #初始化值
errorlist = []
for k in range(steps):
    xk_1 = xk # 上次的xk
    xk = B0*xk+f.T # 本次的xk
    errorlist.append(linalg.norm(xk-xk_1)) # 计算存储误差
    if errorlist[-1] < error: # 判断误差是否小于阈值
        print(k+1) # 输出迭代次数
        break
print(xk) #输出结果

# 误差收敛散点图
plt.plot(range(1,26),errorlist,'o')
plt.show()

输出：

[[3.]
 [2.]
 [1.]]
25
[[3.]
 [2.]
 [1.]]

摘自：《机器学习原理与算法实践》