用深度学习求解高维偏微分方程

用深度学习求解高维偏微分方程简述：我们一般求解的传统的PDE维数也就二维三维，但是在其他领域中，比如说金融学，通过数学建模构建出来的PDE的维数及其之高，维数升高带来的“维数灾难”问题亟待解决。这篇文章，介绍了一种基于深度学习的方法，可以处理一般的高维抛物型方程。这篇文章，用反向随机微分方程构造PDE，并用神经网络近似未知解的梯度，来求解高维的偏微分方程。高维PDE的来源：1、量子多体问...

陆嵩

23911人浏览 · 2018-11-08 20:18:50

陆嵩 · 2018-11-08 20:18:50 发布

用深度学习求解高维偏微分方程

简述：

我们一般求解的传统的PDE维数也就二维三维，但是在其他领域中，比如说金融学，通过数学建模构建出来的PDE的维数及其之高，维数升高带来的“维数灾难”问题亟待解决。
这篇文章，介绍了一种基于深度学习的方法，可以处理一般的高维抛物型方程。这篇文章，用反向随机微分方程构造PDE，并用神经网络近似未知解的梯度，来求解高维的偏微分方程。

高维PDE的来源：

1、量子多体问题的薛定谔方程，在这种情况下，PDE的维数大约是系统中电子或者量子（粒子）的三倍。
2、用于给金融衍生品定价的Black-Scholes方程，其中PDE的维数是所考虑的相关金融资产的数量。
3、动态规划中的Hamilton-Jacobi-Bellman方程。在具有多个代理的博弈论中，维数随着代理的数量线性增加。类似，在资源分配问题中，维数随着设备和资源的数量线性增加。

解决这些问题的计算成本随维度呈指数增长！！！！
为什么深度学习能表示任意函数？万有逼近定理。
基本想法就是随机微分方程中某个梯度算子用神经网络来表示。

方法

我们考虑一类半线性抛物型偏微分方程，长成这个样子：

$∂u∂t(t,x)+12Tr⁡(σσT(t,x)(Hess⁡xu)(t,x))+∇u(t,x)⋅μ(t,x)+f(t,x,u(t,x),σT(t,x)∇u(t,x))=0\frac{\partial u}{\partial t}(t, x)+\frac{1}{2} \operatorname{Tr}\left(\sigma \sigma^{\mathrm{T}}(t, x)\left(\operatorname{Hess}_{x} u\right)(t, x)\right)+\nabla u(t, x) \cdot \mu(t, x)+f\left(t, x, u(t, x), \sigma^{\mathrm{T}}(t, x) \nabla u(t, x)\right)=0$

且有指定的末端条件 $u (T, x) = g (x)$ ，这里的x是d维的， $μ\mu$ 是已知的向量值函数， $σ\sigma$ 是d维矩阵值函数，f是非线性函数。我们感兴趣的是 $t=0,x=ξt=0,x=\xi$ 的解。

让 $X_t$ 是随机变量，满足：
$Xt=ξ+∫0tμ(s,Xs)ds+∫0tσ(s,Xs)dWsX_{t}=\xi+\int_{0}^{t} \mu\left(s, X_{s}\right) d s+\int_{0}^{t} \sigma\left(s, X_{s}\right) d W_{s}$

那么有几篇论文表明了，(1)的解长成下面这样：

$u(t,Xt)=u(0,X0)−∫0tf(s,Xs,u(s,Xs),σT(s,Xs)∇u(s,Xs))ds+∫0t[∇u(s,Xs)]Tσ(s,Xs)dWsu\left(t, X_{t}\right)=u\left(0, X_{0}\right)-\int_{0}^{t} f\left(s, X_{s}, u\left(s, X_{s}\right), \sigma^{\mathrm{T}}\left(s, X_{s}\right) \nabla u\left(s, X_{s}\right)\right) d s+\int_{0}^{t}\left[\nabla u\left(s, X_{s}\right)\right]^{\mathrm{T}} \sigma\left(s, X_{s}\right) d W_{s}$

