remarkRemark

# DeepMartNet - A Martingale Based Deep Neural Network Learning Method for
Dirichlet BVPs and Eigenvalue Problems of Elliptic PDEs in $R^{{d}}$ ††thanks:
December 20, 2023, first version appeared in arxiv preprint arXiv:2311.09456.
2023 Nov 15.

Wei Cai Corresponding author, Department of Mathematics, Southern Methodist
University, Dallas, TX 75275<EMAIL_ADDRESS>Andrew He Department of
Mathematics, Southern Methodist University, Dallas, TX 75275. Daniel Margolis
Department of Mathematics, Southern Methodist University, Dallas, TX 75275.

###### Abstract

In this paper, we propose DeepMartNet - a Martingale based deep neural network
learning method for solving Dirichlet boundary value problems (BVPs) and
eigenvalue problems for elliptic partial differential equations (PDEs) in high
dimensions or domains with complex geometries. The method is based on
Varadhan’s Martingale problem formulation for the BVPs/eigenvalue problems
where a loss function enforcing the Martingale property for the PDE solution
is used for an efficient optimization by sampling the stochastic processes
associated with corresponding elliptic operators. High dimensional numerical
results for BVPs of the linear and nonlinear Poisson-Boltzmann equation and
eigenvalue problems of the Laplace equation and a Fokker-Planck equation
demonstrate the capability of the proposed DeepMartNet learning method in
solving high dimensional PDE problems.

###### keywords:

Martingale problem, Deep neural network, boundary value problems, Eigenvalue
problems.

35Q68, 65N35, 65T99,65K10

##  1 Introduction

Computing eigenvalues and/or eigenfunctions for elliptic operators or solving
boundary value problem of partial differential equations (PDEs) in high
dimensions and optimal stochastic control problems are among the key tasks for
many scientific computing applications, e.g., ground states and band structure
calculation in quantum systems, complex biochemical systems, communications
and manufacturing productions, etc. Deep neural network (DNN) has been
utilized to solve high dimensional PDEs and stochastic control problems due to
its capability of approximating high dimensional functions. The first attempt
to use DNN to solve high dimensional quasi-linear parabolic PDE was carried
out in [7] using Pardoux-Peng theory [19] of nonlinear Feynman-Kac formulas,
which connect the solution of the parabolic PDEs with that of backward
stochastic differential equations (SDEs) [14]. The loss function for the
training of the DeepBSDE uses the terminal condition of the PDEs. The DeepBSDE
framework was also applied to solve stochastic control problems [8]. The
DeepBSDE started a new approach of SDE based DNN approximations to high
dimensional PDEs and stochastic control problems. The work in [21] [24]
extended this idea to use a loss function based on the pathwise comparison of
two stochastic processes, one from the BSDE in the Pardoux-Peng theory and one
from the PDEs solution, this variant of SDE based DNNs has the potential to
find the solution in the whole domain compared with one point solution in the
original DeepBSDE. Recently, a diffusion Monte Carlo DNN method was developed
[9] using the connection between stochastic processes and the solution of
elliptic equations and the backward Kolmogorov equation to build a loss
function for eigen-value calculations. The capability of SDE paths exploring
high dimensional spaces make them a good candidate to be used for DNN
learning. It should be mentioned that other approaches to solving high
dimensional PDEs include the Feynmann-Kac formula based Picard iteration [13]
[6] and stochastic dimension gradient descent method [12].

In addition to the Pardoux-Peng BSDE approach by linking stochastic processes
and the solution of PDEs, another powerful probabilistic method is the
Varadhan’s Martingale problem approach [23, 15], which were used to derive a
probabilistic weak form for the PDE’s solution with a specific Martingale
related to the PDE solution through the Ito formulas [17]. The equivalence of
the classic weak solution of boundary value problems (BVPs) of elliptic PDEs
using bilinear forms and the probabilistic one were established for the
Schrodinger equation for the Neumann BVP [10, 11], and then, for the Robin BVP
[18]. This Martingale weak form for the PDEs solution and in fact also for a
wide class of stochastic control problem [4, 3] provides a new venue to tackle
high dimensional PDE and stochastic control problems [2], and this paper
focuses on the case of high dimensional PDEs only. The Martingale formulation
of the PDE solution is a result of Ito formula and the Martingale nature of
Ito integrals. As a simple conclusion from the broader Martingale property,
the Feynman-Kac formula provides one point solution of the PDE using
expectation of the underlying SDE paths originating from that point. The
Martingale based DNN to be studied in this paper, termed DeepMartNet [2], is
trained based on a loss function enforcing the conditional expectation
definition of the Martingale. It will be shown with extensive numerical tests
that using the same set of SDE paths originating from one single point, the
DeepMartNet can in fact provide approximation to solutions of the BVPs and
eigenvalue problems of elliptic PDEs over the whole solution doamin in high
dimensions. In this sense, the DeepMartNet is able to extract more information
for the PDE solutions from the SDE paths originating from one point than the
(pre-machine learning) traditional use of the one point solution Feynman-Kac
formula.

The rest of this paper is organized as follows. In section 2, a brief review
of existing SDE-based DNN methods for solving PDEs is given and Section 3 will
present the Martingale problem formulation for the BVP problem of a general
elliptic PDE for the case of third kind Robin boundary condition, which
includes both Dirichlet and Neumann BVPs as special limiting cases. Section 4
will present the Martingale based DeepMartNet for solving high dimensional PDE
problems. Numerical results for the Dirichlet BVPs and eigenvalue problems
will be present in Section 5. The implementation of the DeepMartNet for the
Neumann and Robin boundary conditions will be addressed in a follow-up paper,
which will involve reflecting diffusion processes in finite domains and the
computation of local times of the processes. Section 6 will present
conclusions and future work.

##  2 A review of SDE based DNNs for solving PDEs

To set the background for the Martingale based DNNs, we will first briefly
review some existing DNN based on diffusion paths from SDEs.

Let us first consider a terminal value problem for quasi-linear elliptic PDEs

(1) $\partial_{t}u+\mathcal{L}u=\phi,\quad\textbf{x}\in R^{d},$

with a terminal condition $u(T,\mathbf{x})=g(\mathbf{x})$, where the
differential operator $\mathcal{L}$ is given as

(2)
$\mathcal{L}=\mu^{\top}\nabla+\frac{1}{2}Tr(\sigma\sigma^{\top}\nabla\nabla^{\top})=\mu^{\top}\nabla+\frac{1}{2}Tr(A\nabla\nabla^{\top}),$

and $\mu=\mu(t,\textbf{x},u,\nabla u),\sigma=\sigma(t,\textbf{x},u,\nabla u)$
and the diffusion coefficient matrix

(3) $A=(a_{ij})_{d\times d}=\sigma\sigma^{\top}.$

The aim is to find the solution at $\mathbf{x},t=0,$ $u(0,\mathbf{x})$, and
the solution of (1) is related to a coupled FBSDE [19]

(4)
$\displaystyle\begin{split}d\mathbf{X}_{t}&=\mu(t,\textbf{X}_{t},Y_{t},\textbf{Z}_{t})dt+\sigma(t,\textbf{X}_{t},Y_{t})d\textbf{B}_{t},\\\
\textbf{X}_{0}&=\xi,\end{split}$ (5)
$\displaystyle\begin{split}dY_{t}&=\phi(t,\textbf{X}_{t},Y_{t},\textbf{Z}_{t})dt+\textbf{Z}^{T}_{t}\sigma(t,\textbf{X}_{t},Y_{t})d\textbf{B}_{t},\\\
Y_{T}&=g(\textbf{X}_{T}),\end{split}$

where $\textbf{X}_{t}$, $Y_{t}$ and $\textbf{Z}_{t}$ are $d$, 1 and
$d$-dimensional stochastic processes, respectively, that are adapted to
$\\{\mathcal{F}_{t}:0\leq t\leq T\\}$ \- the natural filtration from the
d-dimensional Brownian motion $\textbf{B}_{t}$ . Specifically, we have the
following relations,

(6) $Y_{t}=u(t,\textbf{X}_{t}),\quad\textbf{Z}_{t}=\nabla
u(t,\textbf{X}_{t}).$

DeepBSDE. As the first work of using SDEs to train DNN, the Deep BSDE [7]
trains the network with input $\textbf{X}_{0}=\textbf{x}$ and output
$\textbf{Y}_{0}=u(0,\mathbf{x})$. Applying the Euler–Maruyama scheme (EM) to
the FBSDE (4) and (5), respectively, we have

(7) $\displaystyle\textbf{X}_{n+1}$
$\displaystyle\approx\textbf{X}_{n}+\mu(t_{n},\textbf{X}_{n},\textbf{Y}_{n},\textbf{Z}_{n})\Delta
t_{n}+\sigma(t_{n},\textbf{X}_{n},\textbf{Y}_{n})\Delta\textbf{B}_{n},$ (8)
$\displaystyle Y_{n+1}$ $\displaystyle\approx
Y_{n}+\phi(t_{n},\textbf{X}_{n},Y_{n},\textbf{Z}_{n})\Delta
t_{n}+Z^{T}_{n}\sigma(t_{n},\textbf{X}_{n},Y_{n})\Delta\textbf{B}_{n}.$

