matlab lu分解求线性方程组_线性方程组(14)-直接解法（上）

其实关于迭代解法的预调节器还有很多。前面讲迭代解法讲得特别多，多到快忘记了第一节的直接解法了。我们来复习一下，前面提到LU分解的算法

function

矩阵

存储在返回矩阵的对角线部分和上三角部分；矩阵

的对角线为

，下三角部分存储在返回矩阵的下三角部分。

实际上这个顺序也并不固定，上面是从左上到右下的顺序，也可以是下面这样的从左至右顺序，等等

function

还可以从分块的角度来看

先分解左上角，再分解右下角

最终得到

---

主元选取

实际上这个版本的LU分解对于一般的矩阵来说数值不稳定，因为即使

是可逆的，算法中的主元（pivot）也可能是零或者接近零的小数，一个最简单的例子就是

。如果能有某种排列

、

使得

存在稳定的LU分解，则可以通过解方程

来得到解。在LU分解过程中，不断对

做行变换

在每一步变换前可以从

中选取一个绝对值较大的元素交换到左上角作为主元，即

最终

注意到交换两个元素的排列阵

而

，

，因此

仍然是下三角阵，则有

实际上需要数值稳定性只需要

或者

其中之一即可，即从

的第一列或第一行选绝对值最大的元素作为主元。同时使用

和

可以实现上一节提到的显示秩，即秩显示LU分解（RRLU）。

---

填入

如果

是一个稀疏矩阵，那么

和

有什么样的稀疏特征呢？在LU分解算法中，可以看到这样的更新操作

（

）

如果

原本为零，但

和

都非零，那么更新之后

也非零，这种现象称为填入（fill-in）。从图的角度来说，如果把

看成边

，则填入相当于在

被消去时，因为存在路径

而新增了边

。实际上只要存在路径

，且

均小于

和

则会出现

。因此，

和

的非零结构的计算需要从每个节点出发找到所有可达节点，相当于求解图遍历问题。

为表示方便令

。

非零结构的计算可以用和数值计算同样的方法完成，不过也有一些剪枝方法。这里介绍对称剪枝（symmetric pruning）：如果

非零（

），那么

将会更新到

上，即

的非零位置集合包含了

的非零位置集合。那么对于

，

对第

列的更新完全被

的更新所包含。也就是说，当发现

和

均非零时，从

到

的更新可以跳过不计算。注意这只是在计算非零结构时可以这样，在最终结果的数值计算中还是不能跳过。

前面用到了很多的

之间的大于号小于号这些条件，很容易发现，在一定的重排序下，这些大小关系很容易被改变。也就是说，

的LU分解的填入数和

完全不一样。例如

完全填满

则完全保持稀疏性

---

数据依赖关系图

当

（

）时，容易发现，实际上LU分解的第

列和第

列是互不影响的，可以并行计算的。因此对于稀疏矩阵的LU分解可以挖掘更多的并行性。

从列与列之间的依赖关系来看，这个依赖关系与

有关，即

表示要用第

列更新第

列，所以两列间有依赖，计算需要有一定的顺序来保证正确性。也可以从行与行的关系来看，

表示要用第

行更新第

行。依赖关系可以看成一个有向无环图，和第11节提到的三角求解相似。

L的稀疏结构及其行数据流图

由于依赖关系的传递性，这个图其实是有一定信息冗余的。当图中同时存在

、

和

三条边时，

其实是冗余信息。因为如果

和

的依赖都没有被破坏的话，

也自然不会被破坏。

因此对称剪枝方法又可以用到这里：对于

，如果

和

均非零（这个条件解释了“对称”），那么

或

的依赖关系就可以忽略。这种方法可以有效减少并行计算中调度和同步的开销。

但是要注意，这个依赖关系不是从原始矩阵

直接得到，而是从

或

得到的。这其中就不可忽视填入的作用。在某些应用场景下，填入元可能会使得并行性非常显著的下降。这也是除了内存消耗和计算量以外，填入元带来的另一个问题。对于一些由微分方程所导出的矩阵，如果没有合适的重排，那么填入元将直接破坏掉所有的稀疏并行性，使整个图变成一条长链。