我们对（2）、（3）两个式子用欧拉格式离散，得到：
$Xtn+1−Xtn≈μ(tn,Xtn)Δtn+σ(tn,Xtn)ΔWnX_{t_{n+1}}-X_{t_{n}} \approx \mu\left(t_{n}, X_{t_{n}}\right) \Delta t_{n}+\sigma\left(t_{n}, X_{t_{n}}\right) \Delta W_{n}$
$u(tn+1,Xtn+1)≈u(tn,Xtn)−f(tn,Xtn,u(tn,Xtn),σT(tn,Xtn)∇u(tn,Xtn))(tn+1−tn)+[∇u(tn,Xtn)]Tσ(tn,Xtn)(Wtn+1−Wtn)\begin{aligned} u\left(t_{n+1}, X_{t_{n+1}}\right) \approx & u\left(t_{n}, X_{t_{n}}\right)-f\left(t_{n}, X_{t_{n}}, u\left(t_{n}, X_{t_{n}}\right), \sigma^{\mathrm{T}}\left(t_{n}, X_{t_{n}}\right) \nabla u\left(t_{n}, X_{t_{n}}\right)\right)\left(t_{n+1}-t_{n}\right) \\ &+\left[\nabla u\left(t_{n}, X_{t_{n}}\right)\right]^{\mathrm{T}} \sigma\left(t_{n}, X_{t_{n}}\right)\left(W_{t_{n+1}}-W_{t_{n}}\right) \end{aligned}$
这里，
$Δtn=tn+1−tn,ΔWn=Wtn+1−Wtn\Delta t_{n}=t_{n+1}-t_{n}, \quad \Delta W_{n}=W_{t_{n+1}}-W_{t_{n}}$

这里我有个想法，能不能用其他格式做离散，比如说Milstein格式等。

接下来是关键的步骤，就是在数值格式中，用多层反馈神经网络去逼近 $σ∇u\sigma \nabla u$ 项，即：
$σT(tn,Xtn)∇u(tn,Xtn)=(σT∇u)(tn,Xtn)≈(σT∇u)(tn,Xtn∣θn)\sigma^{\mathrm{T}}\left(t_{n}, X_{t_{n}}\right) \nabla u\left(t_{n}, X_{t_{n}}\right)=\left(\sigma^{\mathrm{T}} \nabla u\right)\left(t_{n}, X_{t_{n}}\right) \approx\left(\sigma^{\mathrm{T}} \nabla u\right)\left(t_{n}, X_{t_{n}} \mid \theta_{n}\right)$

注意到（5）这里每一步都放进去一个神经网络，最后整个堆叠成一个函数表达器。这里边的 $X_{t_n},W_{t_n}$ 是变量。我们可以用一般的深度学习方法来训练这个网络，如随机梯度下降方法等。

还记得我们前面给了一个末端条件，刚好可以用了定义损失函数。
$l(θ)=E[∣g(XtN)−u^({Xtn}0≤n≤N,{Wtn}0≤n≤N)∣2]l(\theta)=\mathbb{E}\left[\left|g\left(X_{t_{N}}\right)-\hat{u}\left(\left\{X_{t_{n}}\right\}_{0 \leq n \leq N},\left\{W_{t_{n}}\right\}_{0 \leq n \leq N}\right)\right|^{2}\right]$

$θn\theta_n$ 是神经网络的参数（权重系数），上述方法被作者称为深度BSDE方法。

例子

带违约风险的Black-Scholes方程

金融衍生品交易的一个关键问题就是确定适当而公平的价格。Black和Scholes说明了衍生品的价格 $u$ 满足一个抛物型的PDE，现在被称为Black-Scholes方程。Black-Scholes模型可以加以考虑实际市场中的几个重要因素，包括可违约证券，借贷利率高于贷款利率，交易成本，模型参数的不确定性等等。这几个效应中的每一个都会在定价模型中产生非线性的贡献。特别是，信贷危机和持续的欧洲主权债务危机凸显了原始的Black-Scholes模型中忽略的最基本风险，即违约风险。

