数值分析：线性方程组迭代求解

前言

常用的线性方程组的直接法求解方法有LU分解、高斯消去、列主元消去法等，但是这些直接法最大的缺点就是对系数矩阵要求太高！导致直接解法的通用性受限。反观近似求解，迭代形式简单、程序方便写、创造收敛格式就可以迭代达到任意精度！看过本文后，相信以后你会选择用迭代法。

本文所有方法对应的matlab程序

迭代方法

常用的迭代方法有3种：雅克比迭代、高斯-赛德尔迭代、(超)松弛迭代。3者的关系是：雅克比迭代 → 赛德尔迭代 → 松弛迭代。所以3者的思路是一样的，这是迭代格式不一样而已。

注意：本文所讨论方程组都是n个方程n个未知数，即系数矩阵为正方形。

雅克比迭代

对于下面的n个方程n个未知数的正方形线性方程组：

$\begin{cases} a_{11}x_1 + a_{12}x_2 + a_{13}x_3 + \cdots + a_{1n}x_n = b_1 & \\ a_{21}x_1 + a_{22}x_2 + a_{23}x_3 + \cdots + a_{2n}x_n = b_2 & \\ \quad \quad \quad\quad\quad\quad\quad\quad\cdots \cdots & \\ a_{n1}x_1 + a_{n2}x_2 + a_{n3}x_3 + \cdots + a_{nn}x_n = b_n & \end{cases}$

如果系数矩阵A可逆，且对角元素 $\color{red}{a_{ii}≠0}$ ，根据克拉默法则该正方形方程组有唯一解！现将方程改写如下，即把第 $i$ 行的 $x_i$ 未知数单提到方程左边：

$\begin{cases} x_1 = \frac{1}{a_{11}}(\quad - a_{12}x_2 - a_{13}x_3 - \cdots - a_{1n}x_n + b_1) & \\ x_2 = \frac{1}{a_{22}}(a_{21}x_1- \quad - a_{23}x_3 - \cdots - a_{2n}x_n + b_2) & \\ \quad \quad \quad\quad\quad\quad\quad\quad\quad \cdots \cdots & \\ x_n = \frac{1}{a_{nn}}(a_{n1}x_1 - a_{n2}x_2 - \cdots - a_{n,n-1}x_{n-1} - \quad + b_n) & \end{cases}$

上式改成迭代格式很简单：

$\begin{cases} x_1^{(k+1)} = \frac{1}{a_{11}}(\quad - a_{12}x_2^{(k)} - a_{13}x_3^{(k)} - \cdots - a_{1n}x_n^{(k)} + b_1) & \\ x_2^{(k+1)} = \frac{1}{a_{22}}(a_{21}x_1^{(k)}- \quad - a_{23}x_3^{(k)} - \cdots - a_{2n}x_n^{(k)} + b_2) & \\ \quad \quad \quad\quad\quad\quad\quad\quad\quad \cdots \cdots & \\ x_n^{(k+1)} = \frac{1}{a_{nn}}(a_{n1}x_1^{(k)} - a_{n2}x_2^{(k)} - \cdots - a_{n,n-1}x_{n-1}^{(k)} - \quad + b_n) & \end{cases}$

或者把第 $i$ 行的 $x_i$ 在右边空的地方填上，前面再补一个可得另一迭代格式：√

$\begin{cases} x_1^{(k+1)} = x_1^{(k)} - \frac{1}{a_{11}}( a_{11}x_1^{(k)} + a_{12}x_2^{(k)} + a_{13}x_3^{(k)} + \cdots + a_{1n}x_n^{(k)} - b_1) & \\ x_2^{(k+1)} = x_2^{(k)} - \frac{1}{\color{blue}{a_{22}}}( a_{21}x_1^{(k)} + a_{22}x_2^{(k)} + a_{23}x_3^{(k)} + \cdots + a_{2n}x_n^{(k)} - b_2) & \\ \quad \quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \cdots \cdots & \\ x_n^{(k+1)} = x_n^{(k)} - \frac{1}{\color{blue}{a_{nn}}}( a_{n1}x_1^{(k)} + a_{n2}x_2^{(k)} + a_{n3}x_3^{(k)} + \cdots + a_{nn}x_n^{(k)} - b_n) & \end{cases}$

