完整推导 XPBD 的约束投影过程

Posted on Edited on Waline:

近日在分析 UE5 物理仿真源码的 XPBD 仿真框架部分，本文将分享作者对 XPBD 的完整推导。

更多相关物理仿真算法的背景知识，包括约束法、PBD，可见之前的一篇文章：基于物理的动画：布料与软体动力学。

概述

继承了 PBD 和 PD 的思想，Extended Position Based Dynamics (XPBD) 是一种较新的约束法，将柔度矩阵 \(\boldsymbol C\) 引入仿真流程，改进了 PBD 中的约束求解步骤，并引入了 PD 的势能的二次形式，使之具有物理意义。

尽管 XPBD 的运行效率不及 PBD，但物理正确性更高，且能稳健地仿真 \(k \to\infty\) 时的弹簧约束。这种约束可用于骨骼动画中关节的连接。UE5 Chaos 在 5.3 版本后加入了基于 XPBD 的布料仿真功能。

XPBD 改进的是 PBD 中的约束投影这一步：

\[ \boldsymbol x^{\text{new}} = \text{Projection}(\boldsymbol x) \]

其他步骤均与 PBD 相同。

推导

XPBD 核心等式

具体地，设系统中的总约束数目为 \(m\)，约束 \(c\) 对应的函数 \(\phi_c(\boldsymbol x)\) 和势能函数 \(E_c(\boldsymbol x)\) 满足

\[ E_c(\boldsymbol x) = \frac{w_c}{2} (\phi_c(\boldsymbol x))^2 \]

那么，系统的总势能为一二次型，将其写成矩阵形式为

\[ E(\boldsymbol x) = \sum_c E_c(\boldsymbol x) = \sum_c \frac{w_c}{2} (\phi_c(\boldsymbol x))^2 = \frac{1}{2} \boldsymbol \phi^T(\boldsymbol x) \boldsymbol C^{-1} \boldsymbol \phi(\boldsymbol x) \]

其中 \(\boldsymbol \phi(\boldsymbol x) \in\mathbb R^{m} = (\phi_{c_1}(\boldsymbol x), \dots, \phi_{c_m}(\boldsymbol x))^T\) 为所有约束对应函数组成的向量，\(\boldsymbol C = \frac{1}{w_c} \boldsymbol I_m \in\mathbb R^{m \times m}\) 为柔度矩阵，其为对角矩阵，且第 \(j\) 个对角元素为第 \(j\) 个约束对应的刚度 \(w_c\) 的倒数。

此时，系统的力为

\[ \boldsymbol f(\boldsymbol x) = - \nabla E(\boldsymbol x) = - \left( \frac{\partial E}{\partial\boldsymbol \phi} \frac{\partial\boldsymbol \phi}{\partial\boldsymbol x} \right)^T \]

记 \(\boldsymbol J(\boldsymbol x) = \frac{\partial\boldsymbol \phi(\boldsymbol x)}{\partial\boldsymbol x} \in\mathbb R^{m \times3n}\) 为 Jacobian 矩阵，那么

\[ \boldsymbol f(\boldsymbol x) = -\boldsymbol J^T \boldsymbol C^{-1} \boldsymbol \phi \]

记

\[ \boldsymbol \lambda = \boldsymbol C^{-1} \boldsymbol \phi(\boldsymbol x) \in \mathbb R^m \]

为拉格朗日乘子，它是一个辅助用的对偶变量。

记

\[ \tilde{\boldsymbol C} = \frac{\boldsymbol C}{\Delta t^2} \]

为时间步长归一化柔度矩阵。

原 XPBD 论文中，用时间步长归一化柔度矩阵重定义了拉格朗日乘子如下：

\[ \tilde{\boldsymbol \lambda} = -\tilde{\boldsymbol C}^{-1} \boldsymbol \phi(\boldsymbol x) = -\boldsymbol \lambda \cdot \Delta t^2. \]

