Sparse QR 论文笔记

⭐⭐⭐Algorithm 9xx, SuiteSparseQR: multifrontal multithreaded rank-revealing sparse QR：

一、重排序与符号分析

1.1寻找单例 (1colamd)

符号分析的第一步是删除单例（singletons）。A的单例指一个列 $j$ 只有一个单独非零项 $a_{ij}$ ，其大小大于给定的阈值 $tol$ ；或者列 $j$ 没有任何非零项（这一情况只会出现在矩阵结构秩亏的情况下）。如果存在使对角无零的置换，则矩阵具有完全结构秩。

如果一个单例存在，则它的列 $j$ 移动至左侧。如果列 $j$ 的非零元单例为 $a_{ij}$ ,则行 $i$ 移动至顶部。然后把列j和行i 从矩阵A中移除，这一过程一直持续到没有单例存在。得到：

$R_{11}$ 这个时候为上三角或上阶梯形矩阵（如果秩亏）， $R_{11}$ 的对角元素都比阈值 $tol$ 更大，因此，R11在数值上预估具有完整的秩，这当然是一个近似值; $R_{11}$ 可能是秩不足，但这将不会被检测到。tol = $20(m+n) \epsilon max_{j}||A_{*j}||_2$ ，其中 $\epsilon$ 是机器舍入误差（双精度下大约是 $10^{-16}$ ）, $max_{j}||A_{*j}||_2$ 为矩阵A最大的列二范数，在实践中，这一数值得到了很好的体现。SuiteSparseQR 则默认为 -2 。

单例在实际应用中经常出现。然而，单例的确定依赖于阈值的数值大小。因此，如果将一个矩阵的符号分析用于另一个具有相同非零元分布模式但不同数值的矩阵分解，则可能不会有相同单例。

$R_{11}$ 的稀疏QR分解不需要做任何数值运算，在R11或R12中不需要进行填充。移动单例的时间仅有 O(n + |R11| + |R12|)。一旦单例被删除，剩余的矩阵 $A_{22}$ 被排序，分析，并使用多波前QR分解方法。

SuiteSparse使用Heath的方法通过“挤出”B的线性相关部分来处理秩亏，需要B分解后才能完成。

1.2 Fill-reducing ordring 填充简化排序 (chol_analyze_ordering 、colamd)

在找到单例之后(如果有的话)，剩下的排序和符号分析阶段只依赖于 剩余矩阵 的非零模式。下文中，统一让 A 表示去掉单例后的剩余矩阵。

如果矩阵 $A$ 的块三角形式（BTF）有单个块，那么因子R的非零模式与 $A^TA$ 的Cholesky分解 $L^T$ 的非零模式相同[Coleman et al. 1986; George and Heath 1980]。否则，LT的非零元分布模式是R模式的上界。

否则， $L^T$ 的模式是 R 的模式的最上部分。符号分析则是在进行完 $A^TA$ 的Cholesky分解后建立的。SuiteSparse使用CHOLMOD 来进行重排序和符号分析：

—— COLAMD on A，给 $A^TA$ 排序但没有形成它。有了这个选项，SuiteSparseQR在不形成 $A^TA$ 的情况下执行它的整个符号分析和数值分解。这是默认值。

——AMD on $A^TA$ 的显式模式。AMD是近似最小度置换的缩写。细节于算法还未弄清楚

Matlab 网站内对于AMD的介绍

——CHOLMOD Cholesky模块。

SuiteSparseQR现在只作用于排序矩阵AP，因此在本文后续内容中，P将被省略，A 表示采用填充-减化排序后的排列矩阵。

这些排序方法的渐近运行时间在数量方面没有已知的紧界，可以事先容易地计算。然而，实验结果表明，除了偶尔的例外，COLAMD和AMD所花费的时间分别与 $A$ 和 $A^TA$ 中的非零数大致成正比。这不是上界，而是基于佛罗里达大学稀疏矩阵集合中实际应用的近2000个矩阵实验的平均值。

对于任何多波前稀疏QR来说，适应任何填充简化排序通常都是非常简单的。

1.3 符号分解（analyze_ordering）

在找到一个减少填充的顺序之后，CHOLMOD通过分析A而不是显式地形成 $A^TA$ ，对排序形式的 $A^TA$ 进行符号Cholesky分解得到 $LL^T$ 。分析从Cholesky因子 $L$ 开始，但为了使讨论更加清晰，本文将只讨论R。

1.3.1 列消去树( etree & postorder)

