#! https://zhuanlan.zhihu.com/p/510654097

论文阅读：XPBD

2022.5.11修订：添加了对引入λ的解释。

本文为对XPBD: Position-Based Simulation of Compliant Constrained Dynamics的阅读笔记。本文并非翻译性质或公式推导的文章，内容与文中会有很大的出入，不过方法还是论文里的那套方法，只是会加入大量我的个人理解，同时还会删去我不感兴趣或不懂的部分。

需要的预备知识：对PBD多少有点了解，熟悉多元微分学，知道牛顿运动方程。

有不正之处欢迎指出。

一、一些约定

有 $n$ 个粒子，粒子位置为 $x_i$ （下文中将忽略 $x_i$ 的维度，看着就像是在处理标量一样，但其实都可以直接迁移到向量上的），设 $x=[x_1, x_2, \cdots, x_n]^T$ 。

每个粒子的质量为 $m_i$ ，质量矩阵：

M=\begin{bmatrix} M_1&&&& \\ &M_2&&& \\ &&M_3&& \\ &&&⋱& \\ &&&&M_n \\ \end{bmatrix}, M_i = \begin{bmatrix} m_i & & \\ & m_i & \\ & & m_i \\ \end{bmatrix}

时步为 $\Delta t$ 。

二、总体

将物体视为只受两类力作用：

势能力：只与粒子位置有关的势能力
- 例如弹性势能的弹力
外力：粒子之间相互独立无关的力
- 例如在任何位置都大小方向不变的重力、只与粒子速度有关的空气阻力
- 通常而言，如果某种力作用到一个粒子上会影响到其他粒子上同种力的大小或方向的话，那么我们就不会把它放进这里的外力

将每一时步分为两个小步：

第一次求解运动方程：计算外力与惯性速度带来的位置变化
第二次求解运动方程：计算势能力带来的位置变化

首先要思考的是“求解运动方程”意味着什么？

我们的最终目标是要求出每一帧起始的粒子位置 $x$ ，求解该位置的方法就是求解常微分方程 $\frac{\partial x}{\partial t} = v$ 。但因为我们通常不可能直接求出解析解，特别是交互式应用里，而且一般也不需要解析解，所以我们一般会把它给离散化。

设第 $n$ 帧起始的粒子位置为 $x^{n}$ ，则 $x^{n+1} = x^n + \int_{t_n}^{t_{n+1}}vdt$ 。各类数值方法（显式欧拉法、隐式欧拉、韦尔莱积分法等等）本质都是在想法子解出来后半那个积分。在我们这边，求解运动方程就是要求这个积分。

那么这里为什么分了两步呢？因为两步考虑的系统是不同的。第一步只考虑到了外力与惯性速度，第二步只考虑到了势能力。通过区分开来这两部分，我们就能对其分别使用不同方法求解其运动方程。这一方法当然有弊端：显然在进行第二步的时候我们会忽略第一步的存在，也就会导致第二步会抵消掉一部分第一步的结果。

三、第一次求解运动方程：外力与惯性速度

XPBD选用的是韦尔莱积分法来进行这第一步：

$x^{n+1} = 2x^n - x^{n-1} + \Delta t^2M^{-1}F_{ext}$

其中 $F_{ext}$ 为外力和。（该式与原论文中的写法不太一样，只是因为我把 $v^n$ 给代入进去了而已）

这一步与PBD没有任何区别。

他之所以选用韦尔莱积分法是因为在第二步求解势能力的位置变化时不会把速度显式更新，所以如果第一步用某种需要显式更新速度的积分法的话，那问题会变得麻烦很多。因此他就选用了韦尔莱积分法这一把速度隐式蕴含在位置里的方法。

[Note: 所以是不是改一改第二步以后，这个第一步也能用别的积分法呢？]

四、第二次求解运动方程：势能力

首先先用隐式欧拉法来求 $x^{n+1} = x^n + \int_{t_n}^{t_{n+1}}vdt$ ：