联立以上各式，整理成关于 \(\boldsymbol x^{\text{new}}\) 和 \(\tilde{\boldsymbol \lambda}^{\text{new}}\) 的方程可得

\[ \begin{cases} \boldsymbol {g}(\boldsymbol {x}^{\text{new}}, \tilde{\boldsymbol \lambda}^{\text{new}}) = \boldsymbol {M}({\boldsymbol {x}^{\text{new}} - \boldsymbol {x}}) - \boldsymbol J^T(\boldsymbol {x}^{\text{new}}) \tilde{\boldsymbol \lambda}^{\text{new}} = \boldsymbol {0} \\ \boldsymbol {h}(\boldsymbol {x}^{\text{new}}, \tilde{\boldsymbol \lambda}^{\text{new}}) = \boldsymbol {\phi}(\boldsymbol {x}^{\text{new}}) + \tilde{\boldsymbol C} \tilde{\boldsymbol \lambda}^{\text{new}} = \boldsymbol {0} \end{cases} \]

符号约定

为简便起见，接下来，我们把需要求解的变量 \(\boldsymbol x^{\text{new}}\) 和 \(\tilde{\boldsymbol \lambda}^{\text{new}}\) 分别记作 \(\boldsymbol x\) 和 \(\tilde{\boldsymbol \lambda}\)，PBD 在第一步仿真下预测的初始位移 \(\boldsymbol x\) 则记作 \(\boldsymbol s\)。即

\[ \begin{cases} \boldsymbol {g}(\boldsymbol {x}, \tilde{\boldsymbol \lambda}) = \boldsymbol {M}({\boldsymbol {x} - \boldsymbol {s}}) - \boldsymbol J^T(\boldsymbol {x}) \tilde{\boldsymbol \lambda} = \boldsymbol {0} \\ \boldsymbol {h}(\boldsymbol {x}, \tilde{\boldsymbol \lambda}) = \boldsymbol {\phi}(\boldsymbol {x}) + \tilde{\boldsymbol C} \tilde{\boldsymbol \lambda} = \boldsymbol {0} \end{cases} \]

离散化

我们将对 \(\boldsymbol {g}\) 和 \(\boldsymbol {h}\) 两个函数等式用泰勒公式线性离散化。设迭代求解次数为 \(l\) 次，那么第 \(k\) 次迭代将产生迭代点对 \((\boldsymbol {x}^{(k)}, \tilde{\boldsymbol \lambda}^{(k)})\)。

对 \(\boldsymbol {g}\) 和 \(\boldsymbol {h}\) 在迭代点 \((\boldsymbol {x}^{(k)}, \tilde{\boldsymbol \lambda}^{(k)})\) 处进行一阶泰勒展开：

\[ \boldsymbol {g}(\boldsymbol {x}^{(k)} + \Delta\boldsymbol {x}, \tilde{\boldsymbol \lambda}_i + \Delta\tilde{\boldsymbol \lambda}) \approx \boldsymbol {g}(\boldsymbol {x}^{(k)}, \tilde{\boldsymbol \lambda}^{(k)}) + \frac{\partial \boldsymbol {g}}{\partial \boldsymbol {x}} \Delta\boldsymbol {x} + \frac{\partial \boldsymbol {g}}{\partial \tilde{\boldsymbol \lambda}} \Delta\tilde{\boldsymbol \lambda} = \boldsymbol {0} \]

\[ \boldsymbol {h}(\boldsymbol {x}^{(k)} + \Delta\boldsymbol {x}, \tilde{\boldsymbol \lambda}_i + \Delta\tilde{\boldsymbol \lambda}) \approx \boldsymbol {h}(\boldsymbol {x}^{(k)}, \tilde{\boldsymbol \lambda}^{(k)}) + \frac{\partial \boldsymbol {h}}{\partial \boldsymbol {x}} \Delta\boldsymbol {x} + \frac{\partial \boldsymbol {h}}{\partial \tilde{\boldsymbol \lambda}} \Delta\tilde{\boldsymbol \lambda} = \boldsymbol {0} \]

