测试函数法推导连续性方程和Fokker-Planck方程
By 苏剑林 | 2023-02-11 | 42163位读者 |在文章《生成扩散模型漫谈(六):一般框架之ODE篇》中,我们推导了SDE的Fokker-Planck方程;而在《生成扩散模型漫谈(十二):“硬刚”扩散ODE》中,我们单独推导了ODE的连续性方程。它们都是描述随机变量沿着SDE/ODE演化的分布变化方程,连续性方程是Fokker-Planck方程的特例。在推导Fokker-Planck方程时,我们将泰勒展开硬套到了狄拉克函数上,虽然结果是对的,但未免有点不伦不类;在推导连续性方程时,我们结合了雅可比行列式和泰勒展开,方法本身比较常规,但没法用来推广到Fokker-Planck方程。
这篇文章我们介绍“测试函数法”,它是推导连续性方程和Fokker-Planck方程的标准方法之一,其分析过程比较正规,并且适用场景也比较广。
分部积分 #
正式推导之前,这里先介绍后面推导会用到的关键结果——分部积分法的高维推广。
一般教程对分部积分的介绍仅限于一维情形,即
∫bauv′dx=uv|ba−∫bavu′dx
这里u,v是x的函数,′表示求函数关于x的导数。下面我们需要它的一个高维版本,为此,我们先回顾一维分部积分的推导过程,它依赖于求导的乘法法则:
(uv)′=uv′+vu′
然后两边对x积分并移项,就得到分部积分公式。对于高维情形,我们考虑类似的公式:
\begin{equation}\nabla\cdot(u\boldsymbol{v}) = \boldsymbol{v}\cdot\nabla u + u\nabla \cdot\boldsymbol{v}\end{equation}
其中u是\boldsymbol{x}的标量函数,\boldsymbol{v}是\boldsymbol{x}的向量函数(维度跟\boldsymbol{v}一致),\nabla表示求函数关于\boldsymbol{x}的梯度。现在我们对两端在区域\Omega积分:
\begin{equation}\int_{\Omega}\nabla\cdot(u\boldsymbol{v})d\boldsymbol{x} = \int_{\Omega}\boldsymbol{v}\cdot\nabla u d\boldsymbol{x} + \int_{\Omega}u\nabla \cdot\boldsymbol{v} d\boldsymbol{x}\end{equation}
根据高斯散度定理,左侧等于\int_{\partial\Omega}u\boldsymbol{v}\cdot\hat{\boldsymbol{n}}dS,\partial\Omega是\Omega的边界,\hat{\boldsymbol{n}}是边界的外向单位法向量,dS是面积微元。所以,移项后有
\begin{equation}\int_{\Omega}\boldsymbol{v}\cdot\nabla u d\boldsymbol{x} = \int_{\partial\Omega}u\boldsymbol{v}\cdot\hat{\boldsymbol{n}}dS - \int_{\Omega}u\nabla \cdot\boldsymbol{v} d\boldsymbol{x}\label{eq:int-by-parts}\end{equation}
这就是我们要推导的高维空间分部积分公式。特别地,对于概率密度函数p,那么由于非负性和积分为1的限制,无穷远处必然有p\to 0和\nabla p\to \boldsymbol{0},所以如果\Omega选为全空间(没有特别注明积分区域的,默认为全空间),那么分别将u=p和\boldsymbol{v}=\nabla p代入上式,得到
\begin{align}\int\boldsymbol{v}\cdot\nabla p d\boldsymbol{x} =&\, - \int p\nabla \cdot\boldsymbol{v} d\boldsymbol{x}\label{eq:int-by-parts-p} \\
\int u\nabla \cdot\nabla p d\boldsymbol{x} = &\,-\int\nabla p\cdot\nabla u d\boldsymbol{x}\label{eq:int-by-parts-gp}\end{align}
如果要进一步严格化上述结论,可以假设p具有紧的支撑集。不过这纯粹是数学上的严格化,事实上对于一般理解来说,直接默认在无穷远处成立p\to 0和\nabla p\to \boldsymbol{0}就够了。
ODE演化 #
测试函数法的原理,是如果对于任意函数\phi(\boldsymbol{x}),都成立
\begin{equation}\int f(\boldsymbol{x})\phi(\boldsymbol{x})d\boldsymbol{x} = \int g(\boldsymbol{x})\phi(\boldsymbol{x})d\boldsymbol{x}\end{equation}
那么就成立f(\boldsymbol{x})=g(\boldsymbol{x}),其中\phi(\boldsymbol{x})就叫做测试函数。更严谨的定义需要声明\phi(\boldsymbol{x})的选取空间,以及等号的具体含义(如严格相等/几乎处处相等/依概率相等之类),这里我们就不引入这些细节了。
对于ODE
\begin{equation}\frac{d\boldsymbol{x}_t}{dt}=\boldsymbol{f}_t(\boldsymbol{x}_t)\label{eq:ode}\end{equation}
我们将它离散化为
\begin{equation}\boldsymbol{x}_{t+\Delta t} = \boldsymbol{x}_t + \boldsymbol{f}_t(\boldsymbol{x}_t)\Delta t\label{eq:ode-diff}\end{equation}
那么就有
\begin{equation}\phi(\boldsymbol{x}_{t+\Delta t}) = \phi(\boldsymbol{x}_t + \boldsymbol{f}_t(\boldsymbol{x}_t)\Delta t)\approx \phi(\boldsymbol{x}_t) + \Delta t\,\,\boldsymbol{f}_t(\boldsymbol{x}_t)\cdot\nabla_{\boldsymbol{x}_t}\phi(\boldsymbol{x}_t)\end{equation}
两边求期望,得到:
\begin{equation}\int p_{t+\Delta t}(\boldsymbol{x}_{t+\Delta t})\phi(\boldsymbol{x}_{t+\Delta t}) d\boldsymbol{x}_{t+\Delta t}\approx \int p_t(\boldsymbol{x}_t)\phi(\boldsymbol{x}_t)d\boldsymbol{x}_t + \Delta t\int p_t(\boldsymbol{x}_t)\boldsymbol{f}_t(\boldsymbol{x}_t)\cdot\nabla_{\boldsymbol{x}_t}\phi(\boldsymbol{x}_t)d\boldsymbol{x}_t\end{equation}
由于积分的结果不依赖于被积自变量的记号,所以左边将\boldsymbol{x}_{t+\Delta t}换成\boldsymbol{x}_t也是等价的:
\begin{equation}\int p_{t+\Delta t}(\boldsymbol{x}_t)\phi(\boldsymbol{x}_t) d\boldsymbol{x}_t\approx \int p_t(\boldsymbol{x}_t)\phi(\boldsymbol{x}_t)d\boldsymbol{x}_t + \Delta t\int p_t(\boldsymbol{x}_t)\boldsymbol{f}_t(\boldsymbol{x}_t)\cdot\nabla_{\boldsymbol{x}_t}\phi(\boldsymbol{x}_t)d\boldsymbol{x}_t\label{eq:change-var}\end{equation}
将右边第一项移到左边,然后取\Delta t\to 0的极限,得到:
\begin{equation}\int \frac{\partial p_t(\boldsymbol{x}_t)}{\partial t}\phi(\boldsymbol{x}_t) d\boldsymbol{x}_t = \int p_t(\boldsymbol{x}_t)\boldsymbol{f}_t(\boldsymbol{x}_t)\cdot\nabla_{\boldsymbol{x}_t}\phi(\boldsymbol{x}_t)d\boldsymbol{x}_t\label{eq:dt-0}\end{equation}
右端利用分部积分公式\eqref{eq:int-by-parts-p}得到
\begin{equation}\int \frac{\partial p_t(\boldsymbol{x}_t)}{\partial t}\phi(\boldsymbol{x}_t) d\boldsymbol{x}_t = -\int \Big[\nabla_{\boldsymbol{x}_t}\cdot\big(p_t(\boldsymbol{x}_t)\boldsymbol{f}_t(\boldsymbol{x}_t)\big)\Big]\phi(\boldsymbol{x}_t)d\boldsymbol{x}_t\end{equation}
根据测试函数法的相等原理,就有
\begin{equation}\frac{\partial p_t(\boldsymbol{x}_t)}{\partial t} = -\nabla_{\boldsymbol{x}_t}\cdot\big(p_t(\boldsymbol{x}_t)\boldsymbol{f}_t(\boldsymbol{x}_t)\big)\end{equation}
这称为“连续性方程”。
SDE演化 #
对于SDE
\begin{equation}d\boldsymbol{x}_t = \boldsymbol{f}_t(\boldsymbol{x}_t) dt + g_t d\boldsymbol{w}\label{eq:sde}\end{equation}
我们离散化为
\begin{equation}\boldsymbol{x}_{t+\Delta t} = \boldsymbol{x}_t + \boldsymbol{f}_t(\boldsymbol{x}_t) \Delta t + g_t \sqrt{\Delta t}\boldsymbol{\varepsilon},\quad \boldsymbol{\varepsilon}\sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{I})\label{eq:sde-diff}\end{equation}
那么
\begin{equation}\begin{aligned}
\phi(\boldsymbol{x}_{t+\Delta t}) =&\, \phi(\boldsymbol{x}_t + \boldsymbol{f}_t(\boldsymbol{x}_t) \Delta t + g_t \sqrt{\Delta t}\boldsymbol{\varepsilon}) \\
\approx&\, \phi(\boldsymbol{x}_t) + \left(\boldsymbol{f}_t(\boldsymbol{x}_t) \Delta t + g_t \sqrt{\Delta t}\boldsymbol{\varepsilon}\right)\cdot \nabla_{\boldsymbol{x}_t}\phi(\boldsymbol{x}_t) + \frac{1}{2} \left(g_t\sqrt{\Delta t}\boldsymbol{\varepsilon}\cdot \nabla_{\boldsymbol{x}_t}\right)^2\phi(\boldsymbol{x}_t)
\end{aligned}\end{equation}
两边求期望,注意右边要同时对\boldsymbol{x}_t和\boldsymbol{\varepsilon}求期望,其中\boldsymbol{\varepsilon}的期望可以事先求出,结果是
\begin{equation}\phi(\boldsymbol{x}_t) + \Delta t\,\,\boldsymbol{f}_t(\boldsymbol{x}_t)\cdot \nabla_{\boldsymbol{x}_t}\phi(\boldsymbol{x}_t) + \frac{1}{2} \Delta t\,g_t^2\nabla_{\boldsymbol{x}_t}\cdot\nabla_{\boldsymbol{x}_t}\phi(\boldsymbol{x}_t)
\end{equation}
于是
\begin{equation}\begin{aligned}
&\,\int p_{t+\Delta t}(\boldsymbol{x}_{t+\Delta t})\phi(\boldsymbol{x}_{t+\Delta t}) d\boldsymbol{x}_{t+\Delta t}\\
\approx&\, \int p_t(\boldsymbol{x}_t)\phi(\boldsymbol{x}_t)d\boldsymbol{x}_t + \Delta t\int p_t(\boldsymbol{x}_t)\boldsymbol{f}_t(\boldsymbol{x}_t)\cdot\nabla_{\boldsymbol{x}_t}\phi(\boldsymbol{x}_t)d\boldsymbol{x}_t + \int\frac{1}{2} \Delta t\,g_t^2 p_t(\boldsymbol{x}_t)\nabla_{\boldsymbol{x}_t}\cdot\nabla_{\boldsymbol{x}_t}\phi(\boldsymbol{x}_t) d\boldsymbol{x}_t
\end{aligned}\end{equation}
跟式\eqref{eq:change-var}、式\eqref{eq:dt-0}类似,取\Delta\to 0的极限,得到
\begin{equation}\int \frac{\partial p_t(\boldsymbol{x}_t)}{\partial t}\phi(\boldsymbol{x}_t) d\boldsymbol{x}_t = \int p_t(\boldsymbol{x}_t)\boldsymbol{f}_t(\boldsymbol{x}_t)\cdot\nabla_{\boldsymbol{x}_t}\phi(\boldsymbol{x}_t)d\boldsymbol{x}_t + \int\frac{1}{2} \,g_t^2 p_t(\boldsymbol{x}_t)\nabla_{\boldsymbol{x}_t}\cdot\nabla_{\boldsymbol{x}_t}\phi(\boldsymbol{x}_t) d\boldsymbol{x}_t\end{equation}
对右边第一项应用式\eqref{eq:int-by-parts-p}、对右边第二项先应用式\eqref{eq:int-by-parts-gp}再应用式\eqref{eq:int-by-parts-p},得到
\begin{equation}\int \frac{\partial p_t(\boldsymbol{x}_t)}{\partial t}\phi(\boldsymbol{x}_t) d\boldsymbol{x}_t = \int \left[-\nabla_{\boldsymbol{x}_t}\cdot\big(p_t(\boldsymbol{x}_t)\boldsymbol{f}_t(\boldsymbol{x}_t)\big)+\frac{1}{2}g_t^2 \nabla_{\boldsymbol{x}}\cdot\nabla_{\boldsymbol{x}}p_t(\boldsymbol{x})\right]\phi(\boldsymbol{x}_t)d\boldsymbol{x}_t\end{equation}
根据测试函数法的相等原理得
\begin{equation}\frac{\partial p_t(\boldsymbol{x}_t)}{\partial t} = -\nabla_{\boldsymbol{x}_t}\cdot\big(p_t(\boldsymbol{x}_t)\boldsymbol{f}_t(\boldsymbol{x}_t)\big)+\frac{1}{2}g_t^2 \nabla_{\boldsymbol{x}}\cdot\nabla_{\boldsymbol{x}}p_t(\boldsymbol{x})\end{equation}
这就是“Fokker-Planck方程”。
文章小结 #
本文介绍了用于推导某些概率方程的测试函数法,主要内容包括分部积分法的高维推广,以及ODE的连续性方程和SDE的Fokker-Planck方程的推导过程。
转载到请包括本文地址:https://spaces.ac.cn/archives/9461
更详细的转载事宜请参考:《科学空间FAQ》
如果您还有什么疑惑或建议,欢迎在下方评论区继续讨论。
如果您觉得本文还不错,欢迎分享/打赏本文。打赏并非要从中获得收益,而是希望知道科学空间获得了多少读者的真心关注。当然,如果你无视它,也不会影响你的阅读。再次表示欢迎和感谢!
如果您需要引用本文,请参考:
苏剑林. (Feb. 11, 2023). 《测试函数法推导连续性方程和Fokker-Planck方程 》[Blog post]. Retrieved from https://spaces.ac.cn/archives/9461
@online{kexuefm-9461,
title={测试函数法推导连续性方程和Fokker-Planck方程},
author={苏剑林},
year={2023},
month={Feb},
url={\url{https://spaces.ac.cn/archives/9461}},
}
March 15th, 2023
您好,我认为应当注明“特别地,对于概率密度函数p,那么由于非负性和积分为1的限制,无穷远处必然有p\rightarrow 0和\nabla p\rightarrow 0是对扩散过程成立。对一般的概率密度函数是不成立的,例如,考虑\href{https://math.stackexchange.com/a/847280}{p的极限不存在的情况}。
严格来说,你是对的,但作为一篇科普文章来说,我们一般只针对形态比较良好的case做介绍,而不考虑刻意构造的特殊case。因为要严格化的话,要补充的内容可就多了,比如测试函数法本身成立的条件,相等是哪种相等,等等,这些对科普是不利的,而且对多数机器学习场景来说也没有太多增益。
感觉可以声明 \phi的支撑集是紧的,这样可以兼顾严谨性和科普性
June 27th, 2023
苏老师,请问第19式最后一项是怎么来的呀?
就是多元泰勒展开,二阶项不就是这样么?
括号的位置是不是写错了呀?
没写错,算子进行平方运算。
错了,不是\phi(x_t),而是\phi(x_t)关于x_t的二阶导数
仔细看了一下,是对的,只是个记法的问题
是的,算子的平方运算,很多时候都是这样记的。
November 16th, 2023
您好,请问从公式5到公式6和7的过程中,公式5的右手边第一项\int_{\partial \Omega}u\textbf{v}\cdot \hat{n} dS积分为什么为0?这个该如何理解?
简单来说,因为\partial\Omega意味着在边界进行积分,这里的边界是无穷远处,而如果u=p是概率密度函数,那么在无穷远处必然为0,因此积分也是0。当然这样说还很不严谨,但严格化也是可以的,这里就不展开了。
January 17th, 2024
苏老师好,请问公式12中,左右求期望时,左边使用 E_{p_{t+\Delta t}},右边使用 E_{p_{t}},这样为什么成立?按道理左右应该使用同一个分布?期待回复!
想了想,是不是因为 \Delta t\to 0, p(t+\Delta t)\approx p(t)+\Delta t \frac{\partial p}{\partial t}= p(t)?
因为两个分布之间是通过一个确定性变换进行转化的,确定性变换之下期望可以直接转化,比如y=x^2,那么\mathbb{E}_y[f(y)] = \mathbb{E}_x[f(x^2)],写成积分形式就是
\int p(y) f(y)dy = \int q(x) f(x^2) dx
其中p,q分别是y,x的概率密度。
嗯嗯,后来意识到了这个事情,还是多谢苏老师解答!
April 7th, 2024
Eq.(21) 最右项少了个 \int 符号
谢谢,已补充
April 7th, 2024
文中在推导(13)和(22)时用 \boldsymbol{x}_t 替换 \boldsymbol{x}_{t+\Delta t},这一步不太理解,因为替换之后意味着:\int p_{t+\Delta t}(\boldsymbol{x}_{t+\Delta t})\phi(\boldsymbol{x}_{t+\Delta t}) d\boldsymbol{x}_{t+\Delta t}\equiv\int p_{t+\Delta t}(\boldsymbol{x}_t)\phi(\boldsymbol{x}_t) d\boldsymbol{x}_t
如果只是符号替换,为何 p_{t+\Delta t}没有替换为 p_{t}?烦请苏神解惑。
你就说你自己写的这个等式成不成立吧?
\boldsymbol{x}_{t+\Delta t}是一个整体记号,p_{t+\Delta t}的t+\Delta t是另一个记号,我要替换的是\boldsymbol{x}_{t+\Delta t},关p_{t+\Delta t}的t+\Delta t什么事?
确实不好意思,当时写问题之后就明白了,但是忘了回复了,又没法删除,谢谢回复。
October 25th, 2024
他还是把公式写的如此清晰。