微分方程速成:ODE 與 SDE
全書用到的微分方程背景都在這一頁。ODE 部分:把 ODE 看成「速度場加上一條處處順著箭頭走的軌跡」;存在唯一性定理保證起點一定就只有一條路,而且兩條路不會交叉;指數積分因子把半線性 ODE 的線性部分精確解掉;Euler、Runge–Kutta、Adams–Bashforth、Picard 迭代分別代表「單步」「多步」「平行於時間」三種取得資訊的方式。SDE 部分:在 Euler 更新上加一個大小為 \(\sqrt{\Delta t}\) 的高斯雜訊,取極限就是 SDE;它要用 Itô 積分(左端點規則)來理解,用 Euler–Maruyama 來數值求解。
讀完這頁你應該能夠
- 透過速度場、解軌跡以及存在唯一性的角色來理解一個 ODE。
- 用積分因子揭露半線性 ODE 的線性結構,並理解為什麼這個結構對 diffusion 的 ODE 解法有用。
- 分辨單步、多步、平行於時間三種數值觀點,認出 Euler、Runge–Kutta、Adams–Bashforth、Picard 迭代是各自的代表。
- 理解 SDE 怎麼用布朗雜訊擴充確定性動態,它的微分記號和 Itô 積分的關係,以及 Euler–Maruyama 怎麼把它離散化。
書在本附錄開頭的提示:這個附錄提供全書用到的微分方程背景。請注意動態、軌跡和數值近似之間的連結;後面的 diffusion 章節會把這些想法特化成 probability flow ODE 和隨機的 diffusion 過程。
前置知識
- 微積分:導數、連鎖律、乘積法則、微積分基本定理;機率:高斯分佈、變異數相加。
- 如果想從更直觀的角度入門,可以先看 補充課 B(微分方程入門)和 補充課 E(布朗運動與 SDE);這一頁照書的順序,並附上用到這些工具的正文頁面。
「書 §」標出原文位置。灰框「轉述者補充」是書裡沒有的內容;藍框「推導」展開書中省略的步驟。本附錄 573–813 行之間被註解掉的舊版解法說明沒有收。
01三種微分方程書 附錄 A 開頭
微分方程是建模動態系統的基本工具,大致分成常微分方程(ODE)、隨機微分方程(SDE)和偏微分方程(PDE)。
- ODE 依一條精確的規則描述系統狀態怎麼隨時間變化,所以知道起點就能精確決定未來的路徑。
- SDE 在這個演化上加入隨機性,描述雜訊或不確定性怎麼影響系統的行為,結果是機率性的,而不是固定的。
- PDE 說明依賴多個變數(例如時間和空間)的函數怎麼一起演化,描述熱的擴散、波的移動,或隨機系統中機率密度隨時間的演化(劇透:Fokker–Planck 方程,附錄 B,B34)。
這幾種微分方程構成理解系統在確定性和隨機影響下怎麼隨時間與空間演化的基本語言。這一章提供微分方程的必要前置知識。
02ODE 的直覺書 §A.1、§A.1.1
這一節介紹 ODE 的基本理論,強調給定初始條件時解的唯一性,也介紹用數值解法實際求解 ODE 的方法。
確定性的過程叫常微分方程(ODE)。在多變數的情形,考慮這種形式的系統:
\(\mathbf x(t)\in\mathbb R^D\) 是表示系統在時間 \(t\) 狀態的向量值函數,\(\mathbf v:\mathbb R^D\times\mathbb R\to\mathbb R^D\) 是一個向量場,指定空間和時間中每一點的變化方向與大小。

