Deep splitting method for parabolic PDEs

TL;DR

This paper presents a novel deep learning-based operator splitting method for solving high-dimensional nonlinear parabolic PDEs efficiently, demonstrated on physics, finance, and control problems.

Contribution

It introduces a new numerical approach combining operator splitting with deep learning to handle extremely high-dimensional PDEs effectively.

Findings

Successfully solves PDEs in up to 10,000 dimensions

Achieves very good accuracy with short computation times

Applicable to physics, stochastic control, and finance problems

Abstract

In this paper we introduce a numerical method for nonlinear parabolic PDEs that combines operator splitting with deep learning. It divides the PDE approximation problem into a sequence of separate learning problems. Since the computational graph for each of the subproblems is comparatively small, the approach can handle extremely high-dimensional PDEs. We test the method on different examples from physics, stochastic control and mathematical finance. In all cases, it yields very good results in up to 10,000 dimensions with short run times.

Tables

Table 1. Table 1: Deep splitting approximations of the solution of the HJB equation ( 3.3 ) for different values of d 𝑑 d , T 𝑇 T and N 𝑁 N .

$d$	$T$	$N$	Expectation	Std. dev.	Ref. value	rel. $L^1$-error	Std. dev. rel. error	avg. runtime
10	$1/3$	8	1.56645	0.00246699	1.56006	0.00410	0.00158134	18.0s
10	$2/3$	16	1.86402	0.00338646	1.85150	0.00677	0.00182904	37.9s
10	1	24	2.07017	0.00634850	2.04629	0.01167	0.00310245	58.2s
50	$1/3$	8	2.39214	0.00151918	2.38654	0.00234	0.00063656	18.0s
50	$2/3$	16	2.84607	0.00140300	2.83647	0.00338	0.00049463	37.9s
50	1	24	3.15098	0.00275839	3.13788	0.00417	0.00087906	58.4s
100	$1/3$	8	2.85090	0.00071267	2.84696	0.00138	0.00025033	18.1s
100	$2/3$	16	3.39109	0.00093368	3.38450	0.00195	0.00027587	38.2s
100	1	24	3.75329	0.00136920	3.74471	0.00229	0.00036564	58.3s
200	$1/3$	8	3.39423	0.00051028	3.39129	0.00087	0.00015047	18.1s
200	$2/3$	16	4.03680	0.00088215	4.03217	0.00115	0.00021878	38.0s
200	1	24	4.46734	0.00079688	4.46172	0.00126	0.00017860	58.2s
300	$1/3$	8	3.75741	0.00063334	3.75530	0.00056	0.00016865	18.3s
300	$2/3$	16	4.46859	0.00049953	4.46514	0.00077	0.00011187	38.5s
300	1	24	4.94586	0.00087736	4.94105	0.00097	0.00017756	58.8s
500	$1/3$	8	4.27079	0.00051256	4.26900	0.00042	0.00012007	18.0s
500	$2/3$	16	5.07900	0.00034792	5.07618	0.00056	0.00006854	38.0s
500	1	24	5.62126	0.00045092	5.61735	0.00070	0.00008027	57.7s
1000	$1/3$	8	5.07989	0.00022764	5.07876	0.00022	0.00004482	20.6s
1000	$2/3$	16	6.04130	0.00030680	6.03933	0.00033	0.00005080	43.7s
1000	1	24	6.68594	0.00040334	6.68335	0.00039	0.00006035	66.5s
5,000	$1/3$	8	7.59772	0.00024745	7.59733	0.00005	0.00003257	120.4s
5,000	$2/3$	16	9.03721	0.00027322	9.03466	0.00028	0.00003024	256.7s
5,000	1	24	9.97266	0.00047098	9.99835	0.00257	0.00004711	393.9s
10,000	$1/3$	8	9.03574	0.00022994	9.03535	0.00004	0.00002545	519.8s
10,000	$2/3$	16	10.74521	0.00026228	10.74478	0.00004	0.00002157	1105.6s
10,000	1	24	11.87860	0.00022705	11.89099	0.00104	0.00001909	1687.7s

Table 2. Table 2: Deep splitting approximations of the solution of the nonlinear Black–Scholes equations ( 3.5 ) for T = 1 / 3 𝑇 1 3 T=1/3 , N = 96 𝑁 96 N=96 and different d 𝑑 d .

$d$	Expectation	Stdev	Ref. value	rel. $L^1$-error	Stdev rel. error	avg. runtime
10	40.6553107	0.1000347132	40.7611353	0.0029624273	0.0019393471	858.3s
50	37.421057	0.0339765334	37.5217732	0.0026842068	0.0009055151	975.4s
100	36.3498646	0.027989905	36.4084035	0.0016078403	0.000768776	1481.5s
200	35.374638	0.035236816	35.4127342	0.0012857962	0.0006625744	951.2s
300	34.8476466	0.0225350305	34.8747946	0.0008818254	0.0004762554	953.3s
500	34.2206181	0.0081072294	34.2357988	0.0004701552	0.0001701012	956.0s
1000	33.4058827	0.0050161752	33.4358163	0.0008952555	0.000150024	1039.6s
5,000	31.7511529	0.0048508218	31.7906594	0.0012427078	0.0001525864	7229.7s
10,000	31.1215014	0.0031131196	31.1569116	0.0011365119	0.00009991746	23,593.2s

Table 3. Table 3: Deep splitting approximations of solutions of Allen–Cahn-type equations of the form ( 3.8 ) for T = 0.3 𝑇 0.3 T=0.3 , N = 10 𝑁 10 N=10 and different d 𝑑 d .

$d$	Expectation	Std. dev.	Ref. value	rel. $L^1$-error	Std. dev. rel. error	avg. runtime
10	0.89327	0.00299962	0.89060	0.00364	0.00258004	22.7s
50	1.01855	0.00073173	1.01830	0.00063	0.00036976	22.2s
100	1.04348	0.00029431	1.04510	0.00156	0.00028161	22.3s
200	1.06119	0.00018821	1.06220	0.00096	0.00017719	22.4s
300	1.06961	0.00017250	1.07217	0.00239	0.00016089	22.6s
500	1.07847	0.00013055	1.08124	0.00256	0.00012074	23.1s
1000	1.08842	0.00005689	1.09100	0.00236	0.00005215	25.9s
5,000	1.10522	0.00005201	1.10691	0.00153	0.00004699	134.6s
10,000	1.11071	0.00004502	1.11402	0.00296	0.00004041	473.6s

Table 4. Table 4: Deep splitting approximations of solutions of the semilinear heat equation ( 3.10 ) for T = 0.3 𝑇 0.3 T=0.3 , N = 20 𝑁 20 N=20 and different d 𝑑 d .

$d$	Expectation	Std. dev.	Ref. value	rel. $L^1$-error	Std. dev. rel. error	avg. runtime
10	0.47138	0.00035606	0.47006	0.00282	0.00075749	46.7s
50	0.34584	0.00018791	0.34425	0.00462	0.00054586	46.7s
100	0.31783	0.00008298	0.31674	0.00343	0.00026198	47.4s
200	0.30210	0.00002238	0.30091	0.00394	0.00007436	48.1s
300	0.29654	0.00001499	0.29534	0.00406	0.00005075	48.3s
500	0.29200	0.00000611	0.29095	0.00361	0.00002099	48.5s
1000	0.28852	0.00000267	0.28753	0.00344	0.00000930	54.1s
5,000	0.28569	0.00000042	0.28469	0.00352	0.00000148	286.3s
10,000	0.28533	0.00000048	0.28433	0.00353	0.00000170	1013.0s

Table 5. Table 5: Deep splitting approximations of solutions of the semilinear heat equation ( 3.10 ) for d = 100 𝑑 100 d=100 , T = 0.3 𝑇 0.3 T=0.3 and different N 𝑁 N .

$N$	Expectation	Std. dev.	Ref. value	rel. $L^1$-error	Std. dev. rel. error	avg. runtime
1	0.33820	0.00005021	0.31674	0.06777	0.00015852	0.5s
2	0.32836	0.00006189	0.31674	0.03669	0.00019539	2.2s
4	0.32259	0.00003834	0.31674	0.01848	0.00012104	5.7s
8	0.31969	0.00006929	0.31674	0.00931	0.00021878	12.8s
16	0.31810	0.00006968	0.31674	0.00430	0.00021998	27.1s
20	0.31783	0.00010565	0.31674	0.00345	0.00033354	34.3s
32	0.31739	0.00006564	0.31674	0.00206	0.00020723	55.6s