得到上面的迭代格式，其实任务已完成！但为了好编程，还是用矩阵表示为好，设：

$x^{(k)} = (x_1^{(k)}, x_2^{(k)}, x_3^{(k)}, \cdots, x_n^{(k)})^T$

$b = (b_1, b_2, b_3, \cdots, b_n)^T$

$D = \left( \begin{matrix} a_{11} & 0 & 0 & \cdots & 0 \\ 0 & a_{22} & 0 & \cdots & 0 \\ 0 & 0 & a_{33} & \cdots & 0 \\ \vdots & \vdots & \vdots & & \vdots \\ 0 & 0 & 0 & \cdots & a_{nn} \end{matrix} \right)$

随意给定一组初始值 $x^{0}$ ，可得到3种完全等价的矩阵迭代格式：

$\begin{cases} x^{(k+1)} = x^{(k)} - D^{-1}(Ax^{(k)} - b) & 用处: 推导后面2种方法 \\ x^{(k+1)} = \color{blue}{(E - D^{-1}A)}x^{(k)} + \color{blue}{D^{-1}b} & 用处: 用此式编程 \\ x^{(k+1)} = \color{blue}{B_1}x^{(k)} + \color{blue}{g_1} & 用处: B_1用来判断收敛性 \end{cases}$

注意：上面两式中的蓝色对应相等。

到此，雅克比迭代格式就得到了，可以看到非常容易编程实现！但是，大家注意到上面"标红"的条件：对角线元素如果等于0了，那么D对角矩阵就不可逆，那迭代不就连格式都写不出来？这里先埋下伏笔，后文会专门说明如何用"预处理"方法完美解决该问题！

高斯-赛德尔迭代

思想：在同一轮迭代中，前面已经更新的结果可以在本轮中给后面的用！而不是像雅克比迭代那样等一轮全结束了再一起更新！因此收敛更快。

迭代格式为：

$\begin{cases} x_1^{(k+1)} = x_1^{(k)} - \frac{1}{a_{11}}( a_{11}x_1^{(k)} + a_{12}x_2^{(k)} + a_{13}x_3^{(k)} + \cdots + a_{1n}x_n^{(k)} - b_1) & \\ x_2^{(k+1)} = x_2^{(k)} - \frac{1}{\color{blue}{a_{22}}}( a_{21}x_1^{\color{blue}{(k+1)}} + a_{22}x_2^{(k)} + a_{23}x_3^{(k)} + \cdots + a_{2n}x_n^{(k)} - b_2) & \\ \quad \quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \cdots \cdots & \\ x_n^{(k+1)} = x_n^{(k)} - \frac{1}{\color{blue}{a_{nn}}}( a_{n1}x_1^{\color{blue}{(k+1)}} + a_{n2}x_2^{\color{blue}{(k+1)}} + a_{n3}x_3^{\color{blue}{(k+1)}} + \cdots + a_{nn}x_n^{(k)} - b_n) & \end{cases}$

同样改为矩阵格式，这里稍微多一步把A系数矩阵单纯拆成上、对角、下3个部分：

$A = L + D + U$

$L = \left( \begin{matrix} 0 & 0 & 0 & \cdots & 0 \\ a_{21} & 0 & 0 & \cdots & 0 \\ a_{31} & a_{32} & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & & \vdots \\ a_{n1} & a_{n2} & a_{n3} & \cdots & 0 \end{matrix} \right) \quad U = \left( \begin{matrix} 0 & a_{12} & a_{13} & \cdots & a_{1n} \\ 0 & 0 & a_{23} & \cdots & a_{2n} \\ 0 & 0 & 0 & \cdots & a_{3n} \\ \vdots & \vdots & \vdots & & \vdots \\ 0 & 0 & 0 & \cdots & 0 \end{matrix} \right)$

注意1："单纯"分解就是单纯、简单的拆分成上、对角、下这3个部分！
注意2：matlab如何快速单纯取上下三角参考本文！

随意给定一组初始值 $x^{0}$ ，可得到2种完全等价的矩阵迭代格式：

$\begin{cases} x^{(k+1)} = \color{blue}{-(D + L)^{-1}U}x^{(k)} + \color{blue}{(D+L)^{-1}b} & 用处: 用此式编程 \\ x^{(k+1)} = \color{blue}{B_2}x^{(k)} + \color{blue}{g_2} & 用处: B_2用来判断收敛性 \end{cases}$

