论文阅读笔记：Parallel Iterative Solvers for Real-time Elastic Deformations (迭代法求解方程组 / 弹性形变仿真)-526互联

材料来源于 Siggraph Asia 2018 的 course note Parallel iterative solvers for real-time elastic deformations, SIGGRAPH Asia 2018 Courses, 2018.

0. 概述

在形变仿真中，许多类方法/模型最终归结为对方程组 \(\boldsymbol{A}\boldsymbol{x} = \boldsymbol{b}\) 的求解，其中 \(\boldsymbol{A} \in \mathbb{R}^{N\times N}\) ， \(\boldsymbol{b} \in \mathbb{R}^{N}\) ， \(\boldsymbol{x} \in \mathbb{R}^{N}\) 是待求解量，即仿真对象的形变位移。由于上述方程组通常维度巨大（几千/几万），直接求解十分耗时。因此，通常可采用迭代法求解。下面将介绍常用的 Jacobi 方法、Gauss-Seidel 方法、以及更多的高级方法。

1. 迭代线性求解器 Iterative Linear Solver

首先，构建迭代形式。 已知线性方程组 \(\boldsymbol{A}\boldsymbol{x} = \boldsymbol{b}\) ，首先，通过某种方式划分矩阵，即 \(\boldsymbol{A} = \boldsymbol{B} - \boldsymbol{C}\) ，由此可得 \(\boldsymbol{B}\boldsymbol{x} - \boldsymbol{C}\boldsymbol{x} = \boldsymbol{b}\) ，即 \(\boldsymbol{B}\boldsymbol{x} = \boldsymbol{C}\boldsymbol{x} + \boldsymbol{b}\) 。那么，标准的迭代法，如 Jacobi 和 Gauss-Seidel ，可以记作如下形式（将上式的左右两端记作前一次迭代和后一次迭代的值）

\[\boldsymbol{B}\boldsymbol{x}^{(k+1)} = \boldsymbol{C}\boldsymbol{x}^{(k)} + \boldsymbol{b} \quad (\text{即} \quad \boldsymbol{x}^{(k+1)} = \boldsymbol{B}^{-1}(\boldsymbol{C}\boldsymbol{x}^{(k)} + \boldsymbol{b})) \]

接下来，对迭代误差进行分析。 在迭代过程中，误差可以记作如下形式

\[\begin{aligned} \boldsymbol{e}^{(k+1)} &= \boldsymbol{x}^{(k+1)} - \boldsymbol{x} \\ &= \boldsymbol{B}^{-1}(\boldsymbol{C}\boldsymbol{x}^{(k)} + \boldsymbol{b}) - \boldsymbol{x} \\ &= \boldsymbol{B}^{-1}\boldsymbol{C}\boldsymbol{x}^{(k)} + \boldsymbol{B}^{-1}\boldsymbol{b} - \boldsymbol{x} \end{aligned}\]

将 \(\boldsymbol{B}\boldsymbol{x} = \boldsymbol{C}\boldsymbol{x} + \boldsymbol{b}\) （即 \(\boldsymbol{x} = \boldsymbol{B}^{-1}\boldsymbol{C}\boldsymbol{x} + \boldsymbol{B}^{-1}\boldsymbol{b}\)）带入上式，可得

\[\boldsymbol{e}^{(k+1)} = \boldsymbol{B}^{-1}\boldsymbol{C}\boldsymbol{x}^{(k)} - \boldsymbol{B}^{-1}\boldsymbol{C}\boldsymbol{x} = \boldsymbol{B}^{-1}\boldsymbol{C}\boldsymbol{e}^{(k)} \]

如果我们将上式中的 \(\boldsymbol{B}^{-1}\boldsymbol{C}\) 做特征值分解 eigenvalue decomposition，即 \(\boldsymbol{B}^{-1}\boldsymbol{C} = \boldsymbol{Q}\boldsymbol{\Lambda}\boldsymbol{Q}^{-1}\)，那么，第 \(k\) 次的迭代误差为