Table 6. Table 6: Deep splitting approximations of the solution of the sine-Gordon-type equation ( 3.12 ) for T = 0.3 𝑇 0.3 T=0.3 , N = 20 𝑁 20 N=20 and different d 𝑑 d .

$d$	Expectation	Std. dev.	Ref. value	rel. $L^1$-error	Std. dev. rel. error	avg. runtime
10	0.3218822	0.00069331	0.3229470	0.0032972	0.0021468	55.6s
50	0.0990598	0.00013433	0.0993633	0.0030544	0.0013519	55.4s
100	0.0526955	0.00005390	0.0528368	0.0026741	0.0010202	55.6s
200	0.0271860	0.00001617	0.0272410	0.0020176	0.0005936	56.2s
300	0.0183162	0.00001010	0.0183617	0.0024765	0.0005502	55.8s
500	0.0110819	0.00000428	0.0111071	0.0022647	0.0003851	54.7s
1000	0.0055775	0.00000090	0.0055896	0.0021697	0.0001613	58.1s
5,000	0.0011209	0.00000012	0.0011231	0.0019464	0.0001070	332.5s
10,000	0.0005608	0.00000004	0.0005621	0.0022495	0.0000720	1083.4s

Equations

$$\frac{\partial }{\partial t}u(t,x)$$

$$F\big{(}x,u(t,x),\nabla_{x}u(t,x)\big{)}=\langle\mu(x),\nabla_{x}u(t,x)\rangle_{{\mathbb{R}}^{d}}+f\big{(}x,u(t,x),\nabla_{x}u(t,x)\big{)}$$

$$\begin{split}\tfrac{\partial}{\partial t}u(t,x)&=f\big{(}x,u(t,x),\nabla_{x}u(t,x)\big{)}+\big{\langle}\mu(x),\nabla_{x}u(t,x)\big{\rangle}_{{\mathbb{R}}^{d}}\\ &\quad+\tfrac{1}{2}\operatorname{Trace}\!\big{(}\sigma(x)\sigma^{*}(x)\operatorname{Hess}_{x}u(t,x)\big{)},\end{split}$$

$$0=t_0<t_1<\dots <t_N=T.$$

$$\begin{split}u(t,x)&=\varphi(x)+\int_{0}^{t}f\big{(}x,u(s,x),\nabla_{x}u(s,x)\big{)}\,ds\\ &\quad+\int_{0}^{t}\Big{[}\tfrac{1}{2}\operatorname{Trace}\!\big{(}\sigma(x)\sigma^{*}(x)\operatorname{Hess}_{x}u(s,x)\big{)}+\big{\langle}\mu(x),\nabla_{x}u(s,x)\big{\rangle}_{{\mathbb{R}}^{d}}\Big{]}\,ds.\end{split}$$

$$\begin{split}u(t,x)&=u(t_{n},x)+\int_{t_{n}}^{t}f\big{(}x,u(s,x),(\nabla_{x}u)(s,x)\big{)}\,ds\\ &\quad+\int_{t_{n}}^{t}\Big{[}\tfrac{1}{2}\operatorname{Trace}\!\big{(}\sigma(x)\sigma^{*}(x)\operatorname{Hess}_{x}u(s,x)\big{)}+\big{\langle}\mu(x),\nabla_{x}u(s,x)\big{\rangle}_{{\mathbb{R}}^{d}}\Big{]}\,ds,\end{split}$$

$$\begin{split}u(t,x)&\approx u(t_{n},x)+\int_{t_{n}}^{t_{n+1}}f\big{(}x,u(t_{n},x),\nabla_{x}u(t_{n},x)\big{)}\,ds\\ &\quad+\int_{t_{n}}^{t}\Big{[}\tfrac{1}{2}\operatorname{Trace}\!\big{(}\sigma(x)\sigma^{*}(x)\operatorname{Hess}_{x}u(s,x)\big{)}+\big{\langle}\mu(x),\nabla_{x}u(s,x)\big{\rangle}_{{\mathbb{R}}^{d}}\Big{]}\,ds,\end{split}$$

$$\begin{split}u(t,x)&\approx u(t_{n},x)+f\bigl{(}x,u(t_{n},x),\nabla_{x}u(t_{n},x)\bigr{)}\,(t_{n+1}-t_{n})\\[4.30554pt] &\quad+\int_{t_{n}}^{t}\Big{[}\tfrac{1}{2}\operatorname{Trace}\!\big{(}\sigma(x)\sigma^{*}(x)\operatorname{Hess}_{x}u(s,x)\big{)}+\big{\langle}\mu(x),\nabla_{x}u(s,x)\big{\rangle}_{{\mathbb{R}}^{d}}\Big{]}\,ds.\end{split}$$

$$\begin{split}v(t,x)&=v(t_{n},x)+f\big{(}x,v(t_{n},x),\nabla_{x}v(t_{n},x)\big{)}\,(t_{n+1}-t_{n})\\[4.30554pt] &\quad+\int_{t_{n}}^{t}\Big{[}\tfrac{1}{2}\operatorname{Trace}\!\big{(}\sigma(x)\sigma^{*}(x)\operatorname{Hess}_{x}v(s,x)\big{)}+\big{\langle}\mu(x),\nabla_{x}v(s,x)\big{\rangle}_{{\mathbb{R}}^{d}}\Big{]}\,ds;\end{split}$$

$$v(t_n,x)\approx u(t_n,x)\text{mbox}foralln\in \{1,\dots ,N\}.$$

$$Y_{t}=\xi+\int_{0}^{t}\mu(Y_{s})\,ds+\int_{0}^{t}\sigma(Y_{s})\,dB_{s}\quad\mbox{$\mathbbm{P}$-a.s.}$$

$$\sup_{t\in[0,T]}{\mathbbm{E}}\big{[}\|Y_{t}\|_{{\mathbb{R}}^{d}}^{p}\big{]}<\infty\quad\mbox{for all }p\in(0,\infty);$$

$$\tfrac{\partial}{\partial t}v(t,x)=\big{\langle}\mu(x),\nabla_{x}v(t,x)\big{\rangle}_{{\mathbb{R}}^{d}}+\tfrac{1}{2}\operatorname{Trace}\!\big{(}\sigma(x)\sigma^{*}(x)\operatorname{Hess}_{x}v(t,x)\big{)}$$

$$\tfrac{\partial}{\partial t}v(T-t,x)+\big{\langle}\mu(x),\nabla_{x}v(T-t,x)\big{\rangle}_{{\mathbb{R}}^{d}}+\tfrac{1}{2}\operatorname{Trace}\big{(}\sigma(x)\sigma^{*}(x)\operatorname{Hess}_{x}v(T-t,x)\big{)}=0.$$

$$\begin{split}&v(T-t,Y_{t})=v(T-t_{n+1},Y_{t_{n+1}})+\int_{T-t_{n+1}}^{t}\big{\langle}\nabla_{x}v(T-s,Y_{s}),\sigma(Y_{s})\,dB_{s}\big{\rangle}_{{\mathbb{R}}^{d}}\\ &+\int_{T-t_{n+1}}^{t}\tfrac{\partial}{\partial s}v(T-s,Y_{s})\,ds+\int_{T-t_{n+1}}^{t}\tfrac{1}{2}\operatorname{Trace}\big{(}\sigma(Y_{s})\sigma^{*}(Y_{s})\operatorname{Hess}_{x}v(T-s,Y_{s})\big{)}\,ds\\ &+\int_{T-t_{n+1}}^{t}\big{\langle}\mu(Y_{s}),\nabla_{x}v(T-s,Y_{s})\big{\rangle}_{{\mathbb{R}}^{d}}\,ds\quad\mbox{$\mathbbm{P}$-a.s.},\end{split}$$

$$v(T-t,Y_{t})=v(T-t_{n+1},Y_{t_{n+1}})+\int_{T-t_{n+1}}^{t}\big{\langle}\nabla_{x}v(T-s,Y_{s}),\sigma(Y_{s})\,dB_{s}\big{\rangle}_{{\mathbb{R}}^{d}}.$$

$$\int_{T-t_{n+1}}^{t}{\mathbbm{E}}\Big{[}\big{\|}\sigma^{*}(Y_{s})\nabla_{x}v(T-s,Y_{s})\big{\|}_{{\mathbb{R}}^{d}}^{2}\Big{]}\,ds<\infty,$$

$${\mathbbm{E}}\bigg{[}\int_{T-t_{n+1}}^{t}\big{\langle}\nabla_{x}v(T-s,Y_{s}),\sigma(Y_{s})\,dB_{s}\big{\rangle}_{{\mathbb{R}}^{d}}\!~{}\Big{|}\!~{}\mathcal{F}_{T-t_{n+1}}\bigg{]}=0\quad\mbox{$\mathbbm{P}$-a.s.}$$