The missing $\textbf{Z}_{n+1}$ at $t_{n+1}$ will be then approximated by a
neural network (NN) with parameters $\theta_{n}$

(9) $\nabla u(t_{n},\textbf{X}_{n}|\theta_{n})\sim\textbf{Z}_{n}=\nabla
u(t_{n},\textbf{X}_{n}).$

Loss function: with an ensemble average approximation, the loss function is
defined as

(10) $Loss_{bsde}(Y_{0},\theta)=E\left\lVert
u(T,\textbf{X}_{T})-g(\textbf{X}_{T})\right\rVert^{2},$

where

$u(T,\textbf{X}_{T})=Y_{N}.$

Trainable parameters are $\\{Y_{0},\theta_{n},n=1,\cdots,N\\}$.

FBSNN. A forward backward neural network (FBSNNs) proposed in [21] uses the
mismatch between two stochastic processes to build the loss function for the
DNN, which aims to train a DNN $u_{\theta}(x,t)$ in the whole domain. The
following is an improved version of the approach in [24].

  * •

Markov chain one. Starting with $\textbf{X}_{0}=\textbf{x}$ ,
$Y_{0}=u_{\theta}(\textbf{x},0)$,

(11)
$\displaystyle\begin{split}\textbf{X}_{n+1}&=\textbf{X}_{n}+\mu(t_{n},\textbf{X}_{n},Y_{n},Z_{n})\Delta
t_{n}+\sigma(t_{n},\textbf{X}_{n},Y_{n})\Delta\textbf{B}_{n},\\\
Y_{n+1}&=Y_{n}+\phi(t_{n},\textbf{X}_{n},Y_{n},\textbf{Z}_{n})\Delta
t_{n}+\textbf{Z}^{T}_{n}\sigma(t_{n},\textbf{X}_{n},Y_{n})\Delta\textbf{B}_{n},\\\
\textbf{Z}_{n+1}&{}={}\nabla u(t_{n+1},\textbf{X}_{n+1}).\end{split}$

  * •

Markov chain two

(12) ${Y^{\star}_{n+1}=u(t_{n+1},\textbf{X}_{n+1}).}$

  * •

The loss function is a Monte Carlo approximation of

(13) $E\left[\frac{1}{N}\sum_{n=1}^{N}\left\lVert
Y_{n}-Y^{\star}_{n}\right\rVert^{2}+0.02\left\lVert
Y^{\star}_{N}-g(\textbf{X}_{N})\right\rVert^{2}+0.02\left\lVert\textbf{Z}_{N}-\nabla
g(\textbf{X}_{N})\right\rVert^{2}\right].$

Numerical results of half-order convergence of $u_{\theta}(\textbf{x},t)$,
similar to the order of the underlying Euler-Maruyama scheme, is observed.

Diffusion Monte Carlo DNN eigensolver. In another approach similar to the
power iteration method, a diffusion-Monte Carlo method was proposed [9]
through a fixed point of semi-group formulation for the eigenvalue problem of
the following linear elliptic operator (i.e.
$\mu=\mu(\textbf{x}),\sigma=\sigma(\textbf{x})$),

(14) $\mathcal{L}\Psi=\lambda\Psi.$

Equation (14) can be reformulated as a virtual time dependent backward
parabolic PDE with the sought-after eigenfunction as the terminal condition,
i.e.,

(15) $\partial_{t}u(t,\textbf{x})+\mathcal{L}u(t,\textbf{x})-\lambda
u(t,\textbf{x})=0,\quad u(T,\textbf{x})=\Psi(\textbf{x}).$

Thus, the following fix-point property holds,

(16) $u(T-t,\cdot)=P^{\lambda}_{t}\Psi,\quad P^{\lambda}_{T}\Psi=\Psi,$

where the semi-group for the evolutionary system is formally defined by

(17) $P^{\lambda}_{t}=e^{-(T-t)\mathcal{L}}.$

The discretized backward in time evolution (with time step $\Delta t$ of the
backward parabolic equation) mimics the power iteration of the semi-group
operator $e^{-n\Delta t\mathcal{L}}$, which will converge to the lowest
eigenfunction for $\mathcal{L}$ as in a power method. The loss function is set
to be $||P^{\lambda}_{T}\Psi-\Psi||^{2}$ , while the evolution of PDE solution
is done by two SDEs, instead of solving the parabolic equation directly,

$\displaystyle\textbf{X}_{n+1}=$
$\displaystyle\textbf{X}_{n}+\sigma\Delta\textbf{B}_{n},$ (18) $\displaystyle
u_{n+1}=$ $\displaystyle
u_{n}+(\lambda\Psi_{\theta}-\mu^{T}\nabla\Psi_{\theta})(\textbf{X}_{n})\Delta
t+\nabla\Psi_{\theta}(\textbf{X}_{n})\Delta\textbf{B}_{n}.$

The original algorithm in [9] uses a separate DNN to approximate the gradient
of $\Psi(\textbf{x})$. And the loss function for the DNN
$\Psi_{\theta}(\textbf{x})$ approximating the eigen-function and the
eigenvalue is then defined by

(19)
$Loss_{semigroup}(\theta,\lambda)=E_{X_{0}\sim\pi_{0}}[|u_{N}-\Psi_{\theta}(\textbf{X}_{N})|^{2}].$

##  3 Martingale problem formulation of elliptic PDEs

In this section, we will present the Martingale problem formulation for the
BVPs and eigenvalue problems of elliptic PDEs, and for this purpose, let us
consider a general PDE for the linear elliptic operator $\mathcal{L}$ with
$\mu=\mu(\mathbf{x}),\sigma=\sigma(\mathbf{x})$,

(20) $\displaystyle\mathcal{L}u+V(\mathbf{x},u,\nabla u)$
$\displaystyle=f(\textbf{x},u),\text{ \ }\mathbf{x}\in D\subset R^{d},$ (21)
$\displaystyle\mbox{or}\quad\mathcal{L}u$
$\displaystyle=f(\textbf{x},u)-V(\mathbf{x},u,\nabla u),$
$\displaystyle\Gamma(u)$ $\displaystyle=g,\text{ \ }\mathbf{x}\in\partial D,$

where $f(\textbf{x},u)=\lambda u,g=0$ for the case of an eigenvalue problem
with an eigenvalue $\lambda$ and eigenfunction $u(x)$, and the boundary
operator $\Gamma$ could be Dirichlet, Neumann, or Robin type or a decay
condition will be given at $\infty$ if $D=R^{d}$. The following shorthand will
be used in the rest of the paper

(22) $v(\textbf{x})=V(\textbf{x},u(\textbf{x}),\nabla u(\textbf{x})).$

The vector $\mu=\mu(\mathbf{x})$, $\sigma_{d\times d}=\sigma_{d\times
d}(\mathbf{x})$ can be associated with the variable drift and diffusion,
respectively, of the following stochastic Ito process $\mathbf{X}_{t}(\omega)$
$\mathbf{\in}$ $R^{d},$ $\omega\in\Omega$ (random sample space) with
$\mathcal{L}$ as its generator,

(23) $\displaystyle d\mathbf{X}_{t}$
$\displaystyle=\mu(\mathbf{X}_{t})dt\mathbf{+}\sigma(\mathbf{X}_{t})\mathbf{\cdot}d\mathbf{B}_{t}$
$\displaystyle\mathbf{X}_{t}$ $\displaystyle=\mathbf{x}_{0}\in D,$

where $\mathbf{B}_{t}=(B_{t}^{1},\cdots,B_{t}^{d})^{\top}\mathbf{\in}R^{d}$ is
the Brownian motion in $R^{d}.$

The transition probability $P(\mathbf{y},t;\mathbf{x},s)$ for the process
$\mathbf{X}_{t}$ will satisfy the following Fokker-Planck equation

(24) $\frac{\partial P(\mathbf{y},t;\mathbf{x},s)}{\partial
t}=\mathcal{L}^{*}_{y}P(\mathbf{y},t;\mathbf{x},s),$

where the adjoint operator

(25)
$\mathcal{L}^{*}_{y}=-\nabla^{\top}_{y}\mu+\frac{1}{2}Tr(\nabla_{y}\nabla_{y}^{\top}A).$

(Robin Problem) Let us consider the BVP of (20) with a Robin type boundary
condition

(26) $\Gamma(u)=\gamma^{\top}\cdot\nabla u+cu=g,$

where the vector

(27) $\gamma(x)=\frac{1}{2}A\cdot\textbf{n},$

and n is the outward normal at $\textbf{x}\in\partial D$. The Dirichlet BC
$u=f$ can be considered as a limiting case of $c\rightarrow N,g\rightarrow
Ng,N\rightarrow\infty$ and the Neumann BC as the case $c=0.$

With the Martingale problem approach [23], the Martingale problem for the BVP
with the third kind boundary condition (26) can formulated using a reflecting
diffusion process $\mathbf{X}^{ref}$ based on the process $\mathbf{X}$ (23)
through the following Skorohod problem.

