47、稀疏矩阵的Cholesky分解及相关算法

稀疏矩阵的Cholesky分解及相关算法

在处理线性方程组时，稀疏矩阵的Cholesky分解是一个重要的问题。本文将详细介绍稀疏矩阵Cholesky分解的相关内容，包括超级节点算法、存储方案以及并行实现等方面。

超级节点Cholesky分解算法

超级节点Cholesky分解算法是一种用于稀疏矩阵Cholesky分解的有效方法。该算法通过将与超级节点相关的计算组合在一起，提高了计算效率。

算法原理

该算法使用 smod() 过程对列进行修改，根据列是否属于超级节点，有两种修改情况：
- 当列 j 属于超级节点 J 时，列 j 仅由 J 中位于其左侧的列修改。
- 当列 j 不属于超级节点 J 时，列 j 由 J 的所有列修改。

以下是超级节点Cholesky分解算法的计算方案，也称为右视超级节点算法：

supernode_cholesky =
for each supernode J do from left to right {
    cdiv(first(J));
    for j = first(J) + 1, ..., last(J) {
        smod(j, J);
        cdiv(j);
    }
    for k ∈ Struct(L*(last(J)))
        smod(k, J);
}

该算法仍然从左到右计算 L 的列。与之前的算法不同的是，它将与超级节点相关的计算组合在一起。在超级节点级别，使用右视方案：对于超级节点 J 的第一列的计算，当与左侧所有列的修改完成后，只需要一次 cdiv() 操作。 J 的列以左视方式计算：在计算超级节点 J 左侧的所有超级节点后，由于右视超级节点方案， J 的列已经用这些超级节点进行了修改，因此列 j 的计算是先将其与 J 中位于其左侧的所有列进行修改，然后执行 cdiv() 操作。在计算完 J 的所有列后， J 右侧依赖于 J 列的所有列 k 都用 J 中的每一列进行修改，即通过 smod(k, J) 过程。

算法优势

超级节点算法的一个优点是增加了内存访问的局部性。因为超级节点 J 的每一列用于修改 J 右侧的几列，并且 J 的所有列用于修改 J 右侧的相同列。

稀疏矩阵的存储方案

由于稀疏矩阵中的大多数元素为零，因此使用特定的存储方案来避免存储零元素。这些压缩存储方案存储非零元素以及关于行和列索引的额外信息，以确定其在完整矩阵中的原始位置。

基本存储方案

稀疏下三角矩阵 L 以大小为$O(n + nz)$的压缩存储方案存储，其中 n 是 L 中的行数（或列数）， nz 是非零元素的数量。存储方案使用两个长度为 nz 的数组 Nonzero 和 Row ，以及三个长度为 n 的数组 StartColumn 、 StartRow 和 Supernode 。
- Nonzero 数组按列主序包含三角矩阵$L = (l_{kj}) {k \geq j}$的所有非零元素的值，即非零元素在一个线性数组中从左到右按列排序。
- StartColumn 数组隐式包含非零元素的对应列索引信息。位置 j 存储 Nonzero 数组中列 j 的第一个非零元素的索引，即 Nonzero[StartColumn[j]] 包含$l {jj}$。
- Row 数组包含 Nonzero 中对应元素的行索引。在不考虑超级节点的简单版本中， Row[r] 包含存储在 Nonzero[r] 中的非零元素的行索引，$r = 0, …, nz - 1$。

考虑超级节点的存储方案

当进一步利用同一超级节点中行的相似稀疏结构时，非零元素的行索引以 Row 和 StartRow 数组的组合方式存储：
- StartRow[j] 存储 Row 中列 j 的第一个非零元素的行索引，即 Row[StartRow[j]] = j ，因为$l_{jj}$是第一个非零元素。
- 对于每一列，行索引仍然存储在 Row 的连续块中。与简单方案不同的是，同一超级节点中不同行的块不是不相交的，而是根据这些列的相似稀疏结构重叠。

通过以下方式可以快速访问集合 Struct(L*j) ：

Struct(L*j) =
{Row[StartRow[j] + i] | 0 ≤ i ≤ StartColumn[j + 1]
 - StartColumn[j - 1]}

Supernode 数组用于管理超级节点：如果列 j 是超级节点 J 的第一列，则 J 的列数存储在 Supernode[j] 中。

共享变量的并行实现

对于稀疏Cholesky分解的并行实现，考虑使用共享内存机器。稀疏Cholesky分解有多种并行源，包括单操作 cmod(j, k) 或 cdiv(j) 内的细粒度并行性，以及左视、右视和超级节点算法中的列导向并行性。

