Linear Algebra Done Left - A Taste of Direct Methods for Linear Systems

最近更新速度比我想的要快，这就出第三篇了。。。也就是稀疏矩阵的直接求解方法。

上一篇讲过了矩阵的一些迭代求解方法，这些方法可以很自然地应用于稀疏矩阵中。然而直接方法就并非如此了。通常，直接方法将矩阵分解为三角阵的乘积（这一步必然是立方复杂度），分别求解。然而，稀疏矩阵的规模巨大，直接套用稠密矩阵的做法并没有充分利用稀疏性，效率必然很低，还有可能导致分解之后退化为稠密矩阵。因此稀疏矩阵求解的直接方法往往比迭代方法复杂很多。

整个话题相当庞杂，因此这篇文章也是长期更新。但是显然我没有那么多工夫去天天折腾数值算法所以说不准哪天本文就永远停更了

国内似乎缺少这方面的专著，专门讲这个的文章也很少，多半没有代码，讲得也较浅。经过一轮搜索和试读后，这里给出$3$本参考书。$\text{Jennifer Scott}$的$\text{Algorithms for Sparse Linear Systems}$这本书，这本书的理论讲解严谨细致而且清晰。另一本经典著作则是$\text{Duff}$的$\text{Direct Methods for Sparse Matrices}$，因为看了一会儿就云里雾里了在此也作为参考。最后一本是$\text{Tim Davis}$的$\text{Direct Methods for Sparse Linear Systems}$，比较入门，主要侧重于实现，有配套的$\text{C}$代码可以学习。

$\text{1. 直接分解算法的基本模式}$

一般来说，直接分解求解稀疏矩阵有$3$个$\text{Pass}$。

第一个$\text{Pass}$被称作$\text{Analyse/Symbolic Factorization}$，顾名思义，这一步对稀疏矩阵的非零元的分布进行分析和追踪，来推算分解之后非零元的位置，以便后续填充。需要注意的是由于这一步不进行实际的数值运算，有可能在后续填充中计算出$0$也被存储进来，一般情况下这问题并不大。这一阶段涉及大量的图论，从篇幅上来看相当复杂，但是实际上每一步都是非常简单显然的，大量简单的小结论一步步堆起来得到了最后的复杂的过程。

第二个$\text{Pass}$是$\text{Numerical Factorization}$，这就相对容易理解了，就是计算数值填入$\text{Analyse}$阶段计算出来的空位里。

最后一个阶段自然是求解了，这个没什么好说的。

$\text{2. Sparse Cholesky Factorization}$

相比于普通矩阵的$\text{LU}$分解来说，对称正定矩阵的$\text{Cholesky}$分解更容易处理一点，各个阶段的性质也更好，因此书上也首先从这里开始讲起。

目前并不打算细讲$\text{Supernodal Cholesky}$和$\text{Multifrontal Cholesky}$，这俩很重要，倒也不难理解，但是实在难写所以必不可能真的去全实现了。。。而且仔细看$\text{Eigen}$的时候发现$\text{Eigen}$自己也没有实现这些更高级的算法，大概说明基本算法对于很多场合来说也够用了但是人家有模板元编程黑魔法编译期计算邪道，你会吗。拿来算斯坦福的兔子这种规模还是差不多的。

不过可能会学习一下怎么并行，$\text{Cholesky}$分解中确实存在适合并行的部分，不过这个可能就需要求助于$\text{CS149}$而非数值线性代数了。

$\text{2.1. Analyse}$

这一步全是图论。。。最麻烦的在于整个过程需要在若干个不同的图或者树之间来回转换。如果对图论的结论感到困惑，结合矩阵消元会比较好理解。

首先按照书上的记号，作出一些符号上的约定以方便后续展开：

$\mathcal{S}(A)\equiv \{(i, j)|a_{ij}\neq 0\}$

这可以简单称为矩阵$A$的稀疏模式，简单来说，它记录了稀疏矩阵$A$的哪些地方是非零的，为了和书上保持一致，下标从$1$开始。

$\mathcal{G}(A)$定义为矩阵$A$定义的图。

$\text{2.1.1. 消元的图论解释}$

