扩散模型的 SDE/ODE 统一：随机微分方程到确定性采样

一、从离散到连续：马尔可夫链的 SDE 极限

DDPM 的前向过程是一个离散的马尔可夫链：每一步从 $x_{t-1}$ 到 $x_t$ 添加一个小的高斯噪声。这个离散过程有一个自然的连续极限——当时间步长 $\Delta t \to 0$ 时，马尔可夫链收敛于一个随机微分方程（SDE）。

1.1 VP-SDE 的推导

回顾 DDPM 前向过程的一步转移：

x_t = \sqrt{1 - \beta_t}\, x_{t-1} + \sqrt{\beta_t}\, \varepsilon, \quad \varepsilon \sim \mathcal{N}(0, I)

当步长很小时， $\beta_t \ll 1$ ，可以近似 $\sqrt{1-\beta_t} \approx 1 - \beta_t/2$ 。令 $\Delta t = 1$ ， $\beta_t = f(t)\Delta t$ ，则：

x_{t+\Delta t} = (1 - f(t)\Delta t/2)\, x_t + \sqrt{f(t)\Delta t}\, \varepsilon

取 $\Delta t \to 0$ 的极限， $\varepsilon \cdot \sqrt{\Delta t}$ 收敛于布朗运动 $dw$ ，得到 VP-SDE（Variance Preserving SDE）：

dx = -\frac{1}{2}\beta(t)\, x\, dt + \sqrt{\beta(t)}\, dw

其中 $\beta(t)$ 是连续化的噪声调度函数。名称"Variance Preserving"来自一个关键性质：如果 $x_0$ 的方差为 1，则 $x_t$ 的方差始终为 1——漂移项的收缩恰好抵消了扩散项的膨胀。

1.2 其他 SDE 变体

除了 VP-SDE，Song et al. (2021) 还提出了 VE-SDE（Variance Exploding）：

dx = \sqrt{\frac{d[\sigma^2(t)]}{dt}}\, dw

其中没有漂移项，噪声方差随时间单调增长。NCSN 系列模型使用 VE-SDE——它们不压缩数据，只加噪声。

VP-SDE 和 VE-SDE 的选择影响的是噪声调度的形态，而非模型的根本能力。从概率流 ODE 的统一视角看，二者等价——都能通过同一框架采样。关键区别在于数值稳定性：VP-SDE 在高维空间中更稳定，因为数据始终被归一化；VE-SDE 在低维空间中更灵活，因为噪声水平可以任意大。

1.3 连续化带来的新视角

从离散到连续的过渡不仅仅是形式上的推广。连续 SDE 框架带来了三个新的能力：

第一，任意时刻的分布。 离散 DDPM 只在整数时间步上定义分布，而 SDE 在连续时间 $t \in [0, T]$ 上定义分布 $p_t(x)$ 。这意味着可以在任意时刻插值采样，无需局限于预定义的时间步。

第二，任意 ODE 求解器。 连续框架允许使用任意高阶 ODE 求解器（Runge-Kutta、Heun 等）进行采样，而非局限于 DDPM/DDIM 的固定格式。这是加速采样的理论基础。

第三，Fokker-Planck 方程。 SDE 描述的是单条轨迹的随机动力学，而 Fokker-Planck 方程描述的是分布 $p_t(x)$ 随时间的确定性演化——这是连接 SDE 和 ODE 的桥梁。

二、反向 SDE：从贝叶斯到逆向随机微分方程

正向 SDE 描述了数据如何被噪声逐渐破坏，而生成模型的核心问题是其逆过程：如何从纯噪声 $x_T$ 出发，逐步恢复数据 $x_0$ 。这个逆过程并非简单地"倒放"正向 SDE——它需要一条全新的随机微分方程，而贝叶斯定理正是推导它的关键工具。

2.1 从离散情形出发

考虑离散前向过程 $x_0 \to x_1 \to \cdots \to x_T$ ，每一步的转移概率 $q(x_t | x_{t-1})$ 已知（高斯）。我们关心的是反向转移 $q(x_{t-1} | x_t)$ ——这不是前向转移的简单反转，因为马尔可夫链的反向转移需要考虑边际分布的信息。

利用贝叶斯定理，以 $x_0$ 为条件：

q(x_{t-1} | x_t, x_0) = \frac{q(x_t | x_{t-1}, x_0) \cdot q(x_{t-1} | x_0)}{q(x_t | x_0)}

由于前向过程是马尔可夫的， $q(x_t | x_{t-1}, x_0) = q(x_t | x_{t-1})$ ，因此：

q(x_{t-1} | x_t, x_0) = \frac{q(x_t | x_{t-1}) \cdot q(x_{t-1} | x_0)}{q(x_t | x_0)}

这个公式的三个分量都可计算： $q(x_t | x_{t-1})$ 是已知的高斯前向转移； $q(x_{t-1} | x_0)$ 和 $q(x_t | x_0)$ 由前向过程的闭式解给出。具体地，对于 VP-SDE：