(Skorohod problem): Assume $D$ is a bounded domain in $R^{d}$ with a $C^{2}$
boundary. Let $\textbf{X}(t)$ be a (continuous) path of (23) in $R^{d}$ with
$\textbf{X}(0)\in\bar{D}$. A pair $(\textbf{X}^{ref}(t),L(t))$ is a solution
to the Skorohod problem $S(\textbf{X};D)$ if the following conditions are
satisfied:

  1. 1.

$\textbf{X}^{ref}$ is a path in $\bar{D}$;

  2. 2.

(local time) $L(t)$ is a non-decreasing function which increases only when
$X^{ref}\in\partial D$, namely,

(28) $L(t)=\int_{0}^{t}I_{\partial D}(\textbf{X}^{ref}(s))L(ds),$

  3. 3.

The Skorohod equation holds:

(29) $S(\textbf{X};D):\qquad\
\textbf{X}^{ref}(t)=X(t)-\int_{0}^{t}\gamma(\textbf{X}^{ref}(s))L(ds).$

Here, $L(t)$ is the local time of the reflecting diffusion process, where an
oblique reflection with respect to the direction $\gamma$ at the boundary
$\partial D$ occurs once the process $\textbf{X}(t)$ hits the boundary [22].
The sampling of reflecting process and the computation of local time $L(t)$
can be found in [5].

As we will only use the reflecting diffusion $\textbf{X}^{ref}(t)$ for the
rest of our discussion, we will keep the same notation

(30) $\textbf{X}(t)\leftarrow\textbf{X}^{ref}(t)$

with the understanding that it now stands for a reflecting diffusion process
within the closed domain $\bar{D}$. Using the Ito formula [17] for the semi-
martingale $\textbf{X}(t)$ (namely, $\textbf{X}^{ref}(t)$) [10] [18],

$du(\textbf{X}(t))={\displaystyle\sum\limits_{i=1}^{d}}\frac{\partial
u}{\partial
x_{i}}(\textbf{X}(t))dX_{i}(t)+\frac{1}{2}{\displaystyle\sum\limits_{i=1}^{d}}{\displaystyle\sum\limits_{j=1}^{d}}a_{ij}(\textbf{X}(t))\frac{\partial^{2}u}{\partial
x_{i}\partial x_{j}}(\textbf{X}(t))dt,$

with the notation of the generator $\mathcal{L}$, we have for the solution
$u(\textbf{x})$ of (20) the following differential