分析首先计算 $A^TA$ 的消去树和 R的每一行的非零数(行计数)。

$A^TA$ 的消去树每一个列都对应一个节点：如果满足 i < k，且 $r_{ik} \neq 0$ 是非对角线上第 i 行的第一个非零元，那么节点 $i$ 的父节点是 $k$ 。 (因为从上往下更新，第 i 行中有非零元的列 j 会影响到之后行i 的更新，即存在依赖)

消去树是符号分析所需要的，它还描述了数值分解的数据依赖关系。 $A^TA$ 的消去树也称为A的列消去树，因为每个节点对应于A的一列。下文将使用后一个术语（列消去树）。这一步复杂度约为O(|A|)，在理论上，“粗略地”在最坏情况下是O(|A| log n)，在实践中，这一步本质上需要O(|A|)时间。

当A的BTF有多个块可以从行合并树中获得时，R的模式具有更简明的描述，以代替列消去树。然而，SuiteSparseQR使用[Heath 1982]的方法处理秩不足的矩阵，它要求R在 $A^TA$ 的Cholesky因数分解中容纳任何非零项。 这需要列消除树，因此SuiteSparseQR不能利用行合并树, 只使用了列消去树。

1.3.2 寻找 supernodes （super_super_symbolic2）

R中的超节点表示一组具有相同或几乎相同非零模式的相邻 t 行。 R中的每一个超节点都产生一个对应的波前矩阵，用于数值QR分解。波前阵是A的稠密子阵，其中执行了A的QR分解的一部分。它包含 A 的行的一个子集。一次分解后，它的前 t 行包含了R对应的超节点的行。它的前 t 列称为波前矩阵的关键列。波前矩阵通常比分解后的矩阵A小得多。

行计数和列消去树决定确切的超级节点。如果 i 的父结点是i+1且| $R_{i*}$ | = | $R_{i+1,*}$ | +1，那么两行 i 和 i+1 在同一个超节点中。这些事实意味着R中的两行具有相同的非零模式。这个过程不需要统计 R矩阵本身的非零模式(需要O(|R|)的时间来查找)，而只需要行数和列消除树(需要O(|A|)的时间来查找)

时间和内存的使用通常通过利用更宽松的超级节点来减少。如果 $R_{i,i+1:n}$ 和 $R_{i+1}$ 的非零模式相同或相似，那么相邻的两行被组合成松弛超节点。松弛超节点只依赖于 R 的每一行的非零数目。

一旦将这些行分组为超级节点，就会计算出每个超级节点的非零模式，所花费的时间与每个超级节点最上面一行的非零数量成比例。这也是来用表示R的supernodal非零分模式所需的整数内存空间。在任何情况下，O(|R|)是时间的一个松散上界;因为有了超级节点，它通常比这个要小得多。实际结果是用R的超节点的非零元分布模式表示。

1.3.3 波前矩阵树(qr_analyze)

列消除树有n个节点，一个节点对应A的一列。将同一超节点中的列分组在一起，就得到了一个波前矩阵树，每个超级节点/每个波前矩阵在波前矩阵树中有一个节点。对于熟悉超节点方法的读者，A的稀疏多波QR分解法的波前矩阵树与 $A^TA$ 的超节点Cholesky分解的超节点消除树相同[Ashcraft and Grimes 1989]。

1.3.4 重排序矩阵A的行(qr_analyze)

在CHOLMOD的超节点分析之后，SuiteSparseQR根据A每一行中最左非零的列索引对A的行进行排序。得到结果排序矩阵 $P_2A$ 。 A的第 i 行在排序矩阵中被组装到了波前矩阵，这个波前矩阵的关键列包含了第 i 行最左列的索引 j 。这一步的时间为O(|A|)。

1.3.5 确定每个波前矩阵的规模(qr_analyze)

每个波前矩阵的列数由R中对应超节点的首行非零元数目给出。为了求每个波前阵的行数，SuiteSparseQR模拟了稀疏多波前阵QR分解，确定了对每个波前矩阵进行HOUSEHOLD QR分解所需的工作量。找到每个波前矩阵的 staircase，即每列中最后一个非零项的行索引。如果矩阵A是秩亏的，则内存使用和工作量可能小于这个估计。这一步所花费的时间与表示R的超节点模式所需的整数数量成比例，时间通常比O(|R|)小得多。

每个波前阵的浮点数工作量估计被计算出来，这一过程也确定了每个波前矩阵的规模，以及后续连续数值分解需要的贡献块所需要的栈大小。

