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

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

---

0. 概述

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

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

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

$$\mathit{B}{\mathit{x}}^{(k+1)}=\mathit{C}{\mathit{x}}^{(k)}+\mathit{b}(\text{即}{\mathit{x}}^{(k+1)}={\mathit{B}}^{-1}(\mathit{C}{\mathit{x}}^{(k)}+\mathit{b}))$$

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

$$\begin{array}{rl}{\mathit{e}}^{(k+1)}& ={\mathit{x}}^{(k+1)}-\mathit{x}\\ & ={\mathit{B}}^{-1}(\mathit{C}{\mathit{x}}^{(k)}+\mathit{b})-\mathit{x}\\ & ={\mathit{B}}^{-1}\mathit{C}{\mathit{x}}^{(k)}+{\mathit{B}}^{-1}\mathit{b}-\mathit{x}\end{array}$$

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

$${\mathit{e}}^{(k+1)}={\mathit{B}}^{-1}\mathit{C}{\mathit{x}}^{(k)}-{\mathit{B}}^{-1}\mathit{C}\mathit{x}={\mathit{B}}^{-1}\mathit{C}{\mathit{e}}^{(k)}$$

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

$${\mathit{e}}^{(k)}=({\mathit{B}}^{-1}\mathit{C}{)}^k{\mathit{e}}^{(0)}=\mathit{Q}{\mathbf{\Lambda }}^k{\mathit{Q}}^{-1}{\mathit{e}}^{(0)}$$

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

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

1.1 Jacobi 迭代求解器

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

$${\mathit{x}}^{(k+1)}={\mathit{D}}^{-1}(\mathit{U}+\mathit{L}){\mathit{x}}^{(k)}+{\mathit{D}}^{-1}\mathit{b}$$

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

$$x_i^{(k+1)}=\frac{1}{d_i}(b_i-\sum _{j=1}^{i-1}{l}_{ij}{x}_j^{(k)}-\sum _{j=i+1}^N{u}_{ij}{x}_j^{(k)})$$

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

1.2 Gauss-Seidel 迭代求解器

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

$${\mathit{x}}^{(k+1)}=（\mathit{D}-\mathit{L}{)}^{-1}\mathit{U}{\mathit{x}}^{(k)}+(\mathit{D}-\mathit{L}{)}^{-1}\mathit{b}$$

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

$$(\mathit{D}-\mathit{L}){\mathit{x}}^{(k+1)}=\mathit{U}{\mathit{x}}^{(k)}+\mathit{b}$$

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

$$\mathit{D}{\mathit{x}}^{(k+1)}=\mathit{L}{\mathit{x}}^{(k+1)}+\mathit{U}{\mathit{x}}^{(k)}+\mathit{b}\text{即}{\mathit{x}}^{(k+1)}={\mathit{D}}^{-1}\mathit{L}{\mathit{x}}^{(k+1)}+{\mathit{D}}^{-1}\mathit{U}{\mathit{x}}^{(k)}+{\mathit{D}}^{-1}\mathit{b}$$

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

$$x_i^{(k+1)}=\frac{1}{d_i}(b_i-\sum _{j=1}^{i-1}{l}_{ij}{x}_j^{(k+1)}-\sum _{j=i+1}^N{u}_{ij}{x}_j^{(k)})$$

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

2. Projective and Position-based Dynamics

（待整理）

3. Parallel Descent Method

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

3.0 问题描述

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

$${\mathit{q}}_{t+1}={\mathit{q}}_t+h{\mathit{v}}_{t+1},{\mathit{v}}_{t+1}={\mathit{v}}_t+h{\mathit{M}}^{-1}\mathit{f}({\mathit{q}}_{t+1})$$

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

$$\mathit{M}({\mathit{q}}_{t+1}-{\mathit{q}}_t-h{\mathit{v}}_{t+1})=h^2\mathit{f}({\mathit{q}}_{t+1})$$

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

$$ϵ(\mathit{q})=\frac{1}{2h^2}(\mathit{q}-{\mathit{q}}_t-h{\mathit{v}}_t{)}^T\mathit{M}(\mathit{q}-{\mathit{q}}_t-h{\mathit{v}}_t)+E(\mathit{q})$$

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

（待添加伪代码）

3.1 常用下降方法

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

3.1.1 梯度下降法 / Gradient Descent Method

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

3.1.2 牛顿法 / Newton's Method

为了提高收敛速度，Newton's method approximates $ϵ({\mathit{q}}^{(k)})$ by a quadratic function. 其搜索方向设定为 $\mathrm{\Delta }{\mathit{q}}^{(k)}=-({\mathit{H}}^{(k)}{)}^{-1}{\mathit{g}}^{(k)}$ ，其中 ${\mathit{H}}^{(k)}$ 为优化函数 $ϵ(\mathit{q})$ 的 Hessian 矩阵（在 ${\mathit{q}}^{(k)}$ 处）。牛顿方法的收敛速度非常快。然而，由于需要求解线性方程组 ${\mathit{H}}^{(k)}\mathrm{\Delta }{\mathit{q}}^{(k)}=-{\mathit{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 定理，搜索方向可计算为：

$$\mathrm{\Delta }{\mathit{q}}^{(k)}=-{\mathit{g}}^{(k)}+\frac{z^k}{z^{k-1}}\mathrm{\Delta }{\mathit{q}}^{(k-1)},z^k={\mathit{g}}^{(k)}\cdot {\mathit{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）:

$$\overline{\mathit{q}}=\text{arg min}ϵ({\mathit{P}}^{-1/2}\overline{\mathit{q}}),\text{for }\overline{\mathit{q}}={\mathit{P}}^{-1/2}\mathit{q}$$

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

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

3.2.2 收敛及性能 / Convergence and Performance

收敛的判别条件为

$$ϵ({\mathit{q}}^{(k)}+{\alpha }^{(k)}\mathrm{\Delta }{\mathit{q}}^{(k)})\le ϵ({\mathit{q}}^{(k)})+{\alpha }^{(k)}\mathrm{\Delta }{\mathit{q}}^{(k)}\cdot {\mathit{g}}^{(k)}+\frac{B}{2}{\left|\left|{\alpha }^{(k)}\mathrm{\Delta }{\mathit{q}}^{(k)}\right|\right|}_2^2$$

3.2.3 步长调节 / Step Length Adjustment

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

3.2.4 加速方法 / Momentum-based Acceleration

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

3.2.5 初始化 / Initialization

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

3.2.6 一些其他实现细节

例如，每八次迭代，更新一次 ${\mathit{H}}^{(k)}$ ，等等。

4. 结语

上述文章主要介绍了形变仿真计算中通过迭代法求解优化问题。特别是给出了一种适合于 GPU 并行计算的下降方法，可实现高效、准确的形变仿真计算。相关源代码及实现细节，可参照随笔 Source Code and Details for the Descent Method