消除树与并行性

稀疏矩阵 L 的稀疏结构可能导致额外的并行源，这在密集分解中是不可用的。当不同的列（以及对它们有影响的列）具有不相交的稀疏结构时，可以避免数据依赖。这种并行性可以通过消除树来描述，消除树使用 parent(j) 关系表达列之间数据依赖的特定情况。
- 对于每一列 j ，$0 \leq j < n$，定义 parent(j) 为：
plaintext parent(j) = min{i | i ∈ Struct(L*j)} if Struct(L*j) ≠ ∅
即 parent(j) 是列 j 的第一个非对角非零元素的行索引。如果 Struct(L*j) = ∅ ，则 parent(j) = j 。
- 对于$0 \leq i < n$，定义 children(i) 为：
plaintext children(i) = {j < i | parent(j) = i}
即 children(i) 包含所有在第 i 行有第一个非对角非零元素的列 j 。

有向图$G = (V, E)$，其中节点集$V = {0, …, n - 1}$，每个列对应一个节点，边集$E$，如果$i = parent(j)$且$i \neq j$，则$(i, j) \in E$。如果矩阵 A 是不可约的，则可以证明$G$是一棵树。

消除树$G$的一个重要性质是它指定了列必须评估的顺序： parent 的定义意味着如果$j = parent(i)$，则列 i 必须在列 j 之前评估。因此，列 j 的所有子节点必须在计算 j 之前完全评估。此外，列 j 不依赖于不在子树$G[j]$中的任何列。因此，如果$G[i]$和$G[j]$是不相交的子树，则列 i 和 j 可以并行计算。特别是，消除树的所有叶子节点可以并行计算，并且计算不需要从列0开始。因此，稀疏结构决定了要利用的并行性。对于给定的矩阵，高度较小的消除树通常比高度较大的树表示更大的并行度。

并行左视算法

并行左视算法（III）的实现基于 n 个列任务 Tcol(0), ..., Tcol(n - 1) ，其中任务 Tcol(j) ，$0 \leq j < n$，包括对所有$k \in Struct(Lj*)$执行 cmod(j, k) 操作以及执行 cdiv(j) 操作。这些任务由于非零元素而相互依赖。并行实现使用任务池来管理任务的执行。任务池有一个中央任务池，用于存储列任务，每个处理器都可以访问。每个处理器负责执行一部分列任务。任务到处理器的分配是动态的，即当处理器空闲时，它从中央任务池获取一个任务。

动态实现的优点是工作负载分布均匀，尽管由于稀疏结构，任务的执行时间可能不同。处理器对中央任务池的并发访问必须无冲突，以确保任务唯一分配给处理器执行。这可以通过锁定机制实现，使得在特定时间只有一个处理器访问任务池。

有三种并行实现变体：
- 变体L1 ：在所有$k \in Struct(Lj )$的列任务 Tcol(k) 完成之前，不将列任务 Tcol(j) 插入任务池。任务池可以初始化为消除树的叶子节点。并行度受树的独立节点数量限制，因为相互依赖的任务按顺序执行。因此，访问任务 Tcol(j) 的处理器可以执行该任务而无需等待其他任务完成。
- 变体L2 ：允许在不需要立即完成执行的情况下开始执行 Tcol(j) 。任务池初始化为所有可用的列任务。列任务由处理器从左到右动态访问，即空闲处理器访问尚未分配给处理器的下一列。在所有$k \in Struct(Lj )$的任务 Tcol(k) 完成之前开始计算列任务 Tcol(j) 。在这种情况下， Tcol(j) 的并非所有 cmod(j, k) 操作都能立即执行，任务只能执行那些对应的任务已经执行的$k \in Struct(Lj*)$的 cmod(j, k) 操作。因此，任务在执行过程中可能需要等待其他任务完成。

以下是变体L2的代码实现：

parallel left cholesky =
c = 0;
while ((j = get_unique_index()) < n) {
    for (i = 0; i < |Struct(Lj*)|; i++) {
        while (Sj empty) wait();
        k = #pop(Sj);
        cmod(j, k);
    }
    cdiv(j);
    for (i ∈ Struct(L*j)) : push(j, Si);
}

