YASPS: A Symbolic Framework for Extensible, High-Performance IPC Simulation

一句话总结

YASPS 用两个关系算子 JOIN 和 UNION，把”谁参与能量、位置如何被参数化、局部导数如何映射回全局自由度”这些结构信息显式写进一个可微符号计算图，从而自动推导导数、全局稀疏结构和装配逻辑，并即时编译成 GPU 核，让 IPC 仿真在保持高性能（媲美手写的 GIPC）的同时获得极强的可扩展性。

研究背景

• 领域现状：Incremental Potential Contact（IPC）把弹性与碰撞处理统一成一个能量最小化问题，已成为接触密集型仿真的稳健范式。要跑得快，主流实现（PolyFEM、GIPC/Stiff GIPC、Stark）都靠针对固定能量、固定图元类型与参数化手写、手调的 GPU 核和装配逻辑。
• 核心痛点：这种高度特化制造了可扩展性瓶颈。加一个新能量、新图元或新参数化（如仿射体 Affine Body Dynamics）往往要重新推导一阶/二阶导数、重写全局梯度与 Hessian 的装配规则。更棘手的是碰撞能量：同一个能量常作用在混合参数化上，会引发参数化组合的指数爆炸——例如点-三角碰撞涉及 4 个顶点，每个顶点若有 3 种参数化，就有 \(3^4 = 81\) 种组合，每种都有不同的导数、Hessian 块布局和装配规则。
• 本文 idea：瓶颈的根源是仿真管线缺少一种统一方式来表示并传播结构信息。已有的自动/符号微分只在具体计算图上工作，遇到分支式的参数化选择就逐个分支各自微分，既无法抽出共享导数项，也无法从同一描述推出统一的稀疏/装配策略。YASPS 的思路是把”关系”和”参数化选择”提升为可微中间表示里的一等公民，用同一份关系描述同时驱动微分与结构感知的装配。

方法

YASPS（Yet Another Symbolic framework for Physical Simulation）提供一个 Python 前端：用户添加网格、创建图元（顶点、三角形、四面体）、声明图元间拓扑关系，然后在”场景/网格/图元”层级上定义符号属性，用 JOIN 与 UNION 表达跨拓扑关系的结构化聚合，最后指定哪些符号属性是能量项、对哪些属性做最小化。后端接手其余全部工作：符号微分器算出符号梯度与 Hessian，索引生成器从同一关系描述推出全局结构与”局部→全局”的块放置方式并做压缩，代码生成器把符号图 JIT 编译成 CUDA 核，最小化器调用这些核组装 \(Hx = g\) 并用 GPU 上的共轭梯度求解。

flowchart LR
  A["Python 前端: 定义能量与最小化目标 (JOIN / UNION)"] --> B["符号微分器: 局部/全局梯度与 Hessian"]
  A --> C["索引生成器: 全局稀疏结构 + 块放置 + 压缩"]
  B --> D["代码生成器: JIT 编译 CUDA 核 (含 PSD 投影)"]
  C --> D
  D --> E["最小化器: 组装 Hx=g, GPU 共轭梯度求解"]
  E --> F["输出 Newton 更新方向"]

关键设计：

• JOIN：跨关系组合依赖量。 用户声明某个量依赖于哪些实体并把它们的属性”拉”进来。比如四面体依赖 4 个顶点、评估体积能量时拉取它们的位置；每个顶点依赖单个仿射体变换、评估体空间形变时拉取该变换；每个排斥能量项拉取两个顶点位置作为输入。JOIN 让能量的构成关系（连接性、实例关系）在符号图里显式可见。
• UNION：表达关系内的可选参数化。 一个顶点位置既可以是直接参数化的自由变量，也可以由仿射体导出（\(p = A_b r + t_b\)，其中 \(A_b \in \mathbb{R}^{3\times 3}\)、\(t_b \in \mathbb{R}^3\) 是体 \(b\) 的变换，\(r\) 是静止位姿位置），不同顶点还能选不同方案。碰撞能量要支持新参数化时，用户唯一要改的就是重新声明哪些属性被 UNION——这正是避开参数化组合指数爆炸的关键：前端与后端代码都不随参数化数量爆炸。
• 在关系算子上做符号微分。 因为 JOIN 与 UNION 在 YASPS 的符号表示里本身可微，系统能对它们做符号微分，得到局部导数，并从同一份描述推导全局梯度与 Hessian 的诱导稀疏与块结构。论文还给出一个高效的二阶过程，复用中间 Jacobian、降低 Hessian 投影成本；局部 Hessian 会通过特征值扰动自动投影到半正定。
• 面向 IPC 负载的 GPU 编译与压缩装配。 IPC 的计算主体是大量局部项的高度并行评估。YASPS 把符号图编译成计算局部能量、导数、以及经关系层做索引抽取的核，并组装块稀疏的 Hessian 与梯度；配套的压缩全局 Hessian 结构（论文附录）显著降低存储与 SpMV 成本，这是它在求解阶段追平乃至超越手写系统的关键。