q(x_t | x_0) = \mathcal{N}(x_t;\, \sqrt{\bar{\alpha}_t}\, x_0,\, (1-\bar{\alpha}_t)I)

将已知的高斯密度代入贝叶斯公式，完成配方法（completing the square）后可以发现：反向转移 $q(x_{t-1} | x_t, x_0)$ 仍然是高斯的。 其均值和方差分别为：

\mu_{t-1} = \frac{\sqrt{\bar{\alpha}_{t-1}}\beta_t}{1-\bar{\alpha}_t}\, x_0 + \frac{\sqrt{\alpha_t}(1-\bar{\alpha}_{t-1})}{1-\bar{\alpha}_t}\, x_t

\sigma_{t-1}^2 = \frac{(1-\bar{\alpha}_{t-1})\beta_t}{1-\bar{\alpha}_t}

这正是 DDPM 反向过程的转移核——DDPM 的训练目标就是学习这个高斯分布的均值（或等价地，预测 $x_0$ 或噪声 $\varepsilon$ ）。

2.2 连续化：反向 SDE

将离散的反向转移取连续极限 $\Delta t \to 0$ ，Anderson (1982) 证明了一个优美的一般性结果：如果正向 SDE 为

dx = f(x,t)\, dt + g(t)\, dw

则反向 SDE 为

dx = \left[f(x,t) - g^2(t)\, \nabla_x \log p_t(x)\right] dt + g(t)\, d\bar{w}

其中 $\bar{w}$ 是反向布朗运动（时间倒流的维纳过程）， $p_t(x)$ 是正向 SDE 在时刻 $t$ 的边际分布。

这个结果包含三个关键观察：

第一，漂移项的修正。 反向 SDE 的漂移项 = 正向漂移项 $-$ 扩散系数的平方 $\times$ 得分函数。额外的 $-g^2(t)\nabla_x \log p_t(x)$ 项是逆向修正：得分函数指向密度增长方向，减去它使得反向过程从高密度区域向低密度区域移动，与正向扩散方向相反。

第二，扩散项不变。 正向 SDE 和反向 SDE 的扩散系数都是 $g(t)$ ，只是时间方向反转（ $dw$ 变为 $d\bar{w}$ ）。这意味着随机扰动的强度在正反两个方向上相同——扩散的"幅度"不因时间反转而改变。

第三，得分函数是唯一的未知量。 在反向 SDE 中， $f(x,t)$ 和 $g(t)$ 由前向 SDE 的设计决定，是已知的。唯一需要学习的是得分函数 $\nabla_x \log p_t(x)$ 。这解释了为什么 SDE 框架天然需要得分匹配——得分函数是反向采样的核心驱动力，必须通过训练来估计。

2.3 DDPM 采样是反向 SDE 的离散化

理解了反向 SDE，再看 DDPM 的采样过程就一目了然了。DDPM 的反向步：

x_{t-1} = \frac{1}{\sqrt{\alpha_t}}\left(x_t - \frac{\beta_t}{\sqrt{1-\bar{\alpha}_t}}\, \varepsilon_\theta(x_t, t)\right) + \sigma_t\, z

其中 $\varepsilon_\theta$ 是噪声预测网络。利用噪声预测与得分函数的关系 $\varepsilon_\theta \approx -\sqrt{1-\bar{\alpha}_t}\, s_\theta(x_t, t)$ ，可以将上式改写为：

x_{t-1} \approx x_t + \Delta t \left[f(x_t, t) - g^2(t)\, s_\theta(x_t, t)\right] + g(t)\sqrt{\Delta t}\, z

这正是反向 SDE 的 Euler-Maruyama 离散化。DDPM 的每一步采样，就是在做反向 SDE 的一阶数值积分。

三、Fokker-Planck 方程

3.1 从 SDE 到分布演化

SDE $dx = f(x,t)\, dt + g(t)\, dw$ 描述了单条样本轨迹如何随时间演化。但我们更关心的是分布 $p_t(x)$ 如何演化。Fokker-Planck 方程（也称 Kolmogorov 前向方程）给出了答案：

\frac{\partial p_t(x)}{\partial t} = -\nabla_x \cdot [f(x,t)\, p_t(x)] + \frac{1}{2}g(t)^2\, \nabla_x^2 p_t(x)

第一项是漂移项——概率质量沿漂移方向流动；第二项是扩散项——概率因随机扰动而散开。对于 VP-SDE， $f(x,t) = -\frac{1}{2}\beta(t)x$ ， $g(t) = \sqrt{\beta(t)}$ 。

3.2 Fokker-Planck 方程的推导

苏剑林（2022）用 Dirac 函数方法给出了 Fokker-Planck 方程的一个简洁推导。核心思路是将边际分布表示为轨迹的期望： $p_t(x) = \mathbb{E}[\delta(x - x_t)]$ ，然后对时间求导，利用 $x_t$ 所满足的 SDE 做 Taylor 展开。