$\displaystyle du(\textbf{X}(t))=$
$\displaystyle\mathcal{L}u(\textbf{X}(t))dt-\gamma^{\intercal}\cdot\nabla
u(u(\textbf{X}_{t}))L(dt)+{\displaystyle\sum\limits_{i=1}^{d}}{\displaystyle\sum\limits_{j=1}^{d}}\sigma_{ij}\frac{\partial
u}{\partial x_{i}}(\textbf{X}(t))dB_{i}(t)$
$\displaystyle=\left[f(\textbf{X}(t),u(\textbf{X}(t))-V(\textbf{X}(t),u(\textbf{X}(t)),\nabla
u(\textbf{X}(t)))\right]dt-\left[g(\textbf{X}(t))-cu(\textbf{X}(t))\right]L(dt)$
$\displaystyle+{\displaystyle\sum\limits_{i=1}^{d}}{\displaystyle\sum\limits_{j=1}^{d}}\sigma_{ij}\frac{\partial
u}{\partial x_{i}}(\textbf{X}(t))dB_{i}(t),$

which in turns gives a Martingale $M_{t}^{u}$ defined by

(31) $\displaystyle M_{t}^{u}\doteq$ (32) $\displaystyle
u(\textbf{X}_{t})-u(\textbf{X}_{0})-\int_{0}^{t}\left[f(\textbf{X}_{s},u(\textbf{X}_{s}))-V(\textbf{X}_{s},u(\textbf{X}_{s}),\nabla
u(\textbf{X}_{s}))\right]ds$
$\displaystyle+\int_{0}^{t}\left[g(\textbf{X}_{s})-cu(\textbf{X}_{s})\right]L(ds)=\int_{0}^{t}{\displaystyle\sum\limits_{i=1}^{d}}{\displaystyle\sum\limits_{j=1}^{d}}\sigma_{ij}\frac{\partial
u}{\partial x_{i}}(\textbf{X}_{s})dB_{i}(s),$

due to the Martingale nature of the Ito integral at the end of the equation
above.

(Dirichlet Problem) For Dirichlet problem of (20) with a boundary condition

(33) $\Gamma[u]=u=g,\quad\textbf{x}\in\partial D,$

the underlying diffusion process is the original diffusion process (20), but
killed at the boundary at the first exit time

(34) $\tau_{D}=\inf\\{t,\textbf{X}_{t}\in\partial D\\},$

and it can be shown that in fact

(35) $\tau_{D}=\inf\\{t>0,L(t)>0\\},$

and also that $M_{t\wedge\tau_{D}}^{u}$ will be still a Martingale [17], which
will not involve the integral with respect to local time $L(t)$, i.e.

(36)
$M_{t\wedge\tau_{D}}^{u}=u(\textbf{X}_{t\wedge\tau_{D}})-u(\textbf{X}_{0})-\int_{0}^{t\wedge\tau_{D}}\left[f(\textbf{X}_{s},u(\textbf{X}_{s}))-V(\textbf{X}_{s},u(\textbf{X}_{s}),\nabla
u(\textbf{X}_{s}))\right]ds.$

For the case of a linear PDE, i.e. $f(\textbf{x},u)=f(\textbf{x}),V=0$ , by
taking expectation of (36) and letting $t\rightarrow\infty$, we will have

$\displaystyle 0=$ $\displaystyle
E[M^{u}_{0}]=lim_{t\rightarrow\infty}E[M^{u}_{t\wedge\tau_{D}}]=E[M^{u}_{\tau_{D}}]$
$\displaystyle=$ $\displaystyle
E[u(\textbf{X}_{\tau_{D}})-u(\textbf{X}_{0})]-\int_{0}^{\tau_{D}}f(\textbf{X}_{s})ds$
(37) $\displaystyle=$ $\displaystyle
E[g(\textbf{X}_{\tau_{D}})]-u(\textbf{x})-E[\int_{0}^{\tau_{D}}f(\textbf{X}_{s})ds],$

resulting in the well-known Feynman-Kac formuula for the Dirichlet boundary
value problem

(38)
$u(\textbf{x})=E[g(\textbf{X}_{\tau_{D}})]-E[\int_{0}^{\tau_{D}}f(\textbf{X}_{s})ds],\quad\textbf{x}\in
D,$

where the diffusion process $\textbf{X}_{t},\textbf{X}_{0}=\textbf{x}$ is
defined by (23).

The Martingale problem of the BVPs states the equivalence of $M_{t}^{u}$ being
a Martingale (i.e. a probabilistic weak form of the BVPs) and the classic weak
form: For every test function $\phi(\mathbf{x})\in C_{\partial
D}^{2}=\\{\phi:\phi\in C^{2}(D)\cap
C^{1}(\overline{D}),\left(\gamma^{\prime}\cdot\nabla-\mu^{\intercal}\textbf{n}\right)\phi=0\\}$,
$\gamma^{\prime}=\gamma-\alpha$,
$\alpha_{j}={\displaystyle\sum\limits_{i=1}^{3}}\frac{\partial
a_{ij}}{\partial x_{i}}$, we have

(39) $\displaystyle\int_{D}u(\textbf{x})\mathcal{L}^{\ast}\phi d\textbf{x}=$
$\displaystyle\int_{D}\left[f(\textbf{x},u(\textbf{x}))-V(\textbf{x},u(\textbf{x}),\nabla
u(\textbf{x}))\right]\phi(\textbf{x})d\textbf{x}$ $\displaystyle+$
$\displaystyle\int_{\partial
D}\phi(\textbf{x})[g(\textbf{x})-cu(\textbf{x})]ds_{x},$

where

(40) $\mathcal{L}^{\ast}\phi=\frac{1}{2}Tr(\nabla\nabla^{\intercal}A)\phi-
div(\mu\phi).$

This equivalence has been proven for the Schrodinger operator
$\mathcal{L}u=\frac{1}{2}\Delta u+qu$ for the Neumann problem [10] and the
Robin problem [18].

##  4 DeepMartNet - a Martingale based neural network

In this section, we will propose a DNN method for solving the BVPs and
eigenvalue problems for elliptic PDEs using the equivalence between its
Martingale problem formulation and classic weak form of the PDEs.

For simplicity of our discussion, let us assume that $s\leq t\leq\tau_{D}$, by
the Martingale property of $M_{t}\equiv M^{u}_{t}$ of (36), we have

(41) $E[M_{t}|\mathcal{F}_{s}]=M_{s},$

which implies for any measurable set $A\in$ $\mathcal{F}_{s},$

(42) $E[M_{t}|A]=M_{s}=E[M_{s}|A],$

thus,

(43) $E[\left(M_{t}-M_{s}\right)|A]=0,$

i,e,

(44) $\int_{A}\left(M_{t}-M_{s}\right)P(d\omega)=0,$

where

$\displaystyle M_{t}-M_{s}=$ $\displaystyle
u(\mathbf{X}_{t})-u(\mathbf{X}_{s})-\int_{s}^{t}\mathcal{L}u(\mathbf{X}_{z})dz$
(45) $\displaystyle=$ $\displaystyle
u(\mathbf{X}_{t})-u(\mathbf{X}_{s})-\int_{s}^{t}(f(z,u(\mathbf{X}_{z}))-v(\mathbf{X}_{z}))dz.$

In particular, if we take $A=\Omega\in\mathcal{F}_{s}$ in (43), we have

(46) $E[M_{t}-M_{s}]=0,$

i.e. the Martingale $M_{t}$ has a constant expectation. However, it should be
noted that constant expectation by itself does not mean a Martingale, for this
we have the following lemma [3].

######  Lemma 4.1.

If $E[M_{S}]=E[M_{T}]$ holds for any two stopping time $S\leq T$, then
$M_{t},t\geq 0$ is a Martingale, i.e. $E[M_{t}|\mathcal{F}_{s}]=M_{s}$ for
$s\leq t.$

For a given time interval $[0,T]$, we define a partition

(47) $0=t_{0}<t_{1}<\cdots<t_{i}<t_{i+1}<\cdots<t_{N}=T,$

and the increment of the $M_{t}$ over $[t_{i},t_{i+k}]$ can be approximated by
using a trapezoidal rule for the integral term

$\displaystyle M_{t_{i+k}}-M_{t_{i}}=$ $\displaystyle
u(\mathbf{X}_{i+k})-u(\mathbf{X}_{i})-\int_{t_{i}}^{t_{i+k}}\mathcal{L}u(\mathbf{X}_{z})dz$
$\displaystyle\doteq$ $\displaystyle
u(\mathbf{X}_{i+k})-u(\mathbf{X}_{i})-\Delta
t\sum_{l=0}^{k}\omega_{l}\mathcal{L}u(\mathbf{X}_{i+l})$ (48) $\displaystyle=$
$\displaystyle u(\mathbf{X}_{i+k})-u(\mathbf{X}_{i})-\Delta
t\sum_{l=0}^{k}\omega_{l}(f(\mathbf{X}_{i+l},u(\mathbf{X}_{i+l}))-v(\mathbf{X}_{i+l})),$

where for $k\geq 1$ $\omega_{0}=\omega_{k}=\frac{1}{2},\quad\omega_{l}=1,\quad
2\leq l\leq k-1$ and for $k=0$, $\omega_{0}=1$.

Adding back the exit time $\tau_{D}$, we note that

$M_{t_{i+k}\wedge\tau_{D}}-M_{t_{i}\wedge\tau_{D}}=u(\mathbf{X}_{t_{i+k}\wedge\tau_{D}})-u(\mathbf{X}_{t_{i}\wedge\tau_{D}})-\int_{t_{i}\wedge\tau_{D}}^{t_{i+k}\wedge\tau_{D}}\mathcal{L}u(\mathbf{X}_{z})dz=0$

if both $t_{i+k},t_{i}\geq\tau_{D}$.

######  Remark 4.2.

We could define a different generator $\mathcal{L}$ by not including
${\mu}^{\top}\nabla$ in (2), then the Martingale in (49) will be changed to

(49) $M_{t}^{\ast}=u(\mathbf{X}_{t})-u(\mathbf{x}_{0})-\int_{0}^{t}(\lambda
u(\mathbf{X}_{s})-\mu^{\top}(\mathbf{X}_{s})\nabla
u(\mathbf{X}_{s})-v(\mathbf{X}_{s}))ds,$

where the process $\mathbf{X}_{t}$ is given simply by
$d\mathbf{X}_{t}=\sigma\mathbf{\cdot}d\mathbf{B}_{t},$ instead.

  * •

DeepMartNet for Dirichlet BVPs

Let $u_{\theta}(\mathbf{x})$ be a neural network approximating the BVP
solution with $\theta$ denoting all the weight and bias parameters of a DNN.
For a given time interval $[0,T]$, we define a partition

$0=t_{0}<t_{1}<\cdots<t_{i}<t_{i+1}<\cdots<t_{N}=T,$

and $M$-discrete realizations

(50) $\Omega^{\prime}=\\{\omega_{m}\\}_{m=1}^{M}\subset\Omega$

of the Ito process using the Euler-Maruyama scheme with $M$-realizations of
the Brownian motions $\mathbf{B}_{i}^{(m)}$, $0\leq m\leq M$,

$\mathbf{X}_{i}^{(m)}(\omega_{m})\sim X(t_{i},\omega_{m}),0\leq i\leq N,$

where

(51) $\displaystyle\mathbf{X}_{i+1}^{(m)}$
$\displaystyle=\mathbf{X}_{i}^{(m)}+\mu\mathbf{(X}_{i}^{(m)}\mathbf{)}\Delta
t_{i}\mathbf{+}\sigma\mathbf{\mathbf{(X}_{i}^{(m)}\mathbf{)}\cdot}\Delta\mathbf{B}_{i}^{(m)},$
(52) $\displaystyle\mathbf{X}_{0}^{(m)}$ $\displaystyle=\mathbf{x}_{0}$

where $\Delta t_{i}=t_{i+1}-t_{i},$

$\Delta\mathbf{B}_{i}^{(m)}=\mathbf{B}_{i+1}^{(m)}-\mathbf{B}_{i}^{(m)}.$

We will build the loss function $Loss(\theta)$ for neural network
$u_{\theta}(\mathbf{x})$ approximation of the BVP solution using the
Martingale property (44) and the M-realization of the Ito diffusion (23).

For each $t_{i}$, we randomly take a subset of $A_{i}\subset\Omega^{\prime}$
with a uniform sampling (without replacement), corresponding to the mini-batch
in computing the stochastic gradient for the empirical training loss. Assuming
that the mini-batch in each $A_{i}$ is large enough such that
$\\{\mathbf{X}_{i+1}^{(m)}\\}$ and $\\{\mathbf{X}_{i}^{(m)}\\}$,
$\omega_{m}\in A_{i}$ sample the distribution
$P(t_{i+1},.,t_{i},A),P(t_{i},.,0,x_{0})$ well, respectively, then equation
(62) gives the following approximate identity for the solution
$u_{\theta}(\mathbf{X}_{t})$ using the $A_{i}-$ensemble average,

As $E[M_{t_{i+k}}^{u_{\theta}}-M_{t_{i}}^{u_{\theta}}]\approx 0$, for a
randomly selected $A_{i}\in\Omega^{\prime}=\mathcal{F}_{t_{i}}$ ( mini-
batches), we can set

$\displaystyle Loss_{mart}(\theta)$
$\displaystyle=\frac{1}{N}\sum_{i=0}^{N-1}(M_{t_{i+k}\wedge\tau_{D}}^{u_{\theta}}-M_{t_{i}\wedge\tau_{D}}^{u_{\theta}})^{2}$
$\displaystyle=\frac{1}{N}\sum_{i=0}^{N-1}\frac{I({t_{i}\leq\tau_{D}})}{|A_{i}|^{2}}\left(\sum_{m=1}^{|A_{i}|}u_{\theta}(\mathbf{X}_{i+k}^{(m)})-u_{\theta}(\mathbf{X}_{i}^{(m)})-\right.$
(53) $\displaystyle\left.\Delta
t\sum_{l=0}^{k}\omega_{l}(f(\mathbf{X}_{i+l}^{(m)},u_{\theta}(\mathbf{X}_{i+l}^{(m)}))-v_{\theta}(\mathbf{X}_{i+l}^{(m)}))\right)^{2},$

where $v_{\theta}(\textbf{x})$ is defined similarly as in (22). Refer to
Remark 4.4 for the discussion of the size of the mini-batch $|A_{i}|$.

Now, we define the total loss for the boundary value problem as

(54) $Loss_{total-
bvp}(\theta)=Loss_{mart}(\theta)+\alpha_{bdry}Loss_{bdry}(\theta),$

and $\alpha_{bdry}$ is a hyper-parameter and
$Loss_{bdry}(\theta)=||u_{\theta}(\textbf{x})-g||_{2}^{2}$, which can be
approximated by evaluations at the boundary.

DeepMartNet solution\- $u_{\theta^{*}}(x)$, here

(55) $\theta^{*}=argminLoss_{total-bvp}(\theta).$

  * •

DeepMartNet for Dirichlet eigenvalue problems

For the eigenvalue problem

(56) $\displaystyle\mathcal{L}u+V(\mathbf{x},u,\nabla u)$
$\displaystyle=\lambda u,\text{ \ }\mathbf{x}\in D\subset R^{d},$
$\displaystyle\Gamma(u)=u$ $\displaystyle=0,\text{ \ }\mathbf{x}\in\partial
D,$

and the Martingale loss becomes

$\displaystyle Loss_{mart}(\lambda,\theta)=\frac{1}{N}\sum_{i=0}^{N-1}$
$\displaystyle\frac{1}{|A_{i}|^{2}}\left(\sum_{m=1}^{|A_{i}|}u_{\theta}(\mathbf{X}_{i+k}^{(m)})-u_{\theta}(\mathbf{X}_{i}^{(m)})-\right.$
(57) $\displaystyle\left.\Delta t\sum_{l=0}^{k}\omega_{l}(\lambda
u_{\theta}(\mathbf{X}_{i+l}^{(m)})-v_{\theta}(\mathbf{X}_{i+l}^{(m)}))\right)^{2}$

and in the case of a bounded domain, the boundary loss $Loss_{bdry}(\theta)$
will be added for the homogeneous boundary condition $g=0$, i.e.,
$Loss_{bdry}(\theta)=||u_{\theta}(\textbf{x})||_{2}^{2}$. For the decay
condition at infinite, a molifier will be used to enforce explicitly the decay
condition there (see (85) ).

Also, in order to prevent the DNN eigenfunction going to a zero solution, we
introduce a simple normalization term using $l_{p}$ (p=1, 2) norm of the
solution at some randomly selected location

(58)
${Loss}_{normal}(\theta)=\left(\frac{1}{m}\sum_{i=1}^{m}|u_{\theta}({\textbf{x}_{i}})|^{p}-c\right)^{2},$

where $\textbf{x}_{i}$ are m arbitrarily selected fixed points and $c$ is a
nonzero constant.

Finally, we have the total loss for the eigenvalue problem as

(59) $Loss_{total-
eig}(\lambda,\theta)=Loss_{mart}(\lambda,\theta)+\alpha_{bdry}Loss_{bdry}(\theta)+\alpha_{normal}{Loss}_{normal}(\theta),$

where $\alpha_{bdry}$ and $\alpha_{normal}$ are hyper-parameters.

DeepMartNet eigen-problem solution \- $(\lambda^{*},u_{\theta^{*}}(x))$, here

(60) $(\lambda^{*},\theta^{*})=argminLoss_{total-eig}(\lambda,\theta).$

######  Remark 4.3.

(Mini-batch in SGD training and Martingale property) Due to the equivalence
between (44) and (41), the loss function defined above ensures that
$M_{t}^{u_{\theta}}$ of (36) for $u_{\theta}(\mathbf{x})$ will be a Martingale
approximately if the mini-batch $A_{i}$ explores all subsets of the sample
space $\Omega^{\prime}$ during the SGD optimization process of the training,
and the sample size $M=|\Omega^{\prime}|$ $\rightarrow$ $\infty,$ the time
step $\max|\Delta t_{i}|\rightarrow 0,$ and the training converges (see Fig.
1).

Also, if we take $A_{i}=\Omega^{\prime}$ for all $i,$ there will be no
stochasticity in the gradient calculation for the loss function, we will have
a traditional full gradient descent method and the full Martingale property
for $u_{\theta}(\mathbf{x})$ is not enforced either. Therefore, the mini-batch
procedure in DNN SGD optimization corresponds perfectly with the Martingale
definition (41).

In summary, the Martingale property implies that for any measurable set $A\in$
$\mathcal{F}_{s},$ we require that

(61) $E[M_{t}|A]=M_{s},$

which then provides a native mechanism for the mini-batch in the SGD.
Therefore, the Martingale based DNN is an ideal fit for deep learning of high-
dimensional PDEs.

Figure 1: DeepMartNet training and Martingale property

######  Remark 4.4.

(Size of mini-batch $A_{i}$) The loss function of the DeepMartNet is based on
the fact that $\int_{A_{i}}\left(M^{u}_{t}-M^{u}_{s}\right)P(d\omega)=0$ for
the exact solution $u(\textbf{x})$, where the expectation will be computed by
an ensemble average of selected paths from the $M$-paths. In theory, using the
transitional probability, using left end point quadrature for the integral
over [s,t] in (45) with $|t-s|<<1$, (44) can be rewritten as

$\displaystyle 0=E[\left(M^{u}_{t}-M^{u}_{s}\right)|A_{i}]$
$\displaystyle\doteq$ $\displaystyle\int_{\textbf{z}\in B}\int_{\textbf{y}\in
R^{d}}\big{[}u(\textbf{y})-u(\textbf{z})-(f(\textbf{z},u(\textbf{z}))-v(\textbf{z}))(t-s)\big{]}P(t,\textbf{y},s,\textbf{z})p(s,\textbf{z},0,\textbf{x}_{0})dzdy$
$\displaystyle=$ $\displaystyle\int_{\textbf{y}\in
R^{d}}u(\textbf{y})P(t,\textbf{y},s,A_{i})dy-\int_{\textbf{z}\in
B}\big{[}u(\textbf{z})+(f(\textbf{z},u(\textbf{z}))-v(\textbf{z}))(t-s)\big{]}P(s,\textbf{z},0,\textbf{x}_{0})dz$
(62) $\displaystyle\sim$
$\displaystyle\frac{1}{|A_{i}|}\sum_{m=1}^{|A_{i}|}\left(u(\textbf{X}_{t}^{(m)})-\big{[}u(\textbf{X}_{s}^{(m)})+(f(\textbf{X}_{s}^{(m)},u(\textbf{X}_{s}^{(m)}))-v(\textbf{X}_{s}^{(m)}))(t-s)\big{]}\right),$

where $B=\mathbf{X}^{-1}_{s}(A_{i})$ and $\textbf{X}_{t}^{(m)},1\leq m\leq M$
are the M-sample paths of the diffusion process $\textbf{X}_{t}$.

Therefore, the size of mini-batch $|A_{i}|$ should be large enough to give an
accurate sampling of the continuous distribution $P(t,\textbf{y},s,A),y\in
R^{d}$ and $P(s,\textbf{z},0,\textbf{x}_{0}),z\in B$, so in our simulation we
select a sufficient large $M$ and set the size of mini-batch $A_{i}$ in the
following range,

(63) $M/m_{1}\leq|A_{i}|\leq M/m_{2},$

where $M$ is the total number of paths used and $m_{1}>m_{2}$ are hyper-
parameters of the training.

The set $A_{i}$ can be the same for each $0\leq i\leq N-1$ per epoch or can be
randomly selected as the subset of $\Omega^{\prime}$ depending how much
stochasticity is put into the calculation of the stochastic gradient in the
SGD optimizations.

For a low memory implementation of the DeepMartNet, at any time we can just
generate enough $|A_{i}|$ number of paths to be used for some epochs of
training, and then regenerate them again without first generating a large
number of paths upfront.

##  5 Numerical Results

The numerical parameters for the DeepMartNet consist of

  * •

M - number of total paths

  * •

T - terminal time of the paths

  * •

N - number of time steps over [0, t]

  * •

$\Delta t=T/N$

  * •

$M_{b}=|A_{i}|$ size of mini-batches of paths selected according to (63).

  * •

Size of networks

The training is carried out on a Nvidia Superpod one GPU node A100. The
Optmizer is Adamax [16].

###  5.1 Dirichet BVPs of the Poisson-Boltzmann equation

We will first apply the DeepMartNet to solve the Dirichlet BVP of the Poisson-
Boltzmann equation (PBE) arising from solvation of biomolecules in inoic
solvents [1].

(64) $\begin{cases}\Delta
u(\textbf{x})+cu(\textbf{x})=f(\textbf{x}),&\textbf{x}\in D\\\
u(\textbf{x})=g(\textbf{x}),&\textbf{x}\in\partial D\end{cases}$

where $c<0$ ($c=-1$ in the numerical tests) with an high-dimensional solution
given by

(65) $u(\textbf{x})=\sum_{i=1}^{d}\cos(\omega x_{i}),\quad\omega=2.$

In this case, the generator for the stochastic process is
$\mathcal{L}=\frac{1}{2}\Delta$, so the corresponding diffusion is simply the
Brownian motion $\textbf{B}(t)$. For the $M-$ Brownian paths
$\textbf{B}^{(j)},j=1,\cdots,M$ originating from $x_{0}$, the Martingale loss
(53) becomes

$\displaystyle Loss_{\text{mart}}(\theta):=$
$\displaystyle\frac{1}{N}\sum_{i=0}^{N-1}\frac{1}{|A_{i}|^{2}}\left(\sum_{j=1}^{|A_{i}|}u_{\theta}(\textbf{B}_{t_{i+1}}^{(j)})-u_{\theta}(\textbf{B}_{t_{i}}^{(j)})\right.$
(66) $\displaystyle-$
$\displaystyle\left.\frac{1}{2}\left(f(\textbf{B}_{t_{i}}^{(j)})-cu(\textbf{B}_{t_{i}}^{(j)})\right)\mathbb{I}(t_{i}\leq\tau_{D}^{(j)})\Delta
t\Big{)}\right)^{2}.$

