使用图着色(graph coloring)加速Gauss-Seidel迭代

Jacobi和Gauss-Seidel都是用于求解线性方程组的迭代方法。Jacobi迭代很容易实现在GPU上并行的版本，但是收敛速度很慢甚至有可能发散。Gauss-Seidel迭代并不能用类似Jacobi的方法并行，因为它是一个依赖于顺序的迭代方法。但是Gauss-Seidel迭代的收敛速度比Jacobi快很多。因此出现了用图论的图着色算法加速Gauss-Seidel迭代的研究。

能用图着色算法加速Gauss-Seidel迭代的原因

这一节翻译自Parallelizing the Gauss-Seidel Method using Graph Coloring。

可以发现对于某些类型的矩阵，实现并行版本的Gauss-Seidel法非常简单。考虑下面这个线性方程组：

\[ \begin{align*} 2x_1 - 1x_2 + 0x_3 + 0x_4 &= 3 \\ -1x_1 + 2x_2 - 1x_3 + 0x_4 &= -4 \\ 0x_1 - 1x_2 + 2x_3 - 1x_4 &= 4 \\ 0x_1 + 0x_2 - 1x_3 + 2x_4 &= -3 \\ \end{align*} \]

假设解的初始值是\(x_1 = x_2 = x_3 = x_4 = 0\)。如果我们要完成一次Gauss-Seidel迭代，我们首先用第一个方程求解\(x_1\)得到：

\[ \begin{align*} x_1 = \frac{1}{2}(3 - 0x_4 - 0x_3 + x_2 - 2x_1) = \frac{3}{2} = 1.5\\ \end{align*} \]

于是得到更优的解\(x_1 = 1.5, x_2 = 0, x_3 = 0, x_4 = 0\)。一般接下来我们会用第二个方程求解\(x_2\)，但是我们不这么做，而是求解线性系统里的第三个方程:

\[ \begin{align*} x_3 = \frac{1}{2}(4 + x_4 + x_2 - 0x_1) = \frac{4}{2} = 2\\ \end{align*} \]

于是得到更优的解\(x_3 = 2\)。注意到，\(x_1\)的值对于更优的\(x_3\)并没有任何影响，因为上面这个式子里\(x_1\)的系数是0。类似的\(x_3\)的值对于更优的\(x_1\)也没有任何影响。这两个计算的顺序无所谓，因此可以并行计算。同时不难看出，更优的解\(x_2\)和\(x_4\)也可以并行计算。因此，对于上面这个线性系统，并行的一次迭代非常容易实现：首先并行地计算更优的\(x_1\)和\(x_3\)，然后类似地，并行地计算更优的\(x_2\)和\(x_4\)。

我们已经发现了如果我们能找到可以独立求解的多个变量的划分，那么Gauss-Seidel算法确实可以并行化。我们的下一个任务是发明一种找到这种划分的算法。首先，以矩阵形式重写线性系统：

\[ \begin{equation*} \begin{bmatrix} 2 & -1 & 0 & 0 \\ -1 & 2 & -1 & 0 \\ 0 & -1 & 2 & -1 \\ 0 & 0 & -1 & 2 \\ \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \\ \end{bmatrix} = \begin{bmatrix} 3 \\ -4 \\ 4 \\ -3 \\ \end{bmatrix} \end{equation*} \]

设这个矩阵为\(M\)。\(x_1\)和\(x_3\)可以独立的求解的原因是\(m_{1,3}\)和\(m_{3,1}\)是0。很明显，\(x_1\)和\(x_2\)不能独立求解，因为\(m_{1,2}\neq 0, m_{2,1}\neq 0\)。也就是说，如果\(m_{i,j} \neq 0\)，那么\(x_i\)和\(x_j\)不能并行地求解。这个情况可以用图来可视化。让变量\(x_1, x_2, x_3, x_4\)对应图中的节点1, 2, 3, 4。仅当\(m_{i,j} \neq 0\)时，节点i和节点j之间有一条边。

在图里可以看到，节点1和节点3之间没有边，因此形成一个划分，可以被独立的求解。找到这些划分的算法即可以找到这些互不相连的节点组的算法。这个工作早已存在，被称为图着色算法。对于上面的图，一种可行的染色为

可以观察到，相同的颜色的节点属于不同的划分。两个相邻的节点不能有相同颜色，这保证互相依赖的变量不能放在同一个划分里。总结一下，找到一个可行的染色方案等价于划分这些变量。

论文Vivace: a practical gauss-seidel method for stable soft body dynamics里使用了图着色算法对Gauss-Seidel迭代进行加速。

图着色算法

太长不看版

1. 为每个顶点\(v\)分配一个\(P_v=\{0, 1, 2, \ldots, \Delta_v/s\}\)。将所有顶点放入集合\(U\)。
2. 重复以下步骤直到\(U\)为空：
   1. 为集合\(U\)每个顶点\(v\)随机选择一个颜色\(c(v)\)。
   2. 检查是否有邻居有相同的颜色，如果没有，接受这个颜色并从邻居的调色板中移除这个颜色。
   3. 如果某顶点的调色板\(P_v\)用完，\(P_v\)加入一个新的颜色。
   4. 从\(U\)中移除已经被染色的顶点。

定义

• \(\Delta\): 一个图的最大度数
• \(V\): 顶点集合
• \(U\): 未染色的顶点集合。
• \(\Delta_v\) : \(v\) 的度数
• \(s\): 调色板收缩因子，\(s > 1\)，该因子对于整个图来说是恒定的。

用图代表animated的对象离散成的三角或四面体网格。虽然这些网格由成百上千的顶点通过约束连起来，但是\(\Delta\)小得足以进行实时染色。

算法输入

无向图\(G\)，对应于矩阵\(Y\)。