具体地，设 $x_t$ 服从 SDE $dx = f(x,t)\, dt + g(t)\, dw$ 。对 $p_t(x) = \mathbb{E}[\delta(x - x_t)]$ 关于 $t$ 求导，将 SDE 代入，对确定性部分 $f\, dt$ 展开，对随机部分 $g\, dw$ 展开至二阶（因为 $\mathbb{E}[(dw)^2] = dt$ ），整理后得到：

\frac{\partial p_t}{\partial t} = -\nabla_x \cdot (f \cdot p_t) + \frac{1}{2} \nabla_x^2 (g^2 \cdot p_t)

物理图像是清晰的：漂移项 $f$ 产生对流效应（ $\nabla_x \cdot (fp_t)$ ），将概率质量沿漂移方向输运；扩散项 $g$ 产生扩散效应（ $\nabla_x^2(g^2 p_t)$ ），使概率质量因随机扰动而散开。Fokker-Planck 方程的本质，就是将对流与扩散这两种微观机制转化为分布演化的确定性方程。

3.3 Fokker-Planck 的直觉理解

考虑一维情况。Fokker-Planck 方程说的是：概率密度 $p_t(x)$ 在每个点 $x$ 处的局部变化率，等于流入的漂移流量减去流出的扩散流量。

当 $\beta(t) = 0$ 时，只有漂移项——概率密度像流体一样沿漂移方向流动。当 $f = 0$ 时，只有扩散项——概率密度像墨水滴入水中一样均匀散开。真实的扩散过程是两种效应的叠加。

四、概率流 ODE：从 SDE 到确定性采样

4.1 核心推导

Fokker-Planck 方程中包含扩散项，这使得 $p_t(x)$ 的演化看起来是随机的。但存在一条确定性的 ODE 路径，使得沿此路径演化的分布与 SDE 完全相同。

推导过程如下。将 Fokker-Planck 方程中的 Laplace 项 $\nabla_x^2 p_t(x)$ 利用恒等式 $\nabla_x^2 p = \nabla_x \cdot (\nabla_x \log p \cdot p)$ 改写：

\frac{\partial p_t}{\partial t} = -\nabla_x \cdot \left[f\, p_t\right] + \frac{1}{2}g^2 \nabla_x \cdot \left[(\nabla_x \log p_t)\, p_t\right]

合并散度项：

\frac{\partial p_t}{\partial t} = -\nabla_x \cdot \left[\left(f - \frac{1}{2}g^2 \nabla_x \log p_t\right) p_t\right]

这恰好是如下 ODE 的 Fokker-Planck 方程（无扩散项）：

dx = \left[f(x,t) - \frac{1}{2}g(t)^2 \nabla_x \log p_t(x)\right] dt

这就是概率流 ODE（Probability Flow ODE）。它是一条确定性的路径，但产生的边缘分布 $p_t(x)$ 与 SDE 完全一致。

4.2 概率流 ODE 的意义

概率流 ODE 意味着：采样不需要随机性。 随机性只是 SDE 参数化的一种选择——概率流 ODE 表明，沿得分函数方向的确定性运动足以复现边缘分布。

这意味着：

DDPM 的随机采样和 DDIM 的确定性采样是同一方程的两种特例——DDPM 对应完整的 SDE，DDIM 对应概率流 ODE
可以使用任意 ODE 求解器加速采样，包括自适应步长方法
概率流 ODE 赋予了扩散模型一个归一化流（Normalizing Flow）的结构——它是一个可逆变换，支持精确的似然计算

4.3 等价变换：同一边际分布的一族 SDE

概率流 ODE 将扩散项完全去除，看似是一个极端的操作。但苏剑林在《一般框架之ODE篇》中证明了一个重要定理，揭示了更深层的事实：固定边际分布 $p_t(x)$ ，存在一族具有不同扩散系数的 SDE，它们产生完全相同的边际分布。

具体地，考虑一般形式的 SDE $dx = f(x,t)\, dt + \hat{g}(t)\, dw$ ，其边际分布为 $p_t(x)$ 。如果将扩散系数 $\hat{g}^2(t)$ 替换为 $\hat{g}^2(t) + h^2(t)$ （其中 $h(t) \geq 0$ 是任意非负函数），同时将漂移项调整为

\tilde{f}(x,t) = f(x,t) - \frac{1}{2}h^2(t)\, \nabla_x \log p_t(x)

则新的 SDE

dx = \tilde{f}(x,t)\, dt + \sqrt{\hat{g}^2(t) + h^2(t)}\, dw

与原 SDE 产生完全相同的边际分布 $p_t(x)$ 。

证明思路是直接的：将新 SDE 代入 Fokker-Planck 方程，利用恒等式 $p_t \nabla_x \log p_t = \nabla_x p_t$ ，可以验证扩散项的增量 $h^2$ 恰好被漂移项的修正量抵消，Fokker-Planck 方程不变。

