稀疏技术

稀疏矩阵的存储

首先何谓稀疏矩阵，就是在矩阵中有众多的零元素。稀疏矩阵可以用稀疏度来进行定量判定。稀疏度的计算公式如下：
$$\rho =\frac{\tau }{m\times n}$$
其中，矩阵的尺寸为 m x n， $\tau$为非零元素的个数 。

稀疏矩阵的存储应该满足如下条件：

1. 不存储 0 元素
2. 能够快速恢复到矩阵形态
3. 存储非零元素的数值和位置。

共有三种存储方式：

1. 散居存储
2. 按列/行存储
3. 三角存储

散居存储

待存储的表格如下：

此种存储方式记录非零元素在矩阵中的行、列位置，以及该元素的数值大小。

按列存储

此种存储方式存储了非零元素的值。并存储每一个元素的列，形成一个列索引；同时存储每一行第一个元素在列索引的位置。

三角存储

对于任何矩阵B，我们都进行 LDU 分解：A= LDU，其中 L 为下三角矩阵、D 为对角线矩阵，U 为上三角矩阵。

于是，对这三个矩阵，特别是 L 和 U，我们只需要用上面两个方法进行存储就行了。比如对 U ，可以采用按列的方式存储：

后面那两个 4 ，代表不存在，你也可以用 0 来表示。

对于 L ，也可以用同样的方法进行存储。而对于 D ，我们只需要存储其值就行。

LDU 分解

对于一个矩阵 A ，我们知道可以对他进行 LDU 分解，从而得到一个下三角、对角线、以及上三角矩阵。在 Python 的 scipy 库中，可以使用 linalg.lu 这个函数进行 PLU 分解。但这个与 LDU 分解不同，他得到的是一个普通的矩阵 P，以及L 和 U。

要进行 LDU 分解，有人可能会觉得：

这种想法大错特错，因为 A =LPU，这三个矩阵的点积明显不是简单拼接起来那么容易。

对于稀疏矩阵而言，这种对称矩阵的 LDU 分解可以采用一种图论的方式求取，下面我们将具体介绍。

对稀疏矩阵的 LDU 分解法

假设 A 是一个对称的稀疏矩阵。由于对称性，我们再对 A 进行 LDU 分解时仅需要讨论上三角部分。为此，我们将 A 的上三角表示成一个带权重的有向无环图，如下所示。

其中，节点和箭头表示非 0 元素的位置，边上的权重表示了数值。比如上面有 节点 1 指向节点 7，代表着有一个非 0 元素在第 1 行的 第 7 列，权重 *代表着这个非 0 元素的大小是 * 。

有了这个概念后，我们再来讨论 A 的 LDU 分解，假设有一个对称矩阵如下，

将其表现为带权重的有向无环图：

首先取任意一个节点，如 p；对 p 发出的每一条边（自边除外），用如下公式进行规则化：
$$a_{pj}=\frac{a_{pj}}{a_{pp}}$$
对其余发出的边，如 pk、pl 均如上计算，更新 $a_{pj},a_{pk},a_{pl}$ 的值。

然后，再消去节点 p，这个消去是“表面消去”，是指消去计算后，节点 p 不再参与后面的计算过程。消去计算算法如下：对于任意两条 p 发出的边，例如 pj 和 pk，修正值 $a_{jk}$，如下：
$$a_{jk}=a_{jk}-a_{pj}{a}_{pk}{a}_{pp}$$
对于 pj 和 pl，我们可以看到图中是没有边与