Meanwhile, we could use the Feynman-Kac formula [17] to compute the solution
at the point $\textbf{x}_{0}$ by

(67)
$u(\textbf{x}_{0})\approx\frac{1}{M}\sum_{j=1}^{M}\left(g\left(\textbf{B}_{\tau_{D}}^{(j)}\right)\mathrm{e}^{\frac{c\tau_{D}}{2}}+\frac{1}{2}\sum_{i=0}^{N-1}f\left(\textbf{B}_{t_{i}}^{(j)}\right)\mathrm{e}^{\frac{ct_{i}}{2}}\mathbb{I}\left(t_{i}\leq\tau_{D}^{(j)}\right)\Delta
t\right),$

and also define a point solution loss (or an integral identity for the
solution for PDE with quasi-linear term $V$ and $f$ in (20) ), termed Feynman-
Kac loss, which is added to the total loss (54)

(68)
$Loss_{\text{F-K}}(\theta)=(u_{\theta}(\textbf{x}_{0})-u(\textbf{x}_{0}))^{2},$

and a boundary loss is approximated as

(69)
$Loss_{bdry}(\theta)=||u_{\theta}-g||^{2}_{2}\approx\frac{1}{N_{bdry}}\sum^{N_{bdry}}_{k=1}|u_{\theta}(\textbf{x}_{k})-g(\textbf{x}_{k})|^{2},$

