论文阅读笔记：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 优化隐式积分方法

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

\[\boldsymbol{x}(t + \Delta t) = \boldsymbol{x}(t) + \boldsymbol{v}(t + \Delta t) \Delta t \]

\[\boldsymbol{M} \boldsymbol{v}(t + \Delta t) = \boldsymbol{M} \boldsymbol{v}(t) + \boldsymbol{f}_{ext}(\Delta t) + \boldsymbol{f}(t + \Delta t) \Delta t \]

进一步可记作

\[\frac{1}{\Delta t^2}\boldsymbol{M}(\boldsymbol{x}(t + \Delta t) - \tilde{\boldsymbol{x}}(t + \Delta t)) = \boldsymbol{f}(t + \Delta t) \]

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

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

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

\[\boldsymbol{x}(t + \Delta t) = \underset{\boldsymbol{x}}{\arg \min}\left(\frac{1}{2 \Delta t^2} ||\boldsymbol{x} - \tilde{\boldsymbol{x}}(t + \Delta t)||^{2}_{\boldsymbol{M}} + U(\boldsymbol{x})\right) \]

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

1.2 形变能函数

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

\[U(\boldsymbol{x}) = \sum_{i=1}^{m} U_{i}(\boldsymbol{D}_{i} \boldsymbol{x}) \]

其中，\(U_{i}\) maps the local coordinates to the corresponding energy. 为了简便起见，将上述公式记作如下形式：

\[U_{*}\left(\begin{bmatrix} \boldsymbol{D}_{1}\boldsymbol{x}\\ \boldsymbol{D}_{2}\boldsymbol{x}\\ \vdots \\ \boldsymbol{D}_{m}\boldsymbol{x} \end{bmatrix}\right) = \sum_{i=1}^{m} U_{i}(\boldsymbol{D}_{i} \boldsymbol{x})\]

因此，即有 \(U(\boldsymbol{x}) = U_{*}(\boldsymbol{D} \boldsymbol{x})\) ，其中 \(\boldsymbol{D} = [\boldsymbol{D}_{1}^{T} \quad \boldsymbol{D}_{2}^{T} \quad \cdots \quad \boldsymbol{D}_{n}^{T}]^{T}\) ，那么，上述最优化问题可记作

\[\underset{\boldsymbol{x}}{\min} \frac{1}{2\Delta t^2}||\boldsymbol{x} - \tilde{\boldsymbol{x}}||^2_{\boldsymbol{M}} + U_{*}(\boldsymbol{D} \boldsymbol{x}) \]

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

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

2.1 ADMM 基础

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

\[\underset{\boldsymbol{x}, \boldsymbol{z}}{\min} f(\boldsymbol{x}) + g(\boldsymbol{z}) \]

\[\text{s.t.} \boldsymbol{Ax} + \boldsymbol{Bz} = \boldsymbol{c} \]

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

\[\boldsymbol{x}^{n+1} = \underset{\boldsymbol{x}}{\arg \min}\left(f(\boldsymbol{x}) + \frac{\rho}{2}||\boldsymbol{Ax} + \boldsymbol{Bz}^{n} - \boldsymbol{c} + \boldsymbol{u}^{n}||^2\right) \]

\[\boldsymbol{z}^{n+1} = \underset{\boldsymbol{z}}{\arg \min}\left(g(\boldsymbol{z}) + \frac{\rho}{2}||\boldsymbol{Ax}^{n+1} + \boldsymbol{Bz} - \boldsymbol{c} + \boldsymbol{u}^{n}||^2\right) \]

\[\boldsymbol{u}^{n+1} = \boldsymbol{u}^{n} + (\boldsymbol{Ax}^{n+1} + \boldsymbol{Bz}^{n+1} - \boldsymbol{c}) \]

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

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

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

\[\underset{\boldsymbol{x}, \boldsymbol{z}}{\min} f(\boldsymbol{x}) + g(\boldsymbol{z}) \]

\[\text{s.t.} \boldsymbol{WAx} + \boldsymbol{WBz} = \boldsymbol{Wc} \]

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

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

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

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

\[\underset{\boldsymbol{x}}{\min} \frac{1}{2\Delta t^2}||\boldsymbol{x} - \tilde{\boldsymbol{x}}||^2_{\boldsymbol{M}} + U_{*}(\boldsymbol{z}) \]