实验结果

主实验是”1–3 层布料落到斯坦福兔子、兔子再落到四角固定的布料上”的接触密集场景（\(dt = 0.01\text{s}\)，在 RTX 4090 + i9-13900F 上跑 200 帧），把 YASPS（未优化 / 手动优化版）与手写的 GIPC、以及 Stark 对比。下表取最大规模（59,997 顶点、3 层布料）一组：

系统	总时间 (s) ↓	CG 单次迭代 (ms) ↓	微分单次 (ms)	备注
YASPS	328.51	0.092	18.99	全自动、无手写核
YASPS 优化版	300.62	0.092	15.21	手动改写能量形式
GIPC	586.46	0.67	11.12	手写手调、含 MAS 预条件子
Stark	141.90	0.91	124.89	CPU、碰撞时线搜索失败仅跑 64 帧

核心结论：GIPC 用解析 Hessian 和 MAS 预条件子，在”微分单次耗时”和”CG 迭代次数”上占优（CG 次数少 2–3×），但因为不做存储压缩，单次 CG 迭代（受 SpMV 主导）慢得多；YASPS 的压缩技术带来接近 \(10\times\) 的单次 CG 加速，最终总时间反而更短。论文另报告：压缩后的约化矩阵相较未优化带来约 \(44\times\) 的相关加速；扩展性上把 1–25 个软体兔子丢进容器，各主要环节随碰撞对数量近似同趋势增长。代码量方面，加入仿射体参数化时 GIPC 仅能量/梯度/Hessian 就要多写至少约 1000 行、Stark 的接触摩擦模块已穷举列出约 1400 行各类接触场景，而 YASPS 后端零改动、仅前端加约 130 行 Python（定义仿射刚性能量 \(E_{\text{affine}}(A) = \tfrac{1}{2}\lVert A^\top A - I \rVert_F^2\) 只需 8 行）。

亮点与局限

• 亮点：
  • 用 JOIN/UNION 这一对关系算子把结构信息显式化，一份描述同时驱动符号微分、稀疏推导与装配，从根上消除了混合参数化的组合爆炸——扩展新参数化在前后端都不产生指数级子核。
  • 在保持高层 Python 接口的前提下，性能追平甚至优于手写手调的 GIPC；存储压缩带来的近 10× CG 加速是性能的主要来源。
  • 加新图元/能量的开发成本极低（约 130 行 Python），把可扩展性落成了可量化的工程收益。
• 局限：
  • 刻意不内置碰撞检测/连续碰撞检测（因参数化任意导致难以通用），实验中直接借用 GIPC 专为线性参数化写的 CCD，这部分未针对 YASPS 优化，占了相当比例的运行时间。
  • 不规定具体的 Newton 步进/线搜索策略，只交付线性系统的解，用户需自行整合步进、阻尼与终止逻辑。
  • 未优化版在微分与 CCD 环节仍慢于 GIPC 的原生集成，总时间的领先主要靠求解阶段的压缩优势换来。

延伸思考

• YASPS 延续了 Ebb、Simit 等关系型网格 DSL 的数据模型，但把”命令式局部计算”换成”声明式能量 + 对关系算子求导”，这条”关系代数 × 可微编程”的路线对更广的基于能量的仿真（FEM-MPM、FEM-SPH 等异构耦合）可能同样适用。
• 把参数化选择建模成 UNION、并能从同一描述推出稀疏结构，思路上和编译器里的联合类型/多态与稀疏张量编译很近；是否能与 Warp、DiffTaichi 这类通用可微编程系统在稀疏装配层面互补，值得追问。
• 由于碰撞检测被解耦到外部，若配上针对任意参数化的通用 CCD，或把 CCD 也纳入符号/关系框架统一优化，有望进一步缩小与原生集成实现的差距。