方程组的直接解法和迭代法 python_线性方程组(12)-多重网格

特征值与收敛性

很多迭代解法的收敛性与矩阵的特征值分布有关。比如静态迭代法

假设精确解为

，误差

，那么

如果初始误差在

的较

大特征值所对应的特征子空间上投影较大，比如接近

，那么迭代收敛较慢。

注意这里比较特征值的大小都比较的是绝对值。

再考虑与

的关系。在迭代中，我们往往是通过残差

来修正当前的解并得到下一步的解。不严格地说，如果当前的误差在

的较

小特征值所对应的特征子空间上投影较大，那么得到的残差也较小，从而不能有效地修正解的误差。

---

几何多重网格（Geometric multigrid, GMG）

对空间的偶数阶偏微分算符的特征函数是正弦函数，且高频正弦函数的特征值更大

带偶数阶空间偏微分的方程如波动方程、热传导方程、拉普拉斯方程等的通解在空间上都有类似正弦函数的分布规律。将这类方程离散化，所导出的线性方程组的系数矩阵的特征向量也往往具有类似正弦函数的分布规律，而且高频特征向量的特征值更大。如一维情况的二阶导（

）和四阶导所导出的系数矩阵分别为

以一维情况的二阶导为例，以

将微分方程离散化，用对称高斯-赛德尔迭代法（在这里保存对称性是为了方便观察特征值），观察

的特征值和特征向量。

A（上）和G（下）的低频特征向量

选取

的规模为

，可以发现

的的低频特征向量所对应的特征值较小，而

的低频特征向量所对应的特征值较大，最低频率的特征值是

。这就使得迭代对于误差有低通滤波特性。如果初值误差的低频分量较大，那么收敛将非常缓慢。再选取更大的

，用对称高斯-赛德尔迭代法求解，并对误差进行傅里叶变换来观察频谱。

低频误差不能有效缩小

不过，如果能选取一个低频分量较小的初值，那么迭代步数将减少。我们注意到低频和高频是相对的。如果一个函数原离散域的频率为

，则离散化后的频率为

。 将它降采样后，只要新的采样频率满足香农采样定理，其频率会升高。

降采样后周期变短，频率升高

如果能用

将微分方程重新离散化，并迭代求解新的方程，那么就能得到一个低频误差被抑制的解。这就是多重网格算法。

将微分方程重新离散化这一过程称为网格粗化（coarsening）。

原始网格下已有迭代结果

，要求方程

相当于在粗化的网格下求解

其中

由

降采样得到，就是相当于左乘矩阵

这里是取三点加权平均的方式，若要直接丢弃一半的数据，可以令

解出低频误差被抑制的

后，插值得到

，将

作为新的解。

插值相当于左乘矩阵

将残差抽样到粗网格，求解后将解插值到细网格，即

，

。

对比可得

。

当然，即使所得的

是

的精确解，

也很可能不是

的精确解。将

作为新的解之后还要在细网格内再进行一定步数的迭代。

网格还可以多层粗化，成为一个递归过程，即

多重网格算法求解

：

1. 如果网格规模足够小，精确或不精确地求解方程，返回解；
2. 否则，用迭代法以
   为初值迭代数步，得到
   ，并计算残差
   ；
3. 将残差粗化
   ；
4. 令
   ，使用多重网格算法求解
   ；
5. 将解插值并累加
   ；
6. 以
   为初值迭代数步，得到
   ，返回
   。

其中，以

为初值的迭代称为前平滑（pre-smoothing），以

为初值的迭代称为后平滑（post-smoothing）。

还是这个一维的例子，选取

，选取

看八倍粗网格的低频特征值，发现从

减小到了

八倍粗化的A（上）和G（下）的低频特征向量

选取

，用三层粗化的多重网格，前平滑和后平滑均为两步高斯-赛德尔迭代，而最粗网格采用四步高斯-赛德尔迭代。通过寥寥数步迭代，就达到了简单高斯-赛德尔迭代法324步才能达到的精度。

粗网格上迭代后显著抑制了低频误差

在一维的例子中，二倍粗化的网格计算开销只有原来的二分之一，而四倍、八倍的粗网格计算开销更小。所以每步迭代的计算量是最细网格层计算量的最多两倍。不过在迭代前需要预先计算矩阵乘

，这个矩阵只需要计算一次，之后每步迭代都可以复用结果。这些额外的矩阵也会占用内存，不过由于规模较小，总的内存占用也是最细网格层的最多两倍。

一维场景只是用来做例子，其实一维情况的多对角方程是很容易直接用LU分解法求解的。对于高维问题，如果问题本身是各向同性的，那么可以多个方向同时粗化，比如在三维空间下可以在每层粗化中每个维度缩小二分之一，总网格点数缩小到八分之一；如果问题本身有各向异性特点，根据每个方向的数值特点，可以选择其中一些方向进行更大幅度粗化，另一些方向进行小幅度粗化或者不粗化。

这里的