理想情况下，定价模型应考虑金融衍生品所依赖的整个底层证券，从而产生了了高维非线性偏微分方程。然而，由于维数的诅咒，现有的定价算法通常无法解决这些问题。为了证明深度BSDE方法的有效性，作者研究了一个具有违约风险的递归估值模型的特例。

我们考虑尚未发生违约的情况下，100个标的资产的欧洲债权的价格是公平的。当违约的情况发生时，索赔的持有人仅收到当前值的分数倍，设为 $δ∈[0,1)\delta \in [0,1)$ 。对于强度 $Q$ (当前值的递减函数)，可能风险由泊松过程的第一个跳跃时间建模得到，换言之，当声明的值低时，风险变得更加可能。然后，可以通过以下方式对价值变化过程进行建模：
$f(t,x,u(t,x),σT(t,x)∇u(t,x))=−(1−δ)Q(u(t,x))u(t,x)−Ru(t,x)f\left(t, x, u(t, x), \sigma^{\mathrm{T}}(t, x) \nabla u(t, x)\right)=-(1-\delta) Q(u(t, x)) u(t, x)-R u(t, x)$

这里的R是无风险资产的利率。我们假设标的资产价格做几何布朗运动，并选择强度函数 $Q$ 为当前值的分段线性函数，有三个区域（ $vh<vl,γh>γlv^h<v^l,\gamma ^h > \gamma ^l$ ），函数形式如下：
$Q(y)=1(−∞,vh)(y)γh+1[vl,∞)(y)γl+1[vh,vl)(y)[(γh−γl)(vh−vI)(y−vh)+γh]Q(y)=1_{\left(-\infty, v^{h}\right)}(y) \gamma^{h}+1_{\left[v^{l}, \infty\right)}(y) \gamma^{l}+1_{\left[v^{h}, v^{l}\right)}(y)\left[\frac{\left(\gamma^{h}-\gamma^{l}\right)}{\left(v^{h}-v^{I}\right)}\left(y-v^{h}\right)+\gamma^{h}\right]$

那么，在 $[0,T]xR100[0,T]\text{x}R^{100}$ 上的非线性Black-Scholes方程就变为了：
$∂u∂t(t,x)+ρx⋅∇u(t,x)+σˉ22∑i=1d∣xi∣2∂2u∂xi2(t,x)−(1−δ)Q(u(t,x))u(t,x)−Ru(t,x)=0\frac{\partial u}{\partial t}(t, x)+\rho x \cdot \nabla u(t, x)+\frac{\bar{\sigma}^{2}}{2} \sum_{i=1}^{d}\left|x_{i}\right|^{2} \frac{\partial^{2} u}{\partial x_{i}^{2}}(t, x)-(1-\delta) Q(u(t, x)) u(t, x)-R u(t, x)=0$

设置如下参数， $1,\delta = 2/3,R=0.02,\bar \mu = 0.02,\bar \sigma = 0.2,v^h = 50,v^l = 70,\gamma ^h =0.2,\gamma ^l =0.02$ ，终止条件 $g(x)=min{x1,x2,...,x100}g(x)=\text{min}\{x_1,x_2,...,x_{100}\}$ 。

下图展示了 $θu0\theta_{u_0}$ （ $u(t=0,x=\{100,...100\})$ 的一个逼近）的平均值和标准差，其最终的相对误差为0.46%。

在这里插入图片描述

方程（11）并没有一个绝对的精确解，它有一个通过多层皮卡方法得到一个近似解为 $u(t=0,x=(100,...100))≈57.300u(t=0,x=(100,...100))\approx 57.300$ 。作为比较，如果我们不考虑违约风险，我们能得到 $u~(t=0,x=(100,...,100))≈60.781\tilde u(t=0,x=(100,...,100))\approx 60.781$ 。