首先给出$\text{Parter's rule}$，这个定理给出了消元的图论解释，也架起了图论和矩阵之间的桥梁，说明了两者的等价性，这是理解接下来的许多结论的基础。如果图论形式难理解，那么就从矩阵去思考，如果矩阵难理解，就从图论去思考如果都难理解那就没办法了。

$\text{Thm.(Parter's rule) 设消元}k\text{步后的矩阵对应的图为}\mathcal{G}^k$

$则从\mathcal{G}^k到\mathcal{G}^{k+1}，删除顶点k，并且把k连接的顶点之间全部加上边。$

这个很容易验证。重复使用这个定理就可以得到如下结论：

$(i,j)\in E(\mathcal{G}(L+L^T)) \Leftrightarrow (i,j)\in\mathcal{G}(A)$ $\text{或}$ $\exists k<i,j \text{ s.t. } (i,k),(k,j) \in E(\mathcal{G}(L+L^T))$

$\text{2.1.2. Elimination Tree}$

在这一步，我们的目的是根据稀疏矩阵构建$\text{Elimination Tree}$（下称“消元树”），这是后续过程的重要基础。

如果我们对$A$作$LL^T$分解，那么我们得到的$L$的稀疏模式一定满足下面的性质：

$$\mathcal{S}(L_{i:n, j})\subseteq \mathcal{S}(L_{i:n, i}),i>j且l_{ij}\neq 0$$

这个性质被称为$\text{column repilcation princle}$，非常容易现。提示：试着对$A$作几次第二类初等行变换，会发现上三角的行向量随着消元被“带到”下三角来。

只要我们能够定义出每个点的父节点，那么自然就能得到树的结构。定义$p(j)=\min\limits_{i>j}\{l_{ij}\neq 0\}$，也就是说，对于每个节点对应的列，我们取$L$矩阵中对应的那一列在对角线以下的第一个元素，就是其在消去树中的父节点。这里的$p$代表$parent$。目前我们并不知道$p(j)$具体是谁，需要等到求出消去树以后反向求解出来。

按照映射合成的记号，$p^k(j)$就代表点$j$的$k$级祖先，当然我们还可以把树论的更多东西都搬过来，定义$anc(j)$为$j$的祖先集合，定义记号$\mathcal{T}(j)$表示消元树中$j$的子树。容易发现，$anc(j)$中的任意标号都比$j$大。

下面的定理说明了消元树的意义——消元树中的祖先关系表达了消元过程中列向量的稀疏模式从左往右（列号从小到大）传播的过程。

$\text{Thm. 设矩阵}A\text{对称正定},A=LL^T,\text{如果}a_{ij}=0,i>j$

$\text{那么}l_{ij}\neq 0\text{的充要条件是存在}k<j,a_{ik}\neq 0,t\ge 1 ,j=p^t(k)。$

这个定理表明了$A$的稀疏模式中的一个非零元$(i,k)$是如何在消元树中传播，被“带到”$(i,j)$来的。这是从矩阵消元的角度来看得到的结果，结合消元的过程理解更为方便。而如果我们从$A$代表的无向图来考虑，在图论的描述下也有一个类似的定理。

$\text{Thm. 对}$ $j < i, i\in anc(j)$ $\text{的充要条件是在}\mathcal{G}(A)$ $\text{中} i,j可以由标号不超过i的点构成的路径联通。$

如果$i$是$j$的祖先，这个结论是非常显然的。反过来的证明则可以使用归纳法。这个定理很显然和上面以矩阵形式表达的定理是等价的，其原因还是$\text{Parter's rule}$。

在上面这些结论的基础上，最后一个结论给出了由$A$构造消元树的方法。

$\text{Thm. }i=p(j) \Leftrightarrow i是\mathcal{G}(A)中与j联通的满足i>j的最小的标号$

根据上面的定理可以知道这是正确的。有了这个结论，构造算法就容易了——由小到大循环标号增量式构造，每次给当前的标号寻找儿子。

$\text{2.1.3. Coding: build elemination tree}$

构造消元树的算法并不难。假设我们已经建好了前$i$个点的消元树，在循环第$i+1$个点的时候，找出与$i+1$在图$\mathcal{G}(A)$中相邻且标号更小的点，这些点所在的子树将成为$i$的子树，找出这些子树的树根，并设其父亲为$i+1$即可。

下面代码表示了这一过程，这份代码将实装在$\text{SpmX}$稀疏计算库中，这个也是这学期的一个大作业。

void BuildFromMatrix(const SparseMatrix &A) {
  int *anc = new int[num_nodes_];
  for (uint i = 1; i < num_nodes_; i++) {
    for (uint j = A.OuterIndex(i);
         j < A.OuterIndex(i + 1) && A.InnerIndex(j) < i; j++) {
      uint x = A.InnerIndex(j);
      while(fa_[x] >= 0) {
        uint t = anc[x];
        anc[x] = i;
        x = t;
      }
      fa_[x] = anc[x] = i;
    }
  }
  delete[] anc;
}

anc用于记录子树树根，可以看到这里使用了路径压缩的技术，免得我们重复跳父亲来找根。

至于复杂度，最坏情况下每个非零元都会被循环一次，而路径压缩的复杂度，均摊至每次合并都是对数的，所以这个算法还是挺快的。

$\text{2.1.4. Sparsity Pattern of L}$

有了消元树，构造$L$的稀疏模式就比较容易了，为了方便不同的存储方式，我们分情况对于行向量和列向量构造稀疏模式。

需要注意的是，为了省去一些不必要的讨论，和书上不同，下面的算法中并未区分对角线元素是否应在稀疏模式中——这并不影响算法的理解。

$\text{2.1.4.1. Sparsity Pattern of Row Vectors}$

定义记号： $row_L(i) = \mathcal{S}(L_{i, 1:i-1})=\{j|j < i, l_{ij}\neq 0\}$

定义节点$i$的行子树$\mathcal{T}_r(i)=\mathcal{T}(i)\cap(\{i\}\cup \{k|a_{ik}\neq 0, k < i\})$

换句话说，我们保留$i$的子树中$A$的下三角中$i$这一行非零元的标号，这就是$i$的行子树了。

因此，$row_L(i)$必定是包含$A$的行子树节点的。然后根据前面的结论，若$a_{ij}=0$而$j\in row_L(i)$，则需要存在$k<j$使得$a_{ik}\neq 0$且$j\in anc(k)$，只需要把这些节点找出来就行。对于算法实现，我们对于第$i$行每个非零的标号$k$，尝试用$k$去更新其不超过$i$的祖先$j$即可。

一个显然的事情是，只有$\mathcal{T}_r(i)$的叶节点$k$会有用，因为在枚举这些$k$的时候就已经将所有的可能的$j$都更新了，不是叶节点的标号必然会在此前的某一轮循环中被它的子孙更新。所以枚举的时候可以不基于$A$的稀疏模式，而是构造$A$的骨架矩阵，其中移除了所有$j < i$且$j$不是$i$的行子树叶节点的非零元$(i,j)$。这样再构造就会快很多。

为了做优化，我们需要求出骨架矩阵，也就需要求出$\mathcal{T}_r(i)$中的非叶节点，这一点可以通过$\mathcal{T}$的后序遍历来做到。如果$\mathcal{T}(i)$中的某个节点$j(j<i)$是$\mathcal{T}_r(i)$的叶子节点，那么必然满足如下几点：

$\text{2.1.4.1. Sparsity Pattern of Col Vectors}$

有了前面的行向量的铺垫，列向量就很简单了。若$i\in col_L(j)$则必有$j\in row_L(i)$，也就是$a_{ij}\neq 0$或者$i,j$满足那么一大串的条件，我们将这个条件改为以$i$为主体叙述，那么也就是说，若$a_{ij}=0$并且$i\in col(j)$，则存在$k, a_{ik}\neq 0$且$k\in T(j)$。

这样写出来，就知道

$\text{Thm. }col_L(j)=\{i>j|\exists k\in T(j), (i,k)\in\mathcal{S}(A)\}$

换句话说，要求第$j$列的稀疏模式，我们只需要把$j$子树中所有点在$\mathcal{G}(A)$中的邻居找出来就行。这一点有很明显的子结构性质，就跟子树和差不多，所以可以自底向上来个小小的树形递推。

不妨记$j$在图$\mathcal{G}(A)$中的邻居为$adj_{\mathcal{G}(A)}(j)$，递推式为

$$col_L(j)=((adj_{\mathcal{G}(A)}(j)\cap (j,n])\cup\bigcup\limits_{k\in son(j)}col_L(k))$$

$ \text{2.1.5. Supernodes} $

$\text{Supernodes}$的优化主要就是找出连续的块，然后可以利用稠密$\text{Cholesky}$分解的良好局部性。

所谓$\text{supernodes}$，其实就是一段极大的连续的列$\{s, s + 1, ..., s+t-1\}$，满足如下性质：

$$col_L(s)\cup\{s\}=col_L(s+t-1)\cup [s, s+t-1]$$

如果把稀疏模式画出来就很明确，这段连续的列的稀疏模式一定是“上面一个稠密的三角阵，下面有一些零散的列数为$t$的稠密矩阵”。既然扯到稠密矩阵了，使用各种优化技巧就方便多了。

如果我们在消元树上把$\text{supernode}$包含的点缩成一个大点，得到的树被称为$\text{assembly tree}$，以下我称为“组装树”。下面的多波前法利用组装树就可以狠狠地优化。

既然$\text{supernodes}$有这么多好处除了难写，我们就希望能很方便地把它找出来，这一点只需要根据$col_L$的大小进行判断。从左往右进行扫描，切分出满足$|col_L(s)|=|col_L(s+t-1)|+t-1$的段即可！

利用$\text{supernodes}$就能整出大量的骚操作。

$\text{2.2. Factorize}$

前面我们构造出了$L$的稀疏模式，接下来的工作就是把数填进去。有许多稀疏$\text{Cholesky}$分解的实践，其主要不同就在于它们如何安排计算任务的顺序，这就影响到能否充分利用稠密矩阵计算的核函数，使用的内存以及并行的可能性。

$\text{2.2.1. Simplicial Cholesky}$

我们首先考虑怎么将稠密的$\text{Cholesky}$分解拓展至稀疏情况。为了简单起见，这里不考虑$\text{supernodes}$，就使用最简单的计算方法，也就是

$$L_{j + 1:n, j}=\frac{A_{j+1:n, j}-\sum\limits_{k=1}^{j-1}L_{j+1:n, k}l_{jk}}{l_{jj}}$$

$$l_{jj}=\sqrt{a_{jj}-\sum\limits_{k=1}^{j-1}l_{jk}^2}$$

于是由这个就可以看出来计算的依赖关系——当且仅当$l_{jk}\neq 0(j>k)$时计算第$j$列依赖于计算第$k$列，考虑到前面已经求出来了稀疏模式，就比较容易安排任务了。这是按列的计算方法，再具体到实现上，有$\text{left-looking}$和$\text{right-looking}$两种方式。

根据存储方式的不同，我们也可以按行来进行计算，这样的实现有时称为$\text{up-looking}$分解。此种方式的计算量对于高度稀疏的矩阵是渐近意义下最优的，但是实际中，这样的方式很难利用上$\text{BLAS}$（基础线性代数子程序库）。

如果计算需要以行的顺序进行，我们有如下关系式：

$$L_{1:i-1,1:i-1}L^T_{i,1:i-1}=A_{1:i-1,i}$$

$$l_{ii}^2=a_{ii}-L_{i,1:i-1}L^T_{i,1:i-1}$$

这里有两种处理方式。

如果我们想要利用此前求出来的行向量的稀疏模式，那么这里直接往里面填数就好。不过假设之前没做这个事情，也有一种替代方案可以现场计算出行稀疏模式，暂时不提。

$\text{2.2.2. Multifrontal Cholesky}$

多波前法倒是还挺流行的，主要原因可能是这个方法对有限元矩阵特攻，它也被用在专业的有限元分析软件中。后面的分析中，我们会看到它的思路和有限元中的“组装”非常相似。

首先有个简单的结论，把前面有关计算任务依赖性的叙述搬到了消元树上，其实也就是在$\text{2.1}$中反复提到的——稀疏模式像“波”一样在消元树上传播。

$\text{Thm. 第k列的计算数值只会影响其在消元树上的所有祖先。}$

这显然是对的。

$\text{2.3. Reordering}$

我们最希望的事情就是排序完成后作分解的时候不会填入过多的数，但是最小化这个量是$\text{NPC}$问题，所以我们只好选一些替代方案。

这就又是一个超级大的话题，这里暂时不加介绍了。