解 ODE 的高層次直覺。把向量場 \(\mathbf v(\mathbf x,t)\) 想成一片會變動的箭頭地景,告訴你一個點 \(\mathbf x\) 在任何時間 \(t\) 該怎麼移動。解微分方程就是在這片場裡描出一條曲線 \(\mathbf x(t)\),使曲線在任一點的切線(也就是瞬時速度)都和 \(\mathbf v(\mathbf x(t),t)\) 給的向量一致。
- 向量場的觀點:\(\mathbf v(\mathbf x,t)\) 定義東西該怎麼動,給出運動或變化的局部「指示」。
- 軌跡的觀點:解 \(\mathbf x(t)\) 是一個粒子每一瞬間都遵守向量場 \(\mathbf v\) 的規則時會走的路徑。
所以解一個 ODE,就像把一個粒子放進流場,觀察它隨時間去了哪裡。
全書的 PF-ODE(B08)、flow matching 的速度場(B12)、flow map(B29)都是這個式子,只是 \(\mathbf v\) 換成各自學到或推出的速度場。
03存在與唯一書 §A.1.2
到目前為止,解 ODE 意味著找一條在每一點都順著向量場方向走的路徑。直觀上,這像是追蹤一個粒子沿速度定義的流動的軌跡。但這個圖像引出一個重要的問題:
回答這些問題很重要,因為它告訴我們系統的行為能不能從起始位置可靠地預測。存在唯一性定理給出向量場的條件,保證從任一個給定的初始點出發恰好有一條路徑。這確保解的行為一致,是 ODE 理論的基石。
局部(時間上)的存在唯一性
下面敘述局部版本的定理:對一個給定的初始條件,在初始時間附近的一段區間上,解存在而且唯一。
證明概要(書)。存在唯一性定理可以用 Picard–Lindelöf 迭代法構造性地證明。這個方法產生一串收斂到解 \(\mathbf x(t)\) 的函數 \(\{\mathbf x_n(t)\}\),迭代定義為
$$\mathbf x_{n+1}(t)=\mathbf x_0+\int_{t_0}^t\mathbf v(\mathbf x_n(s),s)\,ds.$$- 從初始猜測 \(\mathbf x_0(t)=\mathbf x_0\)(常數函數)開始。
- 用積分形式反覆修正 \(\mathbf x_n(t)\)。
- 在 Lipschitz 條件下,套用壓縮映射定理(contraction mapping theorem)保證收斂。
為什麼是壓縮:一行估計
令 \(\mathcal T[\mathbf x](t):=\mathbf x_0+\int_{t_0}^t\mathbf v(\mathbf x(s),s)\,ds\)。對兩條候選路徑 \(\mathbf x,\mathbf y\),在 \(|t-t_0|\le\delta\) 上,
$$\big\|\mathcal T[\mathbf x](t)-\mathcal T[\mathbf y](t)\big\|\le\Big|\int_{t_0}^t L\,\|\mathbf x(s)-\mathbf y(s)\|\,ds\Big|\le L\delta\,\sup_s\|\mathbf x(s)-\mathbf y(s)\|.$$取 \(\delta\lt1/L\),\(\mathcal T\) 就把兩條路徑的距離至少縮小成 \(L\delta\lt1\) 倍:它是壓縮映射,有唯一的不動點,迭代從任何起點收斂到它。不動點滿足 \(\mathbf x=\mathcal T[\mathbf x]\),也就是積分形式的 ODE。「局部」就來自這個 \(\delta\)。(更細的估計 \(\|\mathbf x_n-\mathbf x\|\le\frac{(L|t-t_0|)^n}{n!}\cdot\text{常數}\) 顯示其實在整個 Lipschitz 區域上都收斂,Demo B 會看到這個 \(n!\)。)
證明的核心是 Picard–Lindelöf 迭代法;它的核心想法也在 §9.6(ParaDiGMS,B27)被用來加速 diffusion model 的取樣。
全域(時間上)的存在唯一性
局部定理保證在一小段時間區間上解存在,「全域(時間上)的存在唯一性定理」則在額外的正則性條件下,把這個結果延伸到整個區間 \([t_0,T]\)。這類結果裡著名的是 Carathéodory 定理,它在兩個主要假設下保證 ODE 解的全域存在唯一性:狀態變數上的局部 Lipschitz 連續性,以及線性成長界。
- 對 \(\mathbf x\) 的(局部)Lipschitz 條件:存在一個在 \([0,T]\) 上可積的函數 \(\mathrm{Lip}(t)\),使得對所有 \(\mathbf x_1,\mathbf x_2\in\mathbb R^D\),\(\|\mathbf v(\mathbf x_1,t)-\mathbf v(\mathbf x_2,t)\|\le\mathrm{Lip}(t)\,\|\mathbf x_1-\mathbf x_2\|\)。
- 線性成長條件:存在一個在 \([0,T]\) 上可積的函數 \(M(t)\),使得對所有 \(\mathbf x\in\mathbb R^D\),\(\|\mathbf v(\mathbf x,t)\|\le M(t)\,(1+\|\mathbf x\|)\)。
這些假設的完整討論、正式敘述和詳細證明見 Reid (1971)。
總結:給一個由隨時間變化的速度場定義的 ODE 一個初始條件,粒子流的軌跡就唯一確定。
兩個條件各自防止什麼:兩個反例(轉述者補充)
沒有 Lipschitz,唯一性會壞。\(\frac{dx}{dt}=\sqrt{|x|}\)、\(x(0)=0\)。\(x(t)\equiv0\) 是解,\(x(t)=t^2/4\) 也是解(\(\frac{d}{dt}\frac{t^2}{4}=\frac t2=\sqrt{t^2/4}\)),甚至「先停在 0 一段時間再出發」的都是解。問題出在 \(\sqrt{|x|}\) 在 0 附近斜率無界,不是 Lipschitz。
沒有線性成長,解會在有限時間爆掉。\(\frac{dx}{dt}=x^2\)、\(x(0)=1\),解是 \(x(t)=\frac{1}{1-t}\),在 \(t=1\) 衝到無窮大。\(x^2\) 是局部 Lipschitz 的(所以局部解存在唯一),但它成長得比線性快,解活不過整個區間。線性成長條件 \(\|\mathbf v\|\le M(1+\|\mathbf x\|)\) 保證 \(\|\mathbf x\|\) 最多指數成長,不會在有限時間發散。
04解不相交書 §A.1.2
局部存在唯一性定理保證的唯一性,蘊含一個基本性質:從不同初始條件出發的兩條解軌跡不會相交。這反映 ODE 的確定性:每個狀態沿唯一的路徑演化。下面的推論把這件事形式化。
證明(書,逐步)
用反證法。假設存在某個 \(t^*\in(0,T]\) 使 \(\mathbf x_1(t^*)=\mathbf x_2(t^*)\)。定義兩解第一次相遇的時間
$$t_0:=\inf\{t\in[0,T]:\mathbf x_1(t)=\mathbf x_2(t)\}.$$因為 \(\mathbf x_1(0)\ne\mathbf x_2(0)\) 而 \(t^*\) 在這個集合裡,所以 \(t_0>0\)。由 \(\mathbf x_1\)、\(\mathbf x_2\) 的連續性(相遇點的集合是閉集),\(\mathbf x_1(t_0)=\mathbf x_2(t_0)\)。
考慮初值問題 \(\frac{d\mathbf x}{dt}=\mathbf v(\mathbf x,t)\)、\(\mathbf x(t_0)=\mathbf x_1(t_0)\)。由 ODE 的唯一性定理,\(\mathbf x_1\) 和 \(\mathbf x_2\) 都是它的解,所以在 \([t_0,T]\) 上必須一致。同樣地往時間的反方向套用唯一性,它們在 \([0,t_0]\) 上也一致。所以對所有 \(t\in[0,T]\),\(\mathbf x_1(t)=\mathbf x_2(t)\),和 \(\mathbf x_1(0)\ne\mathbf x_2(0)\) 矛盾。因此兩解在 \([0,T]\) 上不相交。
這個定理保證解路徑不相交,對 flow map 模型(第 10、11 章,B28–B31)提供了隱藏但關鍵的支持。轉述者補充:每個起點都有唯一的終點,不同的起點送到不同的終點,所以 \(\bm\Psi_{s\to t}\) 是一個定義良好的雙射,半群性質(B29)也由此而來。
05指數積分因子書 §A.1.3
即使由一般的、隨時間變化的速度 \(\mathbf v\) 決定的 ODE 沒有閉式解,在某些特殊情形,我們可以解析地求解,或把它的表述化簡成結構更好的形式。
一個示範的例子
考慮線性的純量 ODE
$$\frac{dx(t)}{dt}=L(t)\,x(t),$$\(L(t)\in\mathbb R\) 是連續函數。這個方程有閉式解,答案大家熟悉(對任何 \(s\)、\(t\)):
$$x(t)=x(s)\cdot\exp\Big(\int_s^tL(\tau)\,d\tau\Big).$$這個公式說明解依一個累積了時變係數 \(L(t)\) 效果的指數因子演化。這引出指數積分因子:
它在動態同時含有線性和非線性成分的更一般設定裡特別有用。
半線性 ODE 與指數積分因子
現在考慮更廣的一類 ODE,叫半線性(semilinear)ODE。它把動態分成一個(對狀態變數)線性的部分和一個非線性的剩餘:
\(\mathbf x(t)\in\mathbb R^D\) 是狀態向量,\(L(t)\) 是純量的連續函數,\(\mathbf N:\mathbb R^D\times[0,T]\to\mathbb R^D\) 是非線性的向量場。這種半線性結構自然出現在許多物理和工程系統裡,也出現在 diffusion model 的 probability flow ODE 表述中(B08:\(\frac{d\mathbf x}{dt}=f(t)\mathbf x-\frac12g^2(t)\nabla\log p_t\),\(f(t)\mathbf x\) 是線性部分)。認出這個結構,就能用積分因子表示:已知的線性動態解析地處理,只把剩下的項留給數值近似。這個原則是第 9 章 diffusion ODE 解法(B24–B27)的核心。
第 1 步:用積分因子隔離非線性項。從半線性 ODE 減去線性漂移:
$$\frac{d\mathbf x(t)}{dt}-L(t)\,\mathbf x(t)=\mathbf N(\mathbf x(t),t).$$為了吸收線性項,兩邊乘上反向的積分因子 \(\mathcal E^{-1}(s\to t)=\exp(-\int_s^tL(\tau)\,d\tau)\),再對左邊用乘積法則:
$$\mathcal E^{-1}(s\to t)\Big(\frac{d\mathbf x(t)}{dt}-L(t)\,\mathbf x(t)\Big)=\frac{d}{dt}\Big[\mathcal E^{-1}(s\to t)\,\mathbf x(t)\Big].$$乘積法則那一步
\(\frac{d}{dt}\mathcal E^{-1}(s\to t)=\frac{d}{dt}\exp\big(-\int_s^tL\big)=-L(t)\,\mathcal E^{-1}(s\to t)\)。所以
$$\frac{d}{dt}\big[\mathcal E^{-1}\mathbf x\big]=\mathcal E^{-1}\frac{d\mathbf x}{dt}+\big(-L(t)\mathcal E^{-1}\big)\mathbf x=\mathcal E^{-1}\Big(\frac{d\mathbf x}{dt}-L(t)\mathbf x\Big).$$於是方程變成
$$\frac{d}{dt}\Big[\mathcal E^{-1}(s\to t)\,\mathbf x(t)\Big]=\mathcal E^{-1}(s\to t)\,\mathbf N(\mathbf x(t),t).$$這個轉換隔離出非線性成分,簡化了原方程,讓我們在一個轉換過的座標系裡只專注於非線性動態。
第 2 步:對時間積分。兩邊從 \(s\) 積到 \(t\):
$$\int_s^t\frac{d}{d\tau}\Big[\mathcal E^{-1}(s\to\tau)\,\mathbf x(\tau)\Big]d\tau=\int_s^t\mathcal E^{-1}(s\to\tau)\,\mathbf N(\mathbf x(\tau),\tau)\,d\tau.$$左邊就是轉換後變數在 \(t\) 和 \(s\) 的差:\(\mathcal E^{-1}(s\to t)\mathbf x(t)-\mathbf x(s)\)(因為 \(\mathcal E^{-1}(s\to s)=1\))。所以
$$\mathcal E^{-1}(s\to t)\,\mathbf x(t)=\mathbf x(s)+\int_s^t\mathcal E^{-1}(s\to\tau)\,\mathbf N(\mathbf x(\tau),\tau)\,d\tau.$$第 3 步:解出 \(\mathbf x(t)\)。兩邊乘上指數流 \(\mathcal E(s\to t)\),得到解:
(\(\mathcal E(s\to t)\,\mathcal E^{-1}(s\to\tau)=\exp\big(\int_s^tL-\int_s^\tau L\big)=\mathcal E(\tau\to t)\)。)
解自然地分成線性和非線性兩個成分。指數積分器(exponential integrator)利用這個結構:線性部分以精確的閉式解掉,只離散化非線性的剩餘。這確保步長由非線性動態決定,而不是由可能很大的線性係數決定(例如 §9.1(B24)裡 exponential Euler 和普通 Euler 更新的比較)。Demo A 的「λ = 50」會直接看到這句話的意思。
最簡單的指數積分器:exponential Euler(轉述者補充)
在一步 \([t_n,t_n+h]\) 裡把非線性項凍結成起點的值 \(\mathbf N(\mathbf x_n,t_n)\),常數 \(L\) 時積分可以算出來:\(\int_{t_n}^{t_n+h}e^{L(t_n+h-\tau)}d\tau=\frac{e^{Lh}-1}{L}\)。所以
$$\mathbf x_{n+1}=e^{Lh}\,\mathbf x_n+\frac{e^{Lh}-1}{L}\,\mathbf N(\mathbf x_n,t_n).$$當 \(L=-\lambda\lt0\) 很大時,\(e^{-\lambda h}\) 永遠在 0 和 1 之間,不管 \(h\) 多大都不會放大誤差;普通 Euler 的線性部分是 \(1-\lambda h\),\(h\gt2/\lambda\) 時絕對值超過 1,每一步都把誤差放大。B25 的 DDIM 正是這種「線性部分精確、預測值凍結」的一階指數積分器。
06數值解法的結構書 §A.1.4
考慮 ODE \(\frac{d\mathbf x(t)}{dt}=\mathbf v(\mathbf x(t),t)\) 加上一個初始條件。解 ODE 就是找出一條連續軌跡 \(\mathbf x(t)\),它的瞬時變化在每個時間都和速度場 \(\mathbf v\) 一致。等價地,解滿足積分方程
如果這個積分沒有閉式,就用有限個點上的速度場計算值,數值地近似軌跡。
從連續軌跡到數值更新
數值解法先選一組離散的時間格點 \(t_0,t_1,\dots,t_M\),產生對應的近似 \(\mathbf x_0,\mathbf x_1,\dots,\mathbf x_M\)。間距 \(h_i:=t_{i+1}-t_i\) 叫步長。較小的步長一般提供更多關於軌跡的局部資訊,但需要更多次函數計算。數值解法的目的,是有效率地利用可用的 \(\mathbf v\) 計算值,用相對少的步數得到準確的軌跡。不同的解法主要差在下一步更新用的資訊從哪裡來,由此得到下面這個有用的組織方式。
數值解法的結構觀點
對本書相關的方法,區分兩種大的數值觀點很有用:
- 逐步推進(time-stepping)的方法:從目前的狀態出發,把軌跡從一個時間點推進到下一個。Euler、Runge–Kutta、Adams–Bashforth 是代表。其中再分:
- 單步方法:用目前的狀態,以及在目前這個區間裡產生的額外計算值,構造下一個狀態;
- 多步方法:另外重複使用先前已完成的步裡的資訊。
- 平行於時間(parallel-in-time)的觀點:不是嚴格地一個時間點接一個時間點完成,而是把軌跡的計算寫成可以同時更新好幾個時間點的暫定狀態。Picard 迭代這類不動點構造提供了這個想法的基本機制。
單步方法不要和單次函數計算搞混。例如 Euler 每一步用一次速度計算,Runge–Kutta 方法則可能在中間的階段(stage)計算好幾次速度,但它仍然是單步方法。決定性的區別在於是否需要先前已完成的步的資訊。
逐步推進 → 單步:Euler、Runge–Kutta;多步:Adams–Bashforth。
平行於時間 → Picard 不動點迭代。
同樣的數值想法都可以往時間的正向或反向使用。
這個區分在後面的 diffusion 取樣裡特別有用:Runge–Kutta 類的構造透過看目前這一步裡面取得更多資訊,Adams–Bashforth 類的構造則透過往回看先前算過的值取得資訊(B26、B27)。
正向與反向時間的求解
上面的分類和 ODE 往哪個方向積分無關。一般的正向積分是 \(t_{n+1}\gt t_n\)、\(h_n=t_{n+1}-t_n\gt0\)。diffusion 取樣通常從一個有雜訊的終點狀態出發,用遞減的時間格點:\(t_{n+1}\lt t_n\)、\(h_n=t_{n+1}-t_n\lt0\)。只要把步長當成帶正負號的量,同樣的數值公式在兩個方向都適用。所以反向時間積分一個 ODE 不需要新的一族數值解法;Euler、Runge–Kutta、多步等方法都可以用在遞減的時間格點上。
這種簡單的反向,要和反轉一個 SDE 分開:隨機動態需要額外的分析。書在這裡說下一節會討論,但接下來的 §A.2 只介紹 SDE 本身,沒有談反轉;反向 SDE 在 §4.1(B08)和選讀的 §4.6(B10 第 5 節)。
07代表性的解法書 §A.1.5
上面的分類告訴我們解法從哪裡取得下一步更新的資訊。現在用幾個代表性的例子把這些想法具體化。清單刻意保持簡短:Euler 示範最簡單的單步更新;Runge–Kutta 示範在目前這一步裡使用額外的階段;Adams–Bashforth 示範重複使用過去的歷史;Picard 迭代示範平行於時間方法背後的不動點觀點。以下都令 \(t_{n+1}=t_n+h\),步長 \(h\) 帶正負號:正向積分 \(h\gt0\),反向積分 \(h\lt0\)。
Euler 法
Euler 法是最簡單的單步方法,只用目前區間開頭的速度:
幾何上,Euler 看目前的運動方向,沿這個方向走完整一步。在標準的光滑性假設下,Euler 的局部截斷誤差是 \(\mathcal O(h^2)\),全域誤差是 \(\mathcal O(h)\)。它的簡單使它成為一個有用的一階基準,但當速度在一步之內變化很大時,準確度會變差。
局部 \(\mathcal O(h^2)\) 為什麼變成全域 \(\mathcal O(h)\)(轉述者補充)
局部誤差:精確解的 Taylor 展開 \(\mathbf x(t+h)=\mathbf x(t)+h\mathbf x'(t)+\frac{h^2}{2}\mathbf x''(\xi)\),Euler 只取到第二項,所以一步差 \(\mathcal O(h^2)\)。走完固定長度 \(T\) 需要 \(T/h\) 步,每一步的誤差(在 Lipschitz 條件下最多被放大一個有界的倍數)加起來是 \((T/h)\cdot\mathcal O(h^2)=\mathcal O(h)\)。同樣的道理,二階方法局部 \(\mathcal O(h^3)\)、全域 \(\mathcal O(h^2)\)。所以「階數」指的是全域誤差的冪次:步數加倍,一階方法誤差減半,二階方法誤差變成四分之一(Demo A 的對數圖斜率)。
Runge–Kutta 方法
Euler 只用區間開頭的速度。Runge–Kutta(RK)方法在同一步裡產生的額外階段上計算速度,取得更多資訊。所以它們是單步、多階段的方法:構造 \(\mathbf x_{n+1}\) 用的所有資訊都從目前的起點 \(\mathbf x_n\) 產生,但一步裡可能有好幾次函數計算。下面用兩個二階方法說明基本的 RK 原則。
顯式中點法。先用目前的速度預測區間一半處的狀態,
$$\mathbf x_{\text{mid}}=\mathbf x_n+\frac h2\,\mathbf v(\mathbf x_n,t_n),$$再在這個預測的中點計算速度,用那個方向走完整一步:
$$\mathbf x_{n+1}=\mathbf x_n+h\,\mathbf v\Big(\mathbf x_{\text{mid}},t_n+\frac h2\Big).$$所以中點法不假設起點的速度能代表整個區間,而是先「看進」這一步裡面,再決定怎麼走完它。
Heun 法(改良的 Euler)。Heun 法是另一個兩階段、二階的 RK 方法。它不看中點,而是先用 Euler 法預測終點:
$$\mathbf x_{\text{pred}}=\mathbf x_n+h\,\mathbf v(\mathbf x_n,t_n),$$然後在預測的終點計算速度,把起點和預測終點的速度平均:
$$\mathbf x_{n+1}=\mathbf x_n+\frac h2\Big[\mathbf v(\mathbf x_n,t_n)+\mathbf v(\mathbf x_{\text{pred}},t_n+h)\Big].$$顯式中點法和 Heun 法都用兩個階段,在標準的光滑性假設下局部截斷誤差是 \(\mathcal O(h^3)\)、全域誤差 \(\mathcal O(h^2)\)。它們只差在怎麼收集和組合步內的額外資訊。更高階的 RK 方法用更多階段延伸同一個原則,標準公式這裡省略。在 diffusion 取樣裡,Karras et al. (2022)(EDM)直接用 Heun 法當 ODE 解法。
Heun 為什麼是二階:和 Taylor 展開比對(轉述者補充)
精確解:\(x(t+h)=x+hv+\frac{h^2}{2}\frac{dv}{dt}+\mathcal O(h^3)\),其中 \(\frac{dv}{dt}=\partial_tv+(\partial_xv)v\)(沿軌跡的全導數)。Heun 的第二個速度:\(v(x+hv,t+h)=v+h\,\partial_tv+h(\partial_xv)v+\mathcal O(h^2)=v+h\frac{dv}{dt}+\mathcal O(h^2)\)。代入:\(x+\frac h2\big[v+v+h\frac{dv}{dt}\big]+\mathcal O(h^3)=x+hv+\frac{h^2}{2}\frac{dv}{dt}+\mathcal O(h^3)\),和精確解一致到 \(h^2\),所以局部誤差 \(\mathcal O(h^3)\)。中點法同理:\(v(x+\frac h2v,t+\frac h2)=v+\frac h2\frac{dv}{dt}+\mathcal O(h^2)\),乘上 \(h\) 一樣得到 \(\frac{h^2}{2}\frac{dv}{dt}\)。
Adams–Bashforth 方法
Runge–Kutta 方法透過看目前這一步裡面取得額外資訊。多步方法則是往回看,重複使用先前已完成的步裡的資訊。例如在均勻步長 \(h\) 下,二階 Adams–Bashforth 方法(AB2)是
所以第二份資訊是前一步存下來的速度,而不是目前區間裡的新計算。在標準的光滑性假設和足夠準確的初始化下,AB2 的局部截斷誤差是 \(\mathcal O(h^3)\)、全域誤差 \(\mathcal O(h^2)\)。因為第一步沒有先前的歷史,多步方法必須先用另一個解法初始化。上面的公式假設均勻格點;在非均勻格點上,多步方法的係數依賴相鄰的步長。這個重複使用歷史的原則,後面會在 diffusion 專用的多步解法裡再出現(§9.3–9.5,B26、B27)。
係數 \(\frac32\)、\(-\frac12\) 從哪來(轉述者補充)
精確的一步是 \(\mathbf x_{n+1}=\mathbf x_n+\int_{t_n}^{t_n+h}\mathbf v\,dt\)。用兩個已知點 \((t_{n-1},\mathbf v_{n-1})\)、\((t_n,\mathbf v_n)\) 做直線外插:\(\mathbf v(t)\approx\mathbf v_n+\frac{t-t_n}{h}(\mathbf v_n-\mathbf v_{n-1})\)。在 \([t_n,t_n+h]\) 上積分:\(h\mathbf v_n+\frac h2(\mathbf v_n-\mathbf v_{n-1})=h\big(\frac32\mathbf v_n-\frac12\mathbf v_{n-1}\big)\)。B26 的 DEIS 是同一個做法,只是先把線性部分用指數積分因子解掉、再對非線性部分外插。
Picard 迭代
上面的方法用局部的逐步更新近似軌跡。Picard 迭代走另一條路:透過 ODE 的積分形式,反覆修正整條候選軌跡。從一個滿足 \(\mathbf x^{(0)}(0)=\mathbf x(0)\) 的初始猜測 \(\mathbf x^{(0)}(t)\) 出發,Picard 迭代定義
每次迭代沿前一條候選軌跡計算速度,用這些計算值構造新的一條。這個迭代的不動點滿足原來 ODE 的積分方程。在 Picard–Lindelöf 定理的標準條件下,迭代局部地收斂到 ODE 的解。
對我們後面的目的,重要的觀察是計算上的:在一次 Picard 迭代裡,每個時間的速度只依賴前一條候選軌跡,所以不同時間點的計算可以獨立進行,這開啟了平行於時間計算的可能。這個觀點在 §9.6(ParaDiGMS,B27)再次出現。
Demo A解法實驗室:階數與穩定性轉述者補充
半線性的一維 ODE \(\frac{dx}{dt}=-\lambda x+N(x,t)\),\(N(x,t)=\sin(2x)+2\cos(3t)\),\(x(0)=1\),\(t\in[0,2]\)。精確解用 40,000 步 RK4 算(黑線)。比較五個解法:Euler、中點法、Heun、AB2(第一步用 Heun 起動),以及 exponential Euler(線性部分 \(e^{-\lambda h}\) 精確、\(N\) 在一步裡凍結)。上圖是選定步數 \(M\) 的數值路徑(超出畫面的部分截掉),下圖是 \(t=2\) 的誤差對步數 \(M=5,10,\dots,320\) 的對數圖;某個 \(M\) 下 \(|x|\) 超過 \(10^6\) 就算「發散」,不畫點;誤差超過 10 的點貼在上緣。
讀數(瀏覽器實測,和 Node 原型逐位相同):
- 階數看得到。λ = 1 時,\(M\) 從 80 加倍到 160 再到 320:Euler 的誤差 \(3.9\to1.9\to0.97\)(\(\times10^{-2}\)),每次減半;Heun \(9.7\to2.4\to0.60\)(\(\times10^{-4}\)),每次變成四分之一;中點法、AB2 也是四分之一。對數圖上一階的線斜率 −1,二階的 −2。
- Exponential Euler 在不剛性時沒有階數上的優勢。λ = 1 時它和 Euler 一樣是一階:\(M\) 從 160 到 320,誤差從 \(9.2\times10^{-3}\) 減半成 \(4.6\times10^{-3}\)。常數比 Euler 好一倍(\(M=320\) 時 Euler 是 \(9.7\times10^{-3}\)),但遠不如二階方法。
- 剛性時,穩定性比階數重要。λ = 50 時,四個顯式解法在 \(M=5\)、20(\(\lambda h=20\)、5)全部發散;\(M=80\)(\(\lambda h=1.25\))時 Euler、中點法、Heun 才穩定,AB2 要到 \(M=160\)。Exponential Euler 在 \(M=5\) 就只差 \(3.7\times10^{-2}\),\(M=20\) 時 \(4.8\times10^{-3}\)。這就是正文說的「步長由非線性動態決定,而不是由很大的線性係數決定」。
- 但一旦穩定,exponential Euler 就不是最準的了。λ = 50、\(M=80\)(\(\lambda h=1.25\))時 Euler \(9.0\times10^{-5}\)、中點法 \(6.3\times10^{-5}\)、Heun \(1.3\times10^{-4}\),exponential Euler 反而是 \(7.0\times10^{-4}\);這時還很靠近穩定邊界,Heun 甚至比 Euler 差,二階的好處還沒顯現。到 \(M=320\)(\(\lambda h=0.31\))才看得出階數:Euler \(2.2\times10^{-5}\)、Heun \(4.0\times10^{-6}\)、中點法 \(1.9\times10^{-6}\)、AB2 \(1.5\times10^{-7}\),exponential Euler \(1.5\times10^{-4}\)。所以第 9 章的 DPM-Solver 這類方法把兩者結合:線性部分用積分因子精確處理,非線性部分用高階的外插或多階段。
Demo BPicard 迭代:一次修正整條軌跡轉述者補充
同一個 ODE、λ = 1。在 401 個等距格點上,從常數猜測 \(x^{(0)}(t)\equiv1\) 出發,每次迭代先在所有格點同時算 \(v(x^{(k)}(t_i),t_i)\),再用梯形法累積積分,得到 \(x^{(k+1)}\)。上圖:黑線是精確解,淡藍線是前面的迭代,深藍線是第 \(k\) 次。下圖:和精確解的最大差對迭代次數。
讀數(瀏覽器實測):
- 先慢後快。\(T=2\) 時最大誤差依序約為 1.0、0.41、0.66、0.24、0.23、0.10、0.044、0.012、0.0031、0.00056、0.00010(第 10 次),之後停在 \(3\times10^{-5}\) 左右:那是梯形積分的離散化誤差,不是 Picard 的。前幾次上下起伏、降得不快,然後誤差像 \((LT)^k/k!\) 那樣急遽下降。
- 區間越長,越晚收斂。\(T=4\) 時第 13 次才降到 \(2\times10^{-4}\),之後停在 \(1.6\times10^{-4}\)(同樣是梯形積分的誤差,格點變疏了);\(T=1\) 時第 12 次就到 \(2\times10^{-6}\)。
- 平行的代價與好處。逐步推進要 400 個前後相依的步;Picard 大約 10 次迭代就收斂,每次迭代的 401 個速度計算彼此獨立,可以同時算。總計算量變多(約 10 × 401 次),但「必須一個接一個做」的深度只有約 10。這就是 B27 的 ParaDiGMS 拿多顆 GPU 換時間的原理。
08從 ODE 到 SDE書 §A.2、§A.2.1
隨機微分方程(SDE)是納入隨機性的 ODE 延伸,提供一個建模受不確定性影響之系統的數學框架。這一部分從 ODE 的離散化開始,延伸到 SDE 的離散化,最後討論一般的 SDE,包括 Itô 微積分。
從一個描述狀態變數 \(\mathbf x(t)\in\mathbb R^D\) 確定性演化的 ODE 開始:
$$\frac{d\mathbf x(t)}{dt}=\mathbf f(\mathbf x(t),t),\qquad\mathbf x(0)=\mathbf x_0.$$\(\mathbf f:\mathbb R^D\times[0,T]\to\mathbb R^D\) 是支配 \(\mathbf x(t)\) 動態的時變速度場。這個 ODE 的解是一條光滑的軌跡 \(t\mapsto\mathbf x(t)\),完全由初始條件 \(\mathbf x_0\) 決定。
離散化的觀點
為了建立直覺,考慮這個 ODE 在小時間步 \(\Delta t\) 上的 Euler 離散化:
$$\mathbf x_{t+\Delta t}=\mathbf x_t+\mathbf f(\mathbf x_t,t)\,\Delta t.$$\(\Delta t\to0\) 時這個近似越來越準,(在 \(\mathbf f\) 的標準正則性條件下)收斂到 ODE 的精確解。
引入隨機性:從 ODE 到 SDE
在許多真實系統裡,完全知道動態是不切實際的。雜訊、不確定性或未建模的交互作用都可能影響演化。為了納入這種隨機性,在 ODE 上加一個隨機項:
- \(g:[0,T]\to\mathbb R\) 是擴散係數(一般可以同時依賴狀態和時間,這裡假設只依賴時間);
- \(\bm\epsilon_t\sim\mathcal N(\mathbf 0,\mathbf I_D)\) 是獨立同分佈的標準高斯向量。
這個修改過的更新規則不只反映確定性的漂移,也反映以 \(\sqrt{\Delta t}\) 縮放的隨機擾動。這個縮放確保隨機擾動在 \(\Delta t\to0\) 的極限下保持有限。重要的是,這個表述在 \(\Delta t\to0\) 時產生一個連續時間的隨機過程,帶我們進入 SDE 的框架。
為什麼一定是 \(\sqrt{\Delta t}\)(轉述者補充)
只看雜訊(\(\mathbf f=0\)、\(g\) 為常數),走完 \([0,T]\) 要 \(T/\Delta t\) 步,每步加一個變異數為 \(g^2c^2\) 的獨立高斯(\(c\) 是雜訊的縮放)。獨立的變異數相加,總變異數是 \(\frac{T}{\Delta t}\cdot g^2c^2\)。
- \(c=\sqrt{\Delta t}\):總變異數 \(g^2T\),和 \(\Delta t\) 無關,極限存在。
- \(c=\Delta t\):總變異數 \(g^2T\Delta t\to0\),雜訊在極限中消失,又變回 ODE。
- \(c=1\):總變異數 \(g^2T/\Delta t\to\infty\),極限不存在。
只有 \(\sqrt{\Delta t}\) 給出一個非平凡的極限。代價是路徑在 \(\Delta t\) 時間內移動 \(\sim\sqrt{\Delta t}\),比 \(\Delta t\) 大得多,所以路徑處處不可微。Demo C 可以換這三種縮放。
隨機微分方程
形式上,上面的離散更新在 \(\Delta t\to0\) 時的極限定義了 SDE:
這裡 \(\mathbf w(t)\in\mathbb R^D\) 是一個 Wiener 過程(標準布朗運動),一個由下面性質刻畫的連續時間隨機過程:
- 初始狀態:\(\mathbf w(0)=0\) 幾乎必然成立;
- 獨立增量:對 \(0\le s\lt t\),增量 \(\mathbf w(t)-\mathbf w(s)\) 和過去獨立;
- 高斯增量:\(\mathbf w(t)-\mathbf w(s)\sim\mathcal N\big(\mathbf 0,(t-s)\mathbf I_D\big)\);
- 連續性:樣本路徑 \(t\mapsto\mathbf w(t)\) 幾乎必然連續,但處處不可微。
此外,常用 \(d\mathbf w(t):=\mathbf w(t+dt)-\mathbf w(t)\) 表示 Wiener 過程的無窮小增量。這個記號雖然有啟發性,但不應該解讀成古典的微分(Riemann 或 Lebesgue 意義下),因為布朗路徑幾乎必然處處不可微。它只是表達高斯增量性質的形式簡寫:
$$d\mathbf w(t)\sim\mathcal N(\mathbf 0,dt\,\mathbf I_D),$$意思是在長度為 \(dt\) 的無窮小時間區間上,Wiener 過程的增量像一個平均為零、共變異為 \(dt\,\mathbf I_D\) 的高斯隨機變數。(對照上面:\(\sqrt{\Delta t}\,\bm\epsilon_t\) 正好是 \(\mathcal N(\mathbf 0,\Delta t\,\mathbf I_D)\)。)
09積分形式與 Itô 積分書 §A.2.2
上面的 SDE 應該用它的積分形式理解:
以 Itô 的意義解讀。第一項是古典的(Riemann 或 Lebesgue)積分,代表累積的確定性漂移;第二項是 Itô 隨機積分,對 Wiener 過程 \(\mathbf w(t)\) 積分。這裡不完整嚴謹地構造 Itô 積分,只提供下面的直覺。
Itô 積分的直覺
Itô 積分可以理解成下面這種離散和的極限:
$$\sum_ig(t_i)\big(w(t_{i+1})-w(t_i)\big),$$被積函數 \(g(t)\) 在每個小區間的左端點 \(t_i\) 取值。分割越來越細時,這些和收斂到 Itô 積分。左端點取值的規則是 Itô 積分和古典積分的區別,古典積分可以用中點或其他取值規則。因為布朗路徑連續但幾乎必然處處不可微,古典積分不能直接套用。Itô 積分處理這種不規則性,描述隨機起伏隨時間的累積效應。
書的註腳:對平方可積的被積函數,這個收斂在 \(L^2\)(均方)意義下成立,進而推出依機率收斂。\(L^2\) 的觀點很自然,因為它直接連到 Itô isometry,也就是 Itô 積分標準而嚴謹的構造背後的關鍵恆等式。這個 isometry 把隨機積分和期望中的標準積分連起來:\(\mathbb E\big[\big\|\int_0^T\bm\psi(t)\,d\mathbf w_t\big\|^2\big]=\mathbb E\big[\int_0^T\|\bm\psi(t)\|_F^2\,dt\big]\)(\(\|\cdot\|_F\) 是 Frobenius 範數);附錄 C(B35)的 §C.1.5 會敘述並使用它。
左端點為什麼重要:\(\int_0^Tw\,dw\) 的例子(轉述者補充)
令 \(\Delta w_i=w(t_{i+1})-w(t_i)\)。代數恆等式 \(w_{i+1}^2-w_i^2=2w_i\Delta w_i+(\Delta w_i)^2\) 對 \(i\) 加總:
$$\sum_iw(t_i)\,\Delta w_i=\frac12\Big(w(T)^2-\sum_i(\Delta w_i)^2\Big).$$\(\sum_i(\Delta w_i)^2\) 的期望是 \(\sum_i\Delta t=T\),變異數 \(\sum_i2\Delta t^2\to0\):分割變細時它確定地趨近 \(T\)(二次變分)。所以 Itô 積分 \(\int_0^Tw\,dw=\frac12(w(T)^2-T)\),期望為 0。如果改用左右端點的平均 \(\frac12(w_i+w_{i+1})\Delta w_i\),和式正好是 \(\frac12w(T)^2\),期望 \(T/2\)。對光滑函數取哪個點都一樣,對布朗路徑卻差了一個確定的 \(T/2\):這就是為什麼 SDE 必須說明用哪一種積分,也是附錄 C 的 Itô 公式多出一個二階項的原因。Demo C 的讀數會看到 0 和 0.5。
微分記號的用法
\(d\mathbf x(t)\)、\(dt\)、\(d\mathbf w(t)\) 這些式子不是古典的微分,而是代表各個過程無窮小增量的形式記號。雖然是啟發式的,它們因為方便而被廣泛使用:能用類似 ODE 的方式寫 SDE,也方便在 Itô 微積分裡做形式上的運算。Itô 微積分在 diffusion model 裡怎麼用,在附錄 C(B35)說明。
和 ODE 比較
在 ODE(例如 \(\frac{d\mathbf x}{dt}=\mathbf f(\mathbf x,t)\))裡,積分形式 \(\mathbf x(t)=\mathbf x(0)+\int_0^t\mathbf f(\mathbf x(\tau),\tau)\,d\tau\) 由微積分基本定理保證:可微函數可以從它的導數還原。相反地,在 SDE 裡沒有這個定理的直接類比,因為布朗運動不可微,隨機積分也不遵守古典的連鎖律。Itô 微積分改用別的工具(例如 Itô 引理)來分析和操作隨機動態。所以雖然 SDE 的微分記號精簡又直觀,嚴謹的理解要靠用 Itô 積分解讀它的積分形式。
10SDE 的數值解法書 §A.2.3
和 ODE 一樣,如果 \(\mathbf f(\cdot,t)\) 和 \(g(\cdot)\) 滿足一些光滑性條件,SDE 就有唯一解:\(\mathbf f(\cdot,t)\) 對 \(\mathbf x\) 是 Lipschitz 且線性成長,\(g(\cdot)\) 平方可積。書的註腳:這裡的解是強(strong)意義下的,也就是在一個固定的機率空間上,\(\mathbf x(t)\) 對給定的布朗運動 \(\mathbf w(t)\) 滿足 SDE 的積分形式;詳細的技術定義省略。
一般的 SDE 通常沒有閉式解,所以需要數值方法。常用的是 Euler–Maruyama 法,它推廣了 ODE 的 Euler 法;其實我們已經在第 8 節的更新式裡看過它:在時間步 \(\Delta t\) 上近似漂移項 \(\mathbf f(\mathbf x(t),t)\),並用高斯增量 \(\sqrt{\Delta t}\,\bm\epsilon_t\)(\(\bm\epsilon_t\sim\mathcal N(\mathbf 0,\mathbf I)\))模擬隨機雜訊 \(g(t)\,d\mathbf w(t)\):
之後在附錄 C 的 §C.1.5(B35)會看到,線性的 SDE 有閉式解。
轉述者註:DDPM 的取樣步驟(B04)、score SDE 的反向 SDE 取樣(B09)都是這個更新式的變形;反向時間時 \(\Delta t\) 取負號,雜訊仍然用 \(\sqrt{|\Delta t|}\)。
Demo CEuler–Maruyama 與左端點轉述者補充
Ornstein–Uhlenbeck 過程 \(dx=-x\,dt+\sqrt2\,dw\),\(x(0)=2\)。它是線性 SDE,精確解在時間 \(t\) 是高斯,平均 \(2e^{-t}\)、變異數 \(1-e^{-2t}\)(附錄 C 會推出來)。用 Euler–Maruyama 跑 4,000 條路徑到 \(t=2\)(固定種子)。上圖:12 條路徑,黑線是精確的平均和平均 ± 2 標準差。下圖:\(x(2)\) 的直方圖和精確密度。「雜訊縮放」可以把 \(\sqrt{\Delta t}\) 換成 \(\Delta t\) 或 1,看看會發生什麼。讀數另外用 4,000 條布朗路徑(每條 200 段)算 \(\int_0^1w\,dw\) 的兩種離散和。
讀數(瀏覽器實測,和 Node 原型逐位相同)。\(x(2)\) 的精確平均 0.271、變異數 0.982:
| Δt | √Δt(正確) | Δt | 1 |
|---|---|---|---|
| 0.5 | 0.126 / 1.348 | 0.126 / 0.674 | 0.126 / 2.696 |
| 0.2 | 0.208 / 1.099 | 0.212 / 0.220 | 0.201 / 5.495 |
| 0.05 | 0.240 / 0.998 | 0.253 / 0.050 | 0.181 / 19.97 |
| 0.01 | 0.250 / 0.968 | 0.266 / 0.010 | 0.086 / 96.79 |
每格是 4,000 條路徑的樣本平均 / 樣本變異數。
- 只有 \(\sqrt{\Delta t}\) 收斂到正確的變異數。用 \(\Delta t\) 縮放,變異數隨 \(\Delta t\) 一起趨近 0(0.674 → 0.010),雜訊消失;用 1,變異數像 \(1/\Delta t\) 一樣爆掉(2.7 → 96.8)。三種縮放的平均都差不多,因為雜訊平均為零;出錯的是分佈的寬度。(差最多的是縮放 1、\(\Delta t=0.01\) 的 0.086,但那時變異數 96.8,4,000 條路徑的平均標準誤約 0.16,還在抽樣誤差內。)
- 大步長時 EM 本身有偏差。\(\Delta t=0.5\) 時變異數 1.348,遠大於 0.982,但它和「離散 EM 鏈本身」的精確變異數 1.328 一致(讀數裡另外列出):這是一階方法的離散化誤差,不是抽樣誤差。\(\Delta t\) 變小,EM 鏈的值逐漸接近連續 SDE(0.01 時 0.987)。
- 左端點和左右平均差了 \(T/2\)。\(\int_0^1w\,dw\) 的左端點(Itô)和平均是 0.009,左右平均版本是 0.509:差了 0.5,和第 9 節的推導一致。Itô 和的平均也和 \((w(1)^2-1)/2\) 的平均(0.009)相同。