1.3.6 符号分解的时间和内存占用

用于分析的总时间和内存使用大致是O(|A| + |R|)，不包括排序时间(通常是O(|A|)在实践中)，其中|R|是表征矩阵R的超节点结构所需求的整数数目。这比形成ATA所需的时间和内存要少得多，特别是当m << n时。它也比数字分解步骤所花费的时间要少得多。

1.3.7 样例

图1 稀疏矩阵A 它的因子R，和它的列分解树

图1给出了一个小稀疏矩阵A的多波前稀疏QR分解的例子。A的行已经按照其最左边非零项的列索引排序。因子R在右边表示，R的每个超节点由相邻的行组成，这些行的对角线块是一个完整的上三角矩阵，其他行具有相同的非零模式（ $A->R$ ）。A中的水平线根据各行组合到的波前矩阵把行做了进一步细分。 A中显示的"."表示在对应的波前矩阵中这个位置的0将会变成结构上的非零元素。

列消除树显示在图的右下角，超级节点显示为圆形框。这个树中 $i$ 的父结点是由满足 $r_{ij}\neq0$ 的最小k满足 k >i 给出的（看R矩阵）。这个例子将在第3节继续，它展示了五个正面矩阵中的两个的组装和数字分解。

1.3.8 与之前的算法比较

本文引用的先前的多波前稀疏QR方法都没有使用单纯 $A^TA$ 而不生成它这一思路。只有一种使用ATA的消除树，但必须显式地形成ATA的非零模式。[Edlund 2002]没有形成 $A^TA$ ，但在排序方法中发现从R本身的模式树很像COLAMD;得到的R的树和模式不能适应后续的秩检测，因为它与A的列消除树不相同。

[Amestoy et al. 1996]指出，用这种策略很难预测存储主载体所需的确切内存量。然而，SuiteSparseQR的符号分析阶段精确地模拟了每个前阵在O(|A|+|R|)时间内的装配，以找到这个准确的数量，如果等级检测阈值是无效的(tol < 0)。

四、数值分解

对于其数值分解，SuiteSparseQR从先前的实现多波前稀疏QR分解中使用了大量的理论和算法[Matstoms 1994;Amestoy等，1996年;Lu和Barlow, 1996年;皮尔斯和刘易斯1997;Edlund 2002]。这里描述的是SuiteSparseQR中不同的关键特性，以及足以使本文自洽的足够背景知识。

4.1 波前矩阵组装

对于波前矩阵的组装和分解，SuiteSparse使用了一个MA49中策略3的修正版本[Amestoy et al. 1996] (refer to Figure 7 on page 284 of their paper).

这种方法如图2和图3所示。图2是与图1的列消去树的节点1和节点2对应的一个叶节点在波前矩阵树中的装配。

图2 一个叶子波前矩阵的组装和分解

在本例中，A的第1、2和3行在第1列中具有最左边的非零，第4、5和6行在第2列中具有最左边的非零。由于R的前两行具有相同的非零模式，因此可以在一个波前矩阵中处理。对应的列1和列2是这个波前阵的关键列(pivot)，而其他列(6、8和11)是外部的(exterior)或不重要的(non-pivotal)。

A中的项表示为x。在波前矩阵中，最左边非零项右侧的那些将要变为非零项的0元素用 “.” 来表示。假设列1和列2是线性无关的，对这个波前矩阵进行QR分解得到图2右半部分的结果，其中有两行R（元素标记为r）。

此时，一个3×3的上三角贡献块(右图的c元素)必须组装给父节点，5个household变换向量（记为h）也是一样。一般来说，贡献块是上梯形的。

图3 有三个子节点的波前矩阵的组装和分解

考虑front 1的父元素。父类(front 4)包含front 1的第一个非关键列(在本例中是列6 )作为其关键列之一。在图3中，父节点有三个子节点，其中一个是图2中给出的front 1。图3中只显示了这三个子元素的C块;注意其中一个是上梯形。

C块中的条目根据它们所在的块被赋予下标。父元素front 4的关键列是列5、6和7。A的第16到22行在波前矩阵4 的关键列中有一个最左边的非零项。组装的波前矩阵4的staircase是(4,8,12,14,15,16,17)，它给出了在波前矩阵的每一列中0的起始位置的行索引。