初始化阶段(step 1-3)

给每个顶点\(v\)一个调色盘(有 \(\Delta_v\)/s 个颜色的列表)\(P_v\)。颜色由连续的自然数标识。着色过程开始并重复，直到所有顶点都被着色。每个着色轮包含三个平行步骤。

暂定着色(tentative coloring)回合（step 5-6）

一种 \(c(v)\) 颜色被随机从调色盘\(P_u\)中抽取并赋值给\(v\)。

解除冲突(conflict resolution)阶段（step 7-12）

每个顶点检查是不是没有邻居有选择相同的暂定颜色。如果确实没有，\(v\) 的染色被接受并且 \(c(v)\) 被移除出其邻居的调色板。

为了加速conflict resolution阶段，作者使用了Hungarian heuristic:在冲突的情况下，如果一个节点的索引比邻居的索引都大，这个染色就被认为是合法的。他们发现这种策略可以大大减少算法对图中所有顶点进行着色所需的着色轮数。

喂饱饥饿的人(feed the hungary) 阶段（step 14-16）

一个颜色会被加入到把颜色用完的调色盘。可允许的颜色的最大数量是 \(\Delta_v+1\)，但是在实践中作者从未到达过最大的阈值。

收缩因子的效果是减少染色的数量，但是将其设的太大会导致在减少颜色数量时没有什么收获从而造成更慢的染色。在作者的例子里，他们发现使用图的最小度数作为收缩因子可以得到最好的染色结果。

完整算法

\[ \begin{aligned} &\text { Algorithm } 2 \text { Vivace Graph Coloring Procedure [Grable and Panconesi 2000] }\\ 1:&U \leftarrow V \hspace{8cm} \triangleright \text { Initialization }\\ 2:&\textbf { for all} \text{ vertex } v \in U \textbf { do }\\ 3:&\qquad P_v \leftarrow\left\{0, \ldots, \Delta_v / s\right\}\\ 4:&\text { while }|U|>0 \textbf { do }\\ 5:&\qquad \textbf { for all vertices } v \in U \textbf { do } \hspace{4.3cm} \triangleright \text { Tentative coloring }\\ 6:&\qquad \qquad c(v) \leftarrow \text { random color in } P_v\\ 7:&\qquad I \leftarrow \varnothing\\ 8:&\qquad \textbf { for all} \text{ vertices } v \in U \textbf { do } \hspace{4.5cm} \triangleright \text { Conflict resolution }\\ 9:&\qquad \qquad S \leftarrow\{\text { colors of all the neighbors of } v\}\\ 10:&\qquad \qquad \textbf { if } c(v) \notin S \textbf { then }\\ 11:&\qquad \qquad \qquad I \leftarrow I \cup\{v\}\\ 12:&\qquad \qquad \qquad \text { remove } c(v) \text { from palette of neighbors of } v\\ 13:&\qquad U \leftarrow U-I\\ 14:&\qquad \textbf { for all vertices } v \in U \textbf { do } \hspace{4.3cm} \triangleright \text { Feed the hungry }\\ 15:&\qquad \qquad \textbf { if }\left|P_v\right|=0 \textbf { then }\\ 16:&\qquad \qquad \qquad P_v \leftarrow P_v \cup\left\{\left|P_v\right|+1\right\} \end{aligned} \]

例子

考虑一个更复杂的线性系统：

\[ \begin{equation*} \begin{bmatrix} 3.0& 0.0& 0.0& 0.0& 0.0& 0.0& 0.0& 3.0& 0.0& 0.0\\ 0.0& 15.0& 7.0& 0.0& 0.0& 0.0& 8.0& 0.0& 0.0& 0.0\\ 0.0& 7.0& 28.0& 9.0& 9.0& 0.0& 0.0& 0.0& 3.0& 0.0\\ 0.0& 0.0& 9.0& 18.0& 0.0& 0.0& 0.0& 0.0& 9.0& 0.0\\ 0.0& 0.0& 9.0& 0.0& 20.0& 0.0& 0.0& 8.0& 0.0& 3.0\\ 0.0& 0.0& 0.0& 0.0& 0.0& 8.0& 0.0& 0.0& 0.0& 8.0\\ 0.0& 8.0& 0.0& 0.0& 0.0& 0.0& 14.0& 0.0& 6.0& 0.0\\ 3.0& 0.0& 0.0& 0.0& 8.0& 0.0& 0.0& 11.0& 0.0& 0.0\\ 0.0& 0.0& 3.0& 9.0& 0.0& 0.0& 6.0& 0.0& 18.0& 0.0\\ 0.0& 0.0& 0.0& 0.0& 3.0& 8.0& 0.0& 0.0& 0.0& 11.0\\ \end{bmatrix} \end{equation*} \]

用上面的随机算法得到以下染色结果：

由于矩阵中有许多0，节点之间的边很少，这导致划分的数量较少。具有许多0的矩阵表示为稀疏矩阵，这种矩阵在实际应用中非常常见。最后，算法发现的染色结果使用的颜色很少，这是符合预期的，因为颜色越少，划分就越少。

上图有三个划分，每个划分都可以并行求解（总共只需要串行的执行3个GPU kernel来完成一次Gauss-Seidel迭代）。但是，请注意，如果我们改用Jacobi来求解十个未知数，那么所有十个未知数都可以在执行一次并行的迭代后找到。使用并行高斯-赛德尔方法，我们首先调度一个并行计算来求解第一个划分中的未知数，然后我们必须等待这个计算完成，然后才能为第二个划分调度另一个计算。从中可以看出，与雅可比方法相比，我们的并行化高斯-赛德尔方法表现出较低的并行度。然而，更快的收敛速度将弥补这一点。此外，稀疏矩阵在应用中非常常见，这样的矩阵导致划分数量较少。