$${\mathbbm{E}}\Big{[}v(T-t,Y_{t})\!~{}\big{|}\!~{}\mathcal{F}_{T-t_{n+1}}\Big{]}={\mathbbm{E}}\Big{[}v(t_{n+1},Y_{T-t_{n+1}})\!~{}\big{|}\!~{}\mathcal{F}_{T-t_{n+1}}\Big{]}\quad\mbox{$\mathbbm{P}$-a.s.}$$

$${\mathbbm{E}}\Big{[}v(T-t,Y_{t})\!~{}\big{|}\!~{}Y_{T-t_{n+1}}\Big{]}={\mathbbm{E}}\Big{[}v(t_{n+1},Y_{T-t_{n+1}})\!~{}\big{|}\!~{}Y_{T-t_{n+1}}\Big{]}=v(t_{n+1},Y_{T-t_{n+1}})\quad\mbox{$\mathbbm{P}$-a.s.}$$

$$\begin{split}&\lim_{t\nearrow T-t_{n}}v(T-t,Y_{t}(\omega))=v(t_{n},Y_{T-t_{n}}(\omega))\\ &\qquad+f\big{(}Y_{T-t_{n}}(\omega),v(t_{n},Y_{T-t_{n}}(\omega)),\nabla_{x}v(t_{n},Y_{T-t_{n}}(\omega))\big{)}\,(t_{n+1}-t_{n}).\end{split}$$

$$\sup_{t\in(t_{n},t_{n+1}]}{\mathbbm{E}}\big{[}|v(T-t,Y_{t})|^{p}\big{]}<\infty.$$

$$\begin{split}&v(t_{n+1},Y_{T-t_{n+1}})={\mathbbm{E}}\Big{[}\lim_{t\nearrow T-t_{n}}v(T-t,Y_{t})\!~{}\big{|}\!~{}Y_{T-t_{n+1}}\Big{]}\\ &={\mathbbm{E}}\Bigl{[}v(t_{n},Y_{T-t_{n}})+f\big{(}Y_{T-t_{n}},v(t_{n},Y_{T-t_{n}}),\nabla_{x}v(t_{n},Y_{T-t_{n}})\big{)}\,(t_{n+1}-t_{n})\!~{}\big{|}\!~{}Y_{T-t_{n+1}}\Bigr{]}\quad\mbox{$\mathbbm{P}$-a.s.,}\end{split}$$