将泰勒展开后的方程整理为矩阵形式：

\[ \begin{bmatrix} \frac{\partial \boldsymbol {g}}{\partial \boldsymbol {x}} & \frac{\partial \boldsymbol {g}}{\partial \tilde{\boldsymbol \lambda}} \\ \frac{\partial \boldsymbol {h}}{\partial \boldsymbol {x}} & \frac{\partial \boldsymbol {h}}{\partial \tilde{\boldsymbol \lambda}} \end{bmatrix} \begin{bmatrix} \Delta \boldsymbol x \\ \Delta \tilde{\boldsymbol \lambda} \end{bmatrix} = - \begin{bmatrix} \boldsymbol {g}(\boldsymbol {x}^{(k)}, \tilde{\boldsymbol \lambda}^{(k)}) \\ \boldsymbol {h}(\boldsymbol {x}^{(k)}, \tilde{\boldsymbol \lambda}^{(k)}) \end{bmatrix} \]

计算偏导数

\[ \frac{\partial \boldsymbol {g}}{\partial \boldsymbol {x}} = \boldsymbol {M} - \frac{\partial}{\partial \boldsymbol {x}} \left( \boldsymbol {J}^T(\boldsymbol {x}) \tilde{\boldsymbol \lambda}\right) = \boldsymbol {M} - \frac{\partial^2 \boldsymbol {\phi}(\boldsymbol {x})}{\partial \boldsymbol {x}^2} \tilde{\boldsymbol \lambda} \]

其中

\[ \frac{\partial^2 \boldsymbol {\phi}(\boldsymbol {x})}{\partial \boldsymbol {x}^2} \in \mathbb R^{m \times 3n \times 3n} \]

是各约束函数 \(\boldsymbol {\phi}(\boldsymbol {x})\) 的 Hessian 矩阵组成的张量。

\[ \frac{\partial \boldsymbol {g}}{\partial \tilde{\boldsymbol \lambda}} = -\boldsymbol {J}^T(\boldsymbol {x}) \]

\[ \frac{\partial \boldsymbol {h}}{\partial \boldsymbol {x}} = \boldsymbol {J}(\boldsymbol {x}) \]

\[ \frac{\partial \boldsymbol {h}}{\partial \tilde{\boldsymbol \lambda}} = \tilde{\boldsymbol C} \]

构建线性方程组

代入矩阵形式的线性方程组得到

\[ \begin{bmatrix} \boldsymbol {M} - \frac{\partial \boldsymbol {\phi}^2(\boldsymbol {x}^{(k)})}{\partial \boldsymbol {x}} \tilde{\boldsymbol \lambda_i} & -\boldsymbol {J}^T(\boldsymbol {x}^{(k)}) \\ \boldsymbol {J}(\boldsymbol {x}^{(k)}) & \tilde{\boldsymbol C} \end{bmatrix} \begin{bmatrix} \Delta \boldsymbol x \\ \Delta \tilde{\boldsymbol \lambda} \end{bmatrix} = - \begin{bmatrix} \boldsymbol {g}(\boldsymbol {x}^{(k)}, \tilde{\boldsymbol \lambda}^{(k)}) \\ \boldsymbol {h}(\boldsymbol {x}^{(k)}, \tilde{\boldsymbol \lambda}^{(k)}) \end{bmatrix} \]

简化近似

由于直接计算 Hessian 项 \(\boldsymbol {M} - \frac{\partial^2 \boldsymbol {\phi}(\boldsymbol {x})}{\partial \boldsymbol {x}^2} \tilde{\boldsymbol \lambda}\) 较为复杂，XPBD 作者进行了以下近似：

  • 忽略 Hessian 项。这个近似的分析详见下文。这个近似让 XPBD 节省了大量的计算。
  • \(\boldsymbol {g}(\boldsymbol {x}^{(k)}, \tilde{\boldsymbol \lambda}^{(k)}) \approx \boldsymbol 0\)。这个近似对于 \(i=0\) 时显然正确，并且实验表明是不错的近似。这个近似也让 XPBD 节省了大量的计算。