高斯-赛德尔迭代到此结束，同样非常好编程实现！但是同样有一个潜在问题：如果对角矩阵D中有0，那么(D+L)这个下三角阵的对角元素有0，是不可逆的！也就是说和前面雅克比迭代一样，系数矩阵A若对角元素有0，连迭代格式都写不出来！

(超)松弛迭代：

前面高斯-赛德尔法中提到了将系数矩阵A进行单纯/简单分解： $A = L + D + U$ ，将上式分别带入到雅克比、高斯-赛德尔中可得二者新的迭代公式(实际中不使用)：

新雅克比迭代：

$x^{(k+1)} = x^{(k)} - D^{-1}(Lx^{(k)} + Dx^{(k)} + Ux^{(k)} - b)$

新高斯-赛德尔迭代：

$x^{(k+1)} = x^{(k)} - D^{-1}(Lx^{\color{blue}{(k+1)}} + Dx^{(k)} + Ux^{(k)} - b)$

说明：搞成上面这两种迭代格式，纯粹是为了了解他们真实迭代过程的内涵，不是为了实际使用(实际使用还是用前文说的那些)的哈！在新高斯-塞德斯基础上引入一个常数ω，得到：

(超)松弛迭代：

$x^{(k+1)} = x^{(k)} - \color{blue}{\omega} D^{-1}(Lx^{\color{blue}{(k+1)}} + Dx^{(k)} + Ux^{(k)} - b)$

同样，将上式稍作改写，即把所有的(k+1)次幂项都移到移到等号左边，可得适宜编程和收敛判断的最终的2种形式：

$\begin{cases} x^{(k+1)} = \color{blue}{(D+\omega L)^{-1} [(1-\omega)D - \omega U]} x^{(k)} + \color{blue}{\omega (D+\omega L)^{-1}}b & 用处: 用此式编程 \\ x^{(k+1)} = \color{blue}{B_3}x^{(k)} + \color{blue}{g_3} & 用处: B_3用来判断收敛性 \end{cases}$

同样，(超)松弛迭代也非常容易编程实现！只要参数ω数值给的得当，收敛速度会变高斯-赛德尔快！但至于取多少，很难定夺。

因为(超)松弛迭代就是基于高斯-赛德尔迭代来的，所以它同样存在一个问题：如果系数矩阵A对角元素有0，那么同样连迭代格式也写不出来！下面我们就用一种简单的"预处理"方法来完美解决对角元素有0的情况！目标：通过预处理后，新的等价方程组一定可以写出迭代表达式！

迭代法的收敛性判断

判断所依据的迭代格式：

$x^{(k+1)} = \color{red}{B}x^{(k)} + g$

收敛充要条件： $\rho(B) < 1$ ，或： $\parallel B \parallel < 1$ 。即B的谱半径或任意一种范数要小于1，数值越小收敛越快！

2条重要的特殊推论：

若系数矩阵A为对称正定阵，则高斯-赛德尔迭代对任意初始值都！收！敛！
若系数矩阵A为对称正定阵，且松弛系数 $0 < \omega < 2$ ，则(超)松弛迭代对任意初始值都！收！敛！

线性方程组预处理

预处理解决什么问题：如果n个方程n个未知数的线性方程组，它的系数矩阵A的对角元素有0存在，那么任何一种迭代方法都写不出来迭代表达式，即卡在第一步！因此，预处理就是为了解决这一问题。

开始之前，不给证明的列出3条关于"对角元素有0的方阵"的必备知识：

非奇异方阵A经过一些行之间的调换，一定可以得到一个对角元素非0的新的形式！
非奇异方阵A， $A^TA$ 一定是一个对称阵，且为"对称正定阵"；
(对称)正定阵的主对角线元素一定都是大于0的！

( 注：上面3条信息的更多内容可参看这篇文章。)

有了上面的3条必备知识，下面介绍2条思路来做预处理：

对原对角有0的系数矩阵A做行交换(b也得跟着换)，直到换出一个对角元素无0的新形式！用新形式(新顺序的线性方程组，肯定是等价的)来迭代；缺点：虽思路简单直接，但是编程可能要用到"二分图"来遍历！不好实现。
原方程两边同时左乘原系数矩阵A的转置，形成了等价新形式：