$${\mathbbm{E}}\Big{[}\big{|}v(t_{n-1},Y_{T-t_{n-1}})+f\big{(}Y_{T-t_{n-1}},v(t_{n-1},Y_{T-t_{n-1}}),\nabla_{x}v(t_{n-1},Y_{T-t_{n-1}}\big{)}\,(t_{n}-t_{n-1})\big{|}^{2}\Big{]}<\infty.$$

$$(v(t_{n},x))_{x\in\operatorname{supp}(Y_{T-t_{n}}(\mathbbm{P}))}=\operatornamewithlimits{argmin}_{w\in C(\operatorname{supp}(Y_{T-t_{n}}(\mathbbm{P})),{\mathbb{R}})}{\mathbbm{E}}\Big{[}\big{|}w(Y_{T-t_{n}})-\big{[}v(t_{n-1},Y_{T-t_{n-1}})\\ +f\big{(}Y_{T-t_{n-1}},v(t_{n-1},Y_{T-t_{n-1}}),\nabla_{x}v(t_{n-1},Y_{T-t_{n-1}})\big{)}\,(t_{n}-t_{n-1})\big{]}\big{|}^{2}\Big{]}.$$

$$V_0(\theta ,x)=\phi (x)$$

$$V_n(\theta ,x)=A_{l,1}^{\theta ,(k-1)l(l+1)+l(d+1)}\circ {\mathcal{L}}_l\circ A_{l,l}^{\theta ,(k-2)l(l+1)+l(d+1)}\circ \dots \circ {\mathcal{L}}_l\circ A_{l,l}^{\theta ,l(d+1)}\circ {\mathcal{L}}_l\circ A_{d,l}^{\theta ,0}(x),$$

$$A^{\theta,r}_{i,j}(x)=\left(\begin{array}[]{cccc}\theta_{r+1}&\theta_{r+2}&\dots&\theta_{r+i}\\ \theta_{r+i+1}&\theta_{r+i+2}&\dots&\theta_{r+2i}\\ \theta_{r+2i+1}&\theta_{r+2i+2}&\dots&\theta_{r+3i}\\ \vdots&\vdots&\vdots&\vdots\\ \theta_{r+i(j-1)+1}&\theta_{r+i(j-1)+2}&\dots&\theta_{r+ij}\end{array}\right)\left(\begin{array}[]{c}x_{1}\\ x_{2}\\ x_{3}\\ \vdots\\ x_{i}\end{array}\right)+\left(\begin{array}[]{c}\theta_{r+ij+1}\\ \theta_{r+ij+2}\\ \theta_{r+ij+3}\\ \vdots\\ \theta_{r+ij+j}\end{array}\right),$$

$${\mathcal{L}}_l(x_1,\dots ,x_l)=(\rho (x_1),\dots ,\rho (x_l))$$

$${\mathbb{R}}^{\nu}\ni\theta\mapsto{\mathbbm{E}}\Big{[}\big{|}V_{n}(\theta,Y_{T-t_{n}})-\big{[}V_{n-1}(\theta^{n-1},Y_{T-t_{n-1}})\\ +f\big{(}Y_{T-t_{n-1}},V_{n-1}(\theta^{n-1},Y_{T-t_{n-1}}),\nabla_{x}V_{n-1}(\theta^{n-1},Y_{T-t_{n-1}})\big{)}\,(t_{n}-t_{n-1})\big{]}\big{|}^{2}\Big{]}\in{\mathbb{R}}.$$

Full text

Key words. nonlinear partial differential equations, splitting-up method, neural networks, deep learning

AMS subject classifications. 35K15, 65C05, 65M22, 65M75, 91G20, 93E20

1 Introduction

In this paper we derive a numerical scheme for parabolic partial differential equations (PDEs) of the form

[TABLE]

$(t,x)\in(0,T]\times{\mathbb{R}}^{d}$, with initial condition $u(0,x)=\varphi(x)$, where $F\colon\mathbb{R}^{d}\times\mathbb{R}\times\mathbb{R}^{d}\to\mathbb{R}$, $\sigma\colon\mathbb{R}^{d}\to\mathbb{R}^{d\times d}$, and $\varphi\colon\mathbb{R}^{d}\to\mathbb{R}$ are appropriate continuous functions. Such PDEs describe various phenomena in nature, engineering, economics, and finance. They typically do not admit closed form solutions and, therefore, have to be solved numerically. In some applications, the dimension $d$ can be high. For instance, in physics and engineering applications, $x\in\mathbb{R}^{d}$ typically models the coordinates of all components of a given system, whereas in derivative pricing and optimal investment problems, $d$ usually corresponds to the number of underlying assets. Many classical PDEs, such as the standard heat and Black–Scholes equations are linear. Using the Feynman–Kac representation, their solutions can efficiently be approximated in high dimensions with simple Monte Carlo averages. But if constraints or frictions are taken into account, or the PDE describes a control problem, the function $F$ is no longer linear and equation (1.1) becomes much more challenging to solve for large $d$.

Numerical methods for PDEs have a long history. Classical approaches like finite differences and finite elements (see, e.g., [13, 71, 95]) are deterministic. In their standard form, they work well for $d=1,2$ and $3$, but their complexity grows exponentially in $d$. To tackle higher dimensional problems, different simulation-based approaches have been developed that exploit a stochastic representation of the solution of the PDE. For instance, [1, 7, 8, 12, 15, 16, 17, 18, 19, 20, 21, 22, 24, 26, 36, 37, 38, 39, 40, 41, 54, 70, 73, 74, 77, 78, 79, 84, 85, 87, 89, 90, 91, 96, 98] use BSDE representations of PDEs and study approximation methods based on recursive polynomial regressions, [50, 52, 53, 82, 93, 97] investigate methods based on branching diffusion processes, and  [29, 30, 57, 59, 60] analyze full-history recursive multilevel Picard methods. Recently, numerical methods for high-dimensional PDEs based on the idea to reformulate the PDE as a stochastic learning problem have been proposed in [28, 48]. This opens the door to the application of deep learning; see, e.g., [2, 3, 6, 10, 14, 31, 33, 34, 42, 49, 51, 55, 62, 75, 76, 80, 88, 92] for modifications and extensions. There are also already a few papers studying the convergence of deep learning based approximation methods for PDEs. For instance, [49, 92] derive convergence results without information on the convergence speed, whereas [11, 32, 44, 56, 63, 69] provide convergence and tractability results with dimension-independent convergence rates and error constants depending polynomially on the dimension.

In this paper we develop a new deep learning method for parabolic PDEs that splits the differential operator into a linear and a nonlinear part. More precisely, we write

[TABLE]

for suitable continuous functions $\mu\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$ and $f\colon\mathbb{R}^{d}\times\mathbb{R}\times\mathbb{R}^{d}\to\mathbb{R}$. This decomposition is not unique. But the idea is that $\mu$ is chosen such that the nonlinearity f\big{(}x,u(t,x),\nabla_{x}u(t,x)\big{)} becomes small. Then we solve the PDE iteratively over small time intervals by approximating f\big{(}x,u(t,x),\nabla_{x}u(t,x)\big{)} and using the Feynman–Kac representation locally. This requires a recursive computation of conditional expectations. We approximate them by formulating them as minimization problems that can be approached with deep learning. This decomposes the PDE approximation problem into a sequence of separate learning problems. Since the computational graph for each of the subproblems is comparatively small, the method works for very high-dimensional problems.

The rest of the paper is organized as follows. In Section 2 we introduce the framework and derive the deep splitting method. In Section 3 we test the approach on five different classes of high-dimensional nonlinear PDEs: Hamilton–Jacobi–Bellman (HJB) equations, nonlinear Black–Scholes equations, Allen–Cahn-type equations, nonlinear heat equations, and sine-Gordon-type equations.

2 Derivation of the proposed approximation algorithm

Fix $T\in(0,\infty)$ and $d\in{\mathbb{N}}$. Consider two at most polynomially growing continuous functions $\varphi\colon{\mathbb{R}}^{d}\to{\mathbb{R}}$ and $f\colon{\mathbb{R}}^{d}\times{\mathbb{R}}\times{\mathbb{R}}^{d}\to{\mathbb{R}}$ together with two Lipschitz continuous functions $\mu\colon{\mathbb{R}}^{d}\to{\mathbb{R}}^{d}$ and $\sigma\colon{\mathbb{R}}^{d}\to{\mathbb{R}}^{d\times d}$. Assume $u\colon[0,T]\times{\mathbb{R}}^{d}\to{\mathbb{R}}$ is an at most polynomially growing continuous function that is $C^{1,2}$ on $(0,T]\times{\mathbb{R}}^{d}$ and satisfies the PDE

[TABLE]

$(t,x)\in(0,T]\times{\mathbb{R}}^{d}$, with initial condition $u(0,x)=\varphi(x)$, $x\in{\mathbb{R}}^{d}$.

2.1 Temporal discretization

To approximate the solution $u$ of the PDE (2.1), we discretize the equation in time and use a splitting-up method (see, e.g., [43, 45, 46]) to obtain a semi-discrete approximation problem. To do this, we choose $N\in{\mathbb{N}}$ and let $t_{0},t_{1},\ldots,t_{N}\in[0,T]$ be real numbers such that

[TABLE]

Under appropriate integrability assumptions, it follows from (2.1) that for every $t\in[0,T]$, $x\in{\mathbb{R}}^{d}$ we have

[TABLE]

In particular, for all $n\in\{0,1,\ldots,N-1\}$, $t\in[t_{n},t_{n+1}]$ and $x\in{\mathbb{R}}^{d}$,

[TABLE]

and therefore,

[TABLE]

which can be written as

[TABLE]

To derive the splitting-up approximation, we make a few simplifying assumptions, not all of which are needed to implement the resulting algorithm. First, we suppose that $\varphi$ has an at most polynomially growing gradient $\nabla\varphi\colon{\mathbb{R}}^{d}\to{\mathbb{R}}^{d}$. Moreover, we assume that there exists a function $v\colon[0,T]\times{\mathbb{R}}^{d}\to{\mathbb{R}}$ satisfying $v(0,x)=\varphi(x)$, $x\in{\mathbb{R}}^{d}$, such that for every $n\in\{0,1,\ldots,N-1\}$, $v|_{(t_{n},t_{n+1}]\times{\mathbb{R}}^{d}}$ belongs to $C^{1,2}((t_{n},t_{n+1}]\times{\mathbb{R}}^{d},{\mathbb{R}})$ with at most polynomially growing partial derivatives and for all $n\in\{0,1,\ldots,N-1\}$, $t\in(t_{n},t_{n+1}]$ and $x\in{\mathbb{R}}^{d}$,

[TABLE]

see, e.g., Hairer et al. [47, Section 4.4], Deck & Kruse [23], Krylov [67, Chapter 8], and Krylov [68, Theorem 4.32] for existence, uniqueness, and regularity results for equations of the form (2.7). Comparing (2.7) to (2.6) suggests that

[TABLE]

So $v$ is a specific splitting-up type approximation of the function $u$; see, e.g., [25, 43, 45, 46].

2.2 An approximate Feynman–Kac representation

In the next step we derive a Feynman–Kac representation of $v$; see, e.g., Milstein & Tretyakov [86, Section 2]. Let $B\colon[0,T]\times\Omega\to{\mathbb{R}}^{d}$ be a standard $(\mathcal{F}_{t})_{t\in[0,T]}$-Brownian motion on a filtered probability space $(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t\in[0,T]},\mathbbm{P})$ satisfying the usual conditions. Consider an $\mathcal{F}_{0}$/$\mathcal{B}({\mathbb{R}}^{d})$-measurable random variable $\xi\colon\Omega\to{\mathbb{R}}^{d}$ satisfying ${\mathbbm{E}}[\|\xi\|_{{\mathbb{R}}^{d}}^{p}]<\infty$ for every $p\in(0,\infty)$, and let $Y\colon[0,T]\times\Omega\to{\mathbb{R}}^{d}$ be an $(\mathcal{F}_{t})_{t\in[0,T]}$-adapted process with continuous sample paths satisfying for every $t\in[0,T]$,

[TABLE]

The assumption that ${\mathbbm{E}}[\|\xi\|_{{\mathbb{R}}^{d}}^{p}]<\infty$ for all $p\in(0,\infty)$ ensures that

[TABLE]

see, e.g., Stroock [94, Section 1.2]. Moreover, (2.7) implies that for all $n\in\{0,1,\ldots,N-1\}$, $t\in(t_{n},t_{n+1})$ and $x\in{\mathbb{R}}^{d}$, one has

[TABLE]

from which it follows that for all $n\in\{0,1,\ldots,N-1\}$, $t\in(T-t_{n+1},T-t_{n})$ and $x\in{\mathbb{R}}^{d}$,

[TABLE]

Since for every $n\in\{0,1,\ldots,N-1\}$, $v|_{(t_{n},t_{n+1}]\times{\mathbb{R}}^{d}}$ is in $C^{1,2}((t_{n},t_{n+1}]\times{\mathbb{R}}^{d},{\mathbb{R}})$, we obtain from Itô’s formula that for all $n\in\{0,1,\ldots,N-1\}$ and $t\in[T-t_{n+1},T-t_{n})$ we have

[TABLE]

which by (2.12), gives that $\mathbbm{P}$-a.s.,

[TABLE]

Moreover, since, by assumption, $\sigma\colon{\mathbb{R}}^{d}\to{\mathbb{R}}^{d\times d}$ is Lipschitz continuous and $\nabla_{x}v(t,x)$ at most polynomially growing in $(t,x)\in(t_{n},t_{n+1}]\times{\mathbb{R}}^{d}$, one obtains from (2.10) that for all $n\in\{0,1,\dots,N-1\}$ and $t\in[T-t_{n+1},T-t_{n})$,

[TABLE]

from which it follows that

[TABLE]

Together with (2.14), this shows that for all $n\in\{0,1,\ldots,N-1\}$ and $t\in[T-t_{n+1},T-t_{n})$ we have

[TABLE]

Since $Y_{T-t_{n+1}}$ is $\mathcal{F}_{T-t_{n+1}}$/$\mathcal{B}({\mathbb{R}})$-measurable, one obtains from the tower property of conditional expectations that for all $n\in\{0,1,\ldots,N-1\}$ and $t\in[T-t_{n+1},T-t_{n})$,

[TABLE]