取极端情形 $h^2 = g^2$ （即 $\hat{g} = 0$ ），新 SDE 的扩散项完全消失，漂移项变为 $f - \frac{1}{2}g^2 \nabla_x \log p_t$ ——这正是概率流 ODE。这个定理回答了一个基本问题：“为什么可以去掉随机性？“答案是：不同的随机程度可以产生完全相同的边际分布——零随机性只是这族 SDE 中的一个端点。

4.4 DDIM：概率流 ODE 的一阶求解器

等价变换定理不仅解释了概率流 ODE 的合理性，还揭示了 DDIM 确定性采样的本质。DDIM 的非随机采样公式为：

x_{t-1} = \sqrt{\bar{\alpha}_{t-1}}\, \underbrace{\frac{x_t - \sqrt{1-\bar{\alpha}_t}\, \varepsilon_\theta(x_t, t)}{\sqrt{\bar{\alpha}_t}}}_{\hat{x}_0\text{ 估计}} + \sqrt{1-\bar{\alpha}_{t-1} - \sigma_t^2}\, \varepsilon_\theta(x_t, t)

其中 $\sigma_t$ 控制随机性水平。当 $\sigma_t = 0$ 时，上式退化为：

x_{t-1} = \sqrt{\bar{\alpha}_{t-1}}\, \hat{x}_0 + \sqrt{1-\bar{\alpha}_{t-1}}\, \varepsilon_\theta(x_t, t)

这恰好是概率流 ODE 在 VP-SDE 线性漂移系数下的离散化。从等价变换的视角看，DDIM 的 $\sigma_t$ 参数对应于 $h^2(t)$ 的选择： $\sigma_t$ 越大，随机性越强，趋向 DDPM 采样； $\sigma_t = 0$ 对应零随机性，退化为概率流 ODE。

因此，DDIM 能够实现确定性采样，并非因为设计了一个全新的采样方案，而是因为它就是概率流 ODE 的一阶数值求解（Euler 方法）。DDIM 的跳步能力——从 1000 步压缩到 50 步——也由此而来：确定性 ODE 轨迹允许使用更大的步长和更高阶的求解器。

SDE 轨迹 vs ODE 轨迹

Loading visualization...

同一初始噪声点出发：左侧 SDE 采样产生随机发散的轨迹，右侧 ODE 采样产生确定性收敛的轨迹。

const W = container.clientWidth;
const H = container.clientHeight;
const margin = {top: 25, right: 20, bottom: 40, left: 50};
const iW = (W - margin.left - margin.right) / 2 - 15;
const iH = H - margin.top - margin.bottom;

const svg = d3.select(container).append("svg").attr("width", W).attr("height", H);

const colors = {
  accent: '#3b82f6', accent2: '#ef4444', muted: '#6b7280', grid: '#e5e7eb',
  sde: '#3b82f6', ode: '#ef4444', data: '#10b981'
};

// Left panel: SDE
const gL = svg.append("g").attr("transform", `translate(${margin.left},${margin.top})`);
// Right panel: ODE
const gR = svg.append("g").attr("transform", `translate(${margin.left + iW + 30},${margin.top})`);

// Labels
svg.append("text").attr("x", margin.left + iW/2).attr("y", 14).attr("text-anchor","middle").attr("font-size","13px").attr("font-weight","600").attr("fill","#3b82f6").text("SDE 采样（随机）");
svg.append("text").attr("x", margin.left + iW + 30 + iW/2).attr("y", 14).attr("text-anchor","middle").attr("font-size","13px").attr("font-weight","600").attr("fill","#ef4444").text("ODE 采样（确定性）");

// Common setup
function setupPanel(g) {
  const tScale = d3.scaleLinear().domain([0, 1]).range([0, iW]);
  const xScale = d3.scaleLinear().domain([-3, 3]).range([iH, 0]);
  g.append("g").attr("transform",`translate(0,${iH})`).call(d3.axisBottom(tScale).ticks(5).tickFormat(d => `${Math.round(d*T)}`));
  g.append("g").call(d3.axisLeft(xScale).ticks(5));
  g.append("text").attr("x", iW/2).attr("y", iH+34).attr("text-anchor","middle").attr("fill","#666").attr("font-size","11px").text("时间步 t");
  g.append("text").attr("transform","rotate(-90)").attr("x",-iH/2).attr("y",-38).attr("text-anchor","middle").attr("fill","#666").attr("font-size","11px").text("x");
  return {tScale, xScale};
}

const T = 100;

// VP-SDE parameters
function beta(t) { return 0.1 + t * 0.9; }

// Generate SDE trajectories
function genSDE(nPaths) {
  const paths = [];
  for (let p = 0; p < nPaths; p++) {
    const path = [{t: 0, x: d3.randomNormal(0, 0.1)()}];
    for (let t = 1; t <= T; t++) {
      const s = t / T;
      const dt = 1/T;
      const b = beta(s);
      const drift = -0.5 * b * path[t-1].x * dt;
      const diffusion = Math.sqrt(b * dt) * d3.randomNormal(0, 1)();
      path.push({t: s, x: path[t-1].x + drift + diffusion});
    }
    paths.push(path);
  }
  return paths;
}