\[\boldsymbol{e}^{(k)} = (\boldsymbol{B}^{-1}\boldsymbol{C})^{k}\boldsymbol{e}^{(0)} = \boldsymbol{Q}\boldsymbol{\Lambda}^{k}\boldsymbol{Q}^{-1}\boldsymbol{e}^{(0)} \]

由此可见，迭代方法线性收敛，且其收敛速度取决于特征值 \(\boldsymbol{\Lambda}\) 的大小。

最后，将介绍 Jacobi 和 Gauss-Seidel 两种迭代形式。 将矩阵 \(\boldsymbol{A}\) 的对角元素记作 \(\boldsymbol{D}\) ，将上三角部分和下三角部分分别记作 \(\boldsymbol{U}\) 和 \(\boldsymbol{L}\) ，则有 \(\boldsymbol{A} = \boldsymbol{D} - \boldsymbol{U} - \boldsymbol{L}\) 。通过将上述矩阵组合成不同的 \(\boldsymbol{B}\) 和 \(\boldsymbol{C}\)，迭代方式可转化为如下形式：

1.1 Jacobi 迭代求解器

将上述迭代中的矩阵分别记作 \(\boldsymbol{B} = \boldsymbol{D}\) 以及 \(\boldsymbol{C} = \boldsymbol{U} + \boldsymbol{L}\) 上述迭代形式可转化为如下：

\[\boldsymbol{x}^{(k+1)} = \boldsymbol{D}^{-1}(\boldsymbol{U} +\boldsymbol{L})\boldsymbol{x}^{(k)} + \boldsymbol{D}^{-1}\boldsymbol{b} \]

上述方程右侧全部为已知量，且 \(\boldsymbol{x}^{(k+1)}\) 的各行可以分别求解计算。因此，Jacobi 迭代法非常适合并行计算。将上述公式展开来，则 \(\boldsymbol{x}\) 中的第 \(i\) 个变量迭代计算如下

\[x^{(k+1)}_{i} = \frac{1}{d_{i}}\left(b_{i} - \sum^{i-1}_{j=1}l_{ij}x^{(k)}_{j} - \sum^{N}_{j=i+1}u_{ij}x^{(k)}_{j}\right) \]

其中，\(d_i\) \(b_i\) \(l_{ij}\) \(u_{ij}\) 分别是 \(\boldsymbol{D}\) \(\boldsymbol{b}\) \(\boldsymbol{L}\) \(\boldsymbol{U}\) 中的元素。但是，该方法需要大量的迭代次数，才能得到较小的误差，综合来看，仍然需要很常的计算时间，并不适合实时仿真计算。（后续会介绍加速方法）

1.2 Gauss-Seidel 迭代求解器

将上述迭代中的矩阵分别记作 \(\boldsymbol{B} = \boldsymbol{D} - \boldsymbol{L}\) 以及 \(\boldsymbol{C} = \boldsymbol{U}\) 上述迭代形式可转化为如下：

\[\boldsymbol{x}^{(k+1)} = （\boldsymbol{D} - \boldsymbol{L})^{-1}\boldsymbol{U}\boldsymbol{x}^{(k)} + (\boldsymbol{D} - \boldsymbol{L})^{-1}\boldsymbol{b} \]

对上式进行转化，可得到如下

\[(\boldsymbol{D} - \boldsymbol{L})\boldsymbol{x}^{(k+1)} = \boldsymbol{U}\boldsymbol{x}^{(k)} + \boldsymbol{b} \]

进而有（这样转化的一个原因是，\(\boldsymbol{D}\) 的逆可以直接得到，而其他形式的矩阵求逆十分困难。）