Since $(Y_{t})_{t\in[0,T]}$ has continuous sample paths and $(v(t,x))_{(t,x)\in(t_{n},t_{n+1}]\times{\mathbb{R}}^{d}}$ at most polynomially growing first order partial derivatives, it follows from (2.7) and (2.11) that we have for all $n\in\{0,1,\dots,N-1\}$ and $\omega\in\Omega$,

[TABLE]

In addition, since $\sup_{t\in[0,T]}{\mathbbm{E}}[\|Y_{t}\|_{{\mathbb{R}}^{d}}^{p}]$ $<\infty$ for every $p\in(0,\infty)$, we obtain

[TABLE]

So it follows from (2.18) and (2.19) that for all $n\in\{0,1,\ldots,N-1\}$,

[TABLE]

which is the Feynman–Kac type representation we were aiming for. Note that the nonlinearity $f$ as well as the coefficient functions $\mu$ and $\sigma$ in (2.1) do not depend on $t$. But the above derivation could be extended to time-dependent $f$, $\mu$, and $\sigma$ since the Feynman–Kac formula still holds in this case.

2.3 Formulation as recursive minimization problems

We now reformulate (2.21) as recursive minimization problems. It follows from our assumptions that for every $n\in\{1,2,\ldots,N\}$, $v(t_{n-1},x)+f(x,v(t_{n-1},x),\nabla_{x}v(t_{n-1},x))(t_{n}-t_{n-1})$ is at most polynomially growing in $x\in{\mathbb{R}}^{d}$. Therefore, we obtain from (2.10) that

[TABLE]

Since $v(t_{n},x)$ is continuous in $x\in{\mathbb{R}}^{d}$, it follows from the factorization lemma and the $L^{2}$-minimality property of conditional expectations (see, e.g., Klenke [65, Corollary 8.17]) that for every $n\in\{1,2,\ldots,N\}$ we have

[TABLE]

2.4 Deep artificial neural network approximations

To tackle the minimization problems (2.23) numerically, we approximate the functions $v(t_{n},.)$, $n\in\{1,2,\dots,N\}$, with neural networks $V_{n}$. More precisely, we choose $k\in\{3,4,...\}$ and $l\in{\mathbb{N}}$. Then, we set $\nu=(1+kl-l)(l+1)+l(d+1)$ and consider functions $V_{n}\colon{\mathbb{R}}^{\nu}\times{\mathbb{R}}^{d}\to{\mathbb{R}}$, $n\in\{0,1,\dots,N\}$, such that for every $(\theta,x)\in{\mathbb{R}}^{\nu}\times{\mathbb{R}}^{d}$, $V_{0}$ is given by

[TABLE]

and for $n\in\{1,2,\dots,N\}$, $V_{n}$ is of the form

[TABLE]

where for $r\in{\mathbb{N}}_{0}=\{0\}\cup{\mathbb{N}}$ and $i,j\in{\mathbb{N}}$ with $r+j(i+1)\leq\nu,$ $A^{\theta,r}_{i,j}\colon{\mathbb{R}}^{i}\to{\mathbb{R}}^{j}$ is the affine function defined by

[TABLE]

and ${\cal L}_{l}\colon{\mathbb{R}}^{l}\to{\mathbb{R}}^{l}$ is a mapping of the form

[TABLE]

for a weakly differentiable continuous function $\rho\colon{\mathbb{R}}\to{\mathbb{R}}$.

(2.25) is a feedforward neural network with activation function $\rho$ and $k+1$ layers (an input layer with $d$ neurons, $k-1$ hidden layers with $l$ neurons each, and an output layer with one neuron); see, e.g., [9, 72]. Commonly used activation functions are e.g., the logistic function $x\mapsto e^{x}/(1+e^{x})$ or the ReLU function $x\mapsto\max\{0,x\}$. The logistic function is continuously differentiable. So the corresponding neural networks $V_{n}$ have well-defined $\theta$- and $x$-gradients $\nabla_{\theta}V_{n}$ and $\nabla_{x}V_{n}$. However, in the examples of Section 3 below we use the ReLU function. It is continuously differentiable on $(-\infty,0)\cup(0,\infty)$ and has a left-hand derivative at [math]. This yields weak $\theta$- and $x$-gradients for $V_{n}$, which we also denote by $\nabla_{\theta}V_{n}$ and $\nabla_{x}V_{n}$, respectively.

2.5 Stochastic gradient descent based minimization

Next, we recursively solve quadratic minimization problems to find parameter vectors $\theta^{1},\theta^{2},\ldots,$ $\theta^{N}\in{\mathbb{R}}^{\nu}$ such that $V_{n}(\theta^{n},x)\approx v(t_{n},x)$ for $n\in\{1,2,\dots,N\}$. More specifically, $\theta^{0}\in{\mathbb{R}}^{\nu}$ can be chosen arbitrarily; e.g., $\theta^{0}=(0,\dots,0)\in{\mathbb{R}}^{\nu}$. For given $n\in\{1,2,\ldots,N\}$ and $\theta^{0},\theta^{1},\dots,\theta^{n-1}\in{\mathbb{R}}^{\nu}$, one tries to find an approximate minimizer $\theta^{n}\in{\mathbb{R}}^{\nu}$ of the function

[TABLE]

To do this with a standard stochastic gradient descent, one can initialize $\vartheta^{n}_{0}$ randomly or, e.g., as $\vartheta^{n}_{0}=\theta^{n-1}$. Then one chooses a stepsize $\gamma\in(0,\infty)$ together with a number $M\in{\mathbb{N}}$ and iteratively updates for $m\in\{0,1,\dots,M-1\}$ according to

[TABLE]

where $Y^{m}\colon[0,T]\times\Omega\to{\mathbb{R}}^{d}$, $m\in\{0,1,\dots,M-1\}$, are $(\mathcal{F}_{t})_{t\in[0,T]}$-adapted stochastic process with continuous sample paths satisfying for every $t\in[0,T]$, the SDEs

[TABLE]

corresponding to i.i.d. standard $(\mathcal{F}_{t})_{t\in[0,T]}$-Brownian motions $B^{m}\colon[0,T]\times\Omega\to{\mathbb{R}}^{d}$ and i.i.d. $\mathcal{F}_{0}/\mathcal{B}({\mathbb{R}}^{d})$-measurable functions $\xi^{m}\colon\Omega\to{\mathbb{R}}^{d}$, $m\in\{0,1,\dots,M-1\}$. After $M$ gradient steps one sets $\theta^{n}=\vartheta^{n}_{M}$.

2.6 Discretization of the auxiliary stochastic process $Y$

Equation (2.29) provides an implementable numerical algorithm in the special case where the solutions $Y^{m}$ of the SDEs (2.30) can be simulated exactly. If this is not the case, one can use a numerical approximation method to approximatively simulate $Y^{m}$, $m\in\{0,1,\dots,M-1\}$. In the following we concentrate on the Euler–Maruyama scheme. But it is also possible to use a different approximation method.

Note that it follows from (2.30) that for all $m\in\{0,1,\dots,M-1\}$ and $n\in\{0,1,\ldots,N-1\}$ we have

[TABLE]

or equivalently,

[TABLE]

for $\tau_{n}=T-t_{N-n}$. This suggests that for all $m\in\{0,1,\dots,M-1\}$ and $n\in\{0,1,\dots,N-1\}$,

[TABLE]

Therefore, we introduce the Euler–Maruyama approximations $\mathcal{Y}^{m}\colon\{0,1,\dots,N\}\times\Omega\to{\mathbb{R}}^{d}$ given for every $n\in\{0,1,\dots,N-1\}$ by ${\mathcal{Y}}^{m}_{0}=\xi^{m}$ and

[TABLE]

It can be seen from (2.33) and (2.34) that for every $n\in\{0,1,\dots,N\}$, one has

[TABLE]

which can be used to derive approximations of $(\vartheta^{n}_{m})_{m=0}^{M}$, $n\in\{1,2,\dots,N\}$, from (2.29) that are also implementable if the processes $Y^{m}$ cannot be simulated exactly. More precisely, set $\Theta^{0}_{M}=(0,\dots,0)\in{\mathbb{R}}^{\nu}$. For $n\in\{1,2,\dots,N\}$, initialize $\Theta^{n}_{0}$ randomly or as $\Theta^{n}_{0}=\Theta^{n-1}_{M}$. Then set for every $m\in\{0,1,\dots,M-1\}$,

[TABLE]

Comparing (2.36) to (2.29) suggests that $\Theta^{n}_{m}\approx\vartheta^{n}_{m}$ for $m\in\{0,1,\dots,M\}$ and $n\in\{1,2,\dots,N\}$.