此时线性系统简化为：

\[ \begin{bmatrix} \boldsymbol {M} & -\boldsymbol {J}^T(\boldsymbol {x}^{(k)}) \\ \boldsymbol {J}(\boldsymbol {x}^{(k)}) & \tilde{\boldsymbol C} \end{bmatrix} \begin{bmatrix} \Delta \boldsymbol x \\ \Delta \tilde{\boldsymbol \lambda} \end{bmatrix} = - \begin{bmatrix} 0 \\ \boldsymbol {h}(\boldsymbol {x}^{(k)}, \tilde{\boldsymbol \lambda}^{(k)}) \end{bmatrix} \]

每迭代更新

\[ \begin{cases} \tilde{\boldsymbol \lambda}^{(k+1)} = \tilde{\boldsymbol \lambda}^{(k)} + \Delta \tilde{\boldsymbol \lambda}, \\ \boldsymbol x^{(k+1)} = \boldsymbol x^{(k)} + \Delta \boldsymbol x. \end{cases} \]

迭代完成后，即可从迭代变量 \(\boldsymbol x^{(l)}\) 中得到约束投影结果 \(\boldsymbol x^{\text{new}}\)。由于存在关于 \(\boldsymbol x^{\text{new}}\) 和 \(\tilde{\boldsymbol \lambda}^{\text{new}}\) 的不动点方程，因此迭代可保证收敛。

迭代初值可设为

\[ \boldsymbol x^{(0)} = \boldsymbol s, \tilde{\boldsymbol \lambda}^{(0)} = \boldsymbol 0. \]

计算未知量

为了在每次迭代求解 \(\Delta \boldsymbol x\) 和 \(\Delta \tilde{\boldsymbol \lambda}\)，更新 \(\boldsymbol x\) 和 \(\tilde{\boldsymbol \lambda}\)，可以直接求解整个线性系统，但计算代价太大。论文指出，基于 Schur 补（其实就是二元一次方程），消去 \(\Delta \boldsymbol x\) 可得到一个快速的精确算法。首先，根据以下线性方程组解出 \(\Delta \tilde{\boldsymbol \lambda}\)：

\[ (\boldsymbol J(\boldsymbol x^{(k)}) \boldsymbol M^{-1} \boldsymbol J^T(\boldsymbol x^{(k)}) + \tilde{\boldsymbol C}) \Delta \tilde{\boldsymbol \lambda} = - \boldsymbol {\phi}(\boldsymbol x^{(k)}) - \tilde{\boldsymbol C} \tilde{\boldsymbol \lambda}_i \]

其中系数矩阵 \(\boldsymbol J(\boldsymbol x_i) \boldsymbol M^{-1} \boldsymbol J^T(\boldsymbol x_i) + \tilde{\boldsymbol C}\) 被称为 Schur 补矩阵。

解出 \(\Delta \tilde{\boldsymbol \lambda}\) 后，可直接得到 \(\Delta \boldsymbol x\)：

\[ \Delta \boldsymbol x = \boldsymbol M^{-1} \boldsymbol J^{T}(\boldsymbol x^{(k)}) \Delta \tilde{\boldsymbol \lambda} \]

并行法更新约束

由 Schur 补式子，我们可以发现，对第 \(j\) 个约束，其约束投影过程中，求解 \(\Delta \tilde{\lambda}_j\) 只需

\[ \Delta \tilde{\lambda}_j = -\frac{\phi_j(\boldsymbol x^{(k)}) + \tilde C_j \tilde{\lambda}_{j}^{(k)}}{ \nabla \phi_j(\boldsymbol x^{(k)}) \boldsymbol M^{-1} \nabla \phi_j^T(\boldsymbol x^{(k)}) + \tilde C_j} \]

