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

§1.1FLASH 到底是什么

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

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

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

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

§1.2一帧仿真在算什么

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

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

• $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$ 个单元（一个三角形）的弹性势能 —— 布被拉伸/弯折得越厉害，这一项越大（§2.2 细讲）。

逐块读懂那个"加权范数" $\|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$ 是对角的质量矩阵，所以这一项展开就是

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

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

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

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

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

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

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

核心循环 · 教学伪代码

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

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

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

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

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

§1.4四个痛点

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

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

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

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

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

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

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

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

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

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

§2.1隐式积分 = 优化

我们从最底层讲起：仿真就是反复回答"下一刻物体在哪"。牛顿第二定律 $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 用下一刻的位置算力：

区别看似只是 $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$ 求梯度（高数里的链式法则 + "二次型 $\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$ 的最小值点"是同一件事。

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

§2.2弹性能与形变梯度

§1.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$ 就抵消了：

把"当前两条边"并成一个矩阵 $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$， 于是

这就是形变梯度的闭式公式。注意 $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：静止面积。三角形越大、它的弹性能在总能量里权重越大（§2.3 会用到）。

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

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

于是奇异值 $\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$ 离一个纯旋转有多远"：

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

§2.3Projective Dynamics：local-global

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

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

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

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

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

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

为什么这快得离谱：矩阵 $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)                           # 局部步：投影到最近旋转（§2.4）
    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$）。这一步对每个三角形独立，天然并行（§2.4 细讲）。
• 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$ 次内层迭代依然便宜。

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

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

局部步的核心是一行字：把形变梯度 $F$ 投影到最近的旋转 $R$（这个 $R$ 就是 §2.3 里的投影变量 $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）。

§3.1互补性与 Signorini 条件

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

"非此即彼"：互补性

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

把这一行拆成三个条件读：$\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 法会卡住的根源（§3.3 专门治它）。如果 gap 和 $\lambda$ 对未知量是线性的，整个问题叫 线性互补问题（LCP）；加上摩擦的非线性锥后，叫非线性互补问题（NCP） —— 这就是 FLASH 里 "NCP" 的来历。论文的接触动力学（Eq.2）正是这个形式：

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

§3.2Coulomb 摩擦与摩擦锥

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

• 黏住（stick）：切向力严格在锥内 $\|\lambda_t\| \lt \mu\lambda_n$ —— 不滑动，摩擦力要多少给多少（把物体钉住）；
• 滑动（slip）：切向力达到锥边界 $\|\lambda_t\| = \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)：把法向力截到非负——这一行就是 §3.1 的 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$， 看方块到底是黏住还是滑动，以及当前接触力落在摩擦锥的哪里（空心点 = "想要的力"， 实心点 = 投影后的实际力）。

§3.3非光滑 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 条件——约束优化最优解必须满足的一组方程：

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

用"非光滑 Newton"对付尖角

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

花一分钟验证这个等价（令 $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$ 正好凑齐 §3.1 的 Signorini。于是"一个尖角不等式逻辑"被压成了"一个等式 $\phi_{\text{FB}}=0$"。

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

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

鞍点系统又大又不定（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$ 消掉，只剩 $\Delta\lambda$： $J A^{-1}(g+J^\top h^2\Delta\lambda) + E\,h^2\Delta\lambda = \bar h$，整理同类项：

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

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

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

终于到了 FLASH 真正的核心创新（论文 Eq.10）。前面 §3.4 的 $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}$：

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

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

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

"块对角"到底长什么样

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

每个 $A_i$ 是第 $i$ 个环境的 PD 矩阵（$N\times N$ 量级）。那些零块就是"环境之间没有任何耦合" 的数学体现——这正是 §4.1 那个惯性近似 $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)$，那么

