方程组的直接解法和迭代法 python_线性方程组(16) - 直接解法（下）

最近对上一节又做了一些补充，最后发现过长，于是断开了一下，把减少填入的部分单提一节。

前面讲到，稀疏矩阵的重排序会影响填入数。那么，如何求得一个填入较少的排列呢？很遗憾，求解全局最优是NP完全问题，目前只有遍历几乎所有排列才能得到。不过有不少次优解算法可以使用。

先讨论对称正定问题。在对称问题下，填入元位置的计算从有向图问题变成了无向图问题。即如果存在路径

且

均小于

和

，则会出现填入边

。对称重排序问题又可以这么考虑：在无向图上求一个最佳的节点消去顺序，使得填入边最少，而每消去一个节点，则要产生填入边使得其所有相邻的未消去节点构成一个团。

最小度（minimum degree）算法。每次选取度最小的节点来消去，以此避免产生大团。该算法有近似最小度（approximate minimum degree, AMD），多重最小度（multiple minimum degree, MMD）等变种，均能在不明显影响填入数的前提下显著减少排列的计算时间。

嵌套分割（nested dissection）算法。将图划分成三个部分

、

和

，使得

和

之间没有直接的边相连。随后将矩阵重排成

再递归地对对角线上的子矩阵进行排序。该算法的关键在于图的划分，即在

和

节点数之差的绝对值受限的条件下，找到最小的割集

。而图的最佳划分仍然是一个NP完全问题，所以还需要一个图划分的近似优化算法。图划分作为一个算法问题已经超出了数值计算的范围，不单独成一节，在本节结尾介绍一下。

最小度法更关注图的局部特征，而嵌套分割法更关注整体，因此对于大规模矩阵的重排序，通常用嵌套分割法进行数层划分，到规模足够小时使用最小度法。

在大规模矩阵中，嵌套分割法除了能得到更少的填入元以外，用该算法得到的消去树也通常会有更好的并行性。这是因为得到图划分后，消去树自然地形成了

为根、

和

为子树的结构，而图划分的均衡性使得

和

规模相近，因此消去树更平衡。

最小度算法较容易推广到非对称的情况。但是嵌套分割算法则需要一些变化。一种方法是把有向图看成无向图，即用

来计算重排序；另一种方法是用

来进行计算重排序。

另外，在对称不定和非对称问题中，为了数值稳定性，也需要另一种重排序，也就是主元选取。主元选取和减少填入的兼容也是一个值得考虑的问题。这里介绍几种解决方法：

1. 在最开始可以对矩阵做预先重排，将绝对值大的元素排列在对角线上，同时主元选取有一定阈值，对绝对值达到阈值的对角元可以不用做行/列交换，这样发生行/列交换的概率就会小很多。
2. 在SuperLU中，对于非对称矩阵
   使用
   来计算减少填入的重排，也就是考虑乔列斯基分解
   得到排列
   ，然后做分解
   ，
   是选主元导致的行排列；可以证明对于任意的
   ，都满足
   的非零集合被包含于
   的非零集合，那么减少
   的填入也就相当于减少
   的填入。
3. 在MUMPS中，对结构对称的矩阵
   ，则优先在超节点内部进行选主元，使得这样的行/列交换不改变结果的稀疏性。若超节点内部无法满足主元要求，则将整个超节点中未完成分解的部分交换至消去树上的父超节点之后。由于这部分的非零结构是父超节点非零结构的子集，如果在父超节点选取主元，产生的新填入就会将会使得他们合并成一个新的超节点。如果父超节点上仍然不能满足要求，则逐层向上传递。

---

图划分

早在之前我们就提到过领域分解的概念，对于非结构化的稀疏问题，领域分解也同样能够有效利用并行计算资源。

当时我们提到的问题是边割集问题，找到节点划分

和

，在

和

平衡（即节点权重之差的绝对值不大于某阈值）的条件下，使得连接

和

两部分之间的边权值和最小。而在嵌套分割法中，要解决的问题是点割集问题，找到节点划分

、

和

，在

和

平衡，且没有边连接

和

两部分的条件下，使得

权重最小。两个问题在寻优方法上有相似性，只不过优化的目标函数有不同。这里要注意，无论是点割集问题还是边割集问题都是NP完全问题，实际用的算法都求的是某个局部最优解。

之所以要提到点和边的权重，是因为图划分中最广泛应用的是多层粗化（multilevel coarsening）算法，也就是对原始图

做一系列映射

最终在规模足够小的

上构造一个初始划分，采用一些局部搜索策略来得到局部最优，并逐层映射回原图。在每层后向映射之后，都可以采用一些局部寻优策略来调整当前的解。

都是有权重的。

的多个节点映射到

的一个节点上，所以在

上进行的小范围局部寻优就相当于在

上进行较大范围的局部寻优，从而使得在有限的计算能力下能搜索更大的空间。

在每层粗化中，将权重较大的边作为配对边，其所连接的两个节点映射成一个，权重相加；这条配对边则被消去；这两个节点上与同一个第三节点相连的边也将映射成一条边，权重相加。我们希望在尽可能使每个节点都有且仅有一个配对的条件下（不绝对，可能有无法配对的节点），配对的边权重之和尽可能大。这个优化也是NP完全问题，可以用贪心算法求近似最优，按边权从大到小考虑，尽可能为每个节点都找一个与之配对的节点，使映射到的下一层图规模有效地缩小。

图划分常用的局部寻优方法是FM（Fiduccia-Mattheyses）算法，其主要思想是已知一个解时，在不破坏平衡条件的前提下，将一个节点从划分中的一个子集移动到另一个子集，计算该操作对目标函数的增益，并选择增益最大的操作来执行。

在边割集问题中，当尝试将

中的一个节点

移动到

时，可根据其邻居

计算其增益

而在点割集问题中，尝试将

中的一个点

移动到

，为了维护

和

之间不直连的性质，还要将

中与

相邻的点移动到

，其增益

当所有可能的操作都是负增益时，算法结束。这样将节点一个一个移动，可能会导致陷入某个局部最优解而错过不远处的更优解，例如下面这个例子

如果用FM算法进入到图中的局部最优解，任何的节点移动都会导致负增益

这个问题的主要原因是图中有4个节点构成的簇（cluster），这4个节点内部连接紧密，不应该被拆分，要移动也应该4个节点整体移动。在多层粗化算法框架下，这4个节点有较大概率映射到粗图的同一个节点上，从而在粗图中实现整体移动，但是由于粗化映射算法也不是全局最优的，所以并不能避免这种情况出现。

簇检测迭代划分（cluster-detecting iterative improvement partitioner, CDIP）是一种可以一定程度上避免这个问题的算法，它可以检测到大小不超过

的这样的簇。在FM算法的基础上，当最大增益为负时，将当前的解

保存，然后尝试进行负增益的操作。我们假设存在这样一个簇

需要整体移动，那么在将负增益的节点

从

移动到

之后，把

视为

的一个成员，下一步移动则只考虑

中

的邻居，从中选出节点

进行移动，并把

也加入

。直到得到一个比

更优的解时，可以认为该簇已经被整体移动，回到平常的状态，不再考虑

这个簇。如果当

的节点数达到阈值

时仍未得到更优解，则退回到

（全部木大），避免搜索空间过大。