一篇论文的逐点精读 · 2026-05

读懂 FLASH：从零到 GPU 可变形仿真

想象一下：你要训练一个机器人把一件 T 恤叠整齐。布料会皱、会自己叠到自己上面、 会和桌面摩擦 —— 这是机器人学里最难的一类操作。FLASH 这篇论文的野心是， 在 一块 GPU 上同时跑上千个高保真布料仿真，让强化学习策略 "几分钟"就训练出来，并且 零样本迁移到真实机器人上。

要看懂它，你会撞上一长串名词：隐式积分、Projective Dynamics、ARAP、 Signorini 条件、Coulomb 摩擦锥、非光滑 Newton、Schur 补、惯性近似、 块对角 Cholesky、稀疏直接求解器…… 本文为 只会高数、线代、机器学习的读者，把这些一块一块拆开，从头讲。 不假设你做过任何物理仿真。每一个概念都配一个可以拖动的交互演示， 并且用我们从零复现这篇论文时踩到的真实坑来印证。

25+ 个关键概念 · 6 个交互演示 · 14 篇基础文献 · 复现代码实证 · 最后更新 2026-05-26

阅读路径建议

(1) 完全没碰过仿真 → 从 §1 一路读到底，每一部都建立在前一部之上， 像搭积木。(2) 懂物理仿真、只想看 FLASH 新在哪 → 直接跳到 第四部 §13–§15（惯性近似 Schur + 块对角 + 求解器）。 (3) 只想看复现踩坑的"血泪故事" → 跳到第五部 §16–§19， 那里有"同一个随机种子跑出两种物理"的诡异现象，和一个作者都承认"有点怪"的 workaround。

§1FLASH 到底是什么

FLASH 的全名是 Fast Learning via GPU-Accelerated Simulation for High-Fidelity Deformable Manipulation in Minutes（Luo 等，arXiv:2604.17513，2026 年 4 月）。 拆开看，它其实是两样东西缝在一起：

• 一个算弹性（布怎么被拉伸、弯折）的求解器 —— 用的是 Projective Dynamics（PD，§7）；
• 一个算接触（布和桌面、布和自己怎么不穿透、怎么摩擦）的求解器 —— 用的是把接触写成非线性互补问题（NCP，§9）再用非光滑 Newton（§11）解。

把这两个东西融合成一个 GPU 上能极致并行的迭代循环， 就是 FLASH。它能在单卡（论文用 RTX 5090）上把规模推到 超过 300 万自由度、30 FPS，并且同时跑成百上千个独立环境， 这正是强化学习（RL）疯狂采样所需要的。论文报告：用 FLASH 纯合成数据训练的折布策略， 不看任何真实演示就能迁移到真机叠毛巾、叠衣服。

为什么这件事值得一篇论文？因为现有的主流机器人仿真平台（比如 NVIDIA Isaac Sim） 在刚体上能轻松跑上千并行环境，但碰到可变形物体 （布、绳、软体）就要么慢、要么糙、要么并行不起来。FLASH 想填的， 就是"高保真可变形 × RL 级别并行"这块空白（详见第六部 §20–§21 的全景对比）。

§2一帧仿真在算什么

在讲任何细节之前，先把整个仿真的目标方程摆出来。这是论文的 Eq.1， 也是过去二十年几乎所有快速可变形仿真的共同起点。一帧仿真，本质上是在解一个优化问题：

