这一章节要解的问题和上一章是一样的,依然还是n 元线性方程组的求解问题。
但是,上一章主要是通过矩阵的线性变换转换成可以快速求解的三角阵或者对角阵的方式进行求解,其计算结果是精确的结果。
而本章则是的思路则是将问题Ax=y 转换成x = A'x 的迭代形式,从而,我们就可以给出迭代数组x_{k+1} = A'x_{k}。
此时,如果x_k 满足收敛条件,那么x_{k} 就会收敛到x = A'x 的一组解当中,上述问题同样可以得到解答。
对原始方程进行线性变换,我们显然有:
我们将其写成矩阵形式即为:
我们将其简化为符号表达即为:
基于此,我们就可以给出Jacobi迭代公式如下:
我们给出更加严格的数据表达如下,令:
则有:
进而可以导出:
两边左乘D^{-1} 即有:
记B=D^{-1}(D-A) ,g=D^{-1}y ,即可得到上一小节中一模一样的内容。
迭代矩阵B 即为Jacobi迭代矩阵。
Jacobi迭代收敛的充要条件为:
迭代矩阵B 的谱半径\rho(B) < 1
亦即:
具体而言,有如下定理:
定理 6.1 若方程组Ax=y 的系数矩阵A,满足下列条件之一,则其Jacobi迭代收敛: (1) A为行对角优阵,即|a_{ii}| > \sum_{j \neq i} |a_{ij}| ,i=1,2,...,n ; (2) A为列对角优阵,即
,j=1,2,...,n ;
最后,我们给出Jacobi迭代的python伪代码如下:
def jacobi_iter(A, y, epsilon=1e-6):
n = len(A)
B = [[0 for _ in range(n)] for _ in range(n)]
g = [0 for _ in range(n)]
for i in range(n):
for j in range(n):
if j == i:
continue
B[i][j] = -A[i][j] / A[i][i]
g[i] = y[i] / A[i][i]
x = [0 for _ in range(n)]
for _ in range(10**6):
z = [0 for _ in range(n)]
for i in range(n):
z[i] += g[i]
for j in range(n):
z[i] += B[i][j] * x[j]
if max(abs(z[i] - x[i]) for i in range(n)) < epsilon:
break
x = z
return z
Gauss-Seidel迭代方程和上述Jacobi迭代事实上是非常相似的,唯一的区别在于说Jacobi迭代是以x^{(k)} 为整体每次一起进行迭代更新的,而Guass-Seidel迭代则是在计算每一个x^{(k)}_{i} 的时候就是用当前已经迭代计算完成的所有的x^{(k)}_{j} 的结果。
即是说,以每一个元素为单位不断进行迭代更新。
用公式表达即为:
同样的,我们给出更加严格的数学推导如下:
从而,我们有:
亦即:
两边同乘以(D+l)^{-1} 即有:
令S = -(D+L)^{-1}U ,f = (D+L)^{-1}y ,我们即可写出Gauss-Seidel迭代公式:
称迭代矩阵S 即为Gauss-Seidel迭代矩阵。
当然,如上一小节所述,实际在算法实现当中并没有必要计算出S 和f ,只需要根据前面Jacobi 迭代矩阵进行实时参数更新即可。
同样的,我们给出书中关于Gauss-Seidel迭代的收敛条件如下:
定理6.2 若方程组系数矩阵为行或列对角优时,则Gauss-Seidel迭代收敛。 定理6.3 若方程组系数矩阵A 为对称正定阵,则Gauss-Seidel迭代收敛。
同样的,我们用python给出伪代码如下:
def gauss_seidel_iter(A, y, epsilon=1e-6):
n = len(A)
B = [[0 for _ in range(n)] for _ in range(n)]
g = [0 for _ in range(n)]
for i in range(n):
for j in range(n):
if j == i:
continue
B[i][j] = -A[i][j] / A[i][i]
g[i] = y[i] / A[i][i]
x = [0 for _ in range(n)]
cnt = 0
for _ in range(10**6):
for i in range(n):
z = g[i]
for j in range(n):
z += B[i][j] * x[j]
if abs(z-x[i]) < epsilon:
cnt += 1
else:
cnt = 0
x[i] = z
if cnt >= n:
break
if cnt >= n:
break
return x
松弛迭代的原型依然还是之前的Jacobi迭代,不过,和Gauss-Seidel迭代的实时参数更新不同,松弛迭代在这里是对Jacobi迭代式的批次更新以及Gauss-Seidel迭代式的实时更新取了一个折中,通过一个超参w 来进行参数更新速度的控制。
具体而言,即为:
同样的,我们给出松弛迭代的数学表达如下:
仿照Gauss-Seidel迭代,我们同样可以给出其迭代矩阵格式如下:
其中,
而松弛迭代的收敛判断存在如下两个定理:
定理 6.4 松弛迭代收敛的必要条件为0 < w < 2 ; 定理 6.5 若A 为对称正定矩阵,则当0 < w < 2 时,松弛迭代恒收敛; 特别的,当0 < w < 1 时,称其为亚松弛迭代;当1 < w < 2 时,称其为超松弛迭代;当w=1 时,迭代退化为Gauss-Seidel迭代。
最后,我们同样给出松弛迭代的伪代码实现如下:
def loose_iter(A, y, w=0.5, epsilon=1e-6):
n = len(A)
B = [[0 for _ in range(n)] for _ in range(n)]
g = [0 for _ in range(n)]
for i in range(n):
for j in range(n):
if j == i:
continue
B[i][j] = -A[i][j] / A[i][i]
g[i] = y[i] / A[i][i]
x = [0 for _ in range(n)]
cnt = 0
for _ in range(10**6):
for i in range(n):
z = (1-w) * x[i] + w * g[i]
for j in range(n):
z += w * B[i][j] * x[j]
if abs(z-x[i]) < epsilon:
cnt += 1
else:
cnt = 0
x[i] = z
if cnt >= n:
break
if cnt >= n:
break
return x
最后,我们来考察一下逆矩阵计算。
逆矩阵的计算原则上来说其实算是上述解线性方程组的一个特殊应用,事实上解n 个单元向量然后将其解拼接一下就能得到我们的逆矩阵了。
我们这里就不进行详述了,只给出python伪代码实现如下:
def cal_inverse_matrix(A):
n = len(A)
res = [[0 for _ in range(n)] for _ in range(n)]
for j in range(n):
y = [0 for _ in range(n)]
y[j] = 1
x = gauss_seidel_iter(A, y)
for i in range(n):
res[i][j] = x[i]
return res