\[\boldsymbol{D}\boldsymbol{x}^{(k+1)} = \boldsymbol{L}\boldsymbol{x}^{(k+1)} + \boldsymbol{U}\boldsymbol{x}^{(k)} + \boldsymbol{b} \quad \text{即} \quad \boldsymbol{x}^{(k+1)} = \boldsymbol{D}^{-1}\boldsymbol{L}\boldsymbol{x}^{(k+1)} + \boldsymbol{D}^{-1}\boldsymbol{U}\boldsymbol{x}^{(k)} + \boldsymbol{D}^{-1}\boldsymbol{b} \]

将上述公式展开后，则 \(\boldsymbol{x}\) 中的第 \(i\) 个变量迭代计算如下

\[x^{(k+1)}_{i} = \frac{1}{d_{i}}\left(b_{i} - \sum^{i-1}_{j=1}l_{ij}x^{(k+1)}_{j} - \sum^{N}_{j=i+1}u_{ij}x^{(k)}_{j}\right) \]

其中，\(d_i\) \(b_i\) \(l_{ij}\) \(u_{ij}\) 分别是 \(\boldsymbol{D}\) \(\boldsymbol{b}\) \(\boldsymbol{L}\) \(\boldsymbol{U}\) 中的元素。Gauss-Seidel 迭代的收敛速度要优于 Jacobi 迭代。但是，由上可知，上述计算中，第 \(i\) 个 \(x^{(k+1)}\) 的值依赖于第1至 \(i-1\) 个\(x^{(k+1)}\) 的值，因此，Gauss-Seidel 迭代只能顺序计算，无法实现并行计算。（后续将介绍如何将上述未知的 \(x^{(k+1)}_{i}\) 划分成不同的独立 sets，然后实现并行计算。）

2. Projective and Position-based Dynamics

（待整理）

3. Parallel Descent Method

下面详细介绍非线性优化问题的一类求解方法：下降方法。主要包括优化问题的描述、常见的几类下降方法、以及一种新的下降方法（对并行计算更加友好、计算效率更高）。各类下降方法的基本原理类似，只是在处理搜索方向、步长等各个环节时有所不同。

3.0 问题描述

在形变仿真中，将模型网格节点的位置和速度分别记作 \(\boldsymbol{q} \in \mathbb{R}^{3N}\) 和 \(\boldsymbol{v} \in \mathbb{R}^{3N}\)，通过隐式时间积分方法（implicit time integrationi）计算模型网格从 \(t\) 时刻到 \(t+1\) 时刻的形变，为

\[\boldsymbol{q}_{t+1} = \boldsymbol{q}_{t} + h \boldsymbol{v}_{t+1}, \quad \boldsymbol{v}_{t+1} = \boldsymbol{v}_{t} + h \boldsymbol{M}^{-1}\boldsymbol{f}(\boldsymbol{q}_{t+1}) \]

其中，\(\boldsymbol{M} \in \mathbb{R}^{3N \times 3N}\) 是质量矩阵，\(h\) 是时间积分步长，\(\boldsymbol{f} \in \mathbb{R}^{3N}\) 是作用在网格节点上的合力（包括内部弹性力和外部作用力）。上述方程可合并为

\[\boldsymbol{M}(\boldsymbol{q}_{t+1} - \boldsymbol{q}_{t} - h\boldsymbol{v}_{t+1}) = h^{2}\boldsymbol{f}(\boldsymbol{q}_{t+1}) \]

由于弹性力可由应变能函数计算得到 \(\boldsymbol{f}(\boldsymbol{q}) = -\partial E(\boldsymbol{q})/\partial\boldsymbol{q}\)，所以，上述非线性系统方程可转化为一个无约束非线性优化问题：\(\boldsymbol{q}_{t+1} = \text{arg min} \epsilon(\boldsymbol{q})\)，

\[\epsilon(\boldsymbol{q}) = \frac{1}{2h^{2}}(\boldsymbol{q} - \boldsymbol{q}_{t} - h\boldsymbol{v}_{t})^{T}\boldsymbol{M}(\boldsymbol{q} - \boldsymbol{q}_{t} - h\boldsymbol{v}_{t}) + E(\boldsymbol{q}) \]