其中 \(\tilde{\lambda}_{j}^{(k)}\) 为 \(\tilde{\boldsymbol \lambda}^{(k)}\) 的第 \(j\) 个分量。

对第 \(i\) 个顶点，求解 \(\Delta x_i\) 只需

\[ \Delta x_i = \frac{\sum_{j=1}^m \nabla \phi_{j}(\boldsymbol x^{(k)}) \cdot \Delta \tilde{\lambda}_j}{m_i} \]

这意味着，约束投影过程迭代地进行，其中每次迭代中，可以先并行地求解每个 \(\Delta \tilde{\lambda}_j\)，然后并行地求解每个 \(\Delta x_i\)，即每顶点的约束投影并行执行。这与 PBD 中约束投影的 Jacobian style 匹配。

算法：

  • 迭代若干次：
    • 对 \(j \in [1, m]\) 并行：
      • 求解 \(\Delta \tilde{\lambda}_j\)
    • 对 \(i \in [1, n]\) 并行：
      • 求解 \(\Delta x_i\)
    • 更新 \(\tilde{\boldsymbol \lambda}^{(k+1)}\) 和 \(\boldsymbol x^{(k+1)}\)。

使用该方法在 GPU 上实现并行更新 \(\Delta x_i\) 时，可以考虑用 atomicAdd 让每个与顶点 \(i\) 有关的约束线程自主累加到 \(\Delta x_i\) 上去。

串行法更新约束

算法：

  • 迭代若干次：
    • 对 \(j \in [1, m]\)：
      • 求解 \(\Delta \tilde{\lambda}_j\)
      • 立即更新与约束 \(j\) 有关的顶点的 \(\Delta x_i\)，即 \(\Delta x_i = \frac{ \nabla \phi_{j}(\boldsymbol x^{(k)}) \cdot \Delta \tilde{\lambda}_j}{m_i}\)
      • 更新 \(\tilde{\boldsymbol \lambda}^{(k+1)}\) 和 \(\boldsymbol x^{(k+1)}\)。

该方法直接对应原始 PBD 的 Gauss-Seidel style，一般情况下不方便并行，但易于实现。

但对于在发丝这样的特殊几何中，求解涉及顶点很规则的约束，就很适合在串行法的框架中实现并行求解。UE5 Groom 采用一种被我称之为“红黑边”迭代的方式求解在发丝顶点之间的约束。

XPBD 对势能函数的 Hessian 的近似

势能项

\[ E(\boldsymbol x) = \sum_c E_c(\boldsymbol x) = \sum_c \frac{w_c}{2} (\phi_c(\boldsymbol x))^2 = \frac{1}{2} \boldsymbol \phi^T(\boldsymbol x) \boldsymbol C^{-1} \boldsymbol \phi(\boldsymbol x) \]

的 Hessian 是

\[ \boldsymbol H(\boldsymbol x) = \frac{\partial (\boldsymbol J^T(\boldsymbol x) \boldsymbol C^{-1} \boldsymbol \phi(\boldsymbol x))}{\partial \boldsymbol x} = \frac{\partial^2 \boldsymbol {\phi}(\boldsymbol {x})}{\partial \boldsymbol {x}^2} \boldsymbol C^{-1} \boldsymbol \phi(\boldsymbol x) + \boldsymbol J^T(\boldsymbol x) \boldsymbol C^{-1} \boldsymbol J(\boldsymbol x) = \frac{\partial^2 \boldsymbol {\phi}(\boldsymbol {x})}{\partial \boldsymbol {x}^2} \boldsymbol \lambda + \boldsymbol J^T(\boldsymbol x) \boldsymbol C^{-1} \boldsymbol J(\boldsymbol x). \]

可见，其中的第一项与线性方程组中的 Hessian 项 \(\frac{\partial^2 \boldsymbol {\phi}(\boldsymbol {x})}{\partial \boldsymbol {x}^2} \boldsymbol \lambda\) 有关系。XPBD 采用的近似相当于只取第一项（约束力），忽略第二项（几何刚度）。