随机微分方程 (SDE)
一般的微分方程可以表示为：
$$\frac{d\textbf{x}}{dt} = f(\textbf{x},t)$$几乎所有更高阶的微分方程都可以转化为这种一阶形式。例如，一个二阶系统（如牛顿第二定律 $F=ma$）可以通过引入辅助变量（令速度 $v = \dot{x}$），将其改写为两个一阶方程组成的向量形式：
$$\frac{d}{dt} \begin{bmatrix} x \\ v \end{bmatrix} = \begin{bmatrix} v \\ \frac{F}{m} \end{bmatrix}$$然而，在许多现实世界的系统中，存在随机性或不确定性，这时我们需要引入随机微分方程 (SDE) 来描述这些系统。
随机微分方程在传统微分方程的基础上，加入了一个随机过程项，通常表示为Wiener过程。一个典型的随机微分方程形式为：
$$d\textbf{x} = f(\textbf{x},t)dt + g(\textbf{x},t)d\textbf{W}_t$$其中$f(\mathbf{x}_t, t)$ 是漂移项。$g(t)$ 是扩散项。$\text{d}\mathbf{W}_t$是标准维纳过程。 例如，在slam系统中， imu传感器的bias可以建模为一个随机游走过程：
$$d\textbf{b} = \sigma d\textbf{W}_t$$其离散化形式为：
$$ \mathbf{b}_{k+1} = \mathbf{b}_k + \mathbf{W}_{b,k}, \quad \mathbf{W}_{b,k} \sim \mathcal{N}(0, \Delta t \sigma_{b}^2 ) \\ $$该式描述了bias的方差在每个时间步长 $\Delta t$ 上的变化， 随时间步长线性增长，$\sigma$ 控制了随机游走的强度， 这意味着，随着时间推移，不确定性的范围以时间平方根的速率扩散。
推导福克-普朗克方程 (FPE)
随机微分方程 (SDE) 描述了系统状态的随机演化，而福克-普朗克方程 (FPE) 则描述了系统状态概率密度函数 (PDF) 随时间的变化。通过从 SDE 推导 FPE，我们可以理解系统状态的统计行为。 维纳过程具有以下核心统计特性， 增量 $dW_t = W_{t+dt} - W_t$ 服从正态分布：
$$dW_t \sim \mathcal{N}(0, dt)$$$$\mathbb{E}[dW_t] = 0$$$$\mathbb{E}[(dW_t)^2] = dt$$$dW_t$ 的平方的期望值是 $dt$，而不是 $O(dt^2)$。在极限 $dt \to 0$ 的分析中，$dt$ 级别的项必须保留。
1. 伊藤乘法表（Itô Multiplication Table）
伊藤乘法表列出了在计算 SDE 的乘积 $d(X_t Y_t)$ 时，应该保留哪些二阶项。其核心逻辑是基于 $dW_t \cdot dW_t = dt$。
| $dt$ | $dW_t$ | |
|---|---|---|
| $dt$ | $0$ | $0$ |
| $dW_t$ | $0$ | $dt$ |
2. 伊藤引理（Itô’s Lemma）
假设我们有一个函数 $f(t, X_t)$，其中 $X_t$ 遵循SDE：
$$\text{d}X_t = a \text{d}t + b \text{d}W_t$$对 $f$ 进行泰勒展开：
$$\text{d}f = \frac{\partial f}{\partial t} \text{d}t + \frac{\partial f}{\partial X} \text{d}X_t + \frac{1}{2} \frac{\partial^2 f}{\partial X^2} (\text{d}X_t)^2 + \dots$$现在，我们利用乘法表计算 $(\text{d}X_t)^2$ 的二阶变差项：
$$\begin{aligned}(\text{d}X_t)^2 &= (a \text{d}t + b \text{d}W_t)^2 \\ &= a^2 (\text{d}t)^2 + 2ab (\text{d}t \text{d}W_t) + b^2 (\text{d}W_t)^2 \\ &\stackrel{\text{Itô Rules}}{=} 0 + 0 + b^2 (\text{d}t) \end{aligned}$$将 $(\text{d}X_t)^2 = b^2 \text{d}t$ 代回泰勒展开式，最终得到伊藤引理：
$$\text{d}f = \left(\frac{\partial f}{\partial t} + a \frac{\partial f}{\partial X} + \frac{1}{2} b^2 \frac{\partial^2 f}{\partial X^2} \right) \text{d}t + b \frac{\partial f}{\partial X} \text{d}W_t$$3. 从 SDE 推导 FPE(福克-普朗克方程)
对于任何有界且二次连续可微的测试函数 $f(X_t)$，我们希望其期望值 $\mathbf{E}[f(X_t)]$ 随时间的演化率,
$$\mathbf{E}[f(X_t)] = \int_{-\infty}^{\infty} f(x) p(x, t) \text{d}x $$对时间 $t$ 求导：
$$\frac{\text{d}}{\text{d}t} \mathbf{E}[f(X_t)] = \frac{\text{d}}{\text{d}t} \int f(x) p(x, t) \text{d}x = \int f(x) \frac{\partial p(x, t)}{\partial t} \text{d}x \tag{1}$$根据 SDE 和伊藤引理，我们可以计算$\text{d}f(X_t)$：
$$\text{d}f(X_t) = \left[ a \frac{\partial f}{\partial X} + \frac{1}{2} b^2 \frac{\partial^2 f}{\partial X^2} \right] \text{d}t + \left[ b \frac{\partial f}{\partial X} \right] \text{d}W_t$$计算期望值的微分我们对 $\text{d}f(X_t)$ 求期望 $\mathbf{E}[\dots]$：
$$\mathbf{E}[\text{d}f(X_t)] = \mathbf{E} \left[ \left( a \frac{\partial f}{\partial X} + \frac{1}{2} b^2 \frac{\partial^2 f}{\partial X^2} \right) \text{d}t + \left( b \frac{\partial f}{\partial X} \right) \text{d}W_t \right]$$根据 $\mathbf{E}[\text{d}W_t] = 0$，随机项的期望为零：
$$\mathbf{E}[\text{d}f(X_t)] = \mathbf{E} \left[ a \frac{\partial f}{\partial X} + \frac{1}{2} b^2 \frac{\partial^2 f}{\partial X^2} \right] \text{d}t$$因此，期望值随时间的演化率是：
$$ \begin{aligned} \frac{\text{d}}{\text{d}t} \mathbf{E}[f(X_t)] &= \mathbf{E} \left[ a \frac{\partial f}{\partial X} + \frac{1}{2} b^2 \frac{\partial^2 f}{\partial X^2} \right] \\ &= \int \left( a(x, t) \frac{\partial f}{\partial x} + \frac{1}{2} b^2(x, t) \frac{\partial^2 f}{\partial x^2} \right) p(x, t) \text{d}x \tag{2} \end{aligned} $$现在我们有两个相等的 $(1)$ 和 $(2)$ 来表示 $\frac{\text{d}}{\text{d}t} \mathbf{E}[f(X_t)]$：
$$\int f(x) \frac{\partial p}{\partial t} \text{d}x = \int \left( a \frac{\partial f}{\partial x} p + \frac{1}{2} b^2 \frac{\partial^2 f}{\partial x^2} p \right) \text{d}x$$我们的目标是让右边的表达式也包含 $f(x)$，而不是其导数 $\frac{\partial f}{\partial x}$ 和 $\frac{\partial^2 f}{\partial x^2}$。我们通过分部积分（假设 $f$ 及其导数在边界处为零）：
- 项 1: 漂移项
使用分部积分公式 $\int u \text{d}v = u v | - \int v \text{d}u$：
$$\int a \frac{\partial f}{\partial x} p \text{d}x = \left[ (a p) f(x) \right]_{L}^{R} - \int_{L}^{R} f(x) \frac{\partial}{\partial x} (a p) \text{d}x$$假设$f(x)$ 在边界 $L \to -\infty$ 和 $R \to \infty$ 处收敛到 $0$, 概率密度 $p(x, t)$ 在边界处也收敛到 $0$：
$$\int (a p) \frac{\partial f}{\partial x} \text{d}x = -\int f(x) \frac{\partial}{\partial x} (a p) \text{d}x $$- 项 2: 扩散项 $$\int \frac{1}{2} b^2 \frac{\partial^2 f}{\partial x^2} p \text{d}x = \int \frac{1}{2} (b^2 p) \frac{\partial^2 f}{\partial x^2} \text{d}x$$
需要两次分部积分：
第一次分部积分消去一个 $\frac{\partial}{\partial x}$,
$$-\int \frac{\partial f}{\partial x} \frac{\partial}{\partial x} (\frac{1}{2} b^2 p) \text{d}x$$第二次分部积分消去另一个 $\frac{\partial}{\partial x}$：
$$+\int f(x) \frac{\partial^2}{\partial x^2} (\frac{1}{2} b^2 p) \text{d}x$$- 最终 FPE将分部积分的结果代回等式： $$\int f(x) \frac{\partial p}{\partial t} \text{d}x = \int f(x) \left[ -\frac{\partial}{\partial x} (a p) + \frac{1}{2} \frac{\partial^2}{\partial x^2} (b^2 p) \right] \text{d}x$$
由于这个等式必须对任意测试函数 $f(x)$ 都成立, 被积函数必须相等。
$$\frac{\partial p}{\partial t} = -\frac{\partial}{\partial x} [a(x, t) p] + \frac{1}{2} \frac{\partial^2}{\partial x^2} [b^2(x, t) p]$$这就是福克-普朗克方程（FPE）的形式，它描述了概率密度函数 $p(x, t)$ 随时间的演化，基于随机微分方程 (SDE) 中的漂移项 $a(x, t)$ 和扩散项 $b(x, t)$。
从FPE到FPE -Reverse
修改逆向 FPE 的扩散项符号表示， 我们假设前向过程 SDE 为：
$$\text{d}\mathbf{x}_t = \mathbf{f}(\mathbf{x}_t, t) \text{d}t + g(t) \text{d}\mathbf{w}_t$$对应的前向 FPE 描述了在正向时间 $t$ 增加时，概率密度 $p(\mathbf{x}, t)$ 是如何演化的：
$$\frac{\partial p}{\partial t} = -\nabla \cdot [\mathbf{f} p] + \frac{1}{2} g^2 \Delta p \quad \text{(FPE-Forward)}$$这里利用了散度算子 $\nabla \cdot$ 和拉普拉斯算子 $\Delta$ 来表示多维空间中的漂移和扩散。 其中，梯度（Gradient） $\nabla$：
$$\nabla p_t = \left( \frac{\partial p_t}{\partial x_1}, \frac{\partial p_t}{\partial x_2}, \frac{\partial p_t}{\partial x_3} \right)$$散度（Divergence） $\nabla \cdot (\text{矢量场})$
$$\nabla \cdot (\nabla p_t) = \nabla \cdot \left( \frac{\partial p_t}{\partial x_1}, \frac{\partial p_t}{\partial x_2}, \frac{\partial p_t}{\partial x_3} \right)$$$$\nabla \cdot (\nabla p_t) = \frac{\partial}{\partial x_1} \left( \frac{\partial p_t}{\partial x_1} \right) + \frac{\partial}{\partial x_2} \left( \frac{\partial p_t}{\partial x_2} \right) + \frac{\partial}{\partial x_3} \left( \frac{\partial p_t}{\partial x_3} \right)$$拉普拉斯算子 $\Delta$：
$$\Delta p_t \equiv \nabla \cdot (\nabla p_t)$$$$\Delta p_t = \frac{\partial^2 p_t}{\partial x_1^2} + \frac{\partial^2 p_t}{\partial x_2^2} + \frac{\partial^2 p_t}{\partial x_3^2}$$现在，我们考虑反向过程。
- 定义时间反转（Reverse Time）我们定义反向时间变量 $\tau$。
当 $t$ 从 $0 \to T$ 时（前向），$\tau$ 从 $T \to 0$。当 $t$ 从 $T \to 0$ 时（反向），$\tau$ 从 $0 \to T$。我们需要描述概率密度 $p(\mathbf{x}, t)$ 随正向时间 $t$ 倒退时的变化率 $\frac{\partial p}{\partial t}$。
- 利用链式法则关联时间导数使用链式法则，我们关联对 $t$ 的偏导数和对 $\tau$ 的偏导数：$$\frac{\partial p}{\partial t} = \frac{\partial p}{\partial \tau} \cdot \frac{\partial \tau}{\partial t}$$因为 $\tau = T - t$，所以 $\frac{\partial \tau}{\partial t} = -1$。
- 应用前向 FPE我们知道 $p(\mathbf{x}, t)$ 始终满足前向 FPE 的结构，无论我们是从 $t=0$ 运行到 $T$，还是从 $T$ 运行到 $0$。它只是一个描述当前分布 $p$ 如何被 $\mathbf{f}$ 和 $g$ 驱动的微分方程。现在，我们将前向 FPE 中的 $\frac{\partial p}{\partial t}$ 替换为 $-\frac{\partial p}{\partial \tau}$：
- 得到反向时间的 FPE 结构将整个方程乘以 $-1$，得到概率密度随反向时间 $\tau$ 增加时的演化方程： $$\frac{\partial p}{\partial \tau} = \nabla \cdot [\mathbf{f} p] - \frac{1}{2} g^2 \Delta p$$
这个方程就是描述反向过程（时间 $\tau$ 增加的福克-普朗克方程，我们称之为 FPE-Reverse）。
识别逆向扩散项现在来看 FPE-Reverse 的各项：
$$\frac{\partial p}{\partial \tau} = \underbrace{\nabla \cdot [\mathbf{f} p]}_{\text{逆向漂移和前向漂移的组合}} \underbrace{- \frac{1}{2} g^2 \Delta p}_{\text{逆向扩散项}} \tag{3} $$我们看到，当系统沿着反向时间 $\tau$ 演化时，扩散项的符号翻转成了负号：
$$\mathbf{J}_{\text{逆向扩散}} \propto - \frac{1}{2} g^2 \Delta p$$逆向 SDE 的形式
逆向 SDE 的目标形式: 我们希望找到一个 SDE，使得当 $\tau$ 增加时，其概率密度的演化与前向过程在 $t$ 减小时一致。设其形式为：
$$\text{d}X_\tau = \bar{\mathbf{f}}(X_\tau, \tau) \text{d}\tau + g(\tau) \text{d}\bar{\mathbf{W}}_\tau$$为了凑出标准 FPE 形式，我们在 (3) 右边同时加上和减去一个扩散项 $\frac{1}{2} g^2 \nabla^2 p_t$：
$$\frac{\partial p_\tau}{\partial \tau} = \nabla \cdot (\mathbf{f} p_t) - \frac{1}{2} g^2 \nabla^2 p_t + \underbrace{\left( \frac{1}{2} g^2 \nabla^2 p_t - \frac{1}{2} g^2 \nabla^2 p_t \right)}_{=0}$$整理各项，将一个正的扩散项留在最后：
$$\frac{\partial p_\tau}{\partial \tau} = \nabla \cdot (\mathbf{f} p_t) - g^2 \nabla^2 p_t + \frac{1}{2} g^2 \nabla^2 p_t \quad \tag{4}$$现在我们需要处理中间的 $- g^2 \nabla^2 p_t$ 项。 利用散度定理和恒等式 $\nabla^2 p = \nabla \cdot (\nabla p)$：
$$\nabla \cdot (\mathbf{f} p_t) - g^2 \nabla \cdot (\nabla p_t) = -\nabla \cdot \left[ -\mathbf{f} p_t + g^2 \nabla p_t \right]$$利用 分数函数 $\nabla \log p_t = \frac{\nabla p_t}{p_t} \implies \nabla p_t = p_t \nabla \log p_t$，代入上式：
$$\dots = -\nabla \cdot \left[ (-\mathbf{f} + g^2 \nabla \log p_t) p_t \right]$$将此结果代回 (式 4)，我们得到了 符合标准形式的 FPE-Reverse：
$$\frac{\partial p_\tau}{\partial \tau} = -\nabla \cdot \underbrace{\left[ -\mathbf{f} + g^2 \nabla \log p_t \right]}_{\text{逆向漂移项 } \bar{\mathbf{f}}} p_\tau + \frac{1}{2} g^2 \nabla^2 p_\tau$$根据标准 FPE:
$$\frac{\partial p}{\partial \tau} = -\nabla \cdot (\bar{\mathbf{f}} p) + \frac{1}{2} \bar{g}^2 \nabla^2 p$$与对应 SDE:
$$\text{d}X_\tau = \bar{\mathbf{f}} \text{d}\tau + \bar{g} \text{d}\bar{\mathbf{W}}_\tau$$的对应关系, 我们可以直接读出系数：
逆向扩散系数:
$$\bar{g} = g$$逆向漂移项:
$$\bar{\mathbf{f}} = -\mathbf{f}(x, T-\tau) + g^2(T-\tau) \nabla \log p_{T-\tau}(x)$$因此，逆向 SDE (Reverse SDE) 的形式为：
$$\text{d}X_\tau = \underbrace{\left[ -\mathbf{f}(X_\tau, T-\tau) + g^2(T-\tau) \nabla \log p_{T-\tau}(X_\tau) \right]}_{\text{确定性收拢力}} \text{d}\tau + g(T-\tau) \text{d}\bar{\mathbf{W}}_\tau$$最后将时间变量从 $\tau$（从 $0 \to T$）切换回 $t$（从 $T \to 0$），这个过程非常关键，因为在实际编写代码（如使用 Euler-Maruyama 采样器）时，我们通常习惯于在 $t$ 轴上进行递减迭代。
我们将 $\text{d}\tau = -\text{d}t$ 逆向 SDE 公式中：
$$\text{d}\mathbf{x} = \left[ -\mathbf{f}(\mathbf{x}, t) + g^2(t) \nabla \log p_t(\mathbf{x}) \right] (-\text{d}t) + g(t) \text{d}\bar{\mathbf{W}}$$这里的 $\text{d}\bar{\mathbf{W}}$ 是反向维纳过程。为了让公式看起来更直观，我们将负号 $(-\text{d}t)$ 分配进方括号内：
$$\text{d}\mathbf{x} = \left[ \mathbf{f}(\mathbf{x}, t) - g^2(t) \nabla \log p_t(\mathbf{x}) \right] \text{d}t + g(t) \text{d}\bar{\mathbf{W}}$$在这里，$\text{d}t$ 通常被理解为一个负的时间步长（例如从 $t$ 变到 $t - \Delta t$）。
最终形式以 $t$ 为变量的逆向为 Yang Song et al. 在其论文1给出的标准形式：
$$\text{d}\mathbf{x} = \left[ \mathbf{f}(\mathbf{x}, t) - g^2(t) \nabla_{\mathbf{x}} \log p_t(\mathbf{x}) \right] \text{d}t + g(t) \text{d}\bar{\mathbf{W}}$$其中 $t$ 从 $T$ 运行到 $0$。
# t 从 T 逐渐减小到 0， 令 dt = 0.001 (正数)
for t in reversed(linspace(0, T, steps)):
score = model(x, t)
# 计算漂移项: (f - g^2 * score)
# 注意因为是倒着走，更新时用的是减法或负的 dt
drift = (f(x, t) - g(t)**2 * score) * (-dt)
diffusion = g(t) * sqrt(dt) * noise
x = x + drift + diffusion