该非线性优化问题可通过下降方法（Descent Method）求解得到，其主要迭代步骤如*所示。不同方法之间的区别主要在于如何计算搜索方向，通常涉及到梯度 \(\boldsymbol{g}^{(k)} = \nabla\epsilon(\boldsymbol{q}^{(k)})\) 的使用。

（待添加伪代码）

3.1 常用下降方法

下面将介绍几种常用的下降方法，包括梯度下降法、牛顿法、准牛顿法、非线性共轭梯度法等

3.1.1 梯度下降法 / Gradient Descent Method

在梯度下降法中，搜索方向设定为优化函数的梯度，即 \(\Delta \boldsymbol{q}^{(k)} = - \boldsymbol{g}^{(k)}\)。（优化函数在局部沿梯度方向下降最快。）尽管梯度下降法每个迭代步需要的计算资源少，但其收敛速度是线性的（较慢）。其可视为是通过力的形式来更新位置，本质上与显式时间积分方法类似，同样需要较小的仿真时间步长。

3.1.2 牛顿法 / Newton's Method

为了提高收敛速度，Newton's method approximates \(\epsilon(\boldsymbol{q}^{(k)})\) by a quadratic function. 其搜索方向设定为 \(\Delta \boldsymbol{q}^{(k)} = - (\boldsymbol{H}^{(k)})^{-1}\boldsymbol{g}^{(k)}\) ，其中 \(\boldsymbol{H}^{(k)}\) 为优化函数 \(\epsilon(\boldsymbol{q})\) 的 Hessian 矩阵（在 \(\boldsymbol{q}^{(k)}\) 处）。牛顿方法的收敛速度非常快。然而，由于需要求解线性方程组 \(\boldsymbol{H}^{(k)} \Delta \boldsymbol{q}^{(k)} = -\boldsymbol{g}^{(k)}\) ，因此，其计算量非常大/效率低。（此外，计算 Hessian 矩阵也十分耗时。）（在牛顿法中，并不能保证搜索方向是下降的，除非 Hessian 矩阵是正定的。）

3.1.3 准牛顿法 / Quasi-Newton Method

由于求解线性方程组、计算 Hessian 矩阵十分耗时，因此可以考虑采用某种方法估计 Hessian 矩阵及其逆。例如，在准牛顿法（如 BFGS 方法）中，采用梯度向量（previous gradient vector）估计 Hessian 矩阵的逆。此外，为了避免存储 dense inverse matrix，L-BFGS 方法通过 m 个梯度向量来依次更新 Hessian 矩阵的逆。尽管 L-BFGS 方法的收敛速度比牛顿法慢，但由于其具有更好的计算性能，单步迭代所消耗的计算资源更少。但是，该方法仍然难以实现并行计算。

3.1.4 非线性共轭梯度法 / Nonlinear Conjungate Gradient (CG) Method

在非线性共轭梯度法中，将共轭梯度法应用到非线性优化问题中。基于 Fletcher-Reeves 定理，搜索方向可计算为：

\[\Delta \boldsymbol{q}^{(k)} = -\boldsymbol{g}^{(k)} + \frac{z^{k}}{z^{k-1}} \Delta \boldsymbol{q}^{(k-1)}, \quad z^{k} = \boldsymbol{g}^{(k)}\cdot\boldsymbol{g}^{(k)} \]

非线性共轭梯度法与 \(m=1\) 时的 L-BFGS 十分相似。其收敛速度略快于 L-BFGS，其原因可能是在确定步长时采用了 Hessian 矩阵的精确值（而非估计值）。相对于准牛顿法，非线性共轭梯度法对GPU计算更友好一些，但其仍然需要大量的点乘运算（multiple dot product operation），限制了其在 GPU 上的运算性能。