\[\text{s.t.} \boldsymbol{W}(\boldsymbol{D} \boldsymbol{x} - \boldsymbol{z}) = 0 \]

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

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

\[f(\boldsymbol{x}) = \frac{1}{2\Delta t^2}||\boldsymbol{x} - \tilde{\boldsymbol{x}}||^2_{\boldsymbol{M}} \]

\[g\boldsymbol{z} = U_{*}(\boldsymbol{z}) \]

\[\boldsymbol{A} = \boldsymbol{WD} \]

\[\boldsymbol{B} = -\boldsymbol{W} \]

\[\boldsymbol{c} = \boldsymbol{0} \]

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

\[\boldsymbol{z} = \begin{bmatrix} \boldsymbol{z}_{1} \\ \boldsymbol{z}_{2} \\ \vdots \\ \boldsymbol{z}_{m} \end{bmatrix}, \qquad \boldsymbol{W} = \begin{bmatrix} \boldsymbol{W}_{1} & & & \\ &\boldsymbol{W}_{2} & & \\ & &\ddots &\\ & & &\boldsymbol{W}_{m} \end{bmatrix}\]

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

\[\begin{bmatrix} \boldsymbol{W}_{1} & & & \\ &\boldsymbol{W}_{2} & & \\ & &\ddots &\\ & & &\boldsymbol{W}_{m} \end{bmatrix} \begin{bmatrix} \boldsymbol{D}_{1} \boldsymbol{x} - \boldsymbol{z}_{1} \\ \boldsymbol{D}_{2} \boldsymbol{x} - \boldsymbol{z}_{2} \\ \vdots \\ \boldsymbol{D}_{m} \boldsymbol{x} - \boldsymbol{z}_{m} \end{bmatrix} = \begin{bmatrix} \boldsymbol{0} \\ \boldsymbol{0} \\ \vdots \\ \boldsymbol{0} \end{bmatrix}\]

简记为 \(\boldsymbol{W}_{i}(\boldsymbol{D}_{i} \boldsymbol{x} - \boldsymbol{z}_{i}) = \boldsymbol{0} \quad \text{for all} \quad i=1, \cdots, m.\) 。通常，\(\boldsymbol{W}_{i} = \omega_{i}\boldsymbol{I}\) 。

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

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

2.3 求解计算流程

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

\[\begin{aligned}\boldsymbol{x}^{n+1} = &\underset{\boldsymbol{x}}{\arg \min}\left(\frac{1}{2\Delta t^2}||\boldsymbol{x} - \tilde{\boldsymbol{x}}||^2_{\boldsymbol{M}} + \frac{1}{2}||\boldsymbol{W}(\boldsymbol{Dx} - \boldsymbol{z}^{n} + \bar{\boldsymbol{u}}^{n})||^2\right)\\ = &\left(\boldsymbol{M} + \Delta t^2 \boldsymbol{D}^T \boldsymbol{W}^T \boldsymbol{WD}\right)^{-1} \left(\boldsymbol{M}\tilde{\boldsymbol{x}} + \Delta t^2 \boldsymbol{D}^T \boldsymbol{W}^T \boldsymbol{W}(\boldsymbol{z}^n - \bar{\boldsymbol{u}}^n)\right)\end{aligned}\]

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

\[\boldsymbol{z}^{n+1} = \underset{\boldsymbol{z}}{\arg \min}\left( U_{*}(\boldsymbol{z}) + \frac{1}{2} ||\boldsymbol{W}(\boldsymbol{Dx}^{n+1} - \boldsymbol{z} + \bar{\boldsymbol{u}}^{n})||^2 \right) \]

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

\[\bar{\boldsymbol{u}}^{n+1} = \bar{\boldsymbol{u}}^{n} + \boldsymbol{Dx}^{n+1} - \boldsymbol{z}^{n+1} \]

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

\[\begin{aligned}\boldsymbol{r}^{n+1} &= \boldsymbol{W}(\boldsymbol{Dx}^{n+1} - \boldsymbol{z}^{n+1}) \\ \boldsymbol{s}^{n+1} &= \boldsymbol{D}^T\boldsymbol{W}^T\boldsymbol{W}(\boldsymbol{z}^{n+1} - \boldsymbol{z}^n)\end{aligned}\]

3. ADMM 算法的物理解释

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