论文阅读笔记：ADMM - Projective Dynamics: Fast Simulation of Hyperelastic Models with Dynamic Constraints

材料来源于 ADMM ⊇ Projective Dynamics: Fast Simulation of Hyperelastic Models with Dynamic Constraints, IEEE Transactions on Visualization and Computer Graphics, 2017.

---

1. 形变仿真基础

1.1 优化隐式积分方法

设形变模型网格节点位置 $\mathit{x}\in {\mathbb{R}}^{3n}$ ，速度 $\mathit{v}\in {\mathbb{R}}^{3n}$ ，以及质量 $\mathit{M}\in {\mathbb{R}}^{3n}$。采用隐式时间积分，得到系统方程为

$$\mathit{x}(t+\mathrm{\Delta }t)=\mathit{x}(t)+\mathit{v}(t+\mathrm{\Delta }t)\mathrm{\Delta }t$$

$$\mathit{M}\mathit{v}(t+\mathrm{\Delta }t)=\mathit{M}\mathit{v}(t)+{\mathit{f}}_{ext}(\mathrm{\Delta }t)+\mathit{f}(t+\mathrm{\Delta }t)\mathrm{\Delta }t$$

进一步可记作

$$\frac{1}{\mathrm{\Delta }{t}^2}\mathit{M}(\mathit{x}(t+\mathrm{\Delta }t)-\tilde{\mathit{x}}(t+\mathrm{\Delta }t))=\mathit{f}(t+\mathrm{\Delta }t)$$

其中，$\tilde{\mathit{x}}(t+\mathrm{\Delta }t)=\mathit{x}(t)+\mathit{v}(t)\mathrm{\Delta }t+{\mathit{M}}^{-1}{\mathit{f}}_{ext}(t)\mathrm{\Delta }{t}^2$ ，其为不考虑隐式力（implicit forces $\mathit{f}$）情况下的形变位置预测。

通常情况下，由于弹性力 $\mathit{f}$ 是位置 $\mathit{x}$ 及速度 $\mathit{v}$ 的非线性函数，上述方程是一个大型、非线性系统方程，其求解计算量非常大。在现有的仿真方法中，通常将上述非线性系统线性化后直接求解（linearize the system about the current state, applying one iteration of Newton's method）。但该方法仍有许多局限。

另一种方法是将上述系统转化为优化问题。当 $\mathit{f}$ 是保守力（conservative）时，其可以记作势能的导数，即 $\mathit{f}=-\mathrm{\nabla }U(\mathit{x})$，则上述系统方程可以转化为如下最优化问题

$$\mathit{x}(t+\mathrm{\Delta }t)=\underset{\mathit{x}}{\arg min}(\frac{1}{2\mathrm{\Delta }{t}^2}||\mathit{x}-\tilde{\mathit{x}}(t+\mathrm{\Delta }t)|{|}_{\mathit{M}}^2+U(\mathit{x}))$$

其中 $||\mathit{x}|{|}_{\mathit{M}}=\sqrt{{\mathit{x}}^T\mathit{M}\mathit{x}}$，为了简化描述，可将 $\mathit{x}(t+\mathrm{\Delta }t)$ 和 $\tilde{\mathit{x}}(t+\mathrm{\Delta }t)$ 分别记作 $\mathit{x}$ 和 $\tilde{\mathit{x}}$。

1.2 形变能函数

通常情况下，形变能函数 $U$ 是多个能量子项的和，每个能量子项仅与少量节点有关。在这里，为每个能量子项定义一个“降阶矩阵”（reduction matrix） ${\mathit{D}}_i$，则每个能量子项所对应 ${\mathit{D}}_i\mathit{x}$。那么，形变能函数可以记作：

$$U(\mathit{x})=\sum _{i=1}^m{U}_i({\mathit{D}}_i\mathit{x})$$

其中，$U_i$ maps the local coordinates to the corresponding energy. 为了简便起见，将上述公式记作如下形式：

$$U_{\ast }\left(\left[\begin{array}{c}{\mathit{D}}_1\mathit{x}\\ {\mathit{D}}_2\mathit{x}\\ ⋮\\ {\mathit{D}}_m\mathit{x}\end{array}\right]\right)=\sum _{i=1}^m{U}_i({\mathit{D}}_i\mathit{x})$$