$$ q \;=\; \arg\min_{q'}\left( \frac{1}{2h^2}\,\big\| M^{1/2}(q' - \tilde q) \big\|_F^2 \;+\; \sum_i \psi_i(q') \right) $$

这里每个符号都要认识 —— 它们会贯穿全文：

• $q'$：所有顶点的位置堆成的一个大向量（布上有 $N$ 个点就是 $3N$ 个数，即 $3N$ 个自由度）。我们要找的就是它。
• $h$：时间步长（比如 $0.01$ 秒）。仿真是一帧一帧往前走的。
• $M$：质量矩阵——一个对角矩阵，第 $i$ 个对角元就是第 $i$ 个顶点的质量 （每个顶点只管自己的质量、没有交叉项，所以是对角的）。$M^{1/2}$ 就是把对角元逐个开根号。
• $\tilde q = q_t + h\,v_t + h^2 M^{-1} f_{ext}$："惯性预测位置"（$q_t,v_t$ 是当前帧的位置和速度）—— 意思是"如果没有任何内力，只靠当前速度 $v_t$ 和外力 $f_{ext}$（主要是重力）惯性地往前飘，布会到哪"。
• $\psi_i(q')$：第 $i$ 个单元（一个三角形）的弹性势能 —— 布被拉伸/弯折得越厉害，这一项越大（§6 细讲）。

逐块读懂那个"加权范数" $\|M^{1/2}(q'-\tilde q)\|_F^2$

这一项是初学者最容易一眼跳过、却最该看懂的。我们把记号一层层剥开：

• $q' - \tilde q$ 是两个向量的差：新位置 $q'$ 减去惯性预测位置 $\tilde q$。 它本身是一个"位移向量"——每个顶点都有一个 3 维的小偏移。
• $\|v\|_F$ 是向量的长度（Frobenius 范数；对一个向量它就是最普通的欧氏长度 $\sqrt{\sum_k v_k^2}$）。平方一下就是 $\sum_k v_k^2$。
• 中间夹一个 $M^{1/2}$ 是关键。利用 $\|M^{1/2}u\|^2 = (M^{1/2}u)^\top(M^{1/2}u) = u^\top M u$， 而 $M$ 是对角的质量矩阵，所以这一项展开就是

$$ \frac{1}{2h^2}\,\big\| M^{1/2}(q'-\tilde q) \big\|_F^2 \;=\; \frac{1}{2h^2}\sum_{i} m_i\,\big\| q'_i - \tilde q_i \big\|^2 . $$

读成白话：越重的顶点，越"舍不得"离开它的惯性预测点（前面乘的 $m_i$ 越大， 同样的偏移付出的代价越大）。这正符合直觉——质量大、惯性大、越难被弹性拽走。

把这个 $\arg\min$ 读成一句白话："找一组新位置 $q$，让它既别离惯性预测的地方太远 （第一项），又别把布拉伸/弯折得太狠（第二项）。" 第一项是"惯性 / 动量想让你去哪"，第二项是"弹性想把你拉回什么形状"，两者拔河，平衡点就是这一帧的结果。

为什么"一帧仿真 = 解优化"是对的、而不是随便写的？这就是隐式积分的数学， 也是我们要讲的第一块基础积木（§5）。先记住这个方程，它是后面一切的母体。

§3为什么要"快 + 准 + 并行"

机器学习读者最熟悉的一句话是"数据就是燃料"。强化学习训练一个操作策略，动辄需要 上亿步的环境交互。如果每一步交互都要等一个慢仿真器，训练就是以"周"为单位的。 让仿真又快又能大批量并行，等价于把训练时间从"周"压到"分钟"—— 这就是 FLASH 标题里 "in Minutes" 的来历。

下面这张表把"为什么是这些技术选择"一次说清：

维度	朴素做法	FLASH 的选择	在第几节讲
时间积分	显式 Euler（步长必须极小，否则爆炸）	隐式 Euler（无条件稳定，大步长）	§5
解非线性	每帧完整 Newton（要重新分解大矩阵）	Projective Dynamics（矩阵固定，预分解一次）	§7
接触	加个惩罚弹簧（软、会穿透、要调参）	NCP + 非光滑 Newton（严格不穿透）	§9–§11
线性系统	稠密求逆 / 迭代法（接触一多就病态）	稀疏直接分解（块对角，无迭代悬崖）	§12–§15
并行	多进程跑多个 CPU 仿真	单 GPU 上千环境块对角批处理	§14

这张表也是全文的路线图：每一行的"FLASH 的选择"列，都对应后面一整节的讲解。 把这句话刻在脑子里 —— FLASH 不是发明了一个全新算法，而是把一串经典选择 （隐式 + PD + NCP + 稀疏直接解）重新组装成对 GPU 极度友好的形态。

核心循环 · 教学伪代码

下面这段伪代码就是论文 Algorithm 1 的骨架（多环境、GPU）。现在还不必懂每一行， 先看它的结构：一次性预分解 → 每帧检测一次接触 → 内层 local-global 迭代。 后文每一节都在解释这里的某一行。

# 一次性准备：把所有环境的系统矩阵拼成块对角，做一次 Cholesky 分解
A_bar = block_diagonal(A_1, ..., A_n)        # §7 §14
L_bar = cholesky(A_bar)                       # 预分解，之后只回代

while simulating:
    detect_collisions(all_envs)               # §16 —— 关键：每帧只检测一次！
    q_tilde = q + h*v + h*h * inv_M * f_ext    # §2 惯性预测

    for k in range(k_iters):                  # local-global 内层迭代 §7
        p = project_local(G @ q)              # 局部：每个三角形投影到最近旋转 §8
        # 组装接触 Jacobian J、互补残差，做惯性近似 Schur 补 §12 §13
        Z = J @ inv_M @ J.T + E               # §13 —— FLASH 的核心简化
        dlam = solve(Z, rhs) / (h*h)          # 解接触力增量 §12
        lam = lam + dlam
        q = global_solve(L_bar, b + h*h * J.T @ lam)   # 全局：一次稀疏回代 §7 §15

    v = (q - q_prev) / h                      # 从位移反推速度

§4四个痛点

接下来全文都围绕四个痛点展开。每一个基础概念、每一个 FLASH 的设计， 都是在解决其中某一个。先认识它们，后面读到对应技术时就会"哦，原来是为了治这个"。

痛点 (a)：刚度墙 —— 显式积分需要极小步长

布很"硬"（抗拉伸刚度 $k$ 大；刚度就是"拉它一点点、它回弹多大力"）。 一根又硬又轻的弹簧弹得特别快——它的固有振动频率 $\omega=\sqrt{k/m}$ 很高 （$k$ 越大、质量 $m$ 越小，弹得越快、$\omega$ 越高）。直觉上，东西动得越快，你就得用越小的时间步去"追"它， 否则一步迈太大就会冲过头、越冲越远。

这件事有个精确的门槛：显式积分要稳定，步长必须满足 $h \lesssim 2/\omega$。 布的刚度一高，$\omega$ 动辄上千，于是 $h$ 被逼到 $10^{-4}$ 秒以下——一秒要算上万步，根本跑不动。 这堵翻不过去的"刚度墙"，正是逼着我们改用隐式积分（§5 会证明隐式没有这堵墙）的原因。 下面的 Demo 1 让你亲手把步长推过墙、看显式当场爆炸。

痛点 (b)：接触是"非光滑"的

两个物体要么分离（间隙 $\gt 0$，接触力 $=0$），要么贴住 （间隙 $=0$，接触力 $\gt 0$），绝不会两个都为正。这个"非此即彼"的逻辑 画出来是一个直角拐角，不可导 —— 而 Newton 法需要导数。摩擦的"黏住 / 滑动" 切换同样是个尖角。怎么对一个有尖角的方程做 Newton？这是 §9–§11 的主题。

痛点 (c)：接触会摧毁稀疏性

弹性系统矩阵 $A$ 是稀疏的（一个顶点只和邻居耦合）。但只要把接触力 严格地求解进去，标准做法里会冒出 $A^{-1}$ —— 而稀疏矩阵的逆是稠密的 （§13 会解释为什么）。稠密就意味着没法跨环境分块、没法上千并行。 FLASH 用一个"惯性近似"把稠密的 $A^{-1}$ 换成对角的 $M^{-1}$，奇迹般地保住稀疏 —— 这是它的招牌（§13）。

痛点 (d)：要在一块卡上跑上千个环境

RL 要海量样本。把上千个独立布料环境的矩阵拼成一个块对角巨阵， 对它做一次 Cholesky，就等于同时分解了所有环境（块对角矩阵的 Cholesky 就是各块的 Cholesky）。 痛点 (c) 的稀疏性是这件事成立的前提（§14）。

Demo 1 · 刚度墙：显式 vs 隐式积分

时间步长 h 刚度 k 重置

左右两个弹簧-质点用同样的步长 h 积分，左边显式（Forward Euler）、 右边隐式（Backward Euler）。先把 h 调小：两边都稳。然后慢慢调大 h 或调大刚度 k —— 当读数里 $h\cdot\omega$ 越过 $2$，左边的能量条会冲破上限、当场"爆炸"， 而右边永远稳稳衰减。这就是痛点 (a)：显式有一堵翻不过去的刚度墙，隐式没有。

§5隐式积分 = 优化

我们从最底层讲起：仿真就是反复回答"下一刻物体在哪"。牛顿第二定律 $M\ddot x = f(x)$ 是一个二阶常微分方程，离散化时间后有两种走法。

显式 vs 隐式：在"现在"还是"未来"取力

把加速度写成差分。显式（forward）Euler 用当前位置算力： $x^{n+1} = x^n + h v^n,\; v^{n+1} = v^n + h\,M^{-1}f(x^n)$。 隐式（backward）Euler 用下一刻的位置算力：

$$ x^{n+1} = x^n + h\,v^{n+1}, \qquad v^{n+1} = v^n + h\,M^{-1} f(x^{n+1}). $$

区别看似只是 $f(x^n)$ vs $f(x^{n+1})$，但稳定性天差地别。对一个频率 $\omega$ 的振子， 显式 Euler 每步把振幅放大 $\sqrt{1+(h\omega)^2} \gt 1$ —— 能量只增不减，迟早爆， 而且 $h\omega$ 越大爆得越快。隐式 Euler 每步把振幅乘以 $1/\sqrt{1+(h\omega)^2} \lt 1$ —— 无条件稳定，任意大的 $h$ 都不会发散 （代价是引入一点数值阻尼）。这正是 Demo 1 演示的现象，也是布料仿真必须用隐式的原因（痛点 a）。

关键一步：隐式 Euler 其实是一个优化问题

上面我们说隐式 Euler"用下一刻的力"。把它写成方程：把 $v^{n+1}=(x^{n+1}-x^n)/h$ 代进 $v^{n+1}=v^n+h M^{-1}f(x^{n+1})$，整理得 $\,M(x^{n+1}-\tilde x)/h^2 = f(x^{n+1}) = -\nabla\Psi(x^{n+1})\,$（这里 $\tilde x = x^n+h v^n$ 就是惯性预测， $\Psi$ 是总势能，力等于势能的负梯度 $f=-\nabla\Psi$）。这个方程看着像要解一个非线性方程组—— 但它恰好是下面这个标量函数（叫增量势能，incremental potential）的梯度等于零：

$$ E(x) = \frac{1}{2h^2}\,\| x - \tilde x \|_M^2 \;+\; \Psi(x), \qquad x^{n+1} = \arg\min_x E(x). $$

为什么？直接对 $E$ 求梯度（高数里的链式法则 + "二次型 $\frac12 u^\top M u$ 的梯度是 $Mu$"）： $\nabla E(x) = \frac{1}{h^2}M(x-\tilde x) + \nabla\Psi(x)$。令它为零， $\frac{1}{h^2}M(x-\tilde x) = -\nabla\Psi(x)$ —— 正是上面那个隐式 Euler 方程。 所以"解隐式 Euler"和"找 $E$ 的最小值点"是同一件事。

对照 §2 的 Eq.1 —— 一模一样。"走一步隐式 Euler" 完全等价于 "最小化一个能量"。 这一步转换威力巨大：解方程是脆弱的（容易发散），而最小化能量可以借用优化里成熟的工具 （线搜索、信赖域、预条件），而且"能量下降"天然就是稳定性的保证。FLASH（以及它依赖的 PD） 整套机制，都建立在"把仿真当优化解"这个视角上。

Large Steps — Large Steps in Cloth Simulation

SIGGRAPH 1998 Baraff & Witkin · CMU · PDF

布太硬，显式积分被刚度墙逼到一秒几万步，根本跑不动（痛点 a）。

关键想法：第一次把隐式 Euler 系统地用到布料上 —— 无条件稳定，于是可以"迈大步"。每帧解一个大型稀疏线性系统（用改良共轭梯度）， 这是现代一切快速可变形求解器的祖宗。

vs 之前的显式质点弹簧：把"步长被刚度卡死"换成"每帧解一个线性系统"，从此能跑真实刚度的布。

§6弹性能与形变梯度

§2 的方程里有一项 $\psi_i(q)$ —— 单个三角形的弹性能。这一节就把它彻底拆开。 核心是一个矩阵：形变梯度（deformation gradient）$F$。

先补一个概念：什么是"映射的雅可比"

"变形"就是一个把静止位置 $X$ 映到当前位置 $x$ 的函数 $x = \varphi(X)$。 高数里，一个多元函数在某点附近可以用它的雅可比矩阵（一阶导数排成的矩阵）做线性近似： $\mathrm{d}x \approx F\,\mathrm{d}X$，其中 $F = \partial\varphi/\partial X$。直白说， $F$ 告诉你"在静止状态画的一根无穷小箭头 $\mathrm{d}X$，变形后变成了哪根箭头 $F\,\mathrm{d}X$"。 它就是形变梯度。一个三角形被假设是"线性单元"——内部每个点的形变方式都一样， 所以这个 $F$ 在整片三角形上是常数，能闭式算出来。

从三个顶点闭式算出 $F = D\,D_r^{-1}$

线性单元意味着映射是仿射的：$x = F X + t$（一个线性部分 $F$ 加一个平移 $t$）。 平移 $t$ 很烦人，但有个干净的消除办法——看边向量。取三角形顶点 $X_0,X_1,X_2$， 对两条边作差，平移 $t$ 就抵消了：

$$ x_1 - x_0 = F(X_1 - X_0), \qquad x_2 - x_0 = F(X_2 - X_0). $$

把"当前两条边"并成一个矩阵 $D = [\,x_1-x_0 \ \mid\ x_2-x_0\,]$、 "静止两条边"并成 $D_r = [\,X_1-X_0 \ \mid\ X_2-X_0\,]$，上面两式合写成 $D = F D_r$， 于是

$$ \boxed{\;F = D\, D_r^{-1}\;} $$

这就是形变梯度的闭式公式。注意 $D_r$ 只跟静止形状有关，是个常数， 所以它的逆 $D_r^{-1}$ 可以在仿真开始时预计算一次；之后每一帧只要用当前位置拼出 $D$， 一次矩阵乘法就得到 $F$。

这正是复现代码里 mesh.py 干的事 —— 静止时建局部 2D 帧、预计算 $D_r^{-1}$，每帧用当前位置组 $D$：

# flash/mesh.py · _compute_rest_frame —— 每个三角形预计算静止帧（只跑一次）
p  = rest[faces]                  # (T,3,3)：T 个三角形，各 3 个顶点的 3D 静止坐标
e1 = p[:, 1] - p[:, 0]            # 静止边 1 = 顶点1 - 顶点0
e2 = p[:, 2] - p[:, 0]            # 静止边 2 = 顶点2 - 顶点0

# —— 建一个贴着三角形平面的正交局部 2D 坐标系 (u1, u2) ——
u1 = e1 / norm(e1)                # 第一根轴：沿边 1 的单位向量
n  = cross(u1, e2); n = n/norm(n) # 三角形法向（垂直于平面，只是拿来求 u2）
u2 = cross(n, u1)                 # 第二根轴：平面内、与 u1 垂直的单位向量

# 把两条静止边投影到 (u1,u2)，得到 2x2 的静止边矩阵 Dr，再求逆
Dr = stack([[dot(e1,u1), dot(e2,u1)],
            [dot(e1,u2), dot(e2,u2)]])   # (T,2,2)
Dr_inv = inv(Dr)                  # 预计算的逆 —— 之后每帧用 F = D @ Dr_inv
area   = 0.5 * norm(cross(e1, e2))# 静止面积（叉积模长的一半），给弹性能加权

一行行看它在干嘛：

• e1, e2：从顶点 0 指向另外两个顶点的两条边——前面推导里的 $X_1-X_0$、$X_2-X_0$。
• u1, n, u2：造局部坐标系。u1 沿第一条边；n=cross(u1,e2) 是平面法向（垂直于整片三角形）；u2=cross(n,u1) 于是落在平面内、又垂直于 u1。 $(u_1,u_2)$ 就是三角形自己的"东 / 北"两个方向。
• Dr：把两条 3D 边各自点乘 $u_1$、$u_2$，得到它们的 2D 坐标， 拼成 $2\times2$ 矩阵——这就是 $D_r$。$2\times2$ 才好求逆。
• Dr_inv：预计算的 $D_r^{-1}$。它只依赖静止形状，整场仿真不变， 所以算一次、存起来，是后面每帧算 $F$ 的关键省时点。
• area：静止面积。三角形越大、它的弹性能在总能量里权重越大（§7 会用到）。

为什么"旋转不要钱"：极分解

物理上有一条铁律：纯刚体运动（平移 + 旋转）不应该有任何弹性能 —— 你把一块布整体转个身，它没被拉伸，能量必须是零。所以弹性能只能依赖 $F$ 里 "拉伸"的部分，不能依赖"旋转"的部分。极分解（polar decomposition） 正好把两者分开：

$$ F = R\,S, \qquad R\ \text{是旋转},\quad S = S^\top \succeq 0\ \text{是对称拉伸}. $$

于是奇异值 $\sigma_1, \sigma_2$ 就是两个主方向上的拉伸倍数：$\sigma=1$ 表示该方向 不拉不缩，$\sigma=1.5$ 表示拉长 1.5 倍，$\sigma=0.5$ 表示压扁一半。一个纯旋转（没有任何形变） 意味着 $\sigma_1=\sigma_2=1$，即 $S=I$、$F=R$。最简单也最常用的弹性能 —— As-Rigid-As-Possible（ARAP） —— 就是直接惩罚"$F$ 离一个纯旋转有多远"：

$$ \psi_{\text{ARAP}}(F) = \frac{w}{2}\,\| F - R \|_F^2 = \frac{w}{2}\big[(\sigma_1-1)^2 + (\sigma_2-1)^2\big]. $$

这里 $\|\cdot\|_F$ 是矩阵的 Frobenius 范数（把矩阵所有元素平方求和再开根， 等价于"把矩阵摊平成一个长向量再取欧氏长度"）。右边那个等号是关键：因为旋转不改变 Frobenius 长度， 可以证明 $\|F-R\|_F^2 = \sum_i(\sigma_i-1)^2$ —— 弹性能只数"每个方向拉伸偏离 1 多少"， 完全不关心旋转了多少角度。这把"旋转不要钱"这条物理铁律，干净地写成了一个公式。

Demo 2 · 拖动三角形：形变梯度 F 与最近旋转 R

回到静止

灰色是静止三角形，橙色是可拖动的变形三角形。橙色椭圆 = 把一个圆用 $F$ 映射后的样子 （= 应变椭圆），青色十字 = 最近旋转 $R$ 的朝向。 试试：(1) 整体转一下三角形 → 椭圆还是圆、能量 $\approx 0$（旋转不要钱）； (2) 把一个顶点拉远 → 椭圆被拉长、$\sigma_1\gt 1$、能量上升； (3) 把一个顶点拖过对边把三角形翻面 → $\det F\lt 0$、红色警告 亮起 —— 这就是 §8 要处理的"翻转"陷阱。

§7Projective Dynamics：local-global

现在我们有了优化目标（§5）和弹性能（§6）。但 ARAP 能量里有个 SVD，整个 $\arg\min$ 是非线性的， 直接上 Newton 每帧都要重新组装、分解一个大 Hessian —— 慢。 Projective Dynamics（PD） 是绕过这件事的精妙技巧，也是 FLASH 弹性部分的引擎。

把能量拆成"到约束流形的距离"

PD 的诀窍是把每个弹性能改写成"到一个约束集合的二次距离"，引入一个辅助投影变量 $p_i$：

$$ \psi_i(q) \;=\; \min_{p_i \in \mathcal{C}_i} \frac{w_i}{2}\,\| G_i q - p_i \|_F^2 . $$

对 ARAP 来说，约束集合 $\mathcal{C}_i$ 就是"所有旋转矩阵"，$G_i q = F_i$，最优的 $p_i$ 就是最近旋转 $R_i$。 这一改写让问题变得可分块，于是能用两个超便宜的步骤交替求解：

• 局部步（local）：固定 $q$，每个三角形独立、可完全并行地求自己的 $p_i$ —— 就是把 $F_i$ 投影到最近旋转 $R_i$（§8 的 SVD）。上千万个三角形互不相干，天生适合 GPU。
• 全局步（global）：固定所有 $p_i$，解一个线性系统求 $q$，把所有局部目标和惯性折中：

$$ \underbrace{\Big(M + h^2\textstyle\sum_i w_i G_i^\top G_i\Big)}_{A\ =\ \text{常数}\ =\ LL^\top}\, q \;=\; \underbrace{M\tilde q + h^2\textstyle\sum_i w_i G_i^\top p_i}_{b\ (\text{每次迭代变})}. $$

这里的 $G_i$ 是什么？它是个固定的稀疏矩阵，作用是"从全体顶点位置 $q$ 里， 线性地抽出第 $i$ 个三角形关心的量"——对 ARAP 来说 $G_i q$ 就是那片三角形的形变梯度 $F_i$ （回忆 §6，$F=D D_r^{-1}$ 对顶点位置是线性的，所以能写成一个矩阵 $G_i$ 乘 $q$）。 $G_i^\top G_i$ 因此只把"同一个三角形里的顶点"耦合起来，这正是 $A$ 稀疏的原因 （一个顶点只和它所在三角形的邻居有非零耦合）——这一点到 §13 讲稀疏性时还会回来。

为什么这快得离谱：矩阵 $A$ 是常数

注意全局步里：矩阵 $A$ 只取决于网格拓扑、权重、质量和步长 —— 全是固定的。 每次迭代变的只有右端项 $b$。所以可以在仿真开始时把 $A$ 做一次 Cholesky 分解 $A=LL^\top$， 之后每一步全局求解都只是两次三角回代（前代 + 后代），快到飞起。 "预分解一次，反复回代" 就是 PD（以及 FLASH）速度的核心。 复现代码里这个结构清清楚楚：

# flash/pd_solver.py —— 预分解一次（仅在钉住点变化时重做）
self._L_factor = torch.linalg.cholesky(A_aug)     # A = L L^T，分解一次

# ……每一帧的 local-global 内层迭代：
for _ in range(cfg.k_iters):
    F = compute_F(q, faces, C)                    # 当前形变梯度
    R = project_arap(F)                           # 局部步：投影到最近旋转（§8）
    b = mass * q_tilde + h*h * stretch_rhs(...)   # 组装右端项 b
    q = torch.cholesky_solve(b, self._L_factor)   # 全局步：一次回代，不重新分解

这段代码就是上面整套理论的落地，逐行对应：

• cholesky(A_aug)：把 $A$ 分解成 $LL^\top$，整场仿真只在构造时（和钉住点变化时）做一次。 这是最贵的一步，但只付一次。
• compute_F + project_arap：局部步——算当前形变梯度、投影到最近旋转 $R$（即 $p_i$）。这一步对每个三角形独立，天然并行（§8 细讲）。
• b = mass * q_tilde + h*h * stretch_rhs(...)：组装右端项 $b$。 mass * q_tilde 是惯性项 $M\tilde q$，stretch_rhs 是 $\sum_i w_i G_i^\top p_i$。每次迭代只有它在变。
• cholesky_solve(b, L)：全局步——用已经分解好的 $L$ 做两次三角回代解出 $q$， 不重新分解。这就是为什么 $k$ 次内层迭代依然便宜。

Demo 3 · Projective Dynamics 的 local-global 收敛

单步 local-global 自动迭代 弹簧权重 w 重置

一条两端钉住的链子，在重力下下垂。它从一条直线开始 —— 点"单步 local-global"看它一次次逼近悬链线：每一步先局部把每段投影到静止长度， 再全局解一次三对角线性系统（就是预分解过的 $A$）。看右下读数： 最大位移越来越小，说明在收敛。勾上"自动迭代"让它跑到收敛； 调"弹簧权重"看刚度怎么改变最终形状。这就是 §7 那个"局部并行投影 + 全局一次求解"的循环。

PD 的家谱上有几篇必须认识的论文：

Projective Dynamics — Projective Dynamics: Fusing Constraint Projections for Fast Simulation

SIGGRAPH 2014 Bouaziz, Martin, Liu, Kavan, Pauly · EPFL · project

隐式 Euler 每帧要解非线性系统，完整 Newton 太慢、不稳。

关键想法：把弹性能写成"到约束流形的二次距离"，于是隐式步分解成 并行的局部投影 + 一个常数矩阵的全局解。矩阵固定 → 预分解一次 → 实时、鲁棒、天生并行。这是 FLASH 弹性引擎的直接来源。

vs Baraff–Witkin 的每帧线性系统：把"每帧重新组装求解"换成"矩阵永远不变、只换右端项"。

ADMM ⊇ PD — ADMM ⊇ Projective Dynamics: Fast Simulation of Hyperelastic Models

TVCG 2017 Overby, Brown, Li, Narain · Minnesota · project

原始 PD 只能用特定的二次投影能量，换个本构模型就不灵。

关键想法：证明 PD 其实是ADMM（交替方向乘子法）的一个特例。 一旦看成 ADMM，就能推广到任意非线性材料和硬约束，同时保留 PD 那个常数预分解矩阵。 这给了"PD 的 local-global 到底在优化什么"一个干净的理论解释。

vs 原始 PD：把"只能用固定几种能量"推广成"任意本构 + 硬约束，矩阵照样固定"。

PBD / XPBD — Position Based Dynamics（2007）／ XPBD: Compliant Constrained Dynamics（2016）

JVCIR 2007 · MIG 2016 Müller 等 ／ Macklin, Müller, Chentanez · NVIDIA · PBD · XPBD

游戏引擎要极简、极稳、极快的布料/软体，不想碰力和加速度。

关键想法：PBD 直接在位置层面迭代投影约束（Gauss–Seidel）， 跳过力的计算。XPBD 给每个约束加一个柔度 $\alpha$ 并跟踪 Lagrange 乘子， 让刚度不再依赖迭代次数和步长 —— 这是 FLASH 接触里"柔顺约束 + 乘子"思想的近亲。

vs PD：PBD 更糙更快（无全局解）；XPBD 给了它物理一致的刚度和"接触力估计"，正是接触求解需要的。

§8ARAP 局部投影：SVD 与翻转陷阱

局部步的核心是一行字：把形变梯度 $F$ 投影到最近的旋转 $R$（这个 $R$ 就是 §7 里的投影变量 $p_i$）。 听起来简单，却藏着一个让无数人栽过的坑。先看正确做法：通过 SVD $F = U\Sigma V^\top$， 最近旋转是 $R = UV^\top$。这是经典的 正交 Procrustes 解。

翻转陷阱：$UV^\top$ 可能是镜像而不是旋转

上面的证明有个隐藏漏洞：我们只要求 $R$ 正交，但正交阵分两种—— $\det R=+1$ 的是真旋转，$\det R=-1$ 的是镜像（reflection，照镜子的左右翻转）。 当三角形被压扁、拉过头时，$UV^\top$ 的行列式可能算出 $-1$ —— 物理上这意味着三角形被翻了个面、 "里朝外"了。弹性求解器若把这种镜像当成"目标旋转"，单元会反演（inversion）， 能量曲面在这里有个虚假的低点，数值随即崩溃。

修法很经典：当 $\det(UV^\top)\lt 0$ 时，把 SVD 里最小那个奇异值对应的方向（$U$ 的最后一列） 翻号。几何上这相当于"宁可让那个方向多付一点拉伸能量，也不允许整片三角形翻面"， 强行得到一个真旋转（$\det R=+1$）。复现代码的 project_arap 一字不差地做这件事：

# flash/energy_arap.py —— 投影到最近旋转，并防止"翻面"
def project_arap(F):
    U, S, Vh = torch.linalg.svd(F, full_matrices=False)
    R = U @ Vh
    # 检测是否成了镜像（det < 0 = 翻转），若是则翻最后一列的号
    det_check = torch.linalg.det(R.transpose(-1, -2) @ F)
    sign = torch.where(det_check < 0, -1.0, 1.0).unsqueeze(-1)
    R_fixed = U.clone()
    R_fixed[..., -1] = R_fixed[..., -1] * sign     # 翻号 → 保证是真旋转
    return R_fixed @ Vh

逐行对照前面的数学：

• svd(F)：拿到 $F=U\Sigma V^\top$ 的三块。Vh 是 $V^\top$。
• R = U @ Vh：这就是上面证出来的最近旋转 $UV^\top$（暂未处理翻转）。
• det_check = det(R.T @ F)：检测翻转。$R^\top F = V\Sigma V^\top$ 是对称拉伸 $S$， 它的行列式 $=\sigma_1\sigma_2$；若为负说明发生了镜像（有奇异值"被迫取负"）。
• sign = where(det_check < 0, -1, 1)：负则准备翻号，正则原样。
• R_fixed[..., -1] *= sign：把 $U$ 的最后一列（对应最小奇异值的方向）乘上 $\pm1$， 这正是"翻最小奇异方向的号"。return R_fixed @ Vh 得到保证 $\det=+1$ 的真旋转。

Demo 2 里把一个顶点拖过对边时弹出的红色"⚠ 翻转"警告，就是 $\det F\lt 0$ 的时刻 —— 正是这段代码在守门。这是复现时 [已解决] 的一个经典坑（见仓库 PAPER_QA.md B2）。

ARAP — As-Rigid-As-Possible Surface Modeling

SGP 2007 Sorkine & Alexa · TU Berlin / Darmstadt · PDF

怎么定义一个"只惩罚形变、不惩罚旋转"、又好优化的弹性能？

关键想法：要求每个局部邻域尽量刚体（只许旋转平移）， 能量就是 $\|F-R\|^2$。用 local-global 交替求解：局部步用 SVD 抽出最佳旋转 （含 $\det$ 翻号），全局步解一个固定的 Laplacian 系统。 这个"local 抽旋转、global 解线性"的模板正是后来 PD 推广的原型。

vs 早期线弹性：把"对大旋转敏感、会鼓包"换成"旋转免费、对大变形也鲁棒"。

§9互补性与 Signorini 条件

弹性讲完了，现在进入整篇论文真正难的部分：接触。复现笔记里有句话很到位： "弹性求解器很容易对上论文，接触才是全部的战场。" 先理解接触的逻辑。

"非此即彼"：互补性

考虑布上一个点和桌面。设 $\text{gap}$ 是它到桌面的间隙、$\lambda_n$ 是桌面给它的法向支撑力。 物理只允许两种状态：分离（$\text{gap}\gt 0$，则 $\lambda_n=0$，桌子没碰到它哪来的力） 或接触（$\text{gap}=0$，则 $\lambda_n\ge 0$，可以有支撑力）。 绝不会"既有间隙又有力"。这个条件写成一行，就是大名鼎鼎的 Signorini 条件：

$$ 0 \le \lambda_n \;\perp\; \text{gap}(q) \ge 0 \qquad\Longleftrightarrow\qquad \lambda_n \ge 0,\ \ \text{gap}\ge 0,\ \ \lambda_n\cdot\text{gap}=0. $$

把这一行拆成三个条件读：$\lambda_n\ge 0$（支撑力只能"推"不能"拉"，桌子不会吸住布）、 $\text{gap}\ge 0$（不许穿透）、$\lambda_n\cdot\text{gap}=0$（乘积为零 ⇒ 两个里至少一个是零）。 第三条是灵魂：它把"非此即彼"压成一个等式——要么 $\text{gap}\gt 0$ 逼着 $\lambda_n=0$， 要么 $\lambda_n\gt 0$ 逼着 $\text{gap}=0$。

那个 $\perp$ 读作"互补"（complementarity），就是上面这套关系的缩写。 画在 $\text{gap}$–$\lambda_n$ 平面上，可行状态是一个直角拐角（沿横轴 $\lambda_n=0$ 或沿竖轴 $\text{gap}=0$ 的 L 形），这个尖角不可导——这正是接触"非光滑"、 Newton 法会卡住的根源（§11 专门治它）。如果 gap 和 $\lambda$ 对未知量是线性的，整个问题叫 线性互补问题（LCP）；加上摩擦的非线性锥后，叫非线性互补问题（NCP） —— 这就是 FLASH 里 "NCP" 的来历。论文的接触动力学（Eq.2）正是这个形式：

$$ M(q-\tilde q) - h^2 f_{int}(q) - h^2\!\!\sum_{j\in\mathcal{L}} H_j^\top \lambda = 0, \qquad \forall j:\ \phi_j(q,\lambda)=0, $$

其中 $H$ 是接触雅可比，$\phi_j$ 就编码上面的 Signorini（以及下一节的 Coulomb 摩擦）。

Stewart–Trinkle / Anitescu–Potra — 把接触+摩擦写成可解的时间步进互补问题

IJNME 1996 · Nonlinear Dyn. 1997 Stewart & Trinkle ／ Anitescu & Potra · DOI

接触的"非此即彼"逻辑没法用普通方程写，怎么把它变成能稳定求解的数学问题？

关键想法：把刚体接触 + Coulomb 摩擦写成一个时间步进 LCP， 并证明它对任意接触配置都可解、不穿透。非穿透即 Signorini 互补， 摩擦锥用多面体近似塞进 LCP。这是"接触 = 互补问题"这一整套范式的奠基。

vs 惩罚力法：把"软弹簧 + 调参 + 会穿透"换成"严格互补约束 + 可证不穿透"。

§10Coulomb 摩擦与摩擦锥

法向（不穿透）讲完，还有切向（摩擦）。Coulomb 摩擦的规则也是"非此即彼"：

• 黏住（stick）：切向力严格在锥内 $\|\lambda_t\| \lt \mu\lambda_n$ —— 不滑动，摩擦力要多少给多少（把物体钉住）；
• 滑动（slip）：切向力达到锥边界 $\|\lambda_t\| = \mu\lambda_n$ —— 开始滑动，摩擦力顶满、方向正对着滑动方向（最大耗散原理）。

把切向力约束在一个圆盘里，就是摩擦锥：

$$ \| \lambda_t \| \le \mu\,\lambda_n . $$

注意锥的半径正比于法向力 $\lambda_n$ —— 压得越重，能提供的摩擦越大。这就是为什么你按住一张纸 能搓动它、轻轻一碰却搓不动。为什么叫"锥"？因为在 $(\lambda_n,\ \lambda_{t1},\ \lambda_{t2})$ 三维力空间里，$\|\lambda_t\|\le\mu\lambda_n$ 描出的是一个圆锥（顶点在原点、半张角 $\arctan\mu$）；在固定 $\lambda_n$ 的那一层切片上，它就是一个半径 $\mu\lambda_n$ 的圆盘。

复现里把切向力投影回这个圆盘的代码非常直白，逐行就是上面规则的离散版：

# flash/contact.py —— 把切向冲量投影回 Coulomb 摩擦锥
def project_coulomb_cone(lam_n, lam_t, mu):
    lam_n = lam_n.clamp_min(0.0)               # Signorini：法向力非负（不能拉，只能推）
    lt_norm = norm(lam_t, dim=-1)              # 当前"想要的"切向力大小
    cap = mu * lam_n                           # 锥半径 = μ·λn
    scale = where(lt_norm > cap, cap/lt_norm, 1.0)
    return lam_n, lam_t * scale                # 超出锥 → 径向缩回边界（= 滑动）

逐行读：

• lam_n.clamp_min(0)：把法向力截到非负——这一行就是 §9 的 Signorini "只推不拉"。
• cap = mu * lam_n：算出此刻摩擦锥的半径 $\mu\lambda_n$。法向力越大，圆盘越大。
• scale = where(lt_norm > cap, cap/lt_norm, 1.0)：核心判断。 若"想要的"切向力 lt_norm 还在圆盘内（≤cap）→ scale=1，原样保留 = 黏住；若超出 → 按 cap/lt_norm 比例径向缩回边界， 方向不变、大小压到 $\mu\lambda_n$ = 滑动。
• "沿径向缩回最近的边界点"正是投影到锥的离散实现，也对应最大耗散里"力顶满、方向不变"。

Demo 4 把这套逻辑做成一个可玩的物理场景：你调法向压力、切向拉力和 $\mu$， 看方块到底是黏住还是滑动，以及当前接触力落在摩擦锥的哪里（空心点 = "想要的力"， 实心点 = 投影后的实际力）。

Demo 4 · Signorini 互补 + Coulomb 摩擦锥

摩擦系数 μ 法向压力 λₙ 切向拉力 λₜ

右边大图是摩擦锥（$|\lambda_t|\le\mu\lambda_n$），实心点是实际接触力、空心点是"需求"。 试试：(1) 切向拉力小 → 点在锥内、方块"黏住"不动； (2) 拉力调大越过锥边界 → 点跑到边界上、方块开始"滑动"； (3) 把法向压力调成负（抬起）→ 间隙打开、$\lambda_n=0$、 左上角 Signorini 图的点从竖轴跳到横轴 —— 没接触就没摩擦。

§11非光滑 Newton 与鞍点系统

现在矛盾来了：求解优化要用 Newton（需要导数），但 Signorini / 摩擦的"非此即彼"是个尖角， 在切换点不可导（痛点 b）。怎么办？这就是论文 Eq.5–6 解决的问题。

鞍点系统：位置和接触力一起解

把"弹性平衡"（来自最小化 $E$ 的 $\nabla_q L=0$）和"接触约束"$\phi(q,\lambda)=0$ 写在一起，对未知的位置 $q$ 和接触力（Lagrange 乘子）$\lambda$ 一起做一步 Newton 线性化， 就得到一个鞍点（saddle-point / KKT）系统。这里 KKT 指 Karush–Kuhn–Tucker 条件——约束优化最优解必须满足的一组方程：

$$ \begin{bmatrix} A & -J^\top \\ J & E \end{bmatrix} \begin{bmatrix} q \\ h^2\Delta\lambda \end{bmatrix} = \begin{bmatrix} g \\ \bar h \end{bmatrix}, \qquad J=\frac{\partial\phi}{\partial q},\quad E=\frac{\partial\phi}{\partial\lambda}. $$

上半部分是弹性（$A$ 来自 PD），下半部分是接触约束（$J$ 是接触雅可比）。 叫"鞍点"是因为解在 $q$ 上是极小、在 $\lambda$ 上是极大。

用"非光滑 Newton"对付尖角

对付不可导的尖角，经典武器是 Fischer–Burmeister 函数：

$$ \phi_{\text{FB}}(a,b) = a + b - \sqrt{a^2+b^2}, \qquad \phi_{\text{FB}}(a,b)=0 \iff 0\le a\perp b\ge 0. $$

花一分钟验证这个等价（令 $a=\lambda_n,\ b=\text{gap}$）：要 $a+b=\sqrt{a^2+b^2}$， 两边平方得 $a^2+2ab+b^2=a^2+b^2$，即 $\,2ab=0\Rightarrow ab=0\,$（至少一个为零）； 而 $a+b=\sqrt{a^2+b^2}\ge0$ 又强制 $a,b$ 不能为负。三个条件 $a\ge0,\ b\ge0,\ ab=0$ 正好凑齐 §9 的 Signorini。于是"一个尖角不等式逻辑"被压成了"一个等式 $\phi_{\text{FB}}=0$"。

这个 $\phi_{\text{FB}}$ 在原点 $(0,0)$ 处仍不可微（那个尖角没消失，只是被搬到一个点上）， 但它处处半光滑——可以套用 非光滑 / 半光滑 Newton： 一种允许在不可导点用"广义导数"（次梯度）代替普通导数的 Newton 变体，照样保持快速收敛。 实践中常用带 $\epsilon$ 的平滑版 $\sqrt{a^2+b^2+\epsilon^2}$ 把那个点也磨圆。 FLASH 的接触求解器正是建立在 Macklin 等人这套"NCP 函数 + 半光滑 Newton"框架之上。

Non-Smooth Newton — Non-Smooth Newton Methods for Deformable Multi-Body Dynamics

TOG 2019 Macklin, Erleben, Müller, Chentanez, Jeschke, Makoviychuk · NVIDIA · arXiv:1907.04587

刚体 + 可变形 + 接触 + 摩擦，怎么在一个统一框架里稳定、可并行地解（痛点 b）？

关键想法：把"动力学 + 接触"整体写成 NCP，用 NCP 函数 （Fischer–Burmeister 类）把互补条件变成非光滑方程，再用半光滑 Newton 迭代。 每次迭代归结为解一个对称线性系统，在 GPU 上用共轭残差解， 还配了个"互补预条件子"。FLASH 的接触部分直接长在这篇之上。

vs 把接触当软弹簧的 PBD 系：把"软、近似、调参"换成"严格 NCP + 半光滑 Newton + 每步一个对称解"。

§12Schur 补：把位置消掉，只解接触力

鞍点系统又大又不定（indefinite——它对应的"碗"既有向上弯的方向、又有向下弯的方向， 像个马鞍，所以叫"鞍点"；这种矩阵不能直接 Cholesky）。标准手法是 Schur 补（Schur complement）—— 一个在约束优化里无处不在的消元技巧，它能把这个难缠的鞍点系统化归成一个正定的小系统。

想法很简单：那个 $2\times2$ 分块系统里，位置 $q$ 的块 $A$ 又大又好对付 （它就是 PD 那个预分解过的稀疏矩阵），接触力 $\lambda$ 的块又小又麻烦。 那就先用"好对付"的把"麻烦"的解出来。一步步做：

鞍点系统两行是 $A q - J^\top(h^2\Delta\lambda) = g$ 和 $J q + E(h^2\Delta\lambda)=\bar h$。 第一行解出 $q$：

$$ q = A^{-1}\big(g + J^\top h^2\Delta\lambda\big). $$

把它代入第二行，把 $q$ 消掉，只剩 $\Delta\lambda$： $J A^{-1}(g+J^\top h^2\Delta\lambda) + E\,h^2\Delta\lambda = \bar h$，整理同类项：

$$ \underbrace{\big(J A^{-1} J^\top + E\big)}_{Z\ =\ \text{Schur 补}}\,(h^2\Delta\lambda) \;=\; \bar h - J A^{-1} g , \qquad \text{解出 } \Delta\lambda \text{ 后再回代}\ q = A^{-1}(g + J^\top h^2\Delta\lambda). $$

那个被夹出来的矩阵 $Z = J A^{-1}J^\top + E$ 就叫 $A$ 在这个分块系统里的 Schur 补。 整套手法就是"先消大块、解小块、再回代"——约束优化、电路分析、统计里到处都是它。

这个 $Z$ 的维数等于接触数（通常远小于自由度数），而且通常正定，可以用 Cholesky 或共轭梯度解。 在接触/多体里它有个名字叫 Delassus 算子：直观上回答"如果我在每个接触上施加单位力， 整个弹性体会怎么响应"。论文 Eq.7–9 就是这套：

$$ Z = J A^{-1} J^\top + E, \qquad \Delta\lambda = \tfrac{1}{h^2} Z^{-1}\big(J A^{-1} g - \bar h\big), \qquad q = A^{-1}\big(b + h^2 J^\top(\lambda^{k-1}+\Delta\lambda)\big). $$

§13惯性近似 Schur：FLASH 的招牌

终于到了 FLASH 真正的核心创新（论文 Eq.10）。前面 §12 的 $Z = J A^{-1} J^\top + E$ 有个致命问题，正是痛点 (c)。

为什么 $A^{-1}$ 是稠密的（而 $M^{-1}$ 是对角的）

$A$（弹性系统矩阵）是稀疏的：$A_{ij}\ne 0$ 只当顶点 $i,j$ 共享一个三角形 —— 它是网格连接图上的"图 Laplacian"，只记录"谁直接拉着谁的手"。但它的逆 $A^{-1}$ 是 完全稠密的，因为 $A^{-1}$ 是离散 Green 函数：$(A^{-1})_{ij}$ 表示 "在 $j$ 点戳一下，$i$ 点会动多少" —— 在一张连通的弹性网里，戳任何一点，所有点都会动一点。

于是 $Z = J A^{-1} J^\top$ 会把所有接触全耦合成一块稠密矩阵：跨环境也耦合、没法分块、 没法上千并行。FLASH 的解法干净利落 —— 把稠密的 $A^{-1}$ 直接换成对角的 $M^{-1}$：

$$ Z = J A^{-1} J^\top + E \;\;\approx\;\; J M^{-1} J^\top + E . $$

这个近似为什么物理上站得住？$M^{-1}$ 是惯性响应—— "在一个步长里， 给顶点一个单位冲量它会动多少，先不管邻居的弹性牵扯"。当步长 $h$ 不大时， 接触冲量到位移的映射主要由惯性决定，弹性修正是二阶小量。所以 $J M^{-1} J^\top$ 是真实接触柔度的一个相当好的近似，却便宜得多。论文 §IV-B 把它称作"轻量变体" （lightweight variant），并指出它"在改善 Schur 补稀疏性的同时保持数值鲁棒"。

Demo 5 · 接触度量 Z 的稀疏性：A⁻¹ 稠密 vs M⁻¹ 对角

用 M⁻¹（FLASH Eq.10） 环境数

每个格子是 $Z$ 的一块。勾选框开着（M⁻¹）时，$Z$ 是块对角—— 每个环境一个独立小块，对角线外全空，于是可以"一次分解 = 同时解所有环境"。 取消勾选（切到精确的 A⁻¹）→ 整块涂满：所有接触、所有环境全耦合成一坨稠密， 可并行环境数瞬间掉到 1。拖动"环境数"滑块看块结构怎么长出来。这就是 FLASH 能上千并行的根。

§14块对角多环境与 GPU 批处理

有了 §13 的稀疏性，痛点 (d) 迎刃而解。论文 Eq.11–12 把 $n$ 个环境的系统矩阵拼成一个大的 块对角矩阵：

$$ \bar A = \mathrm{blkdiag}(A_1,\dots,A_n), \qquad \bar q = \mathrm{col}(q_1,\dots,q_n). $$

"块对角"到底长什么样

$n$ 个环境互相独立——第 1 块布永远不会碰到第 2 块布。所以把它们的系统矩阵拼在一起时， 只有对角线上的方块有内容，块与块之间全是零：

$$ \bar A = \begin{bmatrix} A_1 & 0 & \cdots & 0 \\ 0 & A_2 & \cdots & 0 \\ \vdots & & \ddots & \vdots \\ 0 & 0 & \cdots & A_n \end{bmatrix}. $$

每个 $A_i$ 是第 $i$ 个环境的 PD 矩阵（$N\times N$ 量级）。那些零块就是"环境之间没有任何耦合" 的数学体现——这正是 §13 那个惯性近似 $Z\approx JM^{-1}J^\top$ 拼命要保住的结构 （要是用稠密的 $A^{-1}$，接触会把这些零块填上，环境就被搅在一起了）。

为什么块对角让分解"白嫖"：一个两行证明

关键事实：块对角矩阵的 Cholesky 分解，就是各块各自 Cholesky 因子再拼成块对角。 验证一下——设每个块单独分解 $A_i = L_i L_i^\top$，把这些 $L_i$ 拼成块对角 $\bar L=\mathrm{blkdiag}(L_1,\dots,L_n)$，那么

$$ \bar L\,\bar L^\top = \mathrm{blkdiag}(L_1 L_1^\top,\dots,L_n L_n^\top) = \mathrm{blkdiag}(A_1,\dots,A_n) = \bar A. $$

（块对角矩阵相乘，就是对应块相乘——因为零块乘任何东西还是零。）而 Cholesky 因子是唯一的， 所以这个拼出来的 $\bar L$ 就是 $\bar A$ 的 Cholesky 因子。换句话说， "分解那个巨大的 $\bar A$" 从来不需要真的当成一个大矩阵来做——只要分别分解每个小块就行， 块与块之间零填充（fill-in）。

RL 的甜点：同一件布的上千份拷贝 → 只分解一次

强化学习里，那 $n$ 个环境通常是同一件 T 恤的上千份拷贝（同样的网格拓扑、同样的钉住方式）， 只是初始状态/随机扰动不同。这意味着 $A_1=A_2=\dots=A_n$ —— 所有块完全一样！ 于是连"分别分解每个块"都不必，只分解一个 $A$、得到一个 $L$，让所有环境共用。 每一步要做的，只是把各环境不同的右端项 $b_i$ 一起喂给同一个三角回代。复现代码里就是这么写的：

# flash/pd_solver.py —— 一个 L，批量解 B 个环境的右端项
# b 的形状是 (B, N, 3)：B 个环境、N 个顶点、3 个坐标
rhs_flat = b.transpose(0, 1).reshape(N, B * 3)   # 摊成 (N, B*3)：B*3 个右端列
x = torch.cholesky_solve(rhs_flat, self._L_factor)   # 同一个 L，一次回代全部解完
q = x.reshape(N, B, 3).transpose(0, 1)               # 再 reshape 回 (B, N, 3)

• reshape(N, B*3)：把 B 个环境 × 3 个坐标的右端项并排成一个有 $B\!\times\!3$ 列的矩阵。 三角回代天生支持"多列右端"——多一列只是多算一遍前代/后代，矩阵 $L$ 不变。
• cholesky_solve(rhs_flat, L)：一次调用就把所有环境、所有坐标解完。 底层是一个大的批量三角求解 kernel，把 GPU 的几千个核心全喂饱。
• 注意 self._L_factor 只有一个——这就是"分解一次、所有环境复用"。

复现实测的批处理扩展性，和论文 Table II 同款"温和亚线性"：

批大小（环境数）	每步 ms	每环境每步 ms	环境·步/秒
1	4.46	4.46	224
8	7.16	0.89	1118
32	16.81	0.53	1904
64	29.42	0.46	2176
128	54.89	0.43	2332
256	107.4	0.42	2384

这张表怎么读？看"每步 ms"一列：从 1 个环境到 8 个环境，总时间只从 4.46 涨到 7.16 ms —— 环境数 ×8，时间却几乎没动，因为前 8 份订单基本是在"填满闲着的工人"，近乎免费。 这段叫亚线性区。等核心被占满（约 32–64 环境往后），再加环境就只能排队， 总时间转为线性增长（64→128→256 时间大致翻倍）。

再看"每环境每步 ms"一列——这是真正的吞吐指标：从 4.46 ms 一路降到 0.42 ms（快了约 11×）并触底。这个 0.42 ms 就是 GPU 喂饱后的"计算下限"， 再加环境也降不下去了。换算成吞吐（最右列），单卡稳定在约 2400 环境·步/秒。 对 RL 来说，这意味着用一块卡、一个进程，就能持续产出训练所需的海量交互 —— 而这一切的前提，是 §13 那个让矩阵保持稀疏 + 块对角的惯性近似。把链条连起来： 稀疏（§13）→ 块对角可独立分解（§14）→ 一个 $L$ 批量回代喂饱 GPU → 上千环境近乎免费。

Algorithm 1 全图：把所有概念串起来

现在你已经认识每一块了，回看论文主循环（这次每一行都标了对应章节）：

# 论文 Algorithm 1 —— 多环境 GPU 主循环（带章节索引）
A_bar = blkdiag(A_1, ..., A_n)          # §14 块对角
L_bar = cholesky(A_bar);  S_bar = inv(L_bar)   # §7 §12 预分解
while simulating:
    collision_detection(all_envs)        # §16 每步一次（不是每迭代！）
    q_tilde = q + h*v + h*h * inv_M * f_ext      # §2 §5 惯性预测
    for k in range(k_max):               # local-global 内层
        p_bar = project(G_bar @ q)        # §8 局部投影到最近旋转
        # 组装 J, E, g, h，做惯性近似 Schur
        Z_bar = J_bar @ inv_M @ J_bar.T + E_bar  # §13 招牌简化
        dlam  = (1/h**2) * solve(Z_bar, h_bar - J_bar @ S_bar.T @ S_bar @ g_bar)  # §12
        lam   = lam + dlam               # 累积接触力
        q     = S_bar.T @ S_bar @ (b + h*h * J_bar.T @ lam)   # §7 §15 全局回代
    v = (q - q_prev) / h                 # 反推速度

先认一个记号：变量头上的横线（$\bar A,\ \bar q,\ \bar J\dots$，代码里写成 A_bar 等）表示"把所有 $n$ 个环境堆叠在一起"——前面 §14 讲的批处理就藏在这道横线里。 还有一个简写：$\bar S^\top\bar S = \bar A^{-1}$（因为 $A=LL^\top$、$S=L^{-1}$，见 §12）， 所以代码里凡是 S_bar.T @ S_bar @ (...) 都读作"用预分解的 $A$ 回代一次"。 现在逐行走：

准备阶段（循环外，整场只跑一次）

• A_bar = blkdiag(A_1,...,A_n) — 把 $n$ 个环境的 PD 矩阵拼成块对角巨阵（§14）。
• L_bar = cholesky(A_bar); S_bar = inv(L_bar) — 预分解一次（§7）。 $\bar S=\bar L^{-1}$ 就是论文那句"显式逆"的写法；实践中不必真的造出 $\bar S$， 把它当成"用 $\bar L$ 做回代"即可（§12 提醒过这句措辞要小心读）。块对角 ⇒ 这一次分解 = 同时分解所有环境。

每一帧（while 循环）

• collision_detection(all_envs) — 每步只检测一次、缓存接触对（§16 的灵魂）。 把它挪进下面的内层 for 里，就会触发那个让布炸到 1500 mm/s 的能量注入。
• q_tilde = q + h*v + h*h * inv_M * f_ext — 惯性预测 $\tilde q$（§2/§5）： "假如没有任何内力，只有当前速度和重力，布会漂到哪"。它是这一帧优化的"靶心"。

内层 local-global 迭代（for k，§7 的循环，Demo 3 演示的那个收敛过程）

• p_bar = project(G_bar @ q) — 局部步（§8）：G_bar @ q 从当前位置抽出每个三角形的形变梯度 $F$， project 把它投到最近旋转 $R$（即投影变量 $p$）。每个三角形独立、完全并行。
• # 组装 J, E, g, h — 准备接触线性化的零件（§11）：$J=\partial\phi/\partial q$ 是接触雅可比、 $E=\partial\phi/\partial\lambda$ 是小正则、$g$ 和 $\bar h$ 是鞍点系统的右端块。
• Z_bar = J_bar @ inv_M @ J_bar.T + E_bar — FLASH 招牌（§13）： 用对角 inv_M（$M^{-1}$）而不是稠密 $A^{-1}$ 来拼接触度量 $Z$，于是 $Z$ 保持稀疏 + 块对角。
• dlam = (1/h**2) * solve(Z_bar, h_bar - J_bar @ S_bar.T @ S_bar @ g_bar) — 解接触力（§12 的 Schur 补）：括号里 S_bar.T @ S_bar @ g_bar $=\bar A^{-1}\bar g$， 整体就是 $\bar h - \bar J\bar A^{-1}\bar g$，解出接触力增量 $\Delta\lambda$。
• lam = lam + dlam — 把增量累积进总接触力（跨内层迭代热启动，frictional history）。
• q = S_bar.T @ S_bar @ (b + h*h * J_bar.T @ lam) — 全局步（§7/§15）： $=\bar A^{-1}(\bar b + h^2\bar J^\top\bar\lambda)$，把惯性、弹性投影、接触力一起折中解出新位置。 这一步只是对预分解的 $\bar A$ 回代一次（§15 讲它为什么不能落进 PCG 悬崖）。
• 循环回到 p_bar：local 和 global 交替 $k$ 次，一步步逼近这一帧的平衡（§7 Demo 3 里"最大位移越来越小"就是它）。

收尾

• v = (q - q_prev) / h — 用这一帧的位移反推速度（隐式 Euler 的速度更新 $v=(q_{t+h}-q_t)/h$）。 §18 那个 G2D 速度阻尼、以及常规的全局阻尼，就加在这一步之后。

把这张图竖着读一遍，本文前 14 节的每个概念都在里面各就各位： 隐式预测（§2/§5）→ 局部投影（§8）→ 惯性近似 Schur（§13）→ Schur 解接触力（§12）→ 全局回代（§7/§15）， 外面套着块对角批处理（§14）和每步一次的碰撞检测（§16）。这就是 FLASH。

§15把线性求解搬出关键路径

最后一块拼图：那个全局线性解 $A^{-1}b$（或 $Z^{-1}$）到底怎么算最快？这一节也是复现里 最戏剧性的一段，因为我们亲历了所谓的 "PCG 悬崖"。

两条路：迭代法（CG） vs 直接分解（Cholesky）

解 $Ax=b$（$A$ 对称正定）有两种思路，性格完全相反：

• 直接法（Cholesky 分解）：一次性把 $A=LL^\top$ 拆开，再两次回代得到精确解。 它没有"迭代次数"——给定 $A$，要做多少运算是写死的，跟矩阵"好不好解"无关。 代价是分解会产生填充（fill-in，下面解释），大网格上内存/时间会涨。
• 迭代法（共轭梯度 CG）：从一个猜测出发，每步只做稀疏矩阵-向量乘 （便宜、省内存、压根不碰 $L$），一步步逼近真解。但它要迭代多少次取决于矩阵"好不好"—— 由条件数 $\kappa$ 决定。

条件数 κ 与那条 √κ 定律

条件数 $\kappa(A)=\lambda_{\max}/\lambda_{\min}$ 是最大与最小特征值之比。把解 $Ax=b$ 想成在一个碗 $f(x)=\tfrac12 x^\top A x - b^\top x$ 里往最低点滚（碗底正是解）；$A$ 的特征值就是碗在各主方向上的"陡峭度"， $\kappa$ 就是这个碗的长宽比。CG 的误差经 $k$ 步后满足：

$$ \frac{\|x_k-x^\star\|_A}{\|x_0-x^\star\|_A} \;\le\; 2\left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^{k}. $$

要达到给定精度，所需迭代次数约为 $\tfrac12\sqrt{\kappa}\,\ln(2/\varepsilon)$ —— 正比于 $\sqrt{\kappa}$ （注意是根号 κ，这已经是 CG 比朴素梯度下降聪明的地方，后者正比于 κ 本身）。 一句话记住：κ 涨 100 倍，CG 迭代次数约涨 10 倍。

为什么 κ 会爆，而这正是接触的命门

三件事会把 κ 推高：材料越硬、时间步越小、尤其是冗余 / 近乎冲突的约束。 最后一条是接触的命门（§17 细讲）：密集自接触里，相邻接触共享顶点、法向近乎平行， 于是约束矩阵冒出一堆近乎重复的行——两行几乎线性相关，把最小特征值 $\lambda_{\min}\to0$、 $\kappa\to\infty$。此刻 CG 的迭代次数像撞墙一样飙升，性能从悬崖上掉下去。这就是 "PCG 悬崖"。

预条件为什么没把它救回来

对付大 κ 的标准武器是预条件：找一个便宜的近似 $P\approx A$，改解 $P^{-1}Ax=P^{-1}b$， 让 $P^{-1}A$ 的条件数远小于 $A$（几何上 = 把长窄峡谷"重新拉圆"）。最简单的 Jacobi 预条件就是拿 $A$ 的对角线当 $P$。但复现里 Jacobi-PCG 没救回来—— 因为这里的病态是结构性的（约束彼此冗余、近乎线性相关），而 Jacobi 只会缩放对角， 根本消不掉"两行几乎一样"这种问题。实测：约 27 万自由度、开自碰撞时， Jacobi-PCG 撞到 1600 次迭代上限仍不收敛，622 ms / 1.6 fps。

直接法没有"悬崖"：cuDSS 的 42× 一锤定音

换上 NVIDIA cuDSS（GPU 稀疏直接求解器）做 Cholesky 分解后，同一个场景 66.6 ms / 15 fps —— 42× 加速。为什么这么干净？因为直接法根本没有"迭代次数"： 它的代价是 $O(\text{nnz}(L))$（因子 $L$ 里非零元的个数），只跟稀疏结构有关、跟条件数无关。 κ 再大、约束再冗余，分解和回代的运算量一个都不多——悬崖被直接填平。

网格	自由度	self-coll OFF（ms / fps）	self-coll ON（ms / fps）
101×101	30 k	3.54 / 282	5.75 / 174 ← 论文基准点 ×5.8
301×301	272 k	11.4 / 88	66.6 / 15（PCG 时代是 1.6）
601×601	1.08 M	29.7 / 33.7 ★ 命中论文 30fps	744 / 1.3
1001×1001	3.00 M ★ 论文目标	69.4 / 14.4	4847 / 0.2

几个值得记住的数字：关掉自碰撞时，100 万自由度跑到 33.7 fps，命中论文"百万级 @ 30fps"； 300 万自由度（论文目标 N）跑到 14.4 fps，落在论文的 2× 以内。 剩下的差距是内层迭代核函数的常数因子，不是任何算法悬崖 —— 这本身就是个好消息， 说明 FLASH 的算法结构确实没有规模上的硬墙。

§16能量从哪来：每步检测 vs 每迭代检测

从这里开始是复现的血泪故事 —— 它们比任何公式都更能让你真正"懂"接触求解。 第一个，也是最有启发的一个：一块布扔到地面上，开着自碰撞，它永远不肯停， 一直扭动 20 多秒，越扭越欢。用户当时的描述是"好像有鬼在控制"。

烟枪证据：系统在"制造"能量

Projective Dynamics 在弹性部分是无条件稳定、只会耗散能量的。所以布不肯停，只能意味着 接触路径每步在往系统里注入能量。复现里用能量曲线坐实了这一点： 在释放后 1.75 秒，动能飙到 1.93 J，比释放瞬间的峰值 0.97 J 还高 —— 一个耗散系统的动能绝不可能超过初始峰值。积分器变成了一个能量泵。

真凶：接触集在迭代之间"翻动"

根因是个极其微妙、却极其常见的错误。我们最初把碰撞检测放在了 local-global 内层循环里， 每次迭代都重新检测。后果：

一个顶点在 iter 1 距离 9.8 mm 触发接触 → 被推到 iter 2 的 10.2 mm 不再接触 → iter 3 弹性又把它拉回 9.5 mm 再次触发 …… 每一次"重新捕获"都注入一份全量冲量 → 布永远 settle 不下来。

论文 Algorithm 1 第 4 行其实把"碰撞检测"明明白白放在了内层循环之外（每步一次）。 我们照字面改成"每步检测一次、内层只读缓存的接触对"后： silk 布 7 次 lift→release 测试，最大末速从 1500+ mm/s 直接掉到 50 mm/s。 一行代码的位置，决定了系统是稳还是炸。

怎么定位真凶：5 个对照实验

定位这种 bug 不能靠瞎试。复现里先列了 5 个可能机制，再设计 5 组只改一个变量的 对照实验（D1–D5），用"末 2 秒还在动多快"和"动能有没有超过释放峰值"两个指标判读：

实验	改的变量	动能峰值	末速 |v|	结论
D1	关掉自碰撞	0.97 J	0 mm/s	0.5 秒内静止 → 能量来自自碰撞路径
D2	基准（复现用户视频）	1.28 J	856 mm/s	永不静止
D3	关掉自碰撞摩擦	1.24 J	1418 mm/s	更糟 → 摩擦本是耗散，去掉只剩注入
D4	接触带 5mm→1mm	0.90 J	263 mm/s	几乎静止 → 真凶是阈值附近的翻动
D5	迭代 8→32 次	2.30 J	1601 mm/s	更糟 → 迭代越多注入越多

两个最反直觉的结果锁定了机制：D5（多迭代更糟）说明迭代本身在产能； D3（去摩擦更糟）说明摩擦是耗散源、不是问题源。真凶是 D4 指向的—— 顶点在接触阈值附近每次迭代翻动开/关，每翻一次注入一份冲量。把接触带从 5mm 收到 1mm， 阈值附近的边缘接触少了约 200×，注入随之消失。

Demo 6 · 能量注入：每步检测一次 vs 每迭代重检测

每迭代重新检测（错误做法） 重置小球

一个小球落到地面反弹。默认（青色）是"每步检测一次"—— 回弹带恢复系数 < 1，能量曲线衰减、小球 settle，对应复现里的 50 mm/s。 勾选"每迭代重新检测"（橙色）→ 每步的 4 次内层迭代各注入一份冲量， 能量曲线一路爬升直到爆炸，对应复现里的 1500 mm/s。 这就是把碰撞检测放错位置的代价。

§17病态的 Z 与 PCG 悬崖（以及两个吓人的 bug）

§15 讲了悬崖的现象，这里讲它的根因，并附两个复现里"以为是物理 bug、其实是数值 bug"的故事。

为什么 Z 会病态：近乎重复的约束

§15 说"冗余约束让 κ 爆"，这里把它具体到接触上。把一块 51×51 的丝绸抓起来悬空， 自接触对能到 8 万个，$Z = J M^{-1}J^\top + E$（§13）是个 8 万维稀疏矩阵。 回忆 §13 那个矩阵元：两个共享顶点 $v$ 的接触对 $p,q$ 的耦合是

$$ Z_{p,q} = \frac{n_p\cdot n_q}{m_v} + \cdots, \qquad Z_{p,p} = \frac{1}{m_v} + \cdots, $$

其中 $n_p,n_q$ 是两个接触的法向。密集自接触里，相邻的接触法向几乎指同一个方向 （$n_p\cdot n_q\approx 1$），于是 $Z_{p,q}\approx Z_{p,p}$ —— $Z$ 的第 $p$ 行和第 $q$ 行几乎一模一样。

病态的 Z 怎么把布炸上天

接下来就是 §15 那条 $\sqrt\kappa$ 定律的反面教材。解 $Z\,\Delta\lambda=\text{rhs}$ 时， 病态意味着右端项里一丁点噪声，会被 $\kappa$ 倍放大到解 $\Delta\lambda$ 上 （$\|\delta\Delta\lambda\|/\|\Delta\lambda\| \lesssim \kappa\cdot\|\delta\,\text{rhs}\|/\|\text{rhs}\|$）。 于是 CG 要么迭代到上限仍不收敛、要么吐出一个被噪声污染的巨大乱真 $\Delta\lambda$； 这个乱真接触力经 $M^{-1}J^\top\Delta\lambda$ 变成一脚巨大乱真冲量，把布踹飞。 复现实测：纯 CG 让布炸到 24000+ mm/s（正常 settle 应是 ~50 mm/s）， 加 Jacobi 预条件子也救不回来（仍 24000 量级）——因为病态是结构性的（约束彼此冗余）， Jacobi 只缩放对角，消不掉"两行几乎一样"（正应了 §15 的结论）。

真正的修法：active-set（以及复现为什么没走完）

根治办法是 active-set（活跃集）：那些"近乎重复"的约束其实是冗余的—— 保留其中一个就够了。每次迭代只挑出真正在穿透的接触（活跃集），把 $Z$ 缩小到只含这些行/列， 在这个满秩、良态的子空间里求解，κ 就回到正常范围。复现里这条（对应论文 Eq.8/9 的完整 Schur 解） 实现起来工程量大、还没打通，于是暂时禁用，退回到"位置投影 + 一个能量阻尼 workaround"—— 也就是下一节 §18 那个"擦地板"的 G2D。大规模密集自接触下 $Z$ 的良态求解，至今是开放问题。

Bug 一 · float32 让自由落体"偏航"12 cm

早期一块自由下落（没碰任何东西）的布，竟在 0.4 秒内横向漂了 12.5 cm—— 自由落体怎么会拐弯？查到根因是 float32 精度。

为什么这里对精度这么敏感？自由下落时布几乎不变形，形变梯度 $F\approx I$（近单位阵）， 它的最近旋转 $R$ 也近乎单位阵。可"几乎是单位阵"意味着真正有意义的信息藏在小数点后很靠后的位—— 而 float32 只有约 7 位有效数字，这些尾位被截断后，SVD/最近旋转投影会算出一个带方向偏置的微小假旋转。 更糟的是 float32 的 Cholesky 在求解时把这个偏置沿同一方向一帧帧累加（复现实测 $\|A^{-1}\mathbf 1\|$ 的 f32 与 f64 之差最大到 1.144，是灾难级的），于是布稳定地往一边漂。 修法：网格常数与 PD 内部运算一律用 float64，只有对外存储用 float32。 教训刻在脑子里：精度 bug 常常伪装成物理 bug——别一看到"布拐弯"就去找力的方向，先怀疑数值。

Bug 二 · "冰面"：摩擦系数从 1 调到 100 毫无变化

另一次，把 $\mu$ 从 1 一路调到 100，布的行为纹丝不动——摩擦像完全没生效，地面滑得像冰。 原因不在摩擦代码，而是一个量级失配：布被设得太硬，恢复弹性力约 0.84 N， 而那一点接触能提供的摩擦上限（$\mu\lambda_n$）只有约 0.002 N 量级——弹性力比摩擦上限大了约 400 倍。 结果是布刚一沾地，巨大的弹性回弹力立刻把它从地面顶了起来，于是 $\text{gap}\gt0$、 按 §9 的 Signorini 立刻 $\lambda_n=0$ ——没有接触就没有法向力，摩擦锥半径 $\mu\lambda_n=0$，$\mu$ 取多大都是乘以零。 教训：当一个旋钮在好几个数量级上都毫无反应，多半是那个变量压根没进有效代码路径 （这里是"根本没接触"），再怎么调它都没用。

§18G2D：擦地板而不是关水管

这是整个复现里最诚实的一段，也是理解"论文级算法落地有多难"的最佳样本。 §16 把地面场景修好了，但悬空抓着的布（airborne）还是会永远微微扭动 —— 因为位置投影式接触每步仍在注入能量，地面场景靠桌面摩擦把它耗掉了，悬空场景没东西耗。

关键区别：地面有"水池"接着，悬空没有

这里要看清一个能量账。§16 把碰撞检测挪到每步一次后，注入大大减小但没有归零—— 位置投影式接触（只做局部 $\Delta q$、没做 §12 那个全局解）每步仍漏一点能量进系统。 在地面场景，这点漏进来的能量会被桌面的 Coulomb 摩擦不断耗掉， 收支大致平衡，布就静下来了。但把布拎到半空（两个角被夹爪钉住）， 布对布的自接触摩擦极弱、又没有桌面——有进账、没出账，于是那点能量一直攒着， 布就永远微微扭动（airborne wriggle）。所以问题从来不是"悬空时注入更多"，而是"悬空时没有水池接住漏水"。

复现里最终上线的"修法"叫 G2D，代码只有两行：

# flash/pd_solver.py —— G2.D：对"正在自接触"的顶点额外做速度衰减
# ⚠ 这是一个 workaround，不是 paper-faithful 的修法（见 notes/G2D_CAVEATS.md）
in_contact = sc._contact_count > 0
v_new[in_contact] *= (1 - contact_damping)   # 默认 0.05

逐行读，重点是它不碰物理、只擦能量：

• sc._contact_count > 0：标出这一步处于自接触状态的顶点（接触计数大于 0）。 自由下落、没碰到任何东西的顶点不在内。
• v_new[in_contact] *= (1 - contact_damping)：只把这些顶点的速度乘一个略小于 1 的数 （0.95）。这等于人为往接触处加一点黏滞阻力，把 bug 注入的多余能量"擦掉"。
• 注意它没改任何约束、没解任何方程——纯粹是事后给速度打个折。 所以它能压住症状，却没碰到病根（接触投影仍在注入能量）。

它对正在自接触的顶点多做一次速度衰减。效果立竿见影：悬空丝绸 30 秒，平均速度从 90 mm/s 降到 4.4 mm/s，视觉上完全静止。复现里还扫了一遍这个系数（越大越静）： cd=0→90 mm/s、cd=0.02→28、cd=0.05→10、cd=0.20→2，最终选了 0.05。 但它错在哪，复现笔记说得毫不留情：

• 它没有修复 bug，只是加了一个能量吸收槽去抵消另一处 bug 注入的能量 —— 复现笔记的原话是"擦地板，而不是关掉漏水的管子"。
• 那个 0.05 是调出来的，不是物理推导的；换个材料、换个网格分辨率、换个步长都得重调。
• 真实的布对布接触靠 Coulomb 摩擦（切向）耗能，而 G2D 是各向同性黏滞阻尼 —— 连法向速度都一起衰减了，物理上是错的。

三个"它是 hack"的铁证

(1) 它擦地板，不关水管。 G2D 一点没减少接触投影的能量注入——它只是在出口处加了一个 人造的"能量水池"把漏进来的水舀走。病根（缺 §12 的全局解）还在，只是症状被盖住了。

(2) 那个 0.05 的调参曲线是非单调的，这本身就是混沌的告密。 扫系数时：cd=0.02→28 mm/s，cd=0.03→43，cd=0.05→10——注意 0.03 反而比 0.02 更糟。 一个有原理的阻尼/正则项（比如 §13 里 $Z$ 的那个 $E$）应当是单调、可预测的； 而这里"多阻一点反而更动"，说明背后是个混沌/双稳态系统（正是 §19 那回事）， G2D 只是在一片乱流里碰运气找了个还行的常数。

(3) 它过阻了法向、还和已有旋钮打架。 真实 Coulomb 摩擦只作用在切向 （§10），而 v *= 0.95 是各向同性的——连"垂直于接触面"的法向速度也一起衰减了。 后果是慢慢拖动布时会感觉"发黏、发僵"。更乱的是它和 PDConfig 里本来就有的全局 damping 旋钮作用域重叠——两个阻尼旋钮、不同范围，复现笔记自己都嫌这设计"丑"。

正解长什么样、为什么搁置

论文 faithful 的修法（复现代号 G2.C）是：在每次 PD 内层迭代里真正解那个鞍点系统—— 用 §12 的 Schur 补解出接触力 $\Delta\lambda$、再用 §15 的 cuDSS 做 §7 的全局 $q$ 更新、切向上 $\Delta\lambda_t$ 投影到 §10 的 Coulomb 锥。这样 $q$ 和 $\lambda$ 一致地一起解出来， 接触不再"推出去-弹回来"地漏能量，压根没有多余能量需要去擦——G2D 就该彻底消失。

那为什么没上？卡在 §17 那个病态的 $Z$：密集自接触下 $Z$ 极度病态，CG 解出乱真 $\Delta\lambda$ 把布炸到两万多 mm/s，Jacobi 预条件也救不回来。要救它得上 active-set（每迭代把 $Z$ 缩到只含活跃约束的良态子空间）， 工程量大、复现里没打通，于是暂时禁用（g2c_enabled = False），用 G2D 顶着。 所以 G2D 不是"懒"，是它前面横着一道还没解决的数值难关（§17 的开放问题）。

§19当浮点遇上混沌：同种子，不同物理

最后一个故事，对机器学习读者尤其值得警惕。我们发现：同一个随机种子、同一份代码、 同一套配置，跑三次会得到宏观上完全不同的结果：

# seed=42，同一份代码，三次独立进程调用：
trial 0:  |v|_last = 507 mm/s   # 布在"沸腾"
trial 1:  |v|_last =  17 mm/s   # 布已静止
trial 2:  |v|_last = 478 mm/s   # 布在"沸腾"

根因是两件事相乘：一个产生纳米级的微小差异，另一个把它放大成宏观分岔。拆开看。

第一步：浮点加法不满足结合律

数学上 $a+(b+c)=(a+b)+c$，但浮点数不是——因为每次加法都要把结果舍入到有限位数， 舍入顺序不同，结果就不同。一个 float32（约 7 位有效数字）的极端例子：把三个数 $\{10^8,\ 1,\ -10^8\}$ 按不同顺序相加：

# float32 下，同样三个数、不同求和顺序：
(1e8 + 1) - 1e8   # = 0 ：1e8+1 舍入回 1e8（那个 1 比 7 位精度还小，被吞了）→ 再减得 0
(1e8 - 1e8) + 1   # = 1 ：先抵消成 0，再 +1 → 得 1

同样的三个数，仅仅换了加的顺序，答案就从 0 变成 1 —— 100% 的相对误差，纯粹来自舍入。 这就是"加法不结合"在作怪：当数有大有小，谁先加、谁后加，决定了哪些小量被吞掉。

第二步：GPU 的 atomicAdd 让"加的顺序"每次都随机

仿真里"把每个接触/三角形对某顶点的贡献累加起来"这一步，在 GPU 上是几千个线程同时往 同一个内存地址做 atomicAdd。原子操作只保证"不会两个线程同时写、把彼此覆盖掉"， 但完全不保证它们到达的先后顺序——这个顺序取决于当次运行里线程调度的细微时序， 每次跑都可能不同。于是上一段那个"求和顺序"每次都在变，每次累加出的最后几位都略有出入。 代码一模一样，结果却带着一点不可控的随机尾巴。

第三步：混沌 + NCP 开关把纳米差异放大成宏观分岔

纳米级的差异本来无关紧要——除非系统对初值极度敏感。布料接触正是这样一个系统， 而且它带着两个天然的"放大器"，都在前面见过：

• NCP 是个离散开关（§16）：一个顶点位置差 1 个 ulp，就可能把它从"刚好不接触"翻到"刚好接触"， 于是注入 / 不注入一份冲量——一个连续的纳米差被放大成一个离散的"有/无"跳变。
• 病态的 Z 是个误差放大器（§17）：近奇异的接触度量会把输入里的小扰动放大 $\kappa$ 倍。

混沌系统里，微小差异随时间指数级拉开（蝴蝶效应）。布料-接触系统恰好有两个稳定结局 （"attractor"）：settle 盆地（静下来）和 boil 盆地（持续沸腾）。 起点上那点 atomicAdd 带来的纳米噪声，被开关和病态反复放大，最终决定这次轨迹掉进哪个盆地—— 这就是为什么同种子能跑出 17（静止）和 507（沸腾）两种宏观物理。双稳态分岔。

怎么治？复现里没做，但方向清楚：用排序后的确定性 scatter（先按目标顶点把所有贡献用 CUB 基数排序排好，再分段规约）固定求和顺序，结果就逐位可复现；代价是排序比裸 atomicAdd 慢。 退一步也可以升到 float64（缩小但不消除——混沌照样放大）、或干脆接受方差（很多机器人 RL 就是这么做的）、 或加一个"确定性模式"开关按需切换。这是一道速度与可复现性之间的取舍题，没有免费的答案。

§20GPU 并行仿真做 RL：FLASH 的舞台

FLASH 并非凭空出现。它继承的是过去五年"把上千个仿真环境塞进一块 GPU、 让 RL 疯狂采样"这套范式 —— 只不过把它从刚体推向了难得多的可变形。

Isaac Gym / Isaac Lab — GPU 上的大规模并行机器人学习仿真

2021 · 2025 Makoviychuk 等 · NVIDIA · arXiv:2108.10470 · Isaac Lab 2511.04831

CPU 仿真喂不饱 RL 的采样需求，训练以"周"计。

关键想法：物理和策略网络都放在 GPU 上，同时跑成千上万个环境实例， 把 RL 采样吞吐提升一到两个数量级。这是 FLASH 直接继承的范式 —— 但 Isaac 这套主要为刚体/关节体设计，碰到可变形就吃力，FLASH 要填的正是这块。

vs CPU 多进程仿真：把"几十个 CPU 环境"换成"单卡上千 GPU 环境"。

§21可变形仿真引擎全景

FLASH 在和谁赛跑？下面是 2024–2026 这条赛道的地图。一句话总结这张图的张力： 跑得快的（MJX/Brax/Isaac）大多是刚体优先；能算可变形的（FleX/Warp/DiffTaichi/Genesis） 又往往在精度、接触保真或并行吞吐上各有取舍。FLASH 想同时拿下"高保真 + RL 级并行 + 快"。

Warp · FleX · DiffTaichi — GPU 物理 / 可微仿真框架

框架 / SIGGRAPH / ICLR NVIDIA ／ Hu 等 · Warp · FleX 2014 · DiffTaichi 1910.00935

需要能算布/软体、又能上 GPU、最好还可微的底座。

关键想法：FleX 用统一粒子 + PBD 算布/流体/软体（经典 GPU 可变形）； Warp 是 Python→GPU 的可微核框架，含布料 VBD 求解器； DiffTaichi 用源码变换自动微分驱动一大批可微软体仿真。FLASH 与它们的区别在于走的是 严格 NCP（精度优先）这一路，而非 PBD（速度优先、软）。

vs PBD 系：把"软、近似、好并行"换成"严格互补、精确，再想办法保住并行"。

Genesis — "生成式"通用物理平台（2024）

2024-12 · 项目 开源项目（无可验证的同行评审论文） · project

关键想法：把多种求解器（含可微 MPM）、Python 接口、照片级渲染、 语言驱动数据生成统一进一个平台。但要谨慎：它发布时"比 Isaac Gym 快 10–80×" 的宣传受到广泛质疑（对比的是单线程 MuJoCo、基准场景默认关自碰撞）， 后被团队修订。引用它时应注明此争议。

vs 单一求解器引擎：野心是"全都要"，代价是基准与稳定性的可信度问题 —— 这恰是 FLASH 自己的速度数字也需独立复现的警示。

SoftGym · GarmentLab — 可变形 / 衣物操作的 RL 基准

CoRL 2020 · NeurIPS 2024 Lin 等 ／ GarmentLab 团队 · SoftGym 2011.07215 · GarmentLab 2411.01200

研究布料操作得有标准任务和环境（叠布、挂衣、铺平）。

关键想法：SoftGym（建在 FleX 上）定义了绳/布/流体操作的基准词汇； GarmentLab（建在 Isaac Sim 上，PBD + FEM，100+ 件 3D 衣物）是最贴近 FLASH 的同代对手。 FLASH "Isaac Sim 处理可变形吃力"的论断，实质就是对 GarmentLab 底座的挑战。

vs 早期单任务 demo：把"零散演示"换成"统一基准 + 大量 3D 衣物资产"。

IPC · StiffGIPC — GPU 上的高保真接触（增量势能接触）

2024 · 非距离势垒布料 2403.19272 · StiffGIPC 2411.06224

关键想法：图形学这边的高保真前沿 —— IPC 系用对数势垒严格保证不穿透， GPU 优化后比此前快 1–2 个数量级。注意对比：IPC 这类把一个高分辨率场景做到极致， 而 RL 世界（Isaac）是把上千个中等场景并行。FLASH 想跨在两者之间 —— 精度优先的求解器，RL 级别的并行。

vs PBD：把"软接触、可能穿透"换成"严格势垒、保证不穿透"，代价是更贵。

§22sim-to-real 与师生蒸馏：为什么"准+快+并行"缺一不可

把所有线索收束到一个问题：为什么 FLASH 非要同时做到"高保真 + 快 + 并行"？ 因为零样本 sim-to-real 折布同时要求两件历来互相矛盾的事 —— (a) 仿真要足够准，接触/形变行为得贴近现实； (b) 要足够快且能大批量并行，才能做足够多的域随机化（domain randomization） 和大规模数据生成，让策略对"仿真和现实的差异"鲁棒。

Sim-to-Real DOM — Sim-to-Real RL for Deformable Object Manipulation

CoRL 2018 Matas, James, Davison · ICL · arXiv:1806.07851

真实布料数据贵且难标，能不能纯靠仿真训练再迁移？

关键想法：用域随机化在仿真里训练布料操作策略，不看任何真实演示 就迁移到真机。这正是 FLASH 零样本折布所属的传统 —— 而它的卖点是：用更高保真 + 更大规模的仿真， 把这条路推得更稳。师生蒸馏（teacher 用特权信息训练、student 在域随机化图像上模仿）是常配的训练范式。

vs 真机收集数据：把"昂贵的真实演示"换成"海量随机化仿真 + 一次迁移"。

§A仍然没解决的问题

读到这里，你应该能看懂论文的每一块。但诚实地说，无论论文还是复现，都有没解决的硬骨头 —— 这些正是值得做的研究方向：

• 大规模自碰撞。论文实验里自碰撞是关掉的（声称为了跨平台公平比较）。 复现怀疑真实原因之一是自碰撞在密集接触下跑不稳（§17 的病态 + §19 的双稳态）。 在上千并行布料环境里做准确自碰撞，仍是公开难题。
• 病态 Z 的根治。冗余约束让接触度量 $Z$ 病态，普通预条件子救不了。 正解大概率是 active-set（每次迭代只保留真正穿透的约束，把 $Z$ 缩小到良态子空间）， 但复现里这条路还没打通。
• 精度 vs 速度的接触。PBD/XPBD（快、软、近似）与 IPC/NCP（准、硬、贵）的分裂仍在； FLASH 押注在"准的一侧 + 大规模并行"，能否长期站稳还需独立复现。
• 穿过接触的可微性。硬接触梯度不连续、软接触梯度不准，至今没有通解。 FLASH 走的是 RL（绕开可微），但这道题本身仍开放。
• 确定性 / 可复现。GPU 浮点非结合性 + 混沌让"同种子不同物理"（§19）， 对 RL 训练和 sim-to-real 都是隐患。
• 论文未说清的工程细节。bending 能量的具体形式、材料参数到 PD 权重的换算、 friction anchor 的语义、$J$ 的行数（是否含切向）—— 这些复现都得靠猜（见仓库 PAPER_QA.md）。

§B名词表 · Glossary

术语	解释
隐式 Euler	用"下一刻"的力做时间积分，无条件稳定，等价于每步最小化一个增量势能（§5）。
增量势能	$E(x)=\frac{1}{2h^2}\|x-\tilde x\|_M^2+\Psi(x)$；最小化它 = 走一步隐式 Euler。
形变梯度 F	把材料纤维从静止映射到当前的雅可比矩阵 $F=D D_r^{-1}$；编码拉伸 + 旋转（§6）。
极分解 / SVD	$F=RS$ 把旋转 $R$ 和拉伸 $S$ 分开；通过 $F=U\Sigma V^\top$ 得最近旋转 $R=UV^\top$。
ARAP	As-Rigid-As-Possible 弹性能 $\frac{w}{2}\|F-R\|^2$，只惩罚形变、不惩罚旋转（§6 §8）。
Projective Dynamics	把弹性能写成"到约束流形的距离"，用 local（并行投影）+ global（固定矩阵解）交替求解（§7）。
local / global 步	local：每单元投影到最近约束（并行）；global：解一个常数矩阵的线性系统。
预分解 Cholesky	$A=LL^\top$ 分解一次，之后每步只回代；PD/FLASH 速度的核心。
互补性 ⊥	两个非负量至少一个为零；$0\le\lambda\perp \text{gap}\ge0$ 即 Signorini（§9）。
Signorini 条件	不穿透的"非此即彼"：要么有间隙无力，要么贴住有力（§9）。
LCP / NCP	线性 / 非线性互补问题；接触（含摩擦锥）写成 NCP（§9）。
Coulomb 摩擦锥	切向力约束 $\|\lambda_t\|\le\mu\lambda_n$；锥内黏住、锥上滑动（§10）。
Fischer–Burmeister	$\phi_{FB}(a,b)=a+b-\sqrt{a^2+b^2}$，把互补条件写成（半光滑）方程（§11）。
非光滑 Newton	允许在不可导点用广义导数的 Newton，用来解含尖角的接触方程（§11）。
鞍点 / KKT 系统	位置 $q$ 与乘子 $\lambda$ 联立的不定线性系统（§11）。
Schur 补 Z	消掉 $q$ 后只关于 $\lambda$ 的约化系统 $Z=JA^{-1}J^\top+E$；又名 Delassus 算子（§12）。
惯性近似	FLASH 把 Schur 里的稠密 $A^{-1}$ 换成对角 $M^{-1}$，保住稀疏与块对角（§13，Eq.10）。
Green 函数	$A^{-1}$ 的物理意义："在 $j$ 戳一下 $i$ 动多少"；连通弹性网里它是稠密的（§13）。
块对角	多环境矩阵拼成对角块；其 Cholesky = 各块 Cholesky，可同时分解所有环境（§14）。
条件数 κ	矩阵"长宽比"；CG 迭代次数 $\sim\sqrt{\kappa}$，病态时暴涨（§15 §17）。
PCG 悬崖	接触一多 $\kappa$ 飙升，迭代法性能从悬崖跌落；直接分解（cuDSS）无此问题（§15）。
cuDSS	NVIDIA 的 GPU 稀疏直接求解器；把线性求解搬出关键路径（§15）。
自由度 DOF	未知数个数；布有 $N$ 个顶点即 $3N$ 个 DOF。
域随机化	在仿真里随机化物理/视觉参数，让策略对 sim-to-real 差异鲁棒（§22）。

§C推荐阅读顺序

如果你想读原始文献，下面是一个按教学顺序（不是按年份）排的清单 —— 每一篇都为下一篇铺路：

1. Baraff & Witkin 1998《Large Steps in Cloth Simulation》 —— 先建立"为什么布料必须用隐式积分"的直觉，全场的起点。
2. Sifakis FEM 课程笔记（SIGGRAPH 2012） —— 把形变梯度 $F$、应变能、极分解一次讲透；§6 §8 的底层。
3. Sorkine & Alexa 2007《ARAP》 —— 最干净的"local 抽旋转、global 解线性"模板，PD 的原型。
4. Bouaziz 等 2014《Projective Dynamics》 —— FLASH 弹性引擎的直接来源；读懂"矩阵固定、预分解一次"。
5. Overby 等 2017《ADMM ⊇ PD》 —— 给 PD 一个优化理论的解释，知道 local-global 到底在最小化什么。
6. Stewart–Trinkle 1996 / Anitescu–Potra 1997 —— "接触 = 互补问题"的奠基，§9 的源头。
7. Macklin 等 2019《Non-Smooth Newton Methods》 —— FLASH 接触部分的直接父辈；NCP + 半光滑 Newton + 每步一个对称解。
8. Shewchuk《CG Without the Agonizing Pain》 —— 理解条件数、$\sqrt{\kappa}$ 与 PCG 悬崖（§15 §17）的最佳教程。
9. FLASH（arXiv:2604.17513）本体 —— 把以上全部装进 GPU、推到三百万自由度、并上千环境并行。最后读，水到渠成。
10. SIGGRAPH 2022 接触课程笔记 —— 想把接触/摩擦的现代全貌补齐时的最佳一站式参考。

§D参考文献

以下文献的题名、作者、出处均经核对。FLASH 本体的存在、作者与摘要级论断（>3M DOF @ 30 FPS on RTX 5090、 NCP 求解器、零样本 sim-to-real 折布）已在 arXiv 核实；其内部数值与 sim-to-real 结果为作者自报、尚待独立复现。 本文所有交互演示是 Canvas2D 教学玩具，数值不代表论文实际性能。

引用	出处	链接
FLASH（本文主题）	Luo, Zhou, Zhang 等，arXiv 2026	arXiv:2604.17513
Large Steps in Cloth Simulation	Baraff & Witkin，SIGGRAPH 1998	PDF
Projective Dynamics	Bouaziz 等，TOG/SIGGRAPH 2014	project
As-Rigid-As-Possible	Sorkine & Alexa，SGP 2007	PDF
Position Based Dynamics	Müller 等，JVCIR 2007	PDF
XPBD	Macklin, Müller, Chentanez，MIG 2016	PDF
Non-Smooth Newton Methods	Macklin 等，TOG 2019	arXiv:1907.04587
ADMM ⊇ Projective Dynamics	Overby 等，TVCG 2017	project
Quasi-Newton for Hyperelastic	Liu, Bouaziz, Kavan，TOG 2017	arXiv:1604.07378
Implicit Time-Stepping (LCP)	Stewart & Trinkle，IJNME 1996	DOI
Solvable LCP with Friction	Anitescu & Potra，Nonlinear Dyn. 1997	DOI
Exact Coulomb (Nonsmooth Newton)	Bertails-Descoubes 等，TOG 2011	DOI
FEM Simulation of 3D Deformables	Sifakis，SIGGRAPH 2012 课程	notes
CG Without the Agonizing Pain	Shewchuk，CMU 1994	PDF
接触/摩擦课程笔记	Andrews, Erleben，SIGGRAPH 2022	notes
Isaac Gym	Makoviychuk 等，2021	arXiv:2108.10470
DiffTaichi	Hu 等，ICLR 2020	arXiv:1910.00935
SoftGym	Lin 等，CoRL 2020	arXiv:2011.07215
GarmentLab	NeurIPS 2024	arXiv:2411.01200
Sim-to-Real Deformable RL	Matas 等，CoRL 2018	arXiv:1806.07851
cuDSS（GPU 稀疏直接求解）	NVIDIA	developer.nvidia.com/cudss