三个子节点贡献块和A的第16到22行被组装到前面;如图3的左下角所示。HOUSEHOLD QR分解显示在其右侧。一个4乘4的上三角形贡献块仍然需要由front 4的父组件组装。秩不足的情况(图3中最右额矩阵)将在下一节中讨论。

3.2 波前矩阵分解

一旦组装好波前阵，就可以得到波前阵的稠密HOUSEHOLD QR分解。如果波前阵的关键列都是线性无关的，这种分解方法与LAPACK [Anderson et al. 1999]中的DGEQRF方法几乎相同。SuiteSparseQR不调用DGEQRF;而是调用DGEQRF所依赖的功能，即：

(1) DLARFG，构造单个household向量，

(2)DLARF，应用单个household向量，

(3)DLARFT，构造一个household向量块的T矩阵，

(4) DLARFB，它应用了一个块状的household向量。

相比于DGEQRF，稠密的波前矩阵因式分解利用了波前矩阵左下角的零元素。作为一个例子，考虑图3中的波前矩阵。块大小参数b(默认为32)定义了波前矩阵中每个面板的列数。对于这个小示例，假设面板大小为2。第一个Household向量是用DLARFG计算的，并应用于当前面板(即波前矩阵的前两列)。只对前四行进行操作。接下来，可以找到将第6列标记为h的项目归零的Household变换，因为现在面板是完整的，使用DLARFT和DLARFB将它完全应用到波前矩阵的剩余部分。这个Householder更新块仅应用于前矩阵的前8行，因此利用了左下角的大部分零。相反，DGEQRF每次操作每个面板的所有17行。

波前矩阵以后序方式分解，以便贡献块可以容易地放置在堆栈上。为了组装和分解一个front，一个矩形波前矩阵的空间首先被分配在工作空间的上部，这个工作空间将包含R和Household 向量H(假设后者没有被丢弃)。接下来，子节点从堆栈中弹出(位于同一内存块的底部)并组装到波前矩阵。波前阵被分解，然后它的贡献块C被复制出来并推到堆栈上。最后，压缩R和H分量，只保留图3中所示的R和H分量(收回C块所使用的空间和波前矩阵左下角显式的零)。

SuiteSparse使用了[Heath 1982]的一种更简单的方法来处理秩缺陷，在目前的工作中扩展到稀疏多波前QR，并允许Q保留为一组Houehold向量。在Heath的方法中，如果R的对角线低于一个给定的阈值，那么通过Givens的旋转将整行归零，并将这一行从R中删除。结果会得到一个不再是上三角的“压缩”结果R。（论文第10页）。第二个Household变换将会消除第二列对角线以下的条目，将以对角线以下元素的范数替代对角线元素。如果该范数低于阈值tol(与符号分析中使用的“tol”相同)，则跳过此列的Household变换。整个列（对角线及对角线以下）的元素被视为0.第三列的Household变换消去了这一列第二项以下的所有元素。

结果显示在图3的右边。波前矩阵的R比预期少了一行，导致希思方法“压缩”R。在这种情况下，贡献块仍然是4乘4的，但通常它的大小会增加。C中的列数保持不变，即前面的非关键列数。如果C一开始是上梯形的(就像第三个波前阵的情况一样)，那么丢失一列将导致C的大小增加一行。工作空间中没有净增长，因为R比原来小了一行。这些观察导致了以下定理。

定理一:Heath在稀疏QR分解中处理秩亏矩阵的方法需要使用列消去树，且如果取R的非零模式作为 $A^TA$ 的Cholesky因子，那么这个方法不会对R造成填充.

SuiteSparseQR将Household变换应用到波前矩阵，而不是应用到行上的给定旋转，但结构上对秩缺陷R的影响与上述定理相同。

五、并行

在多波前稀疏QR分解中至少存在两个并行性的机会。第一个机会出现在波前矩阵树中。图1中给出的示例中,前三个波前矩阵可以并行分解(一个front计算R的前两行,和接下来两个front用于计算R的行3和行4)。使用这种级别的并行性显式的基于线程的SuiteSparseQR软件。

第二个机会出现在每个波前矩阵中。波前矩阵的QR分解依赖于LAPACK，而LAPACK又使用了level 3 BLAS [Dongarra et al. 1990]。大部分的供应商提供的BLAS,如英特尔R数学内核库(MKL),AMD ACMLBLAS,太阳性能库BLAS,以及Goto BLAS 可以利用共享内存多核架构,且不需要额外的编程工作的开发人员使用BLAS的包。