// Generate ODE trajectories
function genODE(nPaths) {
  const paths = [];
  // Data point at mode
  const dataPoints = [-1.5, 0, 1.5];
  for (let p = 0; p < nPaths; p++) {
    const x1 = dataPoints[p % 3] + d3.randomNormal(0, 0.3)();
    // Reverse ODE: from noise to data
    const path = [{t: 0, x: d3.randomNormal(0, 0.5)()}];
    for (let t = 1; t <= T; t++) {
      const s = t / T;
      const dt = 1/T;
      const b = beta(s);
      // Simplified: drift toward data point with score-like force
      const score = -(path[t-1].x - x1) / (1 - s*0.8 + 0.1);
      const drift = -0.5 * b * path[t-1].x * dt + 0.5 * b * score * dt;
      path.push({t: s, x: path[t-1].x + drift});
    }
    paths.push(path);
  }
  return paths;
}

const sdePaths = genSDE(12);
const odePaths = genODE(12);

const {tScale: tsL, xScale: xsL} = setupPanel(gL);
const {tScale: tsR, xScale: xsR} = setupPanel(gR);

const lineGen = (ts, xs) => d3.line().x(d => ts(d.t)).y(d => xs(d.x)).curve(d3.curveBasis);

// Animation state
let animT = 0;
let playing = true;

// Data distribution markers
function drawDataMarkers(g, xs) {
  [-1.5, 0, 1.5].forEach(x0 => {
    g.append("rect").attr("x", iW - 8).attr("y", xs(x0) - 6).attr("width", 8).attr("height", 12)
      .attr("fill", colors.data).attr("opacity", 0.5).attr("rx", 2);
  });
}

drawDataMarkers(gL, xsL);
drawDataMarkers(gR, xsR);

// Controls
const controls = d3.select(container).append("div")
  .attr("style", "text-align:center; margin-top:8px; font-size:13px;");

const playBtn = controls.append("button").text("暂停")
  .attr("style", "margin-right:12px; padding:2px 12px; cursor:pointer; border:1px solid #999; border-radius:3px;");
const slider = controls.append("input").attr("type","range").attr("min",0).attr("max",100).attr("value",0)
  .attr("style", "width:200px;");

function render(progress) {
  drawPaths(gL, sdePaths, tsL, xsL, colors.sde, progress);
  drawPaths(gR, odePaths, tsR, xsR, colors.ode, progress);
}

render(0);

playBtn.on("click", () => {
  playing = !playing;
  playBtn.text(playing ? "暂停" : "播放");
});
slider.on("input", (e) => {
  animT = +e.target.value / 100;
  render(animT);
});

function animate() {
  if (playing && animT < 1) {
    animT += 0.005;
    if (animT > 1) animT = 1;
    slider.property("value", animT * 100);
    render(animT);
  }
  requestAnimationFrame(animate);
}
animate();

五、得分函数的几何意义

5.1 什么是得分函数？

得分函数定义为对数密度的梯度：

s(x, t) = \nabla_x \log p_t(x)

它的几何意义直观而深刻：得分函数指向概率密度增长最快的方向。 在概率密度高的区域， $p_t(x)$ 大， $\log p_t(x)$ 也大，但梯度小（因为已接近峰值，变化平缓）；在概率密度低的区域， $\log p_t(x)$ 小，但梯度可能很大（因为密度变化剧烈）。

5.2 得分函数与噪声预测的关系

在 DDPM 中，我们训练模型预测噪声 $\varepsilon_\theta(x_t, t)$ 。而得分函数可以通过 Tweedie 公式与噪声预测建立联系：

\nabla_{x_t} \log p_t(x_t) = -\frac{1}{1-\bar{\alpha}_t}\left(x_t - \sqrt{\bar{\alpha}_t}\,\mathbb{E}[x_0|x_t]\right) \approx -\frac{\varepsilon_\theta(x_t, t)}{\sqrt{1-\bar{\alpha}_t}}

这意味着：噪声预测和得分估计是同一对象的两种参数化，只差一个时间相关的缩放因子。 DDPM 训练的 $\varepsilon_\theta$ 就是在学习得分函数。

5.3 得分函数的几何直觉

考虑一个一维混合高斯分布 $p(x) = 0.5\,\mathcal{N}(-2, 0.5) + 0.5\,\mathcal{N}(2, 0.5)$ 。得分函数 $\nabla_x \log p(x)$ 在两个峰值处接近零（局部极值），在峰值之间指向最近的峰（密度增长最快的方向），在远离峰值的区域指向中心（因为尾部密度很低，梯度指向密度更高的区域）。

这种几何结构暗示了一种采样策略：从任意点出发，沿得分函数方向移动，最终会到达高密度区域。 这正是朗之万动力学的思路。

得分函数向量场

Loading visualization...