In the following two Subsections 2.7 and 2.8, we first formalize and then generalize the approximation algorithm derived in Subsections (2.1)–(2.6).

2.7 Description of the algorithm in a special case

In this subsection, we provide a formal description of the algorithm derived in Subsections (2.1)–(2.6) in the case where the standard Euler–Maruyama scheme (cf., e.g., [66, 81, 83]) is used to approximate the solutions $Y^{m}$ of the SDEs (2.30) and optimal parameters $\theta_{1},\dots,\theta_{N}\in{\mathbb{R}}^{\nu}$ are computed with plain vanilla stochastic gradient descent with a constant learning rate $\gamma\in(0,\infty)$ and batch size 1. Note that for the algorithm to be implementable, it is enough if the initial condition $\varphi$ has a weak gradient $\nabla\varphi\colon{\mathbb{R}}^{d}\to{\mathbb{R}}^{d}$ which does not need to satisfy any growth conditions.

In the following Framework 2.1, feedforward neural networks of the form (2.25) are employed to approximate the solution of the PDE. (2.39) describes a stochastic gradient descent scheme with constant learning rate $\gamma$ and (2.40) specifies the quadratic loss functions.

Framework 2.1** (Special case).**

Assume $\varphi$ has a weak gradient $\nabla\varphi\colon{\mathbb{R}}^{d}\to{\mathbb{R}}^{d}$. Consider $N\in{\mathbb{N}}$ and $t_{0},t_{1},\ldots,t_{N}\in[0,T]$ such that

[TABLE]

Set $\tau_{n}=T-t_{N-n}$ for $n\in\{0,1,\dots,N\}$. For a given $M\in{\mathbb{N}}$, let $B^{m}\colon[0,T]\times\Omega\to{\mathbb{R}}^{d}$, $m\in\{0,1,\dots,M-1\}$, be i.i.d. standard Brownian motions on a probability space $(\Omega,{\mathcal{F}},\mathbbm{P})$. Consider i.i.d. random variables $\xi^{m}\colon\Omega\to{\mathbb{R}}^{d}$, $m\in\{0,1,\dots,M-1\}$, that are independent of $B^{m}$, $m\in\{0,1,\dots,M-1\}$, and let the stochastic processes ${\mathcal{Y}}^{m}\colon\{0,1,\ldots,N\}\times\Omega\to{\mathbb{R}}^{d}$, $m\in\{0,1,\dots,M-1\}$, be given by ${\mathcal{Y}}^{m}_{0}=\xi^{m}$ and

[TABLE]

Let $V_{n}\colon{\mathbb{R}}^{\nu}\times{\mathbb{R}}^{d}\to{\mathbb{R}}$, $n\in\{0,1,\dots,N\}$ be the functions given in (2.24)–(2.25). Consider $\gamma\in(0,\infty)$ and let $\Theta^{n}\colon\{0,1,\dots,M\}\times\Omega\to{\mathbb{R}}^{\nu}$, $n\in\{0,1,\ldots,N\}$, be stochastic processes satisfying for all $m\in\{0,1,\dots,M-1\}$ and $n\in\{1,\dots,N\}$,

[TABLE]

where for all $\omega\in\Omega$,

[TABLE]

and $\Phi^{n,m}(\theta,\omega)=\nabla_{\theta}\phi^{n,m}(\theta,\omega)$.

In the setting of Framework 2.1 the random variables $V_{n}(\Theta^{n}_{M},x)\colon\Omega\to{\mathbb{R}}$ provide for all $n\in\{1,2,\dots,N\}$ and $x\in{\mathbb{R}}^{d}$ the approximations

[TABLE]

to the solution $u\colon[0,T]\times{\mathbb{R}}^{d}\to{\mathbb{R}}$ of the PDE (2.1).

2.8 Description of the algorithm in the general case

We now generalize Framework 2.1 so that, besides plain vanilla stochastic gradient descent with constant learning rate and batch size $1$, it also covers more advanced machine learning techniques such as mini-batches, batch normalization and more sophisticated updating rules.

In Framework 2.2 below, functions $V^{j,s}_{n}\colon{\mathbb{R}}^{\nu}\times{\mathbb{R}}^{d}\to{\mathbb{R}}$ parametrized by $(j,s,n)\in{\mathbb{N}}\times{\mathbb{R}}^{\varsigma}\times\{0,1,\dots,N\}$ for some number $\varsigma\in{\mathbb{N}}$, are used to approximate the solution of the PDE. The additional parameters $j$ and $s$ make it possible to describe mini-batches and batch normalization; see Ioffe & Szegedy [61]. As in Framework 2.1, the standard Euler–Maruyama scheme is used to approximate the solutions $Y^{m}$ of the SDEs (2.30). But a different approximation scheme could be employed as well. In (2.43) the quadratic loss functions are given that are used for training the functions $V^{j,s}_{n}$, whereas (2.44) specifies the gradients of the loss functions. The role of the stochastic processes $\mathbb{S}^{n}$ in (2.45) is to describe the variables (running mean and standard deviation) needed for batch normalization. The stochastic processes $\Theta^{n}$ and $\Xi^{n}$ in (2.46) describe the updating rule. Their dynamics are specified by the functions $\Psi^{n}_{m}$ and $\psi^{n}_{m}$. Since we make no assumptions on $\Psi^{n}_{m}$ and $\psi^{n}_{m}$, the framework includes very general stochastic gradient optimization methods. In the examples in Section 3 we use Adam optimization; see Kingma & Ba [64]. The corresponding specification of the functions $\Psi^{n}_{m}$ and $\psi^{n}_{m}$ is given in (3.1)–(3.2) below.

Framework 2.2** (General case).**

Assume $\varphi$ has a weak gradient $\nabla\varphi\colon{\mathbb{R}}^{d}\to{\mathbb{R}}^{d}$. Consider $N\in{\mathbb{N}}$ and $t_{0},t_{1},\dots,t_{N}\in[0,T]$ such that $0=t_{0}<t_{1}<\ldots<t_{N}=T$. Set $\tau_{n}=T-t_{N-n}$ for $n\in\{0,1,\dots,N\}$. For a given $M\in{\mathbb{N}}$, let, for every $n\in\{1,2,\ldots,N\}$, $B^{n,m,j}\colon[0,T]\times\Omega\to{\mathbb{R}}^{d}$, $m\in\{0,1,\dots,M-1\}$, $j\in{\mathbb{N}}$, be i.i.d. standard Brownian motions on a probability space $(\Omega,{\mathcal{F}},\mathbbm{P})$ and $\xi^{n,m,j}\colon\Omega\to{\mathbb{R}}^{d}$, $m\in\{0,1,\dots,M-1\}$, $j\in{\mathbb{N}}$, i.i.d. random variables that are independent of $B^{n,m,j}$, $m\in\{0,1,\dots,M-1\}$, $j\in{\mathbb{N}}$. Let ${\mathcal{Y}}^{n,m,j}\colon\{0,1,\ldots,N\}\times\Omega\to{\mathbb{R}}^{d}$, $n\in\{1,2,\ldots,N\}$, $m\in\{0,1,\dots,M-1\}$, $j\in{\mathbb{N}}$, be stochastic processes given by ${\mathcal{Y}}^{n,m,j}_{0}=\xi^{n,m,j}$ and

[TABLE]

Let $\nu,\varsigma,\varrho,J_{0},\dots,J_{M-1}\in{\mathbb{N}}$. Consider functions $V^{j,s}_{n}\colon{\mathbb{R}}^{\nu}\times{\mathbb{R}}^{d}\to{\mathbb{R}}$, $(j,s,n)\in{\mathbb{N}}\times{\mathbb{R}}^{\varsigma}\times\{0,1,\dots,N\}$, such that $V^{j,s}_{0}(\theta,x)=\varphi(x)$ for all $(j,s,\theta,x)\in{\mathbb{N}}\times{\mathbb{R}}^{\varsigma}\times{\mathbb{R}}^{\nu}\times{\mathbb{R}}^{d}$, and let $\Theta^{n}\colon\{0,1,\dots,M-1\}\times\Omega\to{\mathbb{R}}^{\nu}$, $n\in\{0,1,\dots N\}$, be stochastic processes. For all $n\in\{1,2,\ldots,N\}$, $m\in\{0,1,\dots,M-1\}$ and $s\in{\mathbb{R}}^{\varsigma}$, let the mapping $\phi^{n,m,s}\colon{\mathbb{R}}^{\nu}\times\Omega\to{\mathbb{R}}$ be given by