where uniformly sampled boundary points $\textbf{x}_{k},1\leq k\leq N_{bd}$
are used to compute the boundary integral. Therefore, the total loss for the
boundary value problem of the PB equation is

(70)
$Loss_{bvp}(\theta)=Loss_{\text{mart}}(\theta)+\alpha_{\text{F-K}}Loss_{\text{F-K}}(\theta)+\alpha_{bdry}Loss_{bdry}(\theta),$

where the penalty parameter $\alpha_{\text{F-K}}$ ranges from 10 to 1000 and
$\alpha_{bdry}$ ranges from 1000 to 10,000. An Adamax optimizer with learning
rate 0.05 is applied for training and
$\alpha_{\text{F-K}}=10,\alpha_{\text{bdry}}=10^{3}$ are taken for the
following numerical tests.

  * •

Test 1: PBE in a d=20 dimensional cube $[-1,1]^{d}$. In this test, we solve
the PBE in a 20 dimensional cube with the Dirichlet boundary condition given
by the exact solution (65). The total number of paths is $M=100,000$ starting
from origin, and mini-batch size $M_{b}=|A_{i}|=1000;\Delta t=0.01$ and T=9. A
fully connected network with layers (20, 64, 32, 1) with first hidden layer
Tanh and second hidden layer GeLU activation function is used. A mini-batch of
$N_{\text{bdry}}=2000$ points $\textbf{x}_{k}$ are uniformly sampled points on
the boundary for every epoch of training.

Fig. 2 show the learned solution along diagonal of the cube (top left) as well
as along the first coordinate (top right) and the history of the loss (bottom
left) and relative error L2 (computed by Monte Carlo sample) (bottom right).
The training takes about less than 5 minutes.

Figure 2:  DeepMartNet solution of PBE in $D=[-1,1]^{d},d=20$. (Upper left):
true and predicted value of $u$ along the diagonal of the unite cube; (Upper
right): true and predicted value of $u$ along the first coordinate axis of
$\mathbb{R}^{d}$. (Lower left): The loss $L$ history; (lower right): The
history of relative error $L^{2}$ over the cube.

  * •

Test 2. Effect of path starting point $\textbf{x}_{0}$. In previous test, the
DeepMartNet uses all diffusion paths originating from a fixed point
$\textbf{x}_{0}$ to explore the solution domain. In this test, we consider
paths starting from different $\textbf{x}_{0}=(l,0,\cdots,0),l=0.1,0.3,0.7$ to
investigate the effect of different starting point $\textbf{x}_{0}$ on the
accuracy of the DeepMartNet. A fully connected network with layers (20, 64,
32, 1) with first hidden layer Tanh and second hidden layer GeLU activation
function will be used as in Test 1. The total number of paths $M=100,000$, and
for each epoch, we randomly choose $M_{b}=1000$ from the paths; the time step
of the paths is $\Delta t=0.01$.

Fig. 3 shows that the three different choices of $\textbf{x}_{0}$ produce
similar numerical results with the same numerical parameters as in Test 1.

Figure 3:  DeepMartNet for PBE in $D=[-1,1]^{d},d=20$ with 3 different
starting points for the paths. From up to down:
$\textbf{x}_{0}=(l,0,\cdots,0),l=0.1,0.3,0.7.$ (left):$u(x\bm{e}_{1})$ where
$\bm{e}_{1}=d^{-1/2}(1,1\cdots,1)$,(Right): $u(x\bm{e})$ where
$\bm{e}=(1,0\cdots,0)$.

Moreover, the DeepMartNet can use diffusion paths starting from different
initial position $\textbf{x}_{0}$ in training the DNN as long as the mini-
batch of the paths $A_{i}$ for a given epoch corresponds to paths originating
from a common initial point $\textbf{x}_{0}$. In Fig. 4, we compare the
numerical result using 120,000 total paths starting from
$\textbf{x}_{0}=(0.3,0,\cdots,0,0)$ with that with three choices of
$\textbf{x}_{0}$ as in Fig. 3 with 40,000 paths for each $\textbf{x}_{0}$. In
both cases, for each epoch, we randomly choose mini-batch size $M_{b}=1000$
from the paths; the time step of the paths is $\Delta t=0.01$. An Adamax
optimizer with learning rate 0.05 is applied for training.
$\alpha_{\text{F-K}}=10,\alpha_{\text{bdry}}=10^{3}$. The DeepMartNet produce
similar results for these two cases with same other numerical parameters as in
Test 1.

Figure 4: Comparison of DeepMartNet with paths from one starting point and
from 3 starting points with same number of total paths in both cases. Upper
left: $u(x\bm{e})$ where $\bm{e}=(1,0\cdots,0)$; Upper right: the relative
error $|u-u_{\text{true}}|/\|u_{\text{true}}\|_{\infty}$ for the upper left
plot; Lower left: $u(x\bm{e}_{1})$ where $\bm{e}_{1}=d^{-1/2}(1,1\cdots,1)$;
Lower right: the error $\|u_{\text{true}}\|_{\infty}$ for the lower left plot.

  * •

Test 3: PBE in a d=100 dimensional unit ball. In this test, we solve the PBE
in a 100 dimensional space, We set $T=0.25$ and the time step of the paths is
$\Delta t=0.005$, and the total number of paths is $M=100,000$, and for each
epoch, we randomly choose $M_{b}=1000$ from the paths. A fully connected
network with layers (100, 128, 32, 1) with first hidden layer Tanh and second
hidden layer GeLU activation function is used.

Figure 5:  DeepMartNet solution of the PBE in
$D=\\{\textbf{x}\in\mathbb{R}^{d},|\textbf{x}|<1\\}$, d=100. (Upper left):
true and predicted value of $u$ at the diagonal of $\mathbb{R}^{d}$, i.e. $x$
versus $u(x\bm{e})$ where $\bm{e}=d^{-1/2}(1,1\cdots,1)$; (Upper right): true
and predicted value of $u$ at the first coordinate axis of $\mathbb{R}^{d}$,
i.e. $x$ versus $u(x\bm{e_{1}})$ where $\bm{e}=d^{-1/2}(1,0\cdots,0)$. (
Middle left): $u(x\bm{e}^{\prime})$ where
$\bm{e}^{\prime}=2^{-1/2}(1,1,0,0,\cdots,0)$ (Middle right):
$u(x\bm{e}^{\prime\prime})$ where
$\bm{e}^{\prime\prime}=10^{-1/2}(\underset{ten1^{\prime}s}{1,\cdots,1},0,0,\cdots,0)$
(Lower left): The loss $L$ versus the number of epoch; (lower right): The
relative error $L^{2}$ error versus the number of epoch.

Fig. 5 show the learned solution (Upper left): true and predicted value of $u$
at the diagonal of $\mathbb{R}^{d}$, i.e. $x$ versus $u(x\bm{e})$ where
$\bm{e}=d^{-1/2}(1,1\cdots,1)$; (Upper right): true and predicted value of $u$
at the first coordinate axis of $\mathbb{R}^{d}$, i.e. $x$ versus
$u(x\bm{e_{1}})$ where $\bm{e}=d^{-1/2}(1,0\cdots,0)$. ( Middle left):
$u(x\bm{e}^{\prime})$ where $\bm{e}^{\prime}=2^{-1/2}(1,1,0,0,\cdots,0)$
(Middle right): $u(x\bm{e}^{\prime\prime})$ where
$\bm{e}^{\prime\prime}=10^{-1/2}(\underset{ten1^{\prime}s}{1,\cdots,1},0,0,\cdots,0)$
(Lower left): The loss $L$ versus the number of epoch; (lower right): The
relative error $L^{2}$ error versus the number of epoch. The training takes
about less than 30 minutes.

  * •

Test 4: nonlinear PBE in a d=10 dimensional unit ball. For a 1:1 symmetric two
species ionic solvent, the electrostatic potential based on the Debye-Huckel
theory [1] is given by a nonlinear PBE before linearization, and we will
consider the following model nonlinear PBE to test the capability of the
DeepMartNet for solving nonlinear PDEs,