因此，即有 $U(\mathit{x})=U_{\ast }(\mathit{D}\mathit{x})$ ，其中 $\mathit{D}=[{\mathit{D}}_1^T{\mathit{D}}_2^T\cdots {\mathit{D}}_n^T{]}^T$ ，那么，上述最优化问题可记作

$$\underset{\mathit{x}}{min}\frac{1}{2\mathrm{\Delta }{t}^2}||\mathit{x}-\tilde{\mathit{x}}|{|}_{\mathit{M}}^2+U_{\ast }(\mathit{D}\mathit{x})$$

举例，对于质点弹簧模型，就某根弹簧而言，其所对应的节点为 ${\mathit{x}}_a$ 和 ${\mathit{x}}_b$，弹簧的长度和刚度分别为 $l$ 和 $k$，则选择 ${\mathit{D}}_i$ 使得 ${\mathit{D}}_i\mathit{x}=[{\mathit{x}}_b-{\mathit{x}}_a]$ ，同时，定义 $U_i({\mathit{D}}_i\mathit{x})=\frac{1}{2}k(||{\mathit{D}}_i\mathit{x}-l||{)}^2$ 。这里的 reduction matrix ${\mathit{D}}_i$ 更像是一个选择矩阵，把与能量子项相关的网格节点位置从整体中选择出来，同时将相关的网格节点位置组合计算后得到一个统一变量。（照这样讲，对于四面体单元，${\mathit{D}}_i\mathit{x}$ 可以是相应四面体单元的形变梯度 ${\mathit{F}}_i$）。

2. 基于 ADMM 的形变仿真方法

2.1 ADMM 基础

ADMM (alternating direction method of multipliers) 是一种求解如下最优化问题的方法，最优化问题通常有如下形式

$$\underset{\mathit{x},\mathit{z}}{min}f(\mathit{x})+g(\mathit{z})$$

$$\text{s.t.}\mathit{A}\mathit{x}+\mathit{B}\mathit{z}=\mathit{c}$$

通过引入对偶变量（dual variable） $\mathit{u}$ ，其求解过程如下：

$${\mathit{x}}^{n+1}=\underset{\mathit{x}}{\arg min}(f(\mathit{x})+\frac{\rho }{2}||\mathit{A}\mathit{x}+\mathit{B}{\mathit{z}}^n-\mathit{c}+{\mathit{u}}^n|{|}^2)$$

$${\mathit{z}}^{n+1}=\underset{\mathit{z}}{\arg min}(g(\mathit{z})+\frac{\rho }{2}||\mathit{A}{\mathit{x}}^{n+1}+\mathit{B}\mathit{z}-\mathit{c}+{\mathit{u}}^n|{|}^2)$$

$${\mathit{u}}^{n+1}={\mathit{u}}^n+(\mathit{A}{\mathit{x}}^{n+1}+\mathit{B}{\mathit{z}}^{n+1}-\mathit{c})$$

其中，$\rho$ 是迭代步长，上标 $n$ 是迭代步数。

注1：当 $f$ 和 $g$ 都是凸函数时，ADMM 算法可以保证任意步长 $\rho$ 都能收敛到最优值。特别注意，$f$ 和 $g$ 必须是可微的（differentiable）。但是，实际应用中发现，尽管形变模型中 $f$ 和 $g$ 不可微，ADMM 方法也有很好的表现。

注2：ADMM 有两个特点值得注意。（1）对变量 $\mathit{x}$ 和 $\mathit{z}$ 的缩放（rescaling）并不会影响迭代收敛速度。（we get exactly the same sequence of iterates as we would for the original problem.）；（2）对约束的缩放（rescaling）会极大的影响收敛速率。即，如果将原始的约束最优化问题转化为如下形式：

$$\underset{\mathit{x},\mathit{z}}{min}f(\mathit{x})+g(\mathit{z})$$

$$\text{s.t.}\mathit{W}\mathit{A}\mathit{x}+\mathit{W}\mathit{B}\mathit{z}=\mathit{W}\mathit{c}$$