[TABLE]

and assume $\Phi^{n,m,s}\colon{\mathbb{R}}^{\nu}\times\Omega\to{\mathbb{R}}^{\nu}$ is a function satisfying

[TABLE]

for all $\omega\in\Omega$ and $\theta\in\{\vartheta\in{\mathbb{R}}^{\nu}\colon\phi^{n,m,s}(\cdot,\omega)\colon{\mathbb{R}}^{\nu}\to{\mathbb{R}}~{}\text{is differentiable at}~{}\vartheta\}$. For every $n\in\{1,2,\ldots,N\}$ and $m\in\{0,1,\dots,M-1\}$, let $\mathcal{S}^{n}\colon{\mathbb{R}}^{\varsigma}\times{\mathbb{R}}^{\nu}\times({\mathbb{R}}^{d})^{\{0,1,\ldots,N\}\times{\mathbb{N}}}\to{\mathbb{R}}^{\varsigma}$, $\Psi^{n}_{m}\colon{\mathbb{R}}^{\varrho}\times{\mathbb{R}}^{\nu}\to{\mathbb{R}}^{\varrho}$ and $\psi^{n}_{m}\colon{\mathbb{R}}^{\varrho}\to{\mathbb{R}}^{\nu}$ be functions and $\mathbb{S}^{n}\colon\{0,1,\dots,M-1\}\times\Omega\to{\mathbb{R}}^{\varsigma}$ and $\Xi^{n}\colon\{0,1,\dots,M-1\}\times\Omega\to{\mathbb{R}}^{\varrho}$ stochastic processes such that

[TABLE]

[TABLE]

In the setting of Framework 2.2 the functions $V^{1,\mathbb{S}_{M}^{n}}_{n}(\Theta^{n}_{M},x)\colon\Omega\to{\mathbb{R}}$ yield the approximations

[TABLE]

of the solution $u\colon[0,T]\times{\mathbb{R}}^{d}\to$ to the PDE (2.1).

3 Examples

We now illustrate the performance of the deep splitting method on five concrete example PDEs. In each example we use the general approximation method of Framework 2.2 with approximating functions $V^{j,s}_{n}\colon{\mathbb{R}}^{\nu}\times{\mathbb{R}}^{d}\to{\mathbb{R}}$, $n\in\{1,2,\dots,N\}$, specified as feedforward neural networks with 4 layers (1 input layer, 2 hidden layers, 1 output layer) and ReLU-activation $\rho(x)=\max\{0,x\}$, $x\in{\mathbb{R}}$. We use mini-batches of size $J_{m}=256$ and apply batch normalization before the first affine transformation, before each of the two nonlinear activation functions in front of the hidden layers, and just before the output layer. We use Xavier initialization (see Glorot & Bengio [35]) to initialize all weights in the neural networks together with Adam optimization (see Kingma & Ba [64]) with parameters $\varepsilon=10^{-8}$, $\beta_{1}=0.9$, $\beta_{2}=0.999$ and decreasing learning rates $(\gamma_{m})_{m=0}^{M-1}$ that we choose depending on the form and dimension of the problem. More precisely, we set $\varrho=2\nu$ and denote by $\operatorname{Pow}_{2}\colon{\mathbb{R}}^{\nu}\to{\mathbb{R}}^{\nu}$ the function given by $\operatorname{Pow}_{2}(\eta_{1},\dots,\eta_{\nu})=(\eta_{1}^{2},\dots,\eta^{2}_{\nu})$. Then, Adam optimization corresponds to the following specification of the two functions $\Psi^{n}_{m}\colon{\mathbb{R}}^{3\nu}\to{\mathbb{R}}^{2\nu}$ and $\psi^{n}_{m}\colon{\mathbb{R}}^{2\nu}\to{\mathbb{R}}^{\nu}$ from Framework (2.2):

[TABLE]

and

[TABLE]

In the examples in the following subsections, we approximate $u(T,x)$ for different $T\in(0,\infty)$ and $x\in{\mathbb{R}}^{d}$ with $V^{1,\mathbb{S}^{N}_{M}}_{N}(\Theta^{N}_{M},x)$, which, due to the stochastic gradient optimization method, is a random variable. In each example we report estimates of the expectation and standard deviation of $V^{1,\mathbb{S}^{N}_{M}}_{N}(\Theta^{N}_{M},x)$. We also give relative $L^{1}$-approximation errors with respect to reference values calculated with different alternative methods together with their uncorrected sample standard deviations. The average runtimes needed for calculating one realization of $V^{1,\mathbb{S}^{N}_{M}}_{N}(\Theta^{N}_{M},x)$ are in seconds and were determined as averages over 10 independent runs.

All numerical experiments presented in this paper were implemented in Python using TensorFlow and run on a NVIDIA GeForce GTX 1080 GPU with 1974 MHz core clock and 8 GB GDDR5X memory with 1809.5 MHz clock rate and an underlying system consisting of an Intel Core i7-6800K 3.4 GHz CPU with 64 GB DDR4-2133 memory running TensorFlow 1.5 on Ubuntu 16.04. The Python source codes can be found at https://github.com/seb-becker/deep$\_$pde.

3.1 Hamilton–Jacobi–Bellman (HJB) equations

As a first example, we calculate approximate solutions of the PDE

[TABLE]

with initial condition $u(0,x)=\|x\|^{1/2}_{{\mathbb{R}}^{d}}$ for different $d\in{\mathbb{N}}$. The deep splitting method can be applied to more general HJB equations. But (3.3) has the advantage that it reduces to a linear heat equation under a logarithmic transformation; see, e.g., E et al. [28, Lemma 4.2]. Since solutions of linear equations can be approximated with standard Monte Carlo, this allows us to efficiently compute reference solutions in high dimensions.

Table 1 shows deep splitting approximations of $u(T,0,\dots,0)$ for different values of $d$, $T$ and $N$. We derived it with $\mu(x)=(0,0,\dots,0)\in{\mathbb{R}}^{d}$, $\sigma(x)=\sqrt{2}\mbox{ Id}_{{\mathbb{R}}^{d\times d}}$, and $f(x,y,z)=-\|z\|^{2}_{{\mathbb{R}}^{d}}$ for $(x,y,z)\in{\mathbb{R}}^{d}\times{\mathbb{R}}\times{\mathbb{R}}^{d}$. We used $M=500+100\,\mathbbm{1}_{\{10,000\}}(d)$ and

[TABLE]

We set $\xi^{n,m,j}=(0,0,\dots,0)\in{\mathbb{R}}^{d}$ for every $(n,m,j)\in\{1,2,\dots,N\}\times\{0,1,\dots,M-1\}\times{\mathbb{N}}$, and used feedforward neural networks with a $d$-dimensional input layer, two hidden layers of dimension $d+10$, and a one-dimensional output layer.

The reference values for $u(T,0,0,\dots,0)$ were calculated using a logarithmic transformation and a standard Monte Carlo method; see Han et al. [28, Lemma 4.2].

3.2 Nonlinear Black–Scholes equations

There exist a number of extensions of the classical linear Black–Scholes equation which incorporate nonlinear phenomena such as transaction costs, default risk or Knightian uncertainty. We here consider nonlinear Black–Scholes equations of the form

[TABLE]

$(t,x)\in(0,T]\times{\mathbb{R}}^{d}$, for suitable parameters $\delta$, $R$, $\gamma^{h}$, $\gamma^{l}$, $v^{h}$, $v^{l}$, $\bar{\mu}$, $\bar{\sigma}\in{\mathbb{R}}$. They describe derivative prices under default risk; see, e.g., Han et al. [48, Subsection 3.1] and E et al. [30, Subsection 3.1], from where we adopt the initial condition $u(0,x)=\min_{i\in\{1,2,\dots,d\}}x_{i}$ and the choice of the parameter values $\delta=2/3$, $R=0.02$, $\gamma^{h}=0.2$, $\gamma^{l}=0.02$, $v^{h}=50$, $v^{l}=70$, $\bar{\mu}=0.02$, $\bar{\sigma}=0.2$.

Table 2 reports deep splitting approximations of $u(T,50,\dots,50)$ for $T=1/3$, $N=96$, and different values of $d$. We chose $\mu(x)=\bar{\mu}x$, $\sigma(x)=\bar{\sigma}x$, and

[TABLE]

$(x,y,z)\in{\mathbb{R}}^{d}\times{\mathbb{R}}\times{\mathbb{R}}^{d}$ and used $M=2000+1000\,\mathbbm{1}_{[0,100]}(d)$ together with

[TABLE]

