Force-Dual Modes: Subspace Design from Stochastic Forces

问题与动机

为可变形体的有限元仿真设计降阶模型（Reduced Order Modeling, ROM）子空间是加速图形与工程仿真的关键。困难在于：对于任意动态仿真，很难事先知道哪个子空间是最优的。每个子空间通常只擅长捕捉一小类运动，一旦仿真超出这一范围，就会出现物理上不合理的行为。一个常见的伪影是”远距离的诡异作用”（spooky action at a distance）：由于共享基向量，网格上互不相关的部分会呈现相关运动，例如抬起大象的手臂却导致远处的头和鼻子发生形变。

传统的线性模态分析（Linear Modal Analysis, LMA）之所以产生全局支撑的模态，正是因为它隐式假设了一种力分布：空间上不相关、几乎处处非零的单位方差高斯力。作者指出，现有方法并非缺少统计先验，而是在未言明的情况下假设了一个先验，当真实受力偏离该假设时就会脆弱失效。

核心思路

本文把模型降阶重述为一个概率框架，将子空间构造等价于对物理响应数据做低秩统计模型拟合。关键洞察是：直接指定位移分布不直观（需要预知材料行为与动力学），但对力建模却容易得多。于是作者主张先对力建模，再通过力与位移之间的线性对偶关系隐式地定义位移分布。

把系统抽象为 \(U(\omega) = \mathcal{L}(F(\omega))\)，其中 \(U\) 是关心的量（位移），\(F\) 是系统输入（外力等），二者都视为随机变量。线性降阶把 \(U\) 近似为少量线性分量之和：

\[U \approx BZ + \mu_U\]

给定度量 \(M\) 与 \(M\)-正交约束 \(B^T M B = I\)，最优子空间由如下广义特征值问题的前 \(m\) 个特征向量给出：

\[\Sigma_U M B = B \Lambda\]

因此只要知道关心量的协方差 \(\Sigma_U\)，就能通过求解广义特征值问题得到误差最小的仿真子空间。不同方法的差别，本质上是对该分布的不同假设。

从力分布导出位移协方差

作者从弹性动力学仿真的变分形式出发。将外势做一阶近似暴露外力 \(F = -\nabla C(U)\)，并在静止状态附近对物理能量做二阶泰勒展开，得到二次近似问题，其解为

\[U^* = H^{-1} F\]

其中 \(H\) 是静止状态处的能量 Hessian（刚度矩阵）。若假设力是多元高斯 \(F \sim \mathcal{N}(\mu_F, \Sigma_F)\)，则位移也服从高斯分布 \(U \sim \mathcal{N}(\mu_U, \Sigma_U)\)，且

\[\mu_U = H^{-1} \mu_F\]

\[\Sigma_U = H^{-1} \Sigma_F H^{-1}\]

这就是”Force-Dual”（力对偶）的核心：把易于设计的力分布通过线性化系统推送，得到位移上的对偶分布，再对其拟合低秩高斯模型即得子空间。

统一已有方法

这一框架把两类经典方法解释为力先验的两个极端：

• LMA 作为特例：取力协方差 \(\Sigma_F = M\)（质量矩阵，对应连续白噪声力）时，Force-Dual 恰好还原出 LMA 子空间。这是首次形式化刻画 LMA 隐含假设的力分布，对应”最大不确定性”：处处施加等方差、空间不相关的力。
• 格林函数作为特例：当力分布来自一组秩小于子空间维度的随机激励力时，框架还原格林函数子空间，对应”最小不确定性”：可能的力方向已知且固定。

Force-Dual 模型可以在这两个极端之间取得平滑折中，例如在心脏左半部分施加随机不相关的法向力先验，从而在不产生格林函数式生硬跳变的前提下，比 LMA 更精细地捕捉局部细节。

构造广义力分布

作者给出两类构造策略。

方差指定：只修改 \(\Sigma_F\) 的对角结构。把力方差在空间上局部化即可得到局部化的子空间；由于公式中用到能量 Hessian \(H\)，得到的子空间同时感知材料异质性、光滑性与边界条件。这使得用户可以把用户绘制的绑定权重（如有界双调和权重、双调和坐标或手绘权重）重新解释为载荷方差，得到尊重材料刚性的子空间（例如让鸡翅在不弯折骨头的情况下向后弯曲）。

低秩力分布：假设力由随机激励系数线性参数化 \(F = D A\)，其中 \(D\) 的列张成可能力空间，\(A \sim \mathcal{N}(\mu_A, \Sigma_A)\)，则

\[F \sim \mathcal{N}(D \mu_A, D \Sigma_A D^T)\]

据此覆盖多种交互模式：