和空间有一定的对应关系，所以这种算法称为

几何多重网格。

---

代数多重网格（Algebraic multigrid，AMG）

仍然以微分方程为背景，对于不规则区域，或者区域内各部分的分辨率要求不一致的情况下，可能需要用非结构化网格。下图所示的例子中，每个三角形顶点代表一个变量，每条边代表两个变量之间的关系，也就是系数矩阵的一个非零元。当然，在这种情况下，因为顶点变量只能计算一阶导（即两点值之差除以位移），处理高阶导需要引入额外的边变量，所以通常不会再用上一节提到的简单的有限差分法，而用有限元法或者有限体积法求解，通过对方程积分，来尽量避免高阶导项的直接出现。

非结构化网格（图片来源：Hypre User Manual）

如图的非结构化网格与结构化网格的不同之处在于，每个网格点与哪些网格点相邻是不固定的。从矩阵上来说，矩阵每行的非零元位置是不固定的。

这种情况下我们仍然希望利用粗化网格来抑制低频误差。但是问题在于，如何选取粗网格，也就是如何构造抽样矩阵和插值矩阵呢？这种情况下，如果知道每个点的坐标（比如

），虽然网格点之间的关系不再那么规则（比如

的邻居不再是

），仍然可以利用几何信息构造粗网格，不过这个构造过程稍复杂一些，计算开销会比结构化情况略大一些。也可以不用几何信息，网格点之间的关系只从系数矩阵确定，这种方法就是

代数多重网格。

我们注意到，虽然不再是有限差分法，但是微分项同样是相邻点的变量之差分来代替的，所以和结构化情况一样的是，离得近的两个点

之间对应了一个较大的负数

或

。而偶数阶微分所形成的系数矩阵是对称的，考虑到可能加上其它小项，系数矩阵应该接近对称。

粗化

如果

就可以认为变量

强依赖于变量

，其中

。若

取0，则任意负数都将视为强依赖。

该强依赖关系可以用图表示，对于接近对称的系数矩阵，该图很可能是无向图。于是问题变成：从图

中的节点中选取

，使得某个条件

满足。

是一个启发式条件，有多种选取方法

比如较强的条件

：

，即对任意两个有强依赖关系的变量，要么其中一个变量被选取，要么两个变量共同强依赖的变量被选取。这样一来在往回插值的时候，

的值能把

所强依赖的所有变量都考虑进去。

或者稍弱的条件

：

，即任意一个未被选取的变量强依赖于至少一个被选取的变量。而在插值时部分强依赖变量的影响就考虑不到。

要注意，满足条件的

不止一个，更小的

更能提高粗网格的计算速度。而全局最优是很难计算的。

在很多粗网格点选取算法中都用了最大度优先策略，计算每个变量被多少个变量所强依赖，优先选取值较大的，选取之后更新图，当然并行算法通常不会每次只选一个。满足

的算法有RS、CLJP、Falgout等；满足

的算法有PMIS、HMIS等。

插值

插值公式可以选择简单的插值

其中

，即只考虑粗网格变量的直接影响

或者对满足

的粗网格，使用更精确的插值

其中

（未被选为粗网格的强依赖变量，必然有共同强依赖的粗网格变量），

（弱依赖变量），即考虑所有强依赖变量的影响，包括粗网格变量的直接影响以及两步间接影响。注意该公式对不满足

的粗网格可能会出现除零。

---

由偶数阶偏微分所导出的方程，往往有对低频输入不敏感的特点，其中以椭圆型方程（泊松方程）为典型，在均匀介质中；抛物型方程（热传导方程）和双曲型方程（波动方程）总体上也有这个特点，不过由于加上了时间项，频谱与椭圆形方程有略微的差异，但仍然对误差有一定程度的低通特性。

上面所介绍的迭代模式是从最细网格开始不断粗化、前平滑、粗化……到最粗网格，然后不断插值、后平滑、插值……回到最细网格。这样的一个循环称为V-循环（V-Cycle）。多重网格的迭代模式还有W-循环，完全多重网格（full multigrid, FMG, 也叫F-循环）等。W-循环和F-循环的特点是更多地停留在粗网格层，使得低频误差被消除得更彻底，同时粗网格的计算开销也更小。

V-循环、W-循环和F-循环（图片来源：Gilbert Strang, Computational Science and Engineering）

多重网格可以作为单独的求解方法，反复进行V-循环（或W-循环、F-循环）来求解方程，也可以将一个V-循环作为克雷洛夫子空间法的预调节器。前面的小实验可以看到，即使是用了多重网格，如果最粗网格上还是高斯-赛德尔迭代，误差零频率分量也没能有效消除。所以，作为单独的求解方法时，最粗网格的求解最好是精确求解；在网格足够小时，可以用LU分解法直接求解。作为预调节器时，最粗网格的求解可以是不精确的。

在求解椭圆型方程时，多重网格是目前常用的收敛性最好的预调节器之一，不过对于椭圆形方程以外的情况无法保证。