$\color{blue}{A^TA}x = \color{blue}{A^Tb}$

这个方程左边蓝色看成新的系数矩阵，根据上面的第2、3条知识，这个矩阵一定是对称正定的！它的对角线元素一定是都是大于0的！又因为两边左乘的是一个非奇异矩阵，所以新方程与原方程是完全等价的。优点：一个矩阵的转置很好求。因此我推荐使用这种方法！！

用第2条思路做预处理，就这么简单就完事了！并且这种处理是万能适用的(根据知识2和3)！即：只要原方程确定有解(A非奇异)，不管原系数矩阵A对角有多少个0(哪怕全是0)都能给你处理妥当！！迭代就根据新的方程组进行就行：

$\begin{cases} \color{blue}{B_{4}}x = \color{blue}{g_4} & \\ \color{blue}{B_4} = \color{blue}{A^TA}\quad \color{blue}{g_4} = \color{blue}{A^Tb}& \\ \end{cases}$

更加牛逼的是，根据"收敛性判断"中的2条重要推论可得：预处理后的新方程形式，用高斯-赛德尔迭代一定收敛！！万能！！

第1次更新

本次更新补充预处理中第1种思路的一种解决方法！利用：系数矩阵对角最大化。注意做列对换的时候，x中对应的行要跟着一起变，这样方程才会等价；最后得到的解的顺序还要再做一下调整，恢复到原方程的解的顺序。程序参考本文开头地址；"矩阵行列对换如何保持方程等价"参考本文"第1次"更新。

一个简单实例流程图即可搞清：

图1：对角最大化操作示意图

这种方法虽然既不万能(极端特殊情况导致变换完对角还有0！)，而且还有些麻烦(对换A时x和b要跟着一起！)。但是它有一个优点：原先高斯-赛德尔迭代不收敛的系数阵，这样换完有很大几率变成收敛形式！这里给出程序中的一个例子：

% 原始系数矩阵: 
A =

     1     2     3     1
     1     4     6     2
     2     9     8     3
     3     7     7     2

这个原始的系数矩阵，未预处理时它的高斯-赛德尔迭代格式中的B的谱半径是大于1的！是不收敛的！但是经过对换预处理后变成：

% 对换预处理后:
A =

     9     8     3     2
     7     7     2     3
     4     6     2     1
     2     3     1     1

对它再做高斯-赛德尔迭代，谱半径小于1！可以收敛。

具体内容参考程序中的详细注释。

说明一个问题：为什么系数矩阵做不了对角占优的操作？

因为：对角占优考虑的每一行中的绝对值最大元素在对角位置。在线性方程中的对系数矩阵的等价变换只能是按行或按列来的，不能是元素与元素之间的交换这样的(因为某2个元素改变了，后面x也跟着改变顺序：这样最多就满足了一个方程，剩下的所有方程都不匹配了！所以等价变换只能行与行，列与列整体的换)。举一个很简单的例子，如下面的系数矩阵：

A =

    11    -2     3
     8    -5     2
     1     3     7

第1行不用动，11正好最大而在对角位置。第2行最大的8在第1列，想要换到第2列就要第1列和第2列对换，此时就会把第1列中已经排好的给再次打乱！

A =

    -2    11     3
     -5    8     2
     3     1     7

综上：不是不想把系数矩阵搞成对角占优阵！而是在等价变换的限制下不可能实现！只能被迫选择退而求次的本文的"对角最大化处理"！正因为是替代品，所以无法做到万能。

最后编辑于：2019.04.23 21:29:00

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 216,544评论 6赞 501
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 92,430评论 3赞 392
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 162,764评论 0赞 353
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 58,193评论 1赞 292
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 67,216评论 6赞 388
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 51,182评论 1赞 299
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 40,063评论 3赞 418
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 38,917评论 0赞 274
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 45,329评论 1赞 310
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 37,543评论 2赞 332
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 39,722评论 1赞 348
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 35,425评论 5赞 343
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 41,019评论 3赞 326
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 31,671评论 0赞 22
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 32,825评论 1赞 269
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 47,729评论 2赞 368
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 44,614评论 2赞 353