为了控制单个列任务 Tcol(j) 的执行，为每一列 j 分配一个数据结构 Sj ，其中包含所有$k \in Struct(Lj )$的列，对于这些列， cmod(j, k) 操作已经可以执行。当处理器完成列任务 Tcol(k) 的执行（通过执行 cdiv(k) ）时，它将 k 推送到每个$j \in Struct(L k)$的数据结构 Sj 上。由于不同的处理器可能同时尝试访问同一堆栈，因此必须使用锁定机制来避免访问冲突。执行 Tcol(j) 的处理器从 Sj 中弹出列索引 k 并执行相应的 cmod(j, k) 操作。如果 Sj 为空，处理器等待另一个处理器插入新的列索引。当从 Sj 中检索到$|Struct(Lj*)|$个列索引时，任务 Tcol(j) 可以执行最终的 cdiv(j) 操作。

• 变体L3 ：是L2的一种变体，考虑了消除树的结构。列不是严格从左到右分配给处理器，而是根据它们在消除树中的高度分配，即消除树中列 j 的子节点在其父节点 j 之前分配给处理器。这种变体试图按照完成其他列所需的顺序完成列任务，从而利用矩阵稀疏结构提供的额外并行性。

并行右视算法

并行右视算法（IV）的实现也基于任务池和列任务。列任务的定义与并行左视算法的任务不同：列任务 Tcol(j) ，$0 \leq j < n$，包括执行 cdiv(j) 和对所有$k \in Struct(L j)$执行 cmod(k, j) 操作，即列任务包括列 j 的最终计算以及对列 j 右侧所有依赖于 j 的列 k > j 的修改。任务池初始化为与消除树的叶子节点对应的所有列任务。一个非叶子任务 Tcol(j) 在对所有$k \in Struct(Lj )$执行 cmod(j, k) 操作并且可以执行最终的 cdiv(j) 操作时插入任务池。

以下是并行右视算法的代码实现：

parallel right cholesky =
c = 0;
initialize_task_pool(TP);
while ((j = get_unique_index()) < n) {
    while (!filled_pool(TP, j)) wait();
    get_column(j);
    cdiv(j);
    for (k ∈ Struct(L*j)) {
        lock(k); cmod(k, j); unlock(k);
        if (add_counter(ck, 1) + 1 == |Struct(Lk*)|) add_column(k);
    }
}

任务分配通过为每一列 j 维护一个计数器 cj 来实现。计数器初始化为0，并在每次执行 cmod(j, *) 操作后由相应的处理器使用无冲突的 add_counter() 过程递增。对于任务 Tcol(j) 的 cmod(k, j) 操作的执行，必须锁定列 k ，以防止其他任务同时修改同一列。当计数器 cj 达到值$|Struct(Lj*)|$时，任务 Tcol(j) 插入任务池。

右视实现与左视变体L2的区别在于 cmod() 操作的执行顺序和执行处理器。在L2变体中， cmod(j, k) 操作由计算列 k 的处理器通过将其推送到堆栈 Sj 上发起，但该操作由计算列 j 的处理器执行。此执行不需要在操作发起后立即执行。在右视变体中， cmod(j, k) 操作不仅由计算列 k 的处理器发起，而且由该处理器执行。

综上所述，稀疏矩阵的Cholesky分解及相关算法在处理线性方程组时具有重要的应用价值。通过合理选择存储方案和并行实现方式，可以提高计算效率和性能。不同的算法变体适用于不同的场景，需要根据具体情况进行选择。

稀疏矩阵的Cholesky分解及相关算法

算法对比与总结

为了更清晰地了解各种算法和存储方案的特点，下面通过表格对前面介绍的内容进行总结和对比。

算法/方案	特点	优势	劣势	适用场景
超级节点Cholesky分解算法	将与超级节点相关的计算组合，右视和左视结合计算列	增加内存访问局部性	实现相对复杂	稀疏矩阵且对内存访问效率有要求的场景
基本存储方案（不考虑超级节点）	使用 Nonzero 、 Row 、 StartColumn 数组存储非零元素及索引信息	简单直观	未充分利用超级节点稀疏结构	稀疏结构简单，无明显超级节点特征的矩阵
考虑超级节点的存储方案	结合 Row 和 StartRow 数组存储行索引，利用超级节点相似稀疏结构	存储更紧凑	理解和实现难度较大	具有明显超级节点结构的稀疏矩阵
并行左视算法 - 变体L1	任务按消除树独立节点顺序执行，任务池初始化为叶子节点	执行任务无需等待	并行度受独立节点数量限制	消除树独立节点较多的矩阵
并行左视算法 - 变体L2	任务可提前开始执行，动态从左到右分配任务	工作负载分布均匀	任务可能需等待其他任务完成	任务执行时间差异大的矩阵
并行左视算法 - 变体L3	考虑消除树高度分配任务，利用矩阵稀疏结构并行性	充分利用稀疏结构并行性	实现较复杂	具有明显稀疏结构且可通过消除树高度优化的矩阵
并行右视算法	列任务包含最终计算和右侧列修改，任务池初始化为叶子节点	操作执行顺序和处理器分配明确	需维护计数器和锁定列	对列修改操作顺序有要求的场景