在这种情况下，模型就变成线性的而直接可以用蒙特卡洛方法求解了。但是，如上所述，如果忽略违约风险，那么定价上就会存在相当大的错误。深度BSDE方法允许我们将违约风险纳入定价模型。这反过来又可以使得评估相关各方社会风险较低的金融衍生品成为可能。

HJB方程

我们考虑一个100维的经典的线性二次高斯控制问题，如下：

$[12]dX_t = 2\sqrt{\lambda}m_tdt + \lambda{2}dW_t \space \space \space \space [12]$

这里的 $t∈[0,T],X0=xt\in [0,T],X_0 = x$ ，并且成本函数可以写为：
$J({mt}0≤t≤T)=E[∫0T∥mt∥2dt+g(XT)]J\left(\left\{m_{t}\right\}_{0 \leq t \leq T}\right)=\mathbb{E}\left[\int_{0}^{T}\left\|m_{t}\right\|^{2} d t+g\left(X_{T}\right)\right]$

这里的 ${X_t\}$ 是状态过程， ${m_t\}$ 是控制过程， $λ\lambda$ 表示控制“长度”的一个正的常数。 ${W_t\}$ 是一个标准的布朗运动。我们的目标是通过控制过程最小化损耗函数。

这个问题的HJB方程可以写为：
$∂u∂t(t,x)+Δu(t,x)−λ∥∇u(t,x)∥2=0\frac{\partial u}{\partial t}(t, x)+\Delta u(t, x)-\lambda\|\nabla u(t, x)\|^{2}=0$

方程（13）的解在 $t = 0$ 的值，表示状态从 $x$ 开始时的最佳消耗。应用伊藤公式，我们能证明带有末端条件 $u (T, x) = g (x)$ 的方程（13）的精确解有一个显示的形式如下：
$x)=-\frac{1}{\lambda} \ln \left(\mathbb{E}\left[\exp \left(-\lambda g\left(x+\sqrt{2} W_{T-t}\right)\right)\right]\right)$

由此可以验证所提算法的精确性。

求解（13）方程的一些特殊情况，作图如下：

在这里插入图片描述

Allen-cahn 方程

考虑一个特殊的Allen-Cahn方程如下：
$∂u∂t(t,x)=Δu(t,x)+u(t,x)−[u(t,x)]3\frac{\partial u}{\partial t}(t, x)=\Delta u(t, x)+u(t, x)-[u(t, x)]^{3}$

初值条件为 $u(0,x)=g(x),g(x)=1/(2+0.4∣∣x∣∣2),x∈R100u(0,x)=g(x),g(x)=1/(2+0.4||x||^2),x\in \mathbb{R}^{100}$ ，通过对时间变量做变换 $t→T−t(T>0)t\rightarrow T-t(T>0)$ ，我们能把（15）变成（1）的形式，以便于深度BSDE方法能够使用。

数值实验结果如下图：

在这里插入图片描述

结论

本文提出的算法在许多不同的领域开辟了新的可能性。经济学领域，可以考虑多个不同的“交互代理”；金融领域，可以同时考虑所有参与的因素，而不只是临时假设变量之间的关系；运筹学领域，可以直接处理数百数千个实体参与的案例，而不需要进行临时的近似。

补充材料

BSDE方法构建

在很多文献当中，人们研究了这个（非线性）抛物型方程和倒向随机微分方程的一个关系。