(71) $\begin{cases}-\triangle u+\sinh u=f&x\in D,\\\ u=g&x\in\partial
D,\end{cases}$

where

(72) $D=\\{x\in{\mathbb{R}}^{d}:\|\textbf{x}\|_{2}\leq L\\},L=1.$

We consider the case of a true solution as

(73) $u=\alpha\sum_{i=1}^{d}x_{i}^{2},\quad\alpha=2,$

which gives the boundary data and right hand side of the PBE as

$g\equiv\alpha L^{2},$

and

$f=-2\alpha d+\sinh\left(\alpha\sum_{i=1}^{d}x_{i}^{2}\right).$

This time, we use the loss as

(74)
$Loss(\theta)=Loss_{\text{mart}}(\theta)+\alpha_{\text{bdry}}Loss_{\text{bdry}}(\theta),$

where $\alpha_{\text{bdry}}=10^{-4}$, and the Martingale loss (53) is

(75) $\displaystyle
Loss_{\text{mart}}(\theta)=\frac{1}{N}\sum_{i=0}^{N-1}\frac{\mathbb{I}\left(t_{i}\leq\tau_{\partial
D}^{(j)}\right)}{M_{b}^{2}}$
$\displaystyle\cdot\left(\sum_{j=1}^{M_{b}}u_{\theta}\left(W_{t_{i+1}}^{(j)}\right)-u_{\theta}\left(W_{t_{i}}^{(j)}\right)-\frac{\Delta
t}{2}\left(\sinh
u(W_{t_{i}}^{(j)})-f\left(W_{t_{i}}^{(j)}\right)\right)\right)^{2},$

and the boundary condition loss is again

(76)
$Loss_{\text{bdry}}(\theta)=\frac{1}{N_{\text{bdry}}}\sum_{k=0}^{N_{\text{bdry}}}(u_{\theta}(\textbf{x}_{k})-g(\textbf{x}_{k}))^{2},$

and a mini-batch of $N_{\text{bdry}}=2000$ points $\textbf{x}_{k}$ are
uniformly sampled points on the boundary for every epoch of training. And in
this case, we do not have a Feynman-Kac loss at the starting point
$\textbf{x}_{0}$.

The numerical results in Fig. 6 is done with the DeepMartNet with the fully
connected NN $(10,10,10,1)$. The activation function for the first hidden
layer is Tanh, and the one for the second hidden layer is GELU with Tanh
approximation. $M_{\text{tot}}=10^{6}$ paths are used and the size of mini-
batch $M_{b}=4000$ paths are sampled in each epoch. The total number of epoch
of the training is $3000$. The terminal time for the paths is $T=0.25$ and
time step is $t=0.01$. The learning rate starts at 0.01, and decreased by a
factor of 0.99 every 100 epochs. The training takes less than 20 minutes.

Figure 6: DeepMartNet solution of nonlinear PBE (71) in a 10-dimensional unit
ball. Upper left: $u(x\bm{e_{1}})$ where $\bm{e_{1}}=(1,0,\cdots,0)$; Upper
right: $u(x\bm{e})$ where $\bm{e}=(1,1,\cdots,1)$; Lower left: loss vs. epoch;
right: relative $L^{2}$ error vs epoch.

###  5.2 Eigenvalue problem for elliptic equations

In this section, we will apply the DeepMartNet to solve elliptic eigenvalue
problems in bounded and unbounded domains for both self-adjoint and non-self
adjoint operators.

####  5.2.1 Eigenvalue problem of the Laplace equation in a cube in $R^{10}$

First, we consider a self-adjoint eigenvalue problem

(77) $\displaystyle-\triangle u$ $\displaystyle=\lambda u,\quad\textbf{x}\in
D=[-L,L]^{d},\quad d=10,\quad L=1$ (78) $\displaystyle u|_{\partial D}$
$\displaystyle=0,$

with the following eigenfunction

(79) $u(x)=sin(\frac{\pi x_{1}}{L})\cdots sin(\frac{\pi x_{d}}{L}),$

for the lowest eigenvalue

(80) $\lambda_{1}=d(\frac{\pi}{L})^{2}.$

A DeepMartNet with a fully connected structure $(10,20,10,1)$ and GELU
activation function are used. The total number of paths is $M=100,000$, and
for each epoch, we randomly choose $M_{b}=1000$ from the paths; the time step
of the paths is $\Delta t=0.05$. A mini-batch of $N_{\text{bdry}}=2000$ points
$\textbf{x}_{k}$ are uniformly sampled points on the boundary for every epoch
of training. An Adamax optimizer with learning rate 0.05 is applied for
training. $\alpha_{\text{bdry}}=10^{3}$ in (59). The constant
$\alpha_{normal}=10,p=1,m=1,\textbf{x}_{1}=0,c=1$ in (58).

To accelerate the convergence of the eigenvalue, We included an extra term
into the loss function (59)

(81) $\alpha_{\text{eig}}|\lambda^{(i)}-\lambda^{(i-100)}|^{2},$

where $i$ is the index of the current epoch, which corresponds to the residual
of the equation $\frac{d\lambda}{dt}=0$ as the training time
$t\rightarrow\infty$, and ensures that the eigenvalue $\lambda$ will converge
and the penalty constant $\alpha_{\text{eig}}$ is chosen as $2.5\times
10^{-8}$ in this case. The terminal time for the paths is $T=0.6$ and time
step is $t=0.01$. The learning rate starts at 0.02, and decreased by a factor
of 0.995 every 100 epochs. The training takes less than 20 minutes.

Figure 7: A numerical result for eigenvalue problem of Laplace equation in a
cube $[-L,L]^{10}$ in $R^{10}$. Upper left to middle right: true and predicted
value of $u$ at $x\bm{e}$ where
$\bm{e}=k^{-1/2}(\underset{k1^{\prime}s}{1,\cdots,1},0,0,\cdots,0)$ for
$k=1,2,5,10$ Lower left: The predicted eigenvalue $\lambda$ versus the number
of epoch, (orange) horizontal line shows the true eigenvalue; lower right: The
relative error $L^{2}$ error versus the number of epoch.

####  5.2.2 A non-self adjoint eigenvalue problem for the Fokker-Planck
equations

Here we consider the following a non-self adjoint eigenvalue problem of the
Fokker-Planck equation

(82) $-\Delta\psi-\nabla\cdot\left(\psi\nabla
W\right)+c\psi=-\Delta\psi-\nabla W\cdot\nabla\psi-\Delta
W\psi+c\psi=\lambda\psi,\quad\textbf{x}\in R^{d}$

where the eigenfunction for the eigenvalue $\lambda=c$ with a zero decay
condition at $\infty$ is

(83) $\psi(\textbf{x})=e^{-W(\textbf{x})}.$

Here, we will consider a quadratic potential for our numerical tests

$W(\textbf{x})=||\textbf{x}||^{2},\quad\textbf{x}\in R^{d}.$

Equation (82) can re-written as

(84) $\mathcal{L}\psi=\frac{1}{2}\Delta\psi+\frac{1}{2}\nabla
W\cdot\nabla\psi=-\left(\frac{1}{2}\Delta
W-\frac{1}{2}c+\frac{1}{2}\lambda\right)\psi.$

The generator for the SDE $\mathcal{L}$ will have a drift and a diffusion as

$\mu=\frac{1}{2}\nabla W\quad\text{and}\quad\sigma=I_{d\times d}.$

In order to enforce explicitly the decay condition of the eigenfunction at the
infinite, a mollifier of the following form will be used as a pre-factor for
the DNN solution
$u_{\theta}(\textbf{x})=\rho(\textbf{x})\tilde{u}_{\theta}(\textbf{x})$ to the
eigenfunction,,

(85) $\rho(\textbf{x})=\frac{1}{1+(\frac{||\textbf{x}||}{\alpha})^{2}},$

where $\alpha$ is a constant and $\alpha=\frac{5}{11}$.

The Martingale loss of (57) for this case will be

$\displaystyle{Loss}_{mart}=$ $\displaystyle\frac{1}{\Delta
t}\frac{1}{N}\sum_{i=0}^{N-1}\frac{1}{|A_{i}|^{2}}\left(\sum_{m=1}^{|A_{i}|}u_{\theta}(\textbf{x}_{\textbf{i
+ 1}}^{(m)})-u_{\theta}(\textbf{x}_{\textbf{i}}^{(m)})\right.$ (86)
$\displaystyle+$ $\displaystyle(\frac{1}{2}\Delta
W(\textbf{x}_{\textbf{i}}^{(m)})-\frac{1}{2}c+\frac{1}{2}\lambda)u_{\theta}(\textbf{x}_{\textbf{i}}^{(m)})\Delta
t\Big{)}^{2}.$

For the mini-batches in this case, we take the size of each $A_{i}$ to be
between $\frac{M}{200}$ and $\frac{M}{25}$ as a random assortment of the $M$
total trajectories.

To speed up the convergence of the eigenvalue and eigenfunctions, we found out
that by taking fractional powers of each individual loss term is helpful,
namely, the total loss is now modified as

(87)
$Loss_{total}=\left(\left(Loss_{mart}\right)^{p}+\alpha_{normal}(Loss_{normal})^{q}\right)^{r}.$