通过引入可逆矩阵 $\mathit{W}$ 可以显著的影响 ADMM 算法的收敛速率。

附录 A 中证明，当 $f$ 和 $g$ 是二次凸函数，$\mathit{B}$ 是可逆矩阵时，选择 $\rho =1$ 以及 $\mathit{W}$ 使得 $g(\mathit{z})=\frac{1}{2}||\mathit{W}\mathit{B}\mathit{z}|{|}^2$，则 ADMM 算法每次迭代都使得误差减半（即收敛速度非常快）。已知形变仿真中，$f$ 为二次凸函数，那么，我们就需要恰当的选择 $\mathit{W}$，从而使 $g$ 也是二次凸函数。？？没有特别懂这儿？？

2.2 基于 ADMM 求解优化隐式积分下的形变

（1） 采用 ADMM 求解如 1.2 中所述的最优化问题，则上述最优化问题可以先写作如下形式

$$\underset{\mathit{x}}{min}\frac{1}{2\mathrm{\Delta }{t}^2}||\mathit{x}-\tilde{\mathit{x}}|{|}_{\mathit{M}}^2+U_{\ast }(\mathit{z})$$

$$\text{s.t.}\mathit{W}(\mathit{D}\mathit{x}-\mathit{z})=0$$

注：感觉这里就是将 $\mathit{D}\mathit{x}$ 替换为 $\mathit{z}$，同时，通过添加约束，使得两者相等。

那么，对照 ADMM 方法的最优化函数形式，有

$$f(\mathit{x})=\frac{1}{2\mathrm{\Delta }{t}^2}||\mathit{x}-\tilde{\mathit{x}}|{|}_{\mathit{M}}^2$$

$$g\mathit{z}=U_{\ast }(\mathit{z})$$

$$\mathit{A}=\mathit{W}\mathit{D}$$

$$\mathit{B}=-\mathit{W}$$

$$\mathit{c}=\mathbf{0}$$

（2） 下面，进行详细的分析。如 1.2 所述，物体的形变能被分为多个能量子项 $U_i({\mathit{D}}_i\mathit{x})$ ，因此，在这里也将 $\mathit{z}$ 分成相应的子项 ${\mathit{z}}_i$，并将权重矩阵 $\mathit{W}$ 记作对角块矩阵（block diagonal with blocks），即

$$\mathit{z}=\left[\begin{array}{c}{\mathit{z}}_1\\ {\mathit{z}}_2\\ ⋮\\ {\mathit{z}}_m\end{array}\right],\mathit{W}=\left[\begin{array}{cccc}{\mathit{W}}_1& & & \\ & {\mathit{W}}_2& & \\ & & \ddots & \\ & & & {\mathit{W}}_m\end{array}\right]$$

也就是讲，$\mathit{W}(\mathit{D}\mathit{x}-\mathit{z})=0$ 可以记作

$$\left[\begin{array}{cccc}{\mathit{W}}_1& & & \\ & {\mathit{W}}_2& & \\ & & \ddots & \\ & & & {\mathit{W}}_m\end{array}\right]\left[\begin{array}{c}{\mathit{D}}_1\mathit{x}-{\mathit{z}}_1\\ {\mathit{D}}_2\mathit{x}-{\mathit{z}}_2\\ ⋮\\ {\mathit{D}}_m\mathit{x}-{\mathit{z}}_m\end{array}\right]=\left[\begin{array}{c}\mathbf{0}\\ \mathbf{0}\\ ⋮\\ \mathbf{0}\end{array}\right]$$

简记为 ${\mathit{W}}_i({\mathit{D}}_i\mathit{x}-{\mathit{z}}_i)=\mathbf{0}\text{for all}i=1,\cdots ,m.$ 。通常，${\mathit{W}}_i={\omega }_i\mathit{I}$ 。