We set $\xi^{n,m,j}=(50,50,\dots,50)\in{\mathbb{R}}^{d}$ for every $(n,m,j)\in\{1,2,\dots,N\}\times\{0,1,\dots,M-1\}\times{\mathbb{N}}$, and used feedforward neural networks with a $d$-dimensional input layer, two hidden layers of dimension $d+10+40\,\mathbbm{1}_{[1,100]}(d)$ and a one-dimensional output layer.

The reference values for $u(T,50,50,\dots,50)$ were calculated with the deep BSDE method of E et al. [28]).

3.3 Allen–Cahn-type equations

Next, we approximate solutions of high-dimensional Allen–Cahn-type equations with a cubic nonlinearity of the form

[TABLE]

with initial condition $u(0,x)=\arctan(\max_{i\in\{1,2,\dots,d\}}x_{i})$; see also Beck et al. [3, Section 4.1], E et al. [28, Section 4.2], E et al. [30, Section 3.4], and Han et al. [48, Section 3.3] for further numerical results for equation (3.8).

Table 3 lists deep splitting approximations of $u(T,0,\dots,0)$ for $T=0.3$, $N=10$, and different values of $d$. We chose $\mu(x)=(0,0,\dots,0)\in{\mathbb{R}}^{d}$, $\sigma(x)=\sqrt{2}\mbox{ Id}_{{\mathbb{R}}^{d\times d}}$ and $f(x,y,z)=y-y^{3}$, $(x,y,z)\in{\mathbb{R}}^{d}\times{\mathbb{R}}\times{\mathbb{R}}^{d}$. We used $M=500$ and

[TABLE]

We set $\xi^{n,m,j}=(0,0,\dots,0)\in{\mathbb{R}}^{d}$ for all $(n,m,j)\in\{1,2,\dots,N\}\times\{0,1,\dots,M-1\}\times{\mathbb{N}}$ and used feedforward neural networks with a $d$-dimensional input layer, two hidden layers of dimension $d+10$ and a one-dimensional output layer.

The reference values for $u(T,0,0,\dots,0)$ were calculated with the multilevel Picard method; see, e.g., [5, 29, 30, 57, 59, 60].

3.4 Semilinear heat equations

In this subsection, we calculate approximate solutions of semilinear heat equations of the form

[TABLE]

with initial condition $u(0,x)=\nicefrac{{5}}{{\left(10+{2}\|x\|_{{\mathbb{R}}^{d}}^{2}\right)}}$.

Table 4 shows deep splitting approximations of $u(T,0,\dots,0)$ for $T=0.3$, $N=20$, and different values of $d$. We derived them with $\mu(x)=(0,0,\dots,0)\in{\mathbb{R}}^{d}$, $\sigma(x)=\sqrt{2}\mbox{ Id}_{{\mathbb{R}}^{d\times d}}$, and $f(x,y,z)=-\|z\|^{2}_{{\mathbb{R}}^{d}}$ for $(x,y,z)\in{\mathbb{R}}^{d}\times{\mathbb{R}}\times{\mathbb{R}}^{d}$. We used $M=500$ and

[TABLE]

We set $\xi^{n,m,j}=(0,0,\dots,0)\in{\mathbb{R}}^{d}$ for all $(n,m,j)\in\{1,2,\dots,N\}\times\{0,1,\dots,M-1\}\times{\mathbb{N}}$, and used feedforward neural networks with a $d$-dimensional input layer, two hidden layers of dimension $d+10$, and a one-dimensional output layer. The reference values were computed with the multilevel Picard method; see, e.g., [5, 29, 30, 57, 59, 60]).

To illustrate how the accuracy of the deep splitting method depends on the number of time steps $N$, we report in Table 5, deep splitting approximations of $u(T,0,\dots,0)$ for $d=100$, $T=0.3$, and different values of $N$. It can be seen that the decrease of the relative $L^{1}$-error is approximately linear in $N$. For theoretical convergence results for the deep splitting method, we refer to Beck et al. [4, Theorem 2].

3.5 Sine-Gordon-type equations

As a last example, we approximate solutions of sine-Gordon-type equations of the form

[TABLE]

with initial condition $u(0,x)=\nicefrac{{5}}{{\left(10+{2}\|x\|_{{\mathbb{R}}^{d}}^{2}\right)}}$.

Table 6 shows deep splitting approximations of $u(T,0,\dots,0)$ for $T=0.3$, $N=20$ and different values of $d$. We chose $\mu(x)=(0,0,\dots,0)\in{\mathbb{R}}^{d}$, $\sigma(x)=\sqrt{2}\mbox{ Id}_{{\mathbb{R}}^{d\times d}}$, and $f(x,y,z)=\sin(y)$, $(x,y,z)\in{\mathbb{R}}^{d}\times{\mathbb{R}}\times{\mathbb{R}}^{d}$. We used $M=1000$ and

[TABLE]

We set $\xi^{n,m,j}=(0,0,\dots,0)\in{\mathbb{R}}^{d}$ for every $(n,m,j)\in\{1,2,\dots,N\}\times\{0,1,\dots,M-1\}\times{\mathbb{N}}$ and used feedforward neural networks with a $d$-dimensional input layer, two hidden layers of dimension $d+50$ and a one-dimensional output layer.

The reference values for $u(T,0,0,\dots,0)$ were computed with the multilevel Picard method; see, e.g., [29, 30, 57, 59, 60].

4 Comparison with other methods

In this section, we compare the deep splitting method with the deep BSDE method of E et al. [28] and the multilevel Picard method of [5, 57], which also have been shown to produce good results in approximating solutions of high-dimensional PDEs.

Figure 1 shows estimated relative $L^{1}$-approximation errors as a function of the number of one-dimensional standard normal random variables used by the three methods to approximate $u(T,0,\dots,0)$ for the solution $u\colon[0,T]\times{\mathbb{R}}^{d}\to{\mathbb{R}}$ of the sine-Gordon-type equation (3.12) for $d=10$ and $T=0.3$. It can be seen that for this particular example, the three methods yield comparable results. But it has to be noted that the deep splitting algorithm and the deep BSDE method both involve different hyper-parameters, which, for good results, have to be fine-tuned depending on the form and the parameters of the PDE. Figure 1 just shows results for particular implementations of the three methods.

Generally, since the deep splitting method uses neural networks to approximate the solution $u\colon[0,T]\times{\mathbb{R}}^{d}\to{\mathbb{R}}$ of a PDE on a time grid $0=t_{0}<t_{1}<\dots<t_{N}=T$, it can learn approximations of $u(t_{n},x)$, $n\in\{1,2,\dots,N\}$, simultaneously for all $x\in{\mathbb{R}}^{d}$. Similarly, the deep BSDE method can be implemented so that it approximates $u(T,x)$ directly for all $x\in{\mathbb{R}}^{d}$. However, for a temporal discretization with $N$ subintervals, it then needs to train $N$ neural networks at the same time, whereas the deep splitting method trains one network after the other. So even if the two approaches need a similar total number of one-dimensional standard normal random variables to achieve a given accuracy, the deep splitting method can handle larger problems since it decomposes them into smaller computational tasks which can be solved successively.

The advantage of the multilevel Picard method is that there exist theoretical bounds on the computational effort needed for a given approximation accuracy; see e.g., [29, 57, 59, 60]. However, the method needs to calculate approximations of $u(t,x)$ for different space-time points $(t,x)$ separately and becomes impractical for large $t$.

5 Conclusion

In this paper we have developed a new numerical method to approximate solutions of high-dimensional nonlinear parabolic PDEs. It splits the differential operator into a linear and a nonlinear part and uses deep learning together with the Feynman–Kac formula to iteratively solve linear approximations of the equation over small time intervals. This breaks the PDE approximation task into smaller problems that can be solved successively. As a consequence, the approach can be applied to extremely high-dimensional nonlinear PDEs. We have tested the method on Hamilton–Jacobi–Bellman equations, nonlinear Black–Scholes equations, Allen–Cahn-type equations, semilinear heat equations as well as sine-Gordon-type equations. In all cases, it has produced accurate results in high dimensions with short run times.

Acknowledgments

We are grateful to Adam Andersson for fruitful discussions. This project has been supported by the Swiss National Science Foundation Grant $200020\_175699$ “Higher order numerical approximation methods for stochastic partial differential equations”, by the Deutsche Forschungsgemeinschaft under Germany’s Excellence Strategy EXC 2044-390685587, Mathematics Münster: Dynamics - Geometry - Structure and by the Nanyang Assistant Professorship Grant “Machine Learning based Algorithms in Finance and Insurance”.