2D 混合高斯分布上的得分函数 $\\nabla_x \\log p(x)$。箭头指向密度增长最快的方向。点击任意位置展示朗之万采样轨迹。

const W = container.clientWidth;
const H = container.clientHeight;
const margin = {top: 25, right: 20, bottom: 40, left: 50};
const iW = W - margin.left - margin.right;
const iH = H - margin.top - margin.bottom;

const svg = d3.select(container).append("svg").attr("width", W).attr("height", H);
const g = svg.append("g").attr("transform", `translate(${margin.left},${margin.top})`);

// 2D mixture of Gaussians
const modes = [{mx: -1.5, my: -1, sx: 0.6, sy: 0.6, w: 0.4},
               {mx: 1.5, my: 1, sx: 0.5, sy: 0.5, w: 0.35},
               {mx: 0, my: 1.8, sx: 0.4, sy: 0.4, w: 0.25}];

function pxy(x, y) {
  let sum = 0;
  modes.forEach(m => {
    const ex = -0.5 * ((x-m.mx)/m.sx)**2;
    const ey = -0.5 * ((y-m.my)/m.sy)**2;
    sum += m.w / (2 * Math.PI * m.sx * m.sy) * Math.exp(ex + ey);
  });
  return sum;
}

function score(x, y) {
  // Numerical gradient of log p
  const eps = 0.01;
  const px = pxy(x, y);
  const dpx = (pxy(x+eps, y) - pxy(x-eps, y)) / (2*eps);
  const dpy = (pxy(x, y+eps) - pxy(x, y-eps)) / (2*eps);
  const sx = px > 1e-20 ? dpx / px : 0;
  const sy = px > 1e-20 ? dpy / px : 0;
  return {sx, sy};
}

const xRange = [-4, 4], yRange = [-3, 3];
const xScale = d3.scaleLinear().domain(xRange).range([0, iW]);
const yScale = d3.scaleLinear().domain(yRange).range([iH, 0]);

// Axes
g.append("g").attr("transform",`translate(0,${iH})`).call(d3.axisBottom(xScale).ticks(8));
g.append("g").call(d3.axisLeft(yScale).ticks(6));
g.append("text").attr("x",iW/2).attr("y",iH+32).attr("text-anchor","middle").attr("fill","#666").attr("font-size","11px").text("x₁");
g.append("text").attr("transform","rotate(-90)").attr("x",-iH/2).attr("y",-38).attr("text-anchor","middle").attr("fill","#666").attr("font-size","11px").text("x₂");

// Density heatmap
const nBins = 80;
const xStep = (xRange[1]-xRange[0])/nBins;
const yStep = (yRange[1]-yRange[0])/nBins;
let maxP = 0;
const densityData = [];
for (let i = 0; i < nBins; i++) {
  for (let j = 0; j < nBins; j++) {
    const x = xRange[0] + (i+0.5)*xStep;
    const y = yRange[0] + (j+0.5)*yStep;
    const p = pxy(x, y);
    maxP = Math.max(maxP, p);
    densityData.push({x, y, p});
  }
}

const heatColor = d3.scaleSequential(d3.interpolateBlues).domain([0, maxP]);

g.selectAll(".heat").data(densityData).enter().append("rect")
  .attr("x", d => xScale(d.x - xStep/2)).attr("y", d => yScale(d.y + yStep/2))
  .attr("width", iW/nBins + 0.5).attr("height", iH/nBins + 0.5)
  .attr("fill", d => heatColor(d.p)).attr("opacity", 0.35);

// Score vector field
const nArrows = 14;
const arrowGroup = g.append("g").attr("class", "arrows");
for (let i = 0; i < nArrows; i++) {
  for (let j = 0; j < nArrows; j++) {
    const x = xRange[0] + (i+0.5)*(xRange[1]-xRange[0])/nArrows;
    const y = yRange[0] + (j+0.5)*(yRange[1]-yRange[0])/nArrows;
    const s = score(x, y);
    const mag = Math.sqrt(s.sx*s.sx + s.sy*s.sy);
    if (mag < 0.1 || mag > 10) continue;
    const scale = Math.min(mag, 3) / mag * 12;
    const px = xScale(x), py = yScale(y);
    arrowGroup.append("line")
      .attr("x1", px).attr("y1", py)
      .attr("x2", px + s.sx/mag * scale).attr("y2", py - s.sy/mag * scale)
      .attr("stroke", "#ef4444").attr("stroke-width", 1.2).attr("opacity", 0.6)
      .attr("marker-end", "url(#arrowhead)");
  }
}

// Arrow marker
svg.append("defs").append("marker").attr("id","arrowhead")
  .attr("viewBox","0 0 10 10").attr("refX",8).attr("refY",5)
  .attr("markerWidth",5).attr("markerHeight",5).attr("orient","auto")
  .append("path").attr("d","M 0 0 L 10 5 L 0 10 z").attr("fill","#ef4444").attr("opacity",0.7);

// Langevin trajectory on click
const langGroup = g.append("g").attr("class", "langevin");

