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

这一节会比较杂，列举一下话题：乔列斯基分解、消去树、超节点法与多波前法、填入减少的重排序、图划分。

乔列斯基分解（Cholesky decomposition）

上一节讲直接解法都没有把乔列斯基分解当作一个与LU分解不同的算法，从数学上可以看出LU分解的一种特例，但是在计算上还是有不少差别，所以这里单独提出来。

首先，乔列斯基分解是针对对称阵的算法。

对于对称正定矩阵

，可以求得

或

，其中

为下三角阵，

为对角阵。而且整个过程不需要主元选取也能保证数值稳定性。

对于对称不定矩阵

，可以求得

，其中

为下三角阵，

为由

和

的对称块组成的块对角阵。为了保证数值稳定性，需要一定的主元选取，变为

，例如

。注意这里为了保持对称性，行与列使用同样的排列。

当选取

的块为主元时

从块视角来看

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

最终得到

---

消去树（elimination tree）

上一节提到，对于稀疏矩阵的直接解法，可以根据其稀疏结构得到数据依赖图，而该数据依赖图可以用对称剪枝的方法去除冗余边。对于结构对称的矩阵来说，每个非零元都满足对称剪枝的条件，因此，其数据依赖图中，只有

每列第一个非对角非零元是必需的（同时也是

每行第一个非对角非零元），其它边都可去除掉。这样每个节点只被最多一个节点所依赖，该依赖关系构成树（或森林）。

结构对称矩阵，其中×为原矩阵非零元，□为填入元

消去树（图片来源：IBM WSMP文档，上同）

树比一般的有向无环图更容易调度。在同样的矩阵规模下，消去树越矮越宽，则计算过程中并行性越高。

虽然非对称情况的LU分解并不完全满足树结构，但是也是非常类似的，可以考虑

或

的乔列斯基分解，允许新增一些依赖关系的条件下，构造消去树。从消去树上就可以看出，在靠近根节点附近，并行计算的机会较少，直接解法的高性能并行实现在这里面临较大的挑战。

消去树还有一个重要的意义，那就是可以在得到

或

的非零结构之前就构造出来，在对称问题中，根据填入元产生规则：若

，

中存在路径

，

均小于

，那么分解后一定存在

，

在消去树上就是

的祖先。

取个更简单的形式：若

，

中存在

，那么

在消去树上就是

的祖先。

所以构造算法就是，对

的上三角部分的第

列（或下三角部分的第

行），找到必须是

的子孙的这些节点，将这些节点当前的祖先作为

的子节点。用MATLAB描述如下

function parent = etree_uptri(A)
% 用A的上三角部分计算消去树
% parent(i) 表示节点i的父节点
n = size(A, 1);
T = triu(A, 1); 
parent = 1 : n;                % 初始化每个节点是一棵独立的树
root = 1 : n;                  % 记录每个节点所在的树的根，用于加速查找
for k = 2 : n                  % 第一列没有非对角非零元，跳过
    row_idx = find(T(:, k));   % 第k列的非零元行下标
    for l = 1: length(row_idx)
        r = row_idx(l);        % 必须满足k是row_idx(l)的祖先

        % 这一步查找根的操作也可以在parent上进行
        % 这样就不需要root，但是查找比较慢
        while root(r) != r  
            r_new = root(r);
            root(r) = k;
            r = r_new;
        end

        % 到这里，r是row_idx(l)所在的树的根
        parent(r) = k;
        root(r) = k;
    end
end
end

---

超节点（supernode）

LU分解所产生的矩阵还有一个特点，那就是

的列和

的行有较大概率具有如下所示的非零结构，尤其是右下角比较明显。（SuperLU处理非对称问题，考虑的是第二种结构。MUMP处理对称问题，考虑的是第一种结构）

超节点非零结构，黑色表示稠密，画线代表具有相同非零结构的行/列（图片来源：SuperLU文档）

如之前的例子中例子中的5、6、7和8就满足这个条件。可以将这些节点合并为一个超节点。超节点的出现，主要是因为一行/一列被消去时，将会对相关的行/列进行更新，使得这些行/列的非零结构有较大概率和它相似。

如果按照一般的处理稀疏矩阵的方法来进行分解，那么计算性能仍然会较低下。主要原因是稀疏矩阵的数据表示结构比较复杂，为了压缩而引入了很多间接寻址，计算中也会出现较多分支。而另一方面，对于不需要压缩的稠密矩阵，可以用二维数组表示，利用循环展开、分块、预取等技巧，有BLAS和LAPACK等非常高效的实现。因此在稀疏矩阵分解中，我们还应该注意到这些特点，来帮助我们优化性能。一个超节点的非零元的值可以用一个二维数组表示，这样对超节点的运算便可以使用BLAS和LAPACK中的高效实现。

另外，考虑超节点会使得对矩阵进行符号分析（确定

和

的非零元位置）时处理的图规模下降。