SuiteSparseQR通过使用Intel的Threading Building Blocks（TBB）软件来利用基于树的并行性，这个软件是用于在共享内存多核架构上用c++编写并行应用程序的。通常，基于TBB的应用程序只需要确定任务;TBB本身负责任务调度和同步。这就是SuiteSparseQR的调用办法。尽管TBB为互斥、原子操作、队列等提供了应用程序接口，但SuiteSparseQR不需要这些特性。

请注意，Intel MKL使用的是OpenMP [Chapman et al. 2007]，而不是Intel的 TBB，用于在BLAS中利用并行性。这种设计选择对TBB当前版本的性能有影响，将在第6节中说明。

设 $n_f$ 为波前矩阵的个数。并行分解的符号分析阶段的目标是将 $n_f$ 个波前矩阵分配给TBB tasks，而TBB任务通常比nf任务少。为了简化任务之间的同步，任务之间的关系被保存为一个树。可以将每个波前矩阵放在它自己的任务task中，但这可能导致TBB中的同步开销，特别是对于树底部非常小的波前矩阵。类似的策略也被用于稀疏Cholesky分解[Geist和Ng 1989]。从波前阵到任务的分配花费O(nf)时间。

对于每一个波前矩阵f，工作量以f为根的子树决定。如果一个节点在它子树的工作量大于 $max(\omega /\alpha ,\beta )$ 则称之为一个大节点（big node）；其中 $\omega$ 是整个QR分解的浮点计算量， $\alpha$ 和 $\beta$ 用户定义用来控制任务树粒度的参数。一般情况下， $\alpha$ 应该至少为内核数目的两倍， $\beta =10^6$ 。所有其他节点都是小的(small)。为了确保波前矩阵树是一个树而不是森林，将添加一个占位符 $n_f+1$ 节点，它是所有根节点的父节点，这个节点将直接被标记为大节点（big）。

1. 第一次分配将所有小节点分配给任务。假设front f是一个小节点但是它的父节点 p 是一个大节点。如果f是p的子节点中编号最小的，那么以f为根的子树中，所有的波前矩阵放置到一个新任务中。p的其他孩子中的小节点也被添加到这个任务，直到这个任务达到 $max(\omega /\alpha ,\beta )$ 的计算量。之后，再为p的子树创建一个新的任务。

2. 第二次分配将所有大节点分配给任务，顺序从1到 $n_f+1$ ，这样大节点f的所有子节点在考虑大节点f本身之前都被分配给它们的任务。如果f的所有子结点被分配到相同的任务，那么f也被分配到相同的任务。f 没有子任务，或者如果它的子节点分配给不同任务，则创建一个 f 分配的新任务。

3. 最后，可以找到每个任务的堆栈大小：如果一个任务是另一个任务的祖先，那么两个任务可以共享一个堆栈，因此堆栈的数量与任务树中的叶子数量相同。为树的每个叶分配一个堆栈可以确保不同的TBB任务之间没有内存冲突。

在数值分解阶段，任务使用的所有工作区(包括贡献块堆栈的集合)在TBB调度任务之前被分配。所以在分解的TBB并行阶段，不再需要动态内存分配。每个任务分解波前矩阵树的一个子树中的所有波前矩阵，结果留在该任务使用的堆栈中。不需要同步，除非一个任务可以在它所有的子任务完成后开始;这是通过指定一组任务及其对TBB的依赖关系来处理的，而TBB处理所有的同步和调度。

每个任务中的额外并行通过3级BLAS的多核实现来利用。

并行效果

TBB任务树的叶节点通常可以包含较小的波前矩阵，在BLAS中并行性的作用范围很小。相比之下，树的根节点提供了不需要依赖树结构的并行，但是可以充分利用blas的并行性。这两种并行应该是互补的，但在目前的TBB实现中，它们是相互竞争的。TBB的未来实现希望解决这个问题 [Robison 2008] 。由 TBB 和 OpenMP 创建的线程是独立的，每个 TBB 线程将创建自己的一组 OpenMP 线程(尽管对于小型矩阵，基于OpenMP的BLAS只使用一个线程) 。因此，所使用的最大线程总数是这两组线程的乘积。在未来的实现中，TBB应该使用TBB任务树的叶子和下面部分的所有线程，而BLAS应该使用TBB任务树的根节点的所有线程。中间应该用一种混合物。

实验结果表明，与只在每个波前矩阵中利用基于BLAS的并行性相比，TBB利用基于树的并行性通常能够更有效地利用所有核