考虑如下BSDE：
${Xt=ξ+∫0tμ(s,Xs)ds+∫0tσ(s,Xs)dWsYt=g(XT)+∫tTf(s,Xs,Ys,Zs)ds−∫tT(Zs)TdWs\left\{\begin{array}{l} X_{t}=\xi+\int_{0}^{t} \mu\left(s, X_{s}\right) d s+\int_{0}^{t} \sigma\left(s, X_{s}\right) d W_{s} \\ Y_{t}=g\left(X_{T}\right)+\int_{t}^{T} f\left(s, X_{s}, Y_{s}, Z_{s}\right) d s-\int_{t}^{T}\left(Z_{s}\right)^{\mathrm{T}} d W_{s} \end{array}\right.$

在某种意义下，如下的方程是almost surely成立的：
$Zl−σT(t,Xt)∇u(t,Xt),Y_{i}-u\left(t, X_{i}\right) \quad \text { and } \quad Z_{l}-\sigma^{\mathrm{T}}\left(t, X_{t}\right) \nabla u\left(t, X_{t}\right),$

一般来说，我们可以通过求解（16）、（17）来求解（1），更进一步，我们可以将（18）加入到（1）当中去，重写方程为（3）。

然后我们在时间上，对方程进行离散化，并用神经网络来近似空间梯度和最后的位置函数。

神经网络架构

下图基本上表示了深度BSDE方法的一个神经网络架构：

在这里插入图片描述

这里当前层的 $u$ 值与前一层的 $u$ 值、 $∇u\nabla u$ 值，以及 $X$ 值有关，当前层的 $X$ 值和前一层的 $X$ 值有关。初始层的 $X=ξX=\xi$ ，是我们关注的点。初始层的 $u(t_0,X_{t_0})$ 以及 $∇u(t0,Xt0)\nabla u(t_0,X_{t_0})$ 是未知的，我们作为神经网络的参数来处理，故而当前层的 $u$ 值和 $X$ 值可以从前一层得到，那么当前层只有 $∇u\nabla u$ 值没有，怎么办呢？我们从当前层的 $X$ 搭一个神经网络过来，来表示。每一时间离散步都有一个子神经网络，那么所有的子神经网络就可以堆叠成一个完整的神经网络。这里的 $W$ 是标准的布朗运动。最后Loss Function的构建，可以使用给定的末值条件。

我们所想要的其实是 $u(t_0,X_{t_0})$ 值，当训练好神经网络之后，这个值其实作为参数就已经固定下来了，我们把这个“参数”取出来，就是我们想要的。

实现

作者对文中提到的例子都做了一个实现，每个子网络采用全连接，由一个输入层（d维），一个输出层（d维）和2个隐藏层（d+10维）组成。使用整流函数RELU作为激活函数。在每次线性变换之后和激活之前，都是用了批量归一化（BN，正则化）的技术。该技术通过允许更大的步长和更容易的参数初始化来加速训练。所有参数的初始化都通过正的常数或者均匀分布进行初始化，无需任何预训练。

使用TensorFlow的Adam优化器来优化参数。Adam是随机梯度下降的一个变体，它是基于低阶矩的自适应估计。SGD算法选择批量大小为64，并且选择不同的随机种子5次独立运行，来近似 $L^1$ 误差的平均值和标准差。

隐藏层数目的影响

深度BSDE方法的准确性自然也取决于子网近似中隐藏层的数量。为了测试这种影响，作者也做了一些数值实验，设置不同的隐藏层数来求解反应扩散型PDE。

实验结果表明，当隐藏层数量增加的时候，精度是提高的。

在这里插入图片描述

我的参考文献

用深度学习求解PDE是我的导师给我选定的课题，但是现在我已经不做这个方向了，亦欢迎讨论。能查到的文献并不多，随机微分方程数值解因为之前接触过，相关文献就不列出。参考的主要是鄂维南老师的一篇文章，也就是下面的第一篇文章。

Han J, Jentzen A, Weinan E. Solving high-dimensional partial differential equations using deep learning[J]. Proceedings of the National Academy of Sciences, 2018, 115(34): 8505-8510.

Sirignano J, Spiliopoulos K. DGM: A deep learning algorithm for solving partial differential equations[J]. Journal of Computational Physics, 2018, 375: 1339-1364.

Raissi M, Perdikaris P, Karniadakis G E. Machine learning of linear differential equations using Gaussian processes[J]. Journal of Computational Physics, 2017, 348: 683-693.