• 点控制柄：用二次弹簧目标把顶点拉向目标位置，力为 \(F = \alpha S N A = D A\)。
• 接触：接触事件在表面顶点施加由接触帧法向与切向、用户绘制的接触斑块权重刻画的力，写成 \(F = D A\) 后可得到物理感知的局部响应（如食人魔坐下时床架整体弯曲而非仅床垫局部凹陷）。
• 气动驱动：软体机器人常见的气囊充气，对气囊轮廓顶点沿法向施加力，正确恢复夹爪弯曲行为。
• 肌肉驱动：采用二次肌肉激励能量模型，力为 \(F = -\sum_j v_j F_j d_j d_j^T : \frac{\partial F_j}{\partial U} A_j = D A\)，使先前只能用极粗几何的肌骨仿真得以交互式运行。
• 弹簧驱动：从质量弹簧外势出发，力为 \(F = \frac{\partial l}{\partial U}^T (l_0 - l(U))\)，可写成 \(F = D A\)，用于角色与机器人的弹簧执行器（如实时逆运动学控制龙的运动）。

多模态交互的高斯混合

单一力先验无法覆盖仿真中不同的驱动阶段（如某时刻是局部接触力主导，之后是全局惯性力主导）。作者用高斯混合建模力分布：

\[F = \sum_{k=1}^{K} \chi_k F_k\]

其中 \(\chi_k \in \{0,1\}\) 为掩码变量。混合分布的统计量为

\[\Sigma_F = \sum_{k=1}^{K} \pi_k \left( \Sigma_k + (\mu_k - \mu_F)(\mu_k - \mu_F)^T \right)\]

更重要的是，混合表述提供了运行时自适应选择子空间的概率化依据。给定观测到的力样本 \(\hat{F}\)，可用贝叶斯规则推断各分量的后验责任，取对数展开为

\[\log P(C=k \mid F=\hat{F}) = -\frac{1}{2}(\hat{F} - \mu_k)^T \Sigma_k^{-1} (\hat{F} - \mu_k) - \frac{1}{2}\log\det(\Sigma_k) - \log c\]

据此选出最可能的分量并即时切换对应子空间。由于交互力通常只作用在小范围子区域，可将分布边缘化到该子区域后快速求值，从而实现实时选择。这样即使任意时刻只用两个 skinning eigenmode，也能获得比固定子空间更真实、可交互的仿真。

实现

朴素求解广义特征值问题因需要稠密的 \(H^{-1}\) 而代价高昂。作者按 \(\Sigma_F\) 的结构分两种策略：

• 当 \(\Sigma_F\) 满秩且（块）对角时，求解等价问题 \(H \Sigma_F^{-1} H = M B \Lambda^{-1}\)。由于 \(\Sigma_F^{-1}\) 稀疏、\(H\) 与 \(M\) 也稀疏，可用现成稀疏 GEVP 求解器（scipy 的 eigs）。
• 当 \(\Sigma_F\) 非对角但低秩时，利用类 Cholesky 分解 \(\Sigma_F = L L^T\)、\(M = N N^T\)，通过对 \(N H^{-1} L = U S V^T\) 做奇异值分解，取前 \(m\) 个左奇异向量得到 \(B = N^{-1} U\)、\(\Lambda = S^2\)。在专用力分布下，Cholesky 因子可直接取 \(L = D\)，几乎免费获得。

求解器方面，默认采用带回溯线搜索的牛顿法并配合 cubature 积分加速；对某些凸能量则使用投影动力学（Projective Dynamics）求解器进一步加速。

实验与结论

运行时分析显示子空间仿真比全空间快数个数量级：25,613 顶点、103,307 四面体的龙可在实时下运行（约 448 FPS，加速约 3166 倍）。与数据驱动方法（POD/PCA）相比，在样本充分的极限下两者收敛到同一子空间，但 Force-Dual 直接给出闭式解，无需昂贵、繁琐的离线采样与稠密 SVD。在载荷方差空间衰减的消融中，力先验越接近真实受力，重构误差越低——即使先验只是适度局部化，误差也比朴素 LMA 基线小一个数量级，这在低模态数下尤为重要（因为约减线性系统的求解代价随子空间维度立方增长）。

与 LMA 一样，本方法建立在物理系统变分形式的线性化之上，但可直接复用工程与动画中增强线性降阶的既有工作（如旋转应变坐标 Rotation Strain Coordinates）来处理大非线性形变，无需对这些算法做额外改动。

总体而言，这是首个能够感知任意高斯力分布（涵盖接触、控制柄交互、肌肉、弹簧与气动驱动）来构造仿真子空间的方法，让子空间不仅贴合材料物理属性，也可针对具体应用中预期的力交互进行定制。