svg.on("click", function(event) {
  const [mx, my] = d3.pointer(event, g.node());
  const x0 = xScale.invert(mx);
  const y0 = yScale.invert(my);
  if (x0 < xRange[0] || x0 > xRange[1] || y0 < yRange[0] || y0 > yRange[1]) return;

// Langevin dynamics: x_{k+1} = x_k + eps * score(x_k) + sqrt(2*eps) * z
  const eps = 0.05;
  const nSteps = 150;
  const traj = [{x: x0, y: y0}];
  let cx = x0, cy = y0;
  for (let k = 0; k < nSteps; k++) {
    const s = score(cx, cy);
    const nx = cx + eps * s.sx + Math.sqrt(2*eps) * d3.randomNormal(0,1)();
    const ny = cy + eps * s.sy + Math.sqrt(2*eps) * d3.randomNormal(0,1)();
    cx = nx; cy = ny;
    traj.push({x: cx, y: cy});
  }

langGroup.selectAll("*").remove();
  const line = d3.line().x(d => xScale(d.x)).y(d => yScale(d.y)).curve(d3.curveBasis);
  langGroup.append("path").attr("d", line(traj))
    .attr("fill","none").attr("stroke","#f59e0b").attr("stroke-width",2).attr("opacity",0.8);
  langGroup.append("circle").attr("cx", xScale(x0)).attr("cy", yScale(y0))
    .attr("r",4).attr("fill","#f59e0b");
  langGroup.append("circle").attr("cx", xScale(cx)).attr("cy", yScale(cy))
    .attr("r",4).attr("fill","#10b981");
});

// Instructions
const info = d3.select(container).append("div")
  .attr("style","text-align:center; margin-top:6px; font-size:12px; color:#888;");
info.html("点击图中任意位置，展示从该点出发的朗之万采样轨迹（<span style='color:#f59e0b'>起点</span> → <span style='color:#10b981'>终点</span>）");

六、朗之万动力学与得分匹配

6.1 朗之万动力学

朗之万动力学（Langevin Dynamics）是一种基于得分函数的 MCMC 采样方法。给定目标分布 $p(x)$ ，朗之万迭代为：

x_{k+1} = x_k + \eta\, \nabla_x \log p(x_k) + \sqrt{2\eta}\, z_k, \quad z_k \sim \mathcal{N}(0, I)

其中 $\eta$ 是步长。当 $\eta \to 0$ 且迭代次数 $K \to \infty$ 时， $x_K$ 的分布收敛于 $p(x)$ 。

直觉非常清晰：每一步沿得分函数方向移动（向高密度区域前进），同时添加少量随机扰动以避免陷入局部最优。随机扰动的方差为 $2\eta$ ——步长越大，扰动越大，探索越充分，但收敛越慢。

6.2 得分匹配：如何学习得分函数？

如果已知 $p(x)$ ，可以直接计算 $\nabla_x \log p(x)$ 。但在生成模型中 $p(x)$ 正是我们要学习的目标——这是一个循环依赖。

得分匹配（Hyvärinen, 2005）提供了一个巧妙的绕行方案。核心结果是：以下两个目标关于 $\theta$ 的最优解相同：

\mathbb{E}_{p(x)}\left[\|\nabla_x \log p(x) - s_\theta(x)\|^2\right] = \mathbb{E}_{p(x)}\left[\|\nabla_x \log p(x)\|^2 + 2\,\text{tr}(\nabla_x s_\theta(x)) + \|s_\theta(x)\|^2\right] + C

第一个目标需要知道真实的 $\nabla_x \log p(x)$ ，不可计算。但第二个等价目标中， $\|\nabla_x \log p(x)\|^2$ 是常数（不依赖 $\theta$ ）， $\text{tr}(\nabla_x s_\theta(x))$ 只涉及模型的 Hessian 迹， $\|s_\theta(x)\|^2$ 直接可计算。因此，最小化第二个目标等价于最小化第一个目标，而无需知道 $p(x)$ 的解析形式。

6.3 去噪得分匹配

原始得分匹配在实践中面临高维 Hessian 迹的计算困难。去噪得分匹配（Vincent, 2011）提供了一个更实用的替代方案：

\mathcal{L}_{\text{DSM}} = \mathbb{E}_{q(\tilde{x}|x)}\mathbb{E}_{p(x)}\left[\left\|s_\theta(\tilde{x}) - \nabla_{\tilde{x}} \log q(\tilde{x}|x)\right\|^2\right]

其中 $q(\tilde{x}|x)$ 是加噪分布（如 $\mathcal{N}(\tilde{x}; x, \sigma^2 I)$ ），而 $\nabla_{\tilde{x}} \log q(\tilde{x}|x) = -(\tilde{x}-x)/\sigma^2$ 是解析可知的。

这正是 DDPM 训练目标的得分匹配视角：在噪声水平 $\sigma_t$ 下，用加噪数据训练得分模型，使得模型在每个噪声水平上都能准确估计得分函数。