3.2 基于 Jacobi 和 Chebyshev 的预处理下降方法 / Preconditioning Descent Method by Jacobi+Chebyshev

下面将介绍一种新的方法，参见论文 Descent methods for elastic body simulation on the GPU, ACMTransactions on Graphics (TOG), 2016. 主要包括下降方向计算、步长估计、优化加速、初始化等几个方面。

3.2.1 搜索方向 / Descent Direction

受 preconditioned conjugate gradient 方法的启发，为提高收敛速度，通过 preconditioning 将优化问题转化为 conditioned 状态（a well conditioned one）:

\[\bar{\boldsymbol{q}} = \text{arg min}\epsilon(\boldsymbol{P}^{-1/2} \bar{\boldsymbol{q}}), \quad \text{for } \bar{\boldsymbol{q}} = \boldsymbol{P}^{-1/2} \boldsymbol{q} \]

其中，\(\boldsymbol{P}\) 是 preconditioner 矩阵。（从数学的角度看，该方法相当于在非线性共轭梯度法中用 \(\boldsymbol{P}^{-1} \boldsymbol{g}^{(k)}\) 取代 \(\boldsymbol{g}^{(k)}\) ，同时，\(z^{k} = \boldsymbol{g}^{(k)}\cdot \boldsymbol{P}^{-1}\boldsymbol{g}^{(k)}\)）。

在所有 preconditioners 中，作者更倾向于 Jacobi preconditioner ，这是因为其对 GPU 并行计算更加友好。当优化问题是二次型时，其 Jacobi preconditioner 为常数矩阵 \(\boldsymbol{P} = \text{diag}(\boldsymbol{H})\) ，且在各迭代步中，有 \(\boldsymbol{P}(\boldsymbol{q}^{(k)}) = \text{diag}(\boldsymbol{H}^{(k)})\) 。此时，搜索方向为 \(\Delta \boldsymbol{q}^{(k)} = -\text{diag}^{-1}(\boldsymbol{H}^{(k)})\boldsymbol{g}^{(k)}\) 。

3.2.2 收敛及性能 / Convergence and Performance

收敛的判别条件为

\[\epsilon(\boldsymbol{q}^{(k)} + \alpha^{(k)}\Delta\boldsymbol{q}^{(k)}) \le \epsilon(\boldsymbol{q}^{(k)}) + \alpha^{(k)}\Delta\boldsymbol{q}^{(k)} \cdot \boldsymbol{g}^{(k)} + \frac{B}{2}\left|\left|\alpha^{(k)}\Delta\boldsymbol{q}^{(k)}\right|\right|^{2}_{2} \]

3.2.3 步长调节 / Step Length Adjustment

已知搜索方向 \(\Delta\boldsymbol{q}^{(k)}\) 之后，接下来需要确定搜索步长 \(\alpha^{(k)}\) 。作者采用 backtracking line search 方法来逐渐缩小步长，直至其满足 Wolfe condition 。实际操作中，作者将参数设为 0，即 \(c = 0\) ，因此，只需满足 \(\epsilon(\boldsymbol{q}^{(k)} + \alpha^{(k)}\Delta\boldsymbol{q}^{(k)}) \le \epsilon(\boldsymbol{q}^{(k)})\) 即可。

3.2.4 加速方法 / Momentum-based Acceleration

(待补充，主要是调整了某参数 \(\rho\) )

3.2.5 初始化 / Initialization

在迭代的方式求解优化问题时，初始值的设定也十分重要。在实际操作中，作者采用恒定加速度 constant acceleration 的方式，即 \(\boldsymbol{q}_{t+1} \approx \boldsymbol{q}^{(0)} = \boldsymbol{q}_{t} + h\boldsymbol{v}_{t} + \eta h(\boldsymbol{v}_{t} - \boldsymbol{v}_{t-1})\) ，其中， \(\eta = 0.2\) 。