In our numerical tests, a typical choice is $p=3/8,q=1,r=3/4$, and for the
normalization $Loss_{normal}$, we set $\alpha_{normal}=50$ and $c=30$,
$m=2,\textbf{x}_{1}=0,\textbf{x}_{2}=(1,0,\cdots,0)$ in (58).

Figs. 8 and 9 show the learned eigenvalues ($\lambda=5$ for $d=5$ and
$\lambda=200$ for $d=200$) and eigenfunctions, respectively. The
$k=3$-trapezoidal rule is used in the Martingale loss ${Loss}_{mart}$ and the
other numerical parameters used are listed as follows

  * •

The total number of paths starting from the origin $M=9,000,24,000$ for
$d=5,200$, respectively.

  * •

The number of time steps $N=1350,1300$ and the terminal time $T=9$ for
$d=5,200$, respectively.

  * •

the learning rate is 1/150 and is halved every 500 epochs starting at epoch
500 and halved twice at epoch 7500 for stabilization.

  * •

A In the following numerical tests, a fully connected network with layers (d,
6d, 3d, 1) with a Tanh activation function is used for the eigenfunction while
a fully connected network (1, d, 1) with a $ReLu^{9}$ activation function with
a constant value input is used to represent the eigenvalue.

  * •

An Adamax optimizer is applied for training.

The final relative error in the eigenvalues after 10,000 epochs of training is
$1.3\times 10^{-2},6.7\times 10^{-3}$ for $d=5,200$, respectively. And the
relative $L^{2}$ error of the eigenfunction calculated along the diagonal of
the domain is $2.6\times 10^{-2},2.9\times 10^{-2}$ for $d=5,200$,
respectively. The training takes 25 minutes for the case of $d=200$.

Figure 8:  Eigen-value problem of Fokker-Planck equation in $R^{d}$ $d=5$ for
eigenvalue $\lambda=5$. (Top left) Convergence of eigenvalue, (top right)
history of eigenvalue error, (bottom left) history of loss function, (bottom
right) Learn and exact eigenfunction along the diagonal of the domain.

Figure 9:  Eigenvalue problem of Fokker-Planck equation in $R^{d}$ $d=200$ for
eigenvalue $\lambda=200$. (Top left) Convergence of eigenvalue, (top right)
history of eigenvalue error, (bottom left) history of loss function, (bottom
right) Learn and exact eigenfunction along the diagonal of the domain.

##  6 Conclusion

In this paper, we introduced a Martingale based neural network, DeepMartNet,
for solving the BVPs and eigenvalue problems of elliptic operators with
Dirichlet boundary conditions. The DeepMartNet enforces the Martingale
property for the PDE solutions through the stochastic gradient descent (SGD)
optimization of the Martingale property loss functions. The connection between
the Martingale definition and the mini-batches used in stochastic gradient
computation shows a natural fit between the SGD optimization in DNN training
and the Martingale problem formulation of the PDE solutions. Numerical results
in high dimensions for both BVPs and eigenvalue problems show the promising
potential in approximating high dimensional PDE solutions. The numerical
results show that the DeepMartNet extracts more information for the PDEs
solution over the whole solution domain than the traditional one point
solution Feynman-Kac formula from the same set of diffusion paths originating
from just one point.

Future work on DeepMartNet will include PDEs with Neumann and Robin boundary
conditions where the underlying diffusion will be a reflecting one and the
local time of the reflecting diffusion will be computed and included in the
definition of the Martingale loss function. Another area of application is in
optimal stochastic control where the Martingale optimality condition for the
control and the backward SDE for the value function can be used to construct
the loss function for the control as well as the value function [2]. Another
important area for the application of the DeepMartNet is to solve low-
dimension PDEs in complex geometries as the method uses diffusion paths to
explore the domain and therefore can handle highly complex geometries such as
the interconnects in microchip designs and nano-particles in material sciences
and molecules in biology. Finally, The convergence analysis of the
DeepMartNet, especially to a given eigenvalue, and the choices of mini-batches
of diffusion paths and various hyper-parameters, which can affect the
convergence of the learning strongly, in the loss functions are important
issues to be addressed.

## Acknowledgement

W. C. would like to thank Elton Hsu and V. Papanicolaou for the helpful
discussion about their work on probabilistic solutions of Neumann and Robin
problems.

## References

  * [1] Cai W. Computational Methods for Electromagnetic Phenomena: electrostatics in solvation, scattering, and electron transport. Cambridge University Press; 2013 Jan 3.
  * [2] Cai W. DeepMartNet – A Martingale based Deep Neural Network Learning Algorithm for Eigenvalue/BVP Problems and Optimal Stochastic Controls. arXiv preprint arXiv:2307.11942. 2023 Jul 21.
  * [3] Cohen SN, Elliott RJ. Stochastic calculus and applications. New York: Birkhäuser; 2015 Nov 18.
  * [4] Davis MH. Martingale methods in stochastic control. InStochastic Control Theory and Stochastic Differential Systems: Proceedings of a Workshop of the Sonderforschungsbereich 72 der Deutschen Forschungsgemeinschaft an der Universität Bonn “which took place in January 1979 at Bad Honnef 2005 Oct 6 (pp. 85-117). Berlin, Heidelberg: Springer Berlin Heidelberg.
  * [5] Ding CY, Zhou, YJ, Cai W, Zeng X, and Yan CH. A Path Integral Monte Carlo (PIMC) Method based on Feynman-Kac Formula for Electrical Impedance Tomography, Journal of Computational Physics., 476 (2023) 111862. 121.
  * [6] E W, Han J, Jentzen A. Algorithms for solving high dimensional PDEs: from nonlinear Monte Carlo to machine learning. Nonlinearity. 2021 Dec 9;35(1):278.
  * [7] Han J, Jentzen A, E W. Solving high-dimensional partial differential equations using deep learning. Proceedings of the National Academy of Sciences. 2018 Aug 21;115(34):8505-10.
  * [8] Han J, E W. Deep learning approximation for stochastic control problems, Deep Reinforcement Learning Workshop, NIPS. arXiv preprint arXiv:1611.07422. 2016.
  * [9] Han J, Lu J, Zhou M. Solving high-dimensional eigenvalue problems using deep neural networks: A diffusion Monte Carlo like approach. Journal of Computational Physics. 2020 Dec 15;423:109792.
  * [10] Hsu, P. Reflecting Brownian Motion, Boundary Local Time and the Neumann Problem, Dissertation Abstracts International Part B: Science and Engineering [DISS. ABST. INT. PT. B- SCI. ENG.], 45(6), 1984.
  * [11] Hsu P. Probabilistic approach to the Neumann problem. Communications on pure and applied mathematics. 1985 Jul;38(4):445-72.
  * [12] Hu Z, Shukla K, Karniadakis GE, Kawaguchi K. Tackling the curse of dimensionality with physics-informed neural networks. arXiv preprint arXiv:2307.12306. 2023 Jul 23.
  * [13] Hutzenthaler M, Jentzen A, Kruse, T, Nguyen TA and von Wurstemberger, P. Overcoming the curse of dimensionality in the numerical approximation of semilinear parabolic partial differential equations. Proceedings of the Royal Society A, 476(2244):20190630, 2020.
  * [14] Ji S, Peng, S, Peng Y, Zhang X. Three algorithms for solving high-dimensional fully coupled FBSDEsthrough deep learning, IEEE Intell. Syst. 35(3) (Feb 2020) 71–84.
  * [15] Karatzas I, Shreve S. Brownian motion and stochastic calculus. Springer Science & Business Media; 2012 Dec 6.
  * [16] Kingma DP, Ba J. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980. 2014 Dec 22.
  * [17] Klebaner FC. Introduction to stochastic calculus with applications. World Scientific Publishing Company; 2012 Mar 21.
  * [18] Papanicolaou VG. The probabilistic solution of the third boundary value problem for second order elliptic equations, Probab. Theory Relat. Fields 87 (1990) 27-77.
  * [19] Pardoux E, Peng S. Backward stochastic differential equations and quasilinear parabolic partial differential equations. In Stochastic partial differential equations and their applications 1992 (pp. 200-217). Springer, Berlin, Heidelberg
  * [20] Pfau D, Spencer JS, Matthews AG, Foulkes WM. Ab initio solution of the many-electron SchrÃ¶dinger equation with deep neural networks. Physical Review Research. 2020 Sep 16;2(3):033429.
  * [21] Raissi M. Forward-backward stochastic neural networks: Deep learning of high-dimensional partial differential equations. arXiv preprint arXiv:1804.07010. 2018 Apr 19.
  * [22] Schuss Z. Brownian dynamics at boundaries and interfaces. Springer-Verlag New York; 2015.
  * [23] Stroock DW, Varadhan SRS.Diffusion processes with boundary conditions. Commun. Pure Appl. Math. 24, 147-225 (1971).]
  * [24] Zhang W, Cai W. FBSDE based neural network algorithms for high-dimensional quasilinear parabolic PDEs. Journal of Computational Physics. 2022 Dec 1;470:111557.