6.4 采样 = 沿得分函数走

将以上内容串联起来，生成采样的完整逻辑是：

在多个噪声水平上训练得分函数 $s_\theta(x, t) \approx \nabla_x \log p_t(x)$
从高噪声分布 $p_T(x) \approx \mathcal{N}(0, I)$ 中采样初始点
沿得分函数方向逐步移动，降低噪声水平，最终到达数据分布 $p_0(x)$

这可以用朗之万动力学（随机）或概率流 ODE（确定性）实现。两种方式在分布层面等价，但在单条轨迹上行为迥异：朗之万轨迹是随机的、迂回的，概率流 ODE 轨迹是确定的、直接的。

七、应用：Score-Based Models 与 DDIM 加速

7.1 Score-Based Generative Modeling

Song & Ermon (2019) 的 NCSN（Noise Conditional Score Network）是最早系统使用得分匹配进行生成的模型。其核心思路是：在多个噪声水平上训练一个共享参数的得分网络 $s_\theta(x, \sigma)$ ，然后使用退火朗之万动力学（Annealed Langevin Dynamics）从高噪声逐步退火到低噪声进行采样。

NCSN 的训练目标与 DDPM 的 $L_{\text{simple}}$ 其实相同——只是参数化不同（得分 vs 噪声）。Song et al. (2021) 在 ICLR 论文中明确建立了这一等价关系，并提出了 SDE 框架将二者统一。

7.2 DDIM 加速采样

概率流 ODE 是 DDIM 加速采样的理论基础。由于 ODE 轨迹是确定性的，我们可以使用任意 ODE 求解器在连续时间上积分。高阶求解器（如 Heun’s method、RK45）可以用更少的函数评估达到相同的精度——这正是 DDIM 允许跳步采样的数学根源。

实践中，DDIM 将 1000 步的 DDPM 采样压缩到 20-50 步，质量几乎无损。更激进的 DPM-Solver 使用三阶 ODE 求解器，仅需 10-15 步即可获得高质量样本。这些加速方法都依赖于概率流 ODE 的确定性——随机 SDE 轨迹无法被高阶求解器加速，因为每步的随机性破坏了光滑性假设。

从离散马尔可夫链到连续 SDE，从反向 SDE 到概率流 ODE，我们建立了扩散模型的随机微分方程视角。反向 SDE 揭示了得分的核心地位——它是逆向采样的唯一驱动力；Fokker-Planck 方程是连接随机与确定性的桥梁；等价变换定理则统一了从 DDPM 到 DDIM 的采样谱系。得分函数是贯穿始终的核心对象——它既定义了 ODE 的动力，又指导了朗之万采样的方向，还是反向 SDE 漂移修正的关键。下一篇文章将从 ODE 视角出发，探索 Flow Matching 如何优化传输路径，以及一致性模型如何将 ODE 路径蒸馏为单步映射。

参考文献

苏剑林. 生成扩散模型漫谈（五）：一般框架之SDE篇. https://kexue.fm/archives/9209
苏剑林. 生成扩散模型漫谈（六）：一般框架之ODE篇. https://kexue.fm/archives/9229
Song, Y., et al. (2021). Score-Based Generative Modeling through Stochastic Differential Equations. ICLR 2021.
Anderson, B. D. O. (1982). Reverse-Time Diffusion Equation Models. Stochastic Processes and their Applications.
Hyvärinen, A. (2005). Estimation of Non-Normalized Statistical Models by Score Matching. JMLR.
Vincent, P. (2011). A Connection Between Score Matching and Denoising Autoencoders. Neural Computation.

一、从离散到连续：马尔可夫链的 SDE 极限#

1.1 VP-SDE 的推导#

1.2 其他 SDE 变体#

1.3 连续化带来的新视角#

二、反向 SDE：从贝叶斯到逆向随机微分方程#

2.1 从离散情形出发#

2.2 连续化：反向 SDE#

2.3 DDPM 采样是反向 SDE 的离散化#

三、Fokker-Planck 方程#

3.1 从 SDE 到分布演化#

3.2 Fokker-Planck 方程的推导#

3.3 Fokker-Planck 的直觉理解#

四、概率流 ODE：从 SDE 到确定性采样#

4.1 核心推导#

4.2 概率流 ODE 的意义#

4.3 等价变换：同一边际分布的一族 SDE#

4.4 DDIM：概率流 ODE 的一阶求解器#

五、得分函数的几何意义#

5.1 什么是得分函数？#

5.2 得分函数与噪声预测的关系#

5.3 得分函数的几何直觉#

六、朗之万动力学与得分匹配#

6.1 朗之万动力学#

6.2 得分匹配：如何学习得分函数？#

6.3 去噪得分匹配#

6.4 采样 = 沿得分函数走#

七、应用：Score-Based Models 与 DDIM 加速#

7.1 Score-Based Generative Modeling#

7.2 DDIM 加速采样#

相关概念#

参考文献#

相关文章