流程图展示

下面通过mermaid格式的流程图来展示并行左视算法变体L2和并行右视算法的执行流程。

graph TD;
    A[开始] --> B[c = 0];
    B --> C{j = get_unique_index() < n};
    C -- 是 --> D{Si 非空};
    D -- 是 --> E[k = #pop(Si)];
    E --> F[cmod(j, k)];
    F --> G{是否完成所有 cmod(j, *)};
    G -- 是 --> H[cdiv(j)];
    H --> I{是否完成所有列任务};
    I -- 否 --> C;
    I -- 是 --> J[结束];
    D -- 否 --> K[等待];
    K --> D;
    C -- 否 --> J;

这个流程图展示了并行左视算法变体L2的执行流程，从初始化开始，不断获取列索引，根据堆栈状态执行 cmod 操作，完成后执行 cdiv 操作，直到所有列任务完成。

graph TD;
    A[开始] --> B[c = 0];
    B --> C[initialize_task_pool(TP)];
    C --> D{j = get_unique_index() < n};
    D -- 是 --> E{filled_pool(TP, j)};
    E -- 是 --> F[get_column(j)];
    F --> G[cdiv(j)];
    G --> H{k ∈ Struct(L*j)};
    H --> I[lock(k)];
    I --> J[cmod(k, j)];
    J --> K[unlock(k)];
    K --> L{add_counter(ck, 1) + 1 == |Struct(Lk*)|};
    L -- 是 --> M[add_column(k)];
    L -- 否 --> H;
    M --> N{是否完成所有列任务};
    N -- 否 --> D;
    N -- 是 --> O[结束];
    E -- 否 --> P[等待];
    P --> E;
    D -- 否 --> O;

这个流程图展示了并行右视算法的执行流程，包括任务池初始化、获取列任务、执行 cdiv 和 cmod 操作，以及根据计数器状态插入任务等步骤，直到所有列任务完成。

实际应用中的考虑因素

在实际应用中，选择合适的算法和存储方案需要考虑多个因素。

• 矩阵的稀疏结构 ：如果矩阵具有明显的超级节点结构，那么超级节点Cholesky分解算法和考虑超级节点的存储方案可能更合适。例如，在一些物理模拟场景中，矩阵的稀疏结构可能与物理系统的拓扑结构相关，存在明显的超级节点特征。
• 并行性需求 ：如果需要充分利用多核处理器的并行计算能力，那么并行算法是必要的。根据矩阵的消除树结构和任务执行时间的差异，可以选择不同的并行算法变体。例如，对于消除树独立节点较多的矩阵，并行左视算法变体L1可能更合适；而对于任务执行时间差异较大的矩阵，变体L2可能更能平衡工作负载。
• 内存使用 ：不同的存储方案对内存的使用效率不同。考虑超级节点的存储方案可以更紧凑地存储矩阵，但实现和理解难度较大。在内存资源有限的情况下，需要权衡存储方案的复杂度和内存使用效率。

总结与展望

稀疏矩阵的Cholesky分解及相关算法在处理大规模线性方程组时具有重要的应用价值。通过合理选择存储方案和并行实现方式，可以显著提高计算效率和性能。不同的算法变体适用于不同的场景，需要根据矩阵的具体特征和应用需求进行选择。

随着计算机硬件的不断发展，多核处理器和分布式计算环境的普及，未来对于稀疏矩阵算法的研究将更加注重并行性和可扩展性。同时，如何更好地利用矩阵的稀疏结构和超级节点特征，进一步优化存储方案和算法实现，也是一个值得深入研究的方向。此外，在实际应用中，如何将这些算法与具体的领域知识相结合，提高算法的实用性和性能，也是未来的研究重点之一。

总之，稀疏矩阵的Cholesky分解及相关算法是一个充满挑战和机遇的研究领域，对于推动科学计算和工程应用的发展具有重要意义。