（块对角矩阵相乘，就是对应块相乘——因为零块乘任何东西还是零。）而 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"一列：从 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 来说，这意味着用一块卡、一个进程，就能持续产出训练所需的海量交互 —— 而这一切的前提，是 §4.1 那个让矩阵保持稀疏 + 块对角的惯性近似。把链条连起来： 稀疏（§4.1）→ 块对角可独立分解（§4.2）→ 一个 $L$ 批量回代喂饱 GPU → 上千环境近乎免费。

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

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

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

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

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

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

每一帧（while 循环）

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

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

• p_bar = project(G_bar @ q) — 局部步（§2.4）：G_bar @ q 从当前位置抽出每个三角形的形变梯度 $F$，project 把它投到最近旋转 $R$（即投影变量 $p$）。每个三角形独立、完全并行。
• # 组装 J, E, g, h — 准备接触线性化的零件（§3.3）：$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 招牌（§4.1）：用对角 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) — 解接触力（§3.4 的 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) — 全局步（§2.3/§4.3）：$=\bar A^{-1}(\bar b + h^2\bar J^\top\bar\lambda)$，把惯性、弹性投影、接触力一起折中解出新位置。这一步只是对预分解的 $\bar A$ 回代一次（§4.3 讲它为什么不能落进 PCG 悬崖）。
• 循环回到 p_bar：local 和 global 交替 $k$ 次，一步步逼近这一帧的平衡（§2.3 Demo 3 里"最大位移越来越小"就是它）。

收尾

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

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

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

最后一块拼图：那个全局线性解 $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$ 步后满足：

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

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

三件事会把 κ 推高：材料越硬、时间步越小、尤其是冗余 / 近乎冲突的约束。 最后一条是接触的命门（§5.2 细讲）：密集自接触里，相邻接触共享顶点、法向近乎平行， 于是约束矩阵冒出一堆近乎重复的行——两行几乎线性相关，把最小特征值 $\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$ 里非零元的个数），只跟稀疏结构有关、跟条件数无关。 κ 再大、约束再冗余，分解和回代的运算量一个都不多——悬崖被直接填平。

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

§5.1能量从哪来：每步检测 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 秒还在动多快"和"动能有没有超过释放峰值"两个指标判读：

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

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

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

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

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

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

病态的 Z 怎么把布炸上天

接下来就是 §4.3 那条 $\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 只缩放对角，消不掉"两行几乎一样"（正应了 §4.3 的结论）。

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

根治办法是 active-set（活跃集）：那些"近乎重复"的约束其实是冗余的—— 保留其中一个就够了。每次迭代只挑出真正在穿透的接触（活跃集），把 $Z$ 缩小到只含这些行/列， 在这个满秩、良态的子空间里求解，κ 就回到正常范围。复现里这条（对应论文 Eq.8/9 的完整 Schur 解） 实现起来工程量大、还没打通，于是暂时禁用，退回到"位置投影 + 一个能量阻尼 workaround"—— 也就是下一节 §5.3 那个"擦地板"的 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$、 按 §3.1 的 Signorini 立刻 $\lambda_n=0$ ——没有接触就没有法向力，摩擦锥半径 $\mu\lambda_n=0$，$\mu$ 取多大都是乘以零。 教训：当一个旋钮在好几个数量级上都毫无反应，多半是那个变量压根没进有效代码路径 （这里是"根本没接触"），再怎么调它都没用。

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

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

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

这里要看清一个能量账。§5.1 把碰撞检测挪到每步一次后，注入大大减小但没有归零—— 位置投影式接触（只做局部 $\Delta q$、没做 §3.4 那个全局解）每步仍漏一点能量进系统。 在地面场景，这点漏进来的能量会被桌面的 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 一点没减少接触投影的能量注入——它只是在出口处加了一个 人造的"能量水池"把漏进来的水舀走。病根（缺 §3.4 的全局解）还在，只是症状被盖住了。

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

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

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

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

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

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

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

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

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

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

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

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

§6.2可变形仿真引擎全景

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

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

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

§A仍然没解决的问题

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

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

§B名词表 · Glossary