（3） 如 2.1 所述，设 $\rho =1$，并选择 ${\mathit{W}}_i$ 使得 $U_i$ 近似于二次函数，即选择可逆矩阵 ${\mathit{W}}_i$ 使得 $U_i({\mathit{z}}_i)$ 约等于 $\frac{1}{2}||{\mathit{W}}_i{\mathit{z}}_i|{|}^2$ (原文：we fix $\rho =1$, and choose ${\mathit{W}}_i$ for each energy tem $U_i$ to approximate it by a quadratic function. That is, we seek an invertible matrix ${\mathit{W}}_i$ such that $U_i({\mathit{z}}_i)\approx \frac{1}{2}||{\mathit{W}}_i{\mathit{z}}_i|{|}^2$ over the likely range of values of ${\mathit{z}}_i$.)

具体而言，与论文 Liu et. al 相似，采用 $U_i$ 的 Hessian 矩阵 （using the Hessian of $U_i$ at a single point）。例如，对于弹簧质点模型，弹簧势能为 $U({\mathit{z}}_i)=\frac{1}{2}{k}_i{(||{\mathit{z}}_i||-l_i)}^2$ ，Hessian 矩阵在 ${\mathit{z}}_i=0$ 处无定义，且在 $||{\mathit{z}}_i||=l_i$ 处秩为1（has rank 1），因此，采用 $\frac{1}{2}{k}_i||{\mathit{z}}_i|{|}^2$ 近似形变能 $U_i({\mathit{z}}_i)$ （we may take $U_i({\mathit{z}}_i)\approx \frac{1}{2}{k}_i||{\mathit{z}}_i|{|}^2$ as a reasonable quadratic approximation of the global behavior of the spring energy.）

2.3 求解计算流程

引入 dual variable $\mathit{u}$ 之后，将 ${\mathit{W}}^{-1}\mathit{u}$ 记作 $\overline{\mathit{u}}$，ADMM 算法的计算流程如下：

$$\begin{array}{rl}{\mathit{x}}^{n+1}=& \underset{\mathit{x}}{\arg min}(\frac{1}{2\mathrm{\Delta }{t}^2}||\mathit{x}-\tilde{\mathit{x}}|{|}_{\mathit{M}}^2+\frac{1}{2}||\mathit{W}(\mathit{D}\mathit{x}-{\mathit{z}}^n+{\overline{\mathit{u}}}^n)|{|}^2)\\ =& {(\mathit{M}+\mathrm{\Delta }{t}^2{\mathit{D}}^T{\mathit{W}}^T\mathit{W}\mathit{D})}^{-1}(\mathit{M}\tilde{\mathit{x}}+\mathrm{\Delta }{t}^2{\mathit{D}}^T{\mathit{W}}^T\mathit{W}({\mathit{z}}^n-{\overline{\mathit{u}}}^n))\end{array}$$

注：感觉在这里，是对最优化做了一次求解，导数等于零。

$${\mathit{z}}^{n+1}=\underset{\mathit{z}}{\arg min}(U_{\ast }(\mathit{z})+\frac{1}{2}||\mathit{W}(\mathit{D}{\mathit{x}}^{n+1}-\mathit{z}+{\overline{\mathit{u}}}^n)|{|}^2)$$

注：这里便是优化求解。但为什么不和上面一样，做一次求解呢？

$${\overline{\mathit{u}}}^{n+1}={\overline{\mathit{u}}}^n+\mathit{D}{\mathit{x}}^{n+1}-{\mathit{z}}^{n+1}$$

在迭代求解的过程中，可采用如下方式计算残差，当残差小于给定阈值时，停止迭代。The primal and dual residuals are

$$\begin{array}{rl}{\mathit{r}}^{n+1}& =\mathit{W}(\mathit{D}{\mathit{x}}^{n+1}-{\mathit{z}}^{n+1})\\ {\mathit{s}}^{n+1}& ={\mathit{D}}^T{\mathit{W}}^T\mathit{W}({\mathit{z}}^{n+1}-{\mathit{z}}^n)\end{array}$$

3. ADMM 算法的物理解释

（大致看了一下这段，感觉给出的解释，更像是从广义力，或者说是变换后的空间，解释了优化的流程。并不是从笛卡尔坐标系中的 $\mathit{x}$ 的角度出发，而是变换空间后，从 $\mathit{D}\mathit{x}$ 的视角出发的。感觉好像还挺像那么回事的。后面再具体看吧。）