$\left\{ \begin{aligned} x^{n+1} &= x^n + \Delta t v^{n+1} \\ v^{n+1} &= v^n + \Delta t a^{n+1} \end{aligned} \right.$

变换一下得到 $a^{n+1}$ 的表达式：

$a^{n+1} = \frac{x^{n+1} - 2x^n + x^{n-1}}{\Delta t^2}$

为了方便，下文中将省略 $n+1$ 的上标，并记\tilde{x}≡2x^n-x^{n-1}。

然后我们假设在这一步里考虑的总势能为 $U(x)$ ，那么其对应的势能力为 $F=-\nabla U(x)$ 。因为这一步里只考虑这些势能力的作用效果，所以有 $a=M^{-1}F$ ，代入得：

\begin{equation} -∇U(x) = M\frac{x - \tilde{x}}{\Delta t^2} \tag{1} \end{equation}

上式的复杂程度由 $U$ 的定义式决定，通常上式都是一个关于 $x$ 的非线性方程。

典型的做法是用牛顿法进行求解。设 $g(x)=M(x-\tilde{x})+\Delta t^2 \nabla U$ ，我们就是想求它的零点，使用牛顿法进行迭代（下标 $n$ 表示第 $n$ 次迭代的结果）：

g(x_n) + \left.\frac{∂g}{∂x}\right|_{x_n}(x_{n+1}-x_n) = 0

其中 $g(x_n)$ 很好求，带进去就完事了。但是\left.\frac{∂g}{∂x}\right|_{x_n}就很困难了，把它展开：

\frac{∂g}{∂x} = M + Δt^2H(U)

其中 $H(U)$ 为 $U$ 关于 $x$ 的黑塞矩阵，这一项通常而言都不是那么好求的。

当然也可以直接求 $H(U)$ ，然后就继续用牛顿法进行迭代，找到1式的解，从而也就得到了隐式欧拉法的解。如果这么硬算的话，那就是对只有势能力的系统用牛顿法硬解隐式欧拉法了。

1. 近似求解运动方程

XPBD里求解1式是用了一些技巧来近似计算的。

他首先假设势能可以表示为如下形式：

$U=\frac{1}{2}C^TC$

其中 $C=[C_0, C_1, \cdots, C_m]^T$ ，C_i(x):R^n→R。

形象一点理解这个表示的话，可以把 $C^TC$ 看成是一个关于 $C$ 的二次型。至于他这里的 $\frac{1}{2}$ ，我个人认为是因为他的势能以弹性势能为主，就沿用过来了。

把这一表示代入1式，得到：

\begin{equation} M(x-\tilde{x})+ \Delta t^2∇C^TC=0 \tag{2} \end{equation}

就求解难度而言1式和2式没有区别，只是把 $U$ 展开而已。注意∇U=∇C^TC。

为了找到能简化求解过程的方法，接下来XPBD里又引入了一个新的变量：

λ = -C

把它带入到2式中，将问题转换成一起求解λ, x的非线性方程组：

\begin{equation} \left\{ \begin{aligned} g(x, λ) &= M(x-\tilde{x})-Δt^2∇C^Tλ &= 0 \\ h(x, λ) &=C + λ &= 0 \end{aligned} \right. \tag{3} \end{equation}

求解3式的难度和求解1式的难度是一样的。

这里我们稍作停顿，来考察一下λ是什么，引入它以后的3式与1式的区别。

注意到\frac{∂U}{∂C}=C，所以λ=-C=-\frac{∂U}{∂C}，也就是说如果我们把 $C$ 视为从 $x$ 所在的 $n$ 维空间到某个 $m$ 维空间的R^n→R^m映射的话，那么λ就是 $U$ 在这个 $m$ 维空间里的负梯度，表示的是 $U$ 在这个m维空间里的最速下降方向。（ $U$ 在这个 $m$ 维空间里就只是个二次函数而已）

称这一由 $C$ 规定的 $m$ 维空间为“约束空间”，h(x,λ)=0描述的是λ与 $U$ 在约束空间里的最速下降方向的接近程度。

因为h(x,λ)=0其实并不总是能成立，所以我们设λ所对应的势能为 $U^*$ ，即λ=-\frac{∂U^*}{∂C}。

于是对于g(x,λ)，我们把第二项写开来就很明白了：

\begin{aligned} ∇C^Tλ = (\frac{∂C}{∂x})^T\frac{∂U^*}{∂C} = ∇U^* \end{aligned}

[TODO：这里有个问题，λ所描述的势能也是在约束空间里的吗？还是说是另一个 $m$ 维空间里？还是说其实无所谓在哪个 $m$ 维空间里？]

于是3式可改写为：

\begin{equation} \left\{ \begin{aligned} g(x, U^*) &= M(x-\tilde{x})-Δt^2∇U^* &= 0 \\ h(x, U^*) &=C + \frac{∂U^*}{∂C} &= 0 \end{aligned} \right. \tag{3.5} \end{equation}

如果精确求解该式的话，那么我们得到的解会和1式一样，毕竟本质只是改了个写法么。但如果我们不精确求解的话，那情况就不一样了。

观察3.5式，注意到 $g(x,U^*)$ 描述的是解 $(x, U^*)$ 与精确解的接近程度，而 $h(x,U^*)$ 描述的是 $U^*$ 的形状与原势能场 $U$ 的形状的接近程度。

如果 $h(x,U^*)=0$ 不严格成立的话，那么我们在求解 $g(x,U^*)=0$ 时相当于是换了一个势能场进行求解了。这就是XPBD最核心的部分。

我们先试着求解3.5式……哦，等等， $g$ 和 $h$ 都是关于 $U^*$ 的泛函数，解这方程组大概得用变分法了，我不会求解……

所以我们还是先把λ给取回来，不直接求 $U^*$ ，而是求它在约束空间里的负梯度λ[TODO: 为什么这样能降低求解难度？是因为限定了 $U^*$ 的形状的关系吗？]。用牛顿法求解3式，得到迭代方程组：

\begin{equation} \begin{bmatrix} M-\frac{∂(∇C^T)λ}{∂x} & -∇C^T \\ ∇C & I \end{bmatrix} \begin{bmatrix} Δx \\ Δλ \end{bmatrix} = - \begin{bmatrix} g(x_n, λ_n) \\ h(x_n, λ_n) \end{bmatrix} \tag{4} \end{equation}

其中Δx=x_{n+1}-x_n, Δλ=λ_{n+1}-λ_n， $I$ 为单位矩阵。最左边的系数矩阵里用的值全都是x_n, λ_n。

注意到\frac{∂(∇C^Tλ)}{∂x}=\sum_i λ_iH(C_i)（λ_i为λ的第 $i$ 个分量，注意不要跟λ_n, λ_{n+1}弄混），其中 $H(C_i)$ 为 $C_i$ 的Hessian矩阵。XPBD做的加速优化就集中在这里：

假设\sum_i λ_iH(C_i)=0
假设g(x_n, λ_n)=0[TODO: 我不理解这一近似。按理来说g(x_n, λ_n)并不是多难求的东西。文中说是跟Fast Projection算法相近的做法，尚待调研]

有了这两个近似后，4式就简化成了：

\begin{bmatrix} M & -∇C^T \\ ∇C & I \end{bmatrix} \begin{bmatrix} Δx \\ Δλ \end{bmatrix} = - \begin{bmatrix} 0\\ h(x_n, λ_n) \end{bmatrix}

这一线性方程组中 $M$ 是性质很好的一项，可以用 $M$ 做舒尔补求Δλ，得到计算式：

\begin{equation} \left\{ \begin{aligned} (∇CM^{-1}∇C^T)Δλ = -C-λ \\ Δx = M^{-1}∇C^TΔλ \end{aligned} \right. \tag{5} \end{equation}

显然，5式非常好求。

回过头去观察对4式做的近似与3.5式之间的关系。我们所做的最主要的近似是使\sum_i λ_iH(C_i)=0。我们回忆起刚才的分析λ=-\frac{∂U^*}{∂C}，代入这一近似中即可得到 $H(U^*)=0$ 。这也就是说XPBD所做的这一近似的含义为：

欲求的 $U^*$ 这一势能场，它关于 $x$ 的二阶导数为0.

所以我们迭代计算5式的过程其实可以视为这两步：

先寻找关于 $x$ 的二阶导为零的势能场 $U^*$ ，它在约束空间(由 $C(x)$ 规定的空间)里 $x_n$ 点对应的点上的梯度尽量接近原势能场 $U$ 的梯度。
随后在 $U^*$ 上计算隐式欧拉法。

[Note Q:那为什么是要在约束空间里的梯度相近呢？在原本的 $n$ 维空间里的梯度不是更好吗？A:因为那样的话h(x,λ)就变成h(x,λ)=∇U+λ=0了，在用牛顿法迭代方程组的时候不可避免地要计算 $H(U)$ 。]

五、讨论

原文里还回过头去联系了一下PBD，不过我不打算讨论那个了。阻尼也不讨论。

XPBD的特点我总结为三项：

将运动方程分为“外力与惯性速度”和“势能力”这两部分，然后分别用不同的方法进行求解。
在求解只有势能力的运动方程时使用隐式欧拉法并用牛顿法求解隐式欧拉法的非线性方程组。
在牛顿法迭代时，每一步都先寻找一个关于 $x$ 的二阶导为零的势能场 $U^*$ ，它在由 $C(x)$ 规定的 $m$ 维空间里 $x_n$ 所对应的点上的梯度尽量接近原势能场 $U$ 的梯度（即是说，设 $y_n=C(x_n)$ ，则∇_CU|_{y_n} ≈ ∇_CU^{* }|_{y_n}），随后在求解原问题(对隐式欧拉法应用牛顿法时得到的迭代方程)时使用这一简化的势能场 $U^*$ ，从而加速了求解。

第一点分开求解就相当于是“预测-校准”的方法。第二步求解势能力的运动方程时就相当于是把经过第一步后偏离合法状态的系统投影回由势能函数所规定的势能最低点（n维空间中的一片区域，通常就是指 $U=0$ 的流形）。想象系统在第一步后非常幸运地一直都处在势能最低点，那么这种情况下第二步就是可有可无的。例如无空气阻力的自由下落布料就是一个典型的情况，这种情况下第一步就是对所有粒子的重力，因为粒子之间的相对位置没有发生移动，所以布料内部也就没有弹力发生，第二步也就不会做任何事情。

第二点没什么好说的，就只是选用了隐式欧拉法和牛顿法而已。

第三点是他最主要的提速方法，每一步都寻找一个二阶导为零的势能场来近似原势能场。

这里我的一个思考是 $C$ 所规定的约束空间到底是什么？ $U^*$ 和 $U$ 在 $C$ 上的梯度相近的物理含义又是什么？尚待进一步思考。

除了这三点以外，原文中还作为一个突出亮点的是在构造势能时引入的compliance matrix: U=\frac{1}{2}C^Tα^{-1}C。

这个量的本质是在调控势能力的大小：显然势能的数值范围越大（乘上的系数越大，数值范围也就越大），其梯度也就越大，也就是势能力越大。落实到α的话，就是α越小，势能力越大。当α=0时，势能力就达到无穷大。这一项，咋说呢，我觉得其实不是那么重要，只是用来方便手动调参的而已，如果精巧设计势能函数的话其实也能达到一样的效果。

[TODO 寻找近似势能场这一做法似乎其实和ADMM很类似？还有Projective Dynamics是不是也可以这么解释？]

2022/5/21追记：拉格朗日乘数法视点

在Projective Dynamics的论文里注意到这一解释PBD的视点，觉得也挺有意思，就记录一下。

回顾式1：

\begin{equation} -∇U(x) = M\frac{x - \tilde{x}}{\Delta t^2} \tag{1} \end{equation}

重新排下顺序：

M(x - \tilde{x}) + Δt^2∇U(x) = 0

对它左边关于 $x$ 求个不定积分，得到函数Φ：

Φ(x) = \frac{1}{2}(x-\tilde{x})^TM(x-\tilde{x})+Δt^2U(x)

有∇Φ=M(x - \tilde{x}) + Δt^2∇U(x)。

那么在1式的解 $x^*$ 附近对Φ求极值点等价于直接求解1式。（因为极值点处Φ的梯度为0）（注意现在局部性还不能舍弃，因为 $U$ 是个未知形状的函数）

这就是最近挺流行的最优化视点，在这一视点下看XPBD的话，当然也可以按照本文的看法去处理，不过也可以从拉格朗日乘数法的视点去处理。

我们先把 $U$ 给扔掉，只看约束 $C(x)=0$ ，那么此时要求函数Φ的极值点，就是要求解以下条件极值问题：

\begin{aligned} \min_x\ &Φ(x)=\frac{1}{2}(x-\tilde{x})^TM(x-\tilde{x}) \\ s.t.\ &C(x)=0 \end{aligned}

这应该是很显然的，意思就是在 $C(x)=0$ 的前提下找到Φ(x)的最小值。

使用拉格朗日乘数法的方法将该问题转换为无条件极值问题，引入 $m$ 维列向量λ：

\min_x\ \frac{1}{2}(x-\tilde{x})^TM(x-\tilde{x})+λ^TC(x)

此处做一个近似，在 $\tilde{x}$ 处线性化 $C(x)$ ：C(x)≈C(\tilde{x})+∇C(\tilde{x})(x-\tilde{x})，代入得到：

\min_x\ \frac{1}{2}(x-\tilde{x})^TM(x-\tilde{x})+λ^T[C(\tilde{x})+∇C(\tilde{x})(x-\tilde{x})]

求解该无条件极值问题，就是对它分别求关于 $x$ 和λ的导数，并使导数为0，也就是要求解以下方程组：

\left\{ \begin{aligned} MΔx+∇C(\tilde{x})^Tλ=0 \\ C(\tilde{x})+∇C(\tilde{x})(x-\tilde{x}) = 0 \end{aligned} \right.