#  Solving the inverse potential problem in the parabolic equation by the deep
neural networks method

Mengmeng Zhang<EMAIL_ADDRESS>School of Science, Hebei University of
Technology, Tianjin 300401, China
Nanjing Center for Applied Mathematics
Nanjing, 211135, China
Zhidong Zhang<EMAIL_ADDRESS>School of Mathematics (Zhuhai),
Sun Yat-sen University, Zhuhai 519082, Guangdong, China

###### Abstract

In this work, we consider an inverse potential problem in the parabolic
equation, where the unknown potential is a space-dependent function and the
used measurement is the final time data. The unknown potential in this inverse
problem is parameterized by deep neural networks (DNNs) for the reconstruction
scheme. First, the uniqueness of the inverse problem is proved under some
regularities assumption on the input sources. Then we propose a new loss
function with regularization terms depending on the derivatives of the
residuals for partial differential equations (PDEs) and the measurements.
These extra terms effectively induce higher regularity in solutions so that
the ill-posedness of the inverse problem can be handled. Moreover, we
establish the corresponding generalization error estimates rigorously. Our
proofs exploit the conditional stability of the classical linear inverse
source problems, and the mollification on the noisy measurement data which is
set to reduce the perturbation errors. Finally, the numerical algorithm and
some numerical results are provided.

AMS subject classifications: 34K28, 35R30, 65N15, 62M45.

Keywords: inverse potential problem, deep neural networks, uniqueness,
generalization error estimates, numerical reconstruction.

##  1 Introduction.

###  1.1 Mathematical model.

The following parabolic system is considered in this work:

$\begin{cases}\begin{aligned}
(\partial_{t}-\Delta+q(x))u&=F(x,t),&&(x,t)\in\Omega_{T},\\\
u(x,t)&=b(x,t),&&(x,t)\in\partial\Omega_{T},\\\
u(x,0)&=u_{0}(x),&&x\in\Omega.\end{aligned}\end{cases}$ (1.1)

Here we write $\Omega_{T}=\Omega\times(0,T]$ and
$\partial\Omega_{T}=\partial\Omega\times(0,T]$ for short, and
$\Omega\subset\mathbb{R}^{d}$ is an open bounded domain in $\mathbb{R}^{d}$
with sufficiently smooth boundary. $F(x,t),\ u_{0}(x),\ b(x,t)$ are the source
term, initial status, boundary condition respectively, causing the heat
propagation in the medium. The potential function $q(x)\in
L^{\infty}(\Omega)$, called the heat radiative coefficient of the material, is
a crucial parameter for characterizing the heat conduction process. It
describes the ability of the medium to propagate heat from internal sources or
sinks. For known $(F(x,t),u_{0}(x),b(x,t),q(x))$ with suitable regularities,
the forward problem (1.1) is well-posed in appropriate function space [11]. In
this work, we consider the inverse problem of recovering the unknown $q(x)$,
where the used measurement is the final time data

$\displaystyle u(x,T):=\varphi(x),\quad x\in\Omega.$ (1.2)

In practical applications of inverse problems, the contamination on inverse
problems is unavoidable. So we will be given the noisy data $\varphi^{\delta}$
instead of the exact data $\varphi(x)$ in (1.2), which satisfies

$\displaystyle\|\varphi^{\delta}-\varphi\|_{L^{\infty}(\Omega)}\leq\delta.$
(1.3)

To handle the effect caused by the perturbations, people need to develop
effective methods to improve the accuracy and robustness in applications. In
this study, we choose the deep neural networks (DNNs) to solve the inverse
problem (1.1)-(1.3). Comparing to traditional methods for solving inverse
potential problem, this approach demonstrates the superiority in high-
dimensional space and has the advantage of breaking the curse of
dimensionality.

There are rare works on studying the inverse potential problem for parabolic
equations using deep neural networks, especially the rigorous analysis of its
convergence estimate. In this work, the authors will consider the solution of
the inverse potential problem (1.1)-(1.3) parameterized by DNNs for the
reconstruction scheme. We propose a new loss function with regularization
terms depending on the derivatives of the residuals for PDEs and measurements.
The mollification method has been employed to improve the regularity of the
noisy data. Also, the generalization error estimates are rigorously derived
from the conditional stability of the linear inverse source problem and the
mollification error estimate on noisy data.

###  1.2 Literature.

The reconstructions of $q(x)$ in (1.1) from some inversion input data have
been studied extensively. For zero initial status, the uniqueness for $q(x)$
by (1.1)-(1.2) is established in [24], while the unique reconstruction using
final measurement data is studied in [25]. In the case of non-zero initial
status, the existence and uniqueness of the generalized solution
$(u(x,t),q(x))\in W_{p}^{2,1}(\Omega_{T})\times L^{p}(\Omega)$ with the time-
average temperature measurement are given in [18] for $(u_{0},\varphi)$ with
some regularities. Choulli and Yamamoto [6] prove the generic well-posedness
of the inverse problem in Hölder spaces by final measurement data, and then
the conditional stability result in a Hilbert space setting for sufficiently
small $T$ is studied in [7]. Chen et al [3] consider the inverse potential
problem from a partial measurements over $[T_{0},T_{1}]\times\Omega$ with
$[T_{0},T_{1}]\subset[0,T]$, where the conditional stability estimates of the
inverse problem in some Sobolev space and the reasonable convergence rates of
the Tikhonov regularization are derived. Recently, Jin et al [16] uses the
same observational data and shows a weighted $L^{2}$ stability in the standard
$L^{2}$ norm under a positivity condition. They provide an error analysis of
reconstruction scheme based on the standard output least-squares formulation
with Tikhonov regularization (by an $H^{1}$-seminorm penalty). Zhang et al
[39] prove the uniqueness of the identification from final time data for
(sub)diffusion equation and show the conditional stability in Hilbert spaces
under some suitable conditions on the problem data. The convergence and error
analysis of the reconstruction discrete scheme are rigorously analyzed. The
investigations in the inverse non-smooth potential problem are given in [38],
where the uniqueness for this nonlinear inverse problem is proved.
Numerically, an iterative process called two-point gradient method is proposed
by minimizing the data-fit term and the penalty term alternatively, with a
convergence analysis in terms of the tangential condition. There also exists
some works involving multiple coefficient identification. For example,
Yamamoto and Zou [35] investigate the simultaneous reconstruction of the
initial temperature and heat radiative coefficient in a heat conductive
system, with stability of the inverse problem and the reconstruction scheme.
Kaltenbacher and Rundell [17] consider the inverse problem of simultaneously
recovering two unknowns, spatially dependent conductivity and the potential
function from overposed data consisting of $u(x,T)$. The uniqueness result and
the convergence of an iteration scheme are established. We also refer to [2,
4, 8, 14, 28, 34] and the references therein for the inverse potential
problems in diffusion models from different types of observational data.

Recently, deep learning methods for solving PDEs have been realized as an
effective approach, especially in high dimensional PDEs. Such methods have the
advantage of breaking the curse of dimensionality. The basic idea is to use
neural networks (nonlinear functions) to approximate the unknown solutions of
PDEs by learning the parameters. For the forward problems, there exists many
numerical works with deep neural networks involving the depth Ritz method
(DRM) [10], the depth Galerkin method (DGM) [33], the DeepXDE method [20],
depth operator network method (DeepONet) [19], physical information neural
networks (PINNs) [26], the weak adversary neural network (WAN) [1, 36] and so
on. Theoretically, there are some rigorous analysis works investigating the
convergence and error estimates for the solution of PDEs via neural networks,
but the result are still far from complete. For example, the convergence rate
of DRM with two layer networks and deep networks are studied in [9, 13, 21,
22]; the convergence of PINNs is given in [15, 23, 29, 30, 31]. For the
inverse problems, the PINNs frameworks can be employed to solve the so-called
data assimilation or unique continuation problems, and rigorous estimates on
the generalization error of PINNs are established in [23]. Bao et al [36]
develop the WAN to solve electrical impedance tomography (EIT) problem. In
[37], the authors study a classical linear inverse source problem using the
final time data under the frameworks of neural networks, where a rigorous
generalization error estimate is proposed with a novel loss function including
the Sobolev norm of some residuals. For more specific inverse problems applied
in engineering and science, we refer to [5, 27, 32].

###  1.3 Outline.

The rest of this article is organized as follows. In Section 2 we introduce
the knowledge of neural networks and the setting of mollification. In Section
3, we introduce a conditional stability of the linear inverse source problem
first. Then the uniqueness theorem (Theorem 3.1) of this inverse potential
problem can be proved followed from the conditional stability. In Section 4, a
novel loss function with specific regularization terms is introduced. Then we
prove the generalization error estimates of data-driven solution of inverse
problems, which is stated in Theorem 4.1. In Section 5, we propose the
reconstruction algorithm and provide several experiments to show the validity
of the proposed algorithm.

##  2 Preliminaries.

###  2.1 Neural network architecture.

First we introduce the basic knowledge of neural network briefly. Note that
$u_{\theta}$ and $q_{\eta}$ are two separate networks with different variables
$(x,t)$ and $x$. Thus, we use $\xi$ to denote collectively the network
parameters for a parametric function $s_{\xi}(z)$ such that a general scheme
can be applied for either $u_{\theta}(x,t)$ (with $z=(x,t),\ \xi=\theta$) or
$q_{\eta}(x)$ (with $z=x,\ \xi=\eta$). For a positive integer
$K\in\mathbb{N}$, a $K$-layer feed-forward neural network of $s_{\xi}(z)$ for
$z\in\mathbb{R}^{d_{0}}$ is a function $s_{\xi}(z)$ defined by

$\displaystyle s_{\xi}(z):=W_{K}l_{K-1}\circ\cdots\circ l_{1}(z)+b_{K},$ (2.1)

where the $k$-th layer
$l_{k}:\mathbb{R}^{d_{k-1}}\rightarrow\mathbb{R}^{d_{k}}$ is given by
$l_{k}(z)=\sigma\left(W_{k}z+b_{k}\right)$ with weights
$W_{k}\in\mathbb{R}^{d_{k}\times d_{k-1}}$ and biases
$b_{k}\in\mathbb{R}^{d_{k}}$ for $k=2,\cdots,K$. The activation function
$\sigma(\cdot)$ includes sigmoid, tanh, ReLU (Rectified Linear Unit), softmax
and so on [12]. These activation functions introduce non-linearities and
enable the network to learn complex patterns and relationships in the data.
The neural network (2.1) consists of an input layer with argument $z$, where
$d_{0}=d$ is the problem dimension (also known as the size of input layer), an
output layer which has the weights $W_{K}\in\mathbb{R}^{d_{K}\times d_{K-1}}$
and biases $b_{K}\in\mathbb{R}^{d_{K}}$, and $K-1$ hidden layers for some
$K\in\mathbb{N}$. The network parameters of all layers are collectively
denoted by

$\displaystyle\xi:=\left(W_{K},b_{K},W_{K-1},b_{K-1},\cdots,W_{1},b_{1}\right).$

In Figure 1, we give a simple architectures of fully connected neural
networks, where $z=(x_{1},x_{2},\cdots,x_{d})$ is d-dimensional input
variables, and the neural networks function is given as $s_{\xi}(z)=y_{NN}$.

Figure 1: The fully connected neural networks.

###  2.2 Mollification.

In the practical applications of inverse problems, the noise of the
measurements is unavoidable. The noisy data will make the residuals
uncontrollable, which can be seen in the next section. Hence, we choose to
mollify the measured data beforehand. The next is the introduction of
mollification.

Fix one function $\rho\in C^{2}(\mathbb{R})$ as

$\text{supp}\rho=(0,1),\;\rho(0)=\rho(1)=\rho^{\prime}(0)=\rho^{\prime}(1)=0,$

and

$\displaystyle\int_{0}^{\infty}\rho(t)t^{d-1}\ dt=\frac{1}{\pi_{d}},$

with $\pi_{d}$ is the surface area of unit sphere $B(0,1)$ in $R^{d}$. Set
$\rho_{\epsilon}(x):=\epsilon^{-d}\rho(x/\epsilon)$, and define the mollifier
as

$G_{\epsilon}\psi=\int_{|x-y|\leq\epsilon}\rho_{\epsilon}(|x-y|)\psi(y)\ dy.$
(2.2)

Then we have

$\int_{\mathbb{R}^{d}}\rho_{\epsilon}(|x-y|)\ dy=1.$

In the next lemma, we concern with the estimate of $\Delta\varphi-\Delta
G_{\epsilon}(\varphi^{\delta})$.

######  Lemma 2.1.

Assume that the noisy data $\varphi^{\delta}\in L^{\infty}(\Omega)$ and the
exact data $u(x,T):=\varphi(x)\in H^{2}(\Omega)$ satisfy

$\displaystyle\|\varphi-\varphi^{\delta}\|_{L^{\infty}(\Omega)}\leq\delta.$

Also, the exact data imposes the high-order Lipschitz continuous condition.
More precisely, we can find a positive constant $C_{\varphi}$ such that

$\displaystyle|\varphi(x)-\varphi(y)|$ $\displaystyle\leq C_{\varphi}|y-x|,$
$\displaystyle|\Delta\varphi(x)-\Delta\varphi(y)|$ $\displaystyle\leq
C_{\varphi}|y-x|,$

for $x,y\in\overline{\Omega}$ uniformly. For the mollification operator (2.2),
if we pick $\epsilon=O(\delta^{1/3})$, then we can achieve the following
optimal error bound

$\displaystyle\|\Delta\varphi-\Delta
G_{\epsilon}(\varphi^{\delta})\|_{L^{\infty}(\Omega)}\leq C\delta^{1/3}.$

###### Proof.

We split the subtraction $\Delta\varphi-\Delta G_{\epsilon}(\varphi^{\delta})$
as following:

$\Delta\varphi-\Delta G_{\epsilon}(\varphi^{\delta})=(\Delta\varphi-
G_{\epsilon}(\Delta\varphi))+(G_{\epsilon}(\Delta\varphi)-\Delta
G_{\epsilon}(\varphi^{\delta}))=:I_{1}+I_{2}.$

For $I_{1}$, we have that

$|I_{1}|\leq\int_{|x-y|\leq\epsilon}\rho_{\epsilon}(|x-y|)~{}|\Delta\varphi(x)-\Delta\varphi(y)|\
dy\leq C\epsilon.$

For $I_{2}$, Green’s identities and the properties of the kernel function
$\rho$ give that $\Delta G_{\epsilon}(\varphi)=G_{\epsilon}(\Delta\varphi)$.
Hence,

$\displaystyle I_{2}$
$\displaystyle=\Delta\Big{[}\int_{R^{d}}\rho_{\epsilon}(|x-y|)~{}(\varphi(y)-\varphi^{\delta}(y))\
dy\Big{]}$
$\displaystyle=\int_{R^{d}}\Delta\rho_{\epsilon}(|x-y|)~{}(\varphi(y)-\varphi^{\delta}(y))\
dy.$

From the straightforward calculation, we can deduce that

$\Delta\rho_{\epsilon}(|x-y|)=\epsilon^{-d-2}\rho^{\prime\prime}(|x-y|/\epsilon)+(d-1)\epsilon^{-d-1}|x-y|^{-1}\rho^{\prime}(|x-y|/\epsilon),$

which gives

$|I_{2}|\leq\delta\int_{|x-y|\leq\epsilon}|\Delta\rho_{\epsilon}(|x-y|)|\
dy\leq C\delta\epsilon^{-2}.$

So we have

$|\Delta\varphi-\Delta G_{\epsilon}(\varphi^{\delta})|\leq
C\epsilon(1+\delta\epsilon^{-3}).$

By picking $\epsilon=O(\delta^{1/3})$, we can achieve the desired estimate and
complete the proof. ∎

##  3 Uniqueness.

The uniqueness of this inverse potential problem is one of our main results.
In this section, we will prove the uniqueness and the proof relies on the
conditional stability of the inverse source problem of equation (1.1). The
conditional stability will be stated in the next subsection.

###  3.1 Conditional stability of the inverse source problem.

Under the framework of DNNs, the total error of the reconstructed solution
depends on the training error and the measurement error. This connection
relies on the conditional stability of linear inverse source problem, i.e.,
the quantitative dependence of the unknown source on the measurement data.
Sequentially, here we will introduce some known results for the linear inverse
source problem.

The mathematical statement of inverse source problem in parabolic equations
with final time data is given below. For the parabolic equation

$\begin{cases}\begin{aligned}
(\partial_{t}-\Delta+\overline{q}(x))v(x,t)&=p(x)h(x,t),&&(x,t)\in\Omega_{T},\\\
v(x,t)&=0,&&(x,t)\in\partial\Omega_{T},\\\
v(x,0)&=0,&&x\in\Omega,\end{aligned}\end{cases}$ (3.1)

we set $\overline{q}\geq 0$ and $\overline{q}\in L^{\infty}(\Omega)$, and
$h(x,t)$ is given. Then the inverse source problem is to use the measurement

$\displaystyle\overline{\varphi}(x):=v[p](x,T)$ (3.2)

to recover the unknown $p(x)$ in the source term.

Recalling the norm of the classical Sobolev space $W^{2,1}_{2}(\Omega_{T})$ as

$\displaystyle\|u\|_{W^{2,1}_{2}(\Omega_{T})}=\sqrt{\sum\limits_{|\alpha|\leq
2}\|D^{\alpha}u\|^{2}_{L^{2}(\Omega_{T})}+\|u_{t}\|^{2}_{L^{2}(\Omega_{T})}},$

the following classical result on the inverse source problem (3.1)-(3.2) can
be found in [24].

######  Lemma 3.1.

For equation (3.1), we assume that

$h\in L^{\infty}(\Omega_{T}),\ h_{t}\in L^{\infty}(\Omega_{T}),\ p(x)\in
L^{2}(\Omega),\ ph\in L^{2}(\Omega_{T}),\ ph_{t}\in L^{2}(\Omega_{T}),$

and

$h(x,t)\geq 0,\ h_{t}(x,t)\geq 0\ \text{on}\ \Omega_{T},\ |h(x,T)|\geq\nu>0\
\text{on}\ \Omega.$

Here $\nu$ is a fixed positive number. Then, for known $\overline{q}\in
L^{\infty}(\Omega)$ and input data $\overline{\varphi}\in H^{2}(\Omega)$,
there exists a unique solution $(v(x,t),p(x))\in W^{2,1}_{2}(\Omega_{T})\times
L^{2}(\Omega)$ to (3.1)-(3.2), following the estimate

$\displaystyle\|p\|_{L^{2}(\Omega)}+\|v\|_{W^{2,1}_{2}(\Omega_{T})}\leq
C\|(-\Delta+\overline{q})\overline{\varphi}\|_{L^{2}(\Omega)}.$

The constant $C$ depends on
$\left\|\overline{q}\right\|_{L^{\infty}(\Omega)}$, $\nu,$ $\Omega$ and $T$.

###  3.2 Uniqueness theorem.

Now it is time to show the uniqueness theorem. First we introduce the
admissible set for the unknown potential $q(x)$ as

$\mathcal{A}:=\\{\psi\in L^{\infty}(\Omega):0\leq\psi(x)\leq M\text{ a.e. on
}\Omega\\}\subset L^{2}(\Omega).$

The constant $M$ is the given upper bound of the admissible set. Next,
recalling equation (1.1), we collect some restrictions on the controllable
source $F(x,t)$, initial status $u_{0}(x)$ and boundary condition $b(x,t)$.

######  Assumption 1.

The assumptions on $F(x,t)$, $u_{0}(x)$ and $b(x,t)$ are given as follows.

  * •

$u_{0}(x)\in H^{2}(\Omega)$, $u_{0}(x)=b(x,0)$ on $\partial\Omega$,
$\exists\nu>0$ such that $u_{0}(x)\geq\nu>0$ on $\Omega$;

  * •

$b\in H^{2}(\partial\Omega)$, $b\geq\nu>0$ on $\partial\Omega$, $b_{t}\geq 0$
on $\partial\Omega$;

  * •

$F\in L^{2}(\Omega_{T})$, $F_{t}\in L^{2}(\Omega_{T})$, $F\geq 0$ on
$\Omega_{T}$, $F_{t}\geq 0$ on $\Omega_{T}$;

  * •

$\Delta u_{0}(x)-Mu_{0}(x)+F(x,0)\geq 0$ on $\Omega$.

######  Theorem 3.1.

Under Assumption 1, the inverse problem (1.1)-(1.3) has at most one solution
in $W^{2,1}_{2}(\Omega_{T})\times\mathcal{A}$.

###### Proof.

Assume that there are two distinct pairs $(u[q_{1}],q_{1})$ and
$(u[q_{2}],q_{2})$ satisfying (1.1)-(1.3) with same data

$u[q_{1}](x,T)=u[q_{2}](x,T)=\varphi(x).$

Setting

$w(x,t):=u[q_{1}](x,t)-u[q_{2}](x,t),\quad\overline{q}(x):=q_{2}(x)-q_{1}(x),$

then $w(x,t)$ meets the system

$\begin{cases}\begin{aligned}
(\partial_{t}-\Delta+q_{1}(x))w(x,t)&=\overline{q}(x)u[q_{2}](x,t),&&(x,t)\in\Omega_{T},\\\
w(x,t)&=0,&&(x,t)\in\partial\Omega_{T},\\\
w(x,0)&=0,&&x\in\Omega,\end{aligned}\end{cases}$ (3.3)

with

$w(x,T)=0,\quad x\in\Omega.$ (3.4)

We need to prove that

$(w(x,t),\overline{q}(x))=(0,0)\ \text{in}\ W^{2,1}_{2}(\Omega_{T})\times
L^{\infty}(\Omega).$

Obviously $\overline{q}(x)\in L^{2}(\Omega)$. Also there holds $u[q_{2}]\in
L^{\infty}(\Omega_{T})$ and $u_{t}[q_{2}]\in L^{\infty}(\Omega_{T})$ by [38,
Lemma 2.1]. Then we have

$\overline{q}u[q_{2}]\in L^{2}(\Omega_{T}),\quad\overline{q}u_{t}[q_{2}]\in
L^{2}(\Omega_{T}).$

Under Assumption 1 and the maximum principle, we can see that $u[q_{2}]\geq 0$
on $\Omega_{T}$. For $u_{t}[q_{2}]$, with Assumption 1 and equation (1.1), it
satisfies

$\begin{cases}\begin{aligned}
(\partial_{t}-\Delta+q_{2}(x))(u_{t}[q_{2}])&=F_{t}(x,t)\geq
0,&&(x,t)\in\Omega_{T},\\\ u_{t}[q_{2}](x,t)&=b_{t}(x,t)\geq
0,&&(x,t)\in\partial\Omega_{T},\\\ u_{t}[q_{2}](x,0)&=\Delta
u_{0}(x)-q_{2}u_{0}(x)+F(x,0)\geq 0,&&x\in\Omega.\end{aligned}\end{cases}$

Then the maximum principle leads to $u_{t}[q_{2}]\geq 0$ straightforwardly.
With the positivity of $u_{t}[q_{2}]$, we derive that

$u[q_{2}](x,t)=u_{0}(x)+\int_{0}^{t}\partial_{s}u[q_{2}](x,s)ds\geq
u[q_{2}](x,0)\geq\nu>0,\ (x,t)\in\Omega_{T},$

which yields $u[q_{2}](x,T)\geq\nu>0$. Now the conditions of Lemma 3.1 are
satisfied, and we conclude $(w(x,t),\overline{q}(x))=(0,0)$ by applying Lemma
3.1 on (3.3)-(3.4). The proof is complete. ∎

##  4 Generalization error estimates.

In this section, we will discuss the error estimate of our approach for the
inverse potential problem. Firstly, we introduce the corresponding residuals
and define the loss function.

###  4.1 Loss function and training errors.

We propose a formulation of loss function for data-driven solutions of inverse
problems, which can ensure the accuracy with the conditional stability of the
given linear inverse source problem. To achieve it, we define suitable
residuals that measure the errors of the governed system and the input data.

Assume that the activation function is of $C^{2}$ regularity for the neural
network $u_{\theta}$ defined by (2.1), which leads to $u_{\theta}\in
H^{2}(\overline{\Omega}\times[0,T])$. For the network parameters

$\theta\in\Theta:=\\{(W_{k},b_{k})\\}_{k=1}^{K}:W_{k}\in\mathbb{R}^{d_{k}\times
d_{k-1}},b_{k}\in\mathbb{R}^{d_{k}}\\},$

the set of all possible trainable parameters $u_{\theta}(x,t)$ up to its
second order weak derivatives are bounded in $\overline{\Omega}\times[0,T]$
for any specific $\theta$. Similarly, noticing that $q_{\eta}(x)$ is the
parametric neural network to approximate the potential function $q(x)$, we
assume the activation function for the neural network $q_{\eta}(x)$ is of
$L^{\infty}$ regularity such that $q_{\eta}(x)\in L^{\infty}(\Omega)$. We
define

  * •

Interior PDE residual

$\displaystyle\mathcal{R}_{{int},\theta,\eta}(x,t):=\partial_{t}u_{\theta}(x,t)-\Delta
u_{\theta}(x,t)+q_{\eta}(x)u_{\theta}(x,t)-F(x,t),\quad(x,t)\in\Omega_{T}.$

  * •

Spatial boundary residual

$\displaystyle\mathcal{R}_{sb,\theta}(x,t):=u_{\theta}(x,t)-b(x,t),\quad(x,t)\in\partial\Omega_{T}.$

  * •

Initial status residual

$\displaystyle\mathcal{R}_{tb,\theta}(x):=u_{\theta}(x,0)-u_{0}(x),\quad
x\in\Omega.$

  * •

Data residual

$\displaystyle\mathcal{R}_{d,\theta}(x):=u_{\theta}(x,T)-G_{\epsilon}\varphi^{\delta}(x),\quad
x\in\Omega.$ (4.1)

Note that in the data residual (4.1), we use the mollified data
$G_{\epsilon}\varphi^{\delta}(x)$ instead of the noisy data
$\varphi^{\delta}(x)$. A loss function minimization scheme for data-driven
inverse problems seeks to minimize these residuals comprehensively with some
weights balancing different residuals. The loss function is defined as
follows:

$\displaystyle J_{\lambda}(\theta,\eta)=$
$\displaystyle\|q_{\eta}\mathcal{R}_{d,\theta}\|^{2}_{L^{2}(\Omega)}+\|\Delta\mathcal{R}_{d,\theta}\|^{2}_{L^{2}(\Omega)}+\lambda\|\mathcal{R}_{{int},\theta,\eta}\|^{2}_{H^{1}(0,T;L^{2}(\Omega))}$
(4.2)
$\displaystyle+\|\mathcal{R}_{{tb},\theta}\|^{2}_{L^{2}(\Omega)}+\|q_{\eta}\mathcal{R}_{{tb},\theta}\|^{2}_{L^{2}(\Omega)}+\|\Delta\mathcal{R}_{{tb},\theta}\|^{2}_{L^{2}(\Omega)}+\|\mathcal{R}_{{sb},\theta}\|^{2}_{H^{2}(0,T;L^{2}(\partial\Omega))},$

where $\lambda$ is a hyper-parameter to balance the residuals between the
knowledge of PDE and the measurements. The proposed loss function (4.2)
includes derivative penalties on the residuals. This is motivated by the
conditional stability result for linear inverse source problem, which requires
higher regularity on the measurement data $u(\cdot,T)$ (see Lemma 3.1). To
improve the regularity of the noisy measurement data, we employ the
mollification method by applying the mollification operator $G_{\varepsilon}$
on the noisy data $\varphi^{\delta}$. The design of the loss function for
inverse problems distinguishes itself from that for forward problems such as
physics-informed neural networks. The smoothness requirements not only ensure
the existence of forward problem solutions, but also ensure the well-posedness
of the inverse problem within the optimization framework.

######  Remark 1.

The following standard loss function

$\displaystyle
J^{s}(\theta,\eta)=\|\mathcal{R}_{d,\theta}\|_{L^{2}(\Omega)}^{2}+\lambda\|\mathcal{R}_{{int},\theta,\eta}\|_{L^{2}(\Omega_{T})}^{2}+\|\mathcal{R}_{{tb},\theta}\|_{L^{2}(\Omega)}^{2}+\|\mathcal{R}_{{sb},\theta}\|_{L^{2}(\partial\Omega_{T})}^{2}$
(4.3)

has often been used in the literature. For example, the DGM workflow adopts
this form of loss function and minimizes it by least squares scheme [33].

To determine $(\theta,\eta)$ from the discrete training set, accurate
numerical evaluation of the integrals in (4.2) is essential. We introduce the
following training sets that facilitate efficient computation of the
integrals, leading to better performance:

$\displaystyle\mathcal{S}_{d}$
$\displaystyle:=\left\\{(x_{n},T):x_{n}\in\Omega,\quad
n=1,2,\cdots,N_{d}\right\\},$ $\displaystyle\mathcal{S}_{int}$
$\displaystyle:=\left\\{(\widetilde{x}_{n},\widetilde{t}_{n}):(\widetilde{x}_{n},\widetilde{t}_{n})\in\Omega_{T},\quad
n=1,2,\cdots,N_{int}\right\\},$ $\displaystyle\mathcal{S}_{tb}$
$\displaystyle:=\left\\{(\overline{x}_{n},0):\overline{x}_{n}\in\Omega,\quad
n=1,2,\cdots,N_{tb}\right\\},$ $\displaystyle\mathcal{S}_{sb}$
$\displaystyle:=\left\\{(\widehat{x}_{n},\widehat{t}_{n}):(\widehat{x}_{n},\widehat{t}_{n})\in\partial\Omega_{T},\quad
n=1,2,\cdots,N_{sb}\right\\}.$

Applying these sets and the numerical quadrature rules [23], we get the
following empirical loss function

$\displaystyle J_{\lambda}^{N}(\theta,\eta)$
$\displaystyle=\sum_{n=1}^{N_{d}}\omega_{n}^{d,0}|q_{\eta}(x_{n})\mathcal{R}_{d,\theta}(x_{n})|^{2}+\sum_{n=1}^{N_{d}}\omega_{n}^{d,1}|\Delta\mathcal{R}_{d,\theta}(x_{n})|^{2}$
(4.4)
$\displaystyle\quad+\lambda\sum_{n=1}^{N_{int}}\omega_{n}^{int,0}|\mathcal{R}_{{int},\theta,\eta}(\widetilde{x}_{n},\widetilde{t}_{n})|^{2}+\lambda\sum_{n=1}^{N_{int}}\omega_{n}^{int,1}|\partial_{t}\mathcal{R}_{{int},\theta,\eta}(\widetilde{x}_{n},\widetilde{t}_{n})|^{2}$
$\displaystyle\quad+\sum_{n=1}^{N_{tb}}\omega_{n}^{tb,0}|\mathcal{R}_{tb,\theta}(\overline{x}_{n})|^{2}+\sum_{n=1}^{N_{tb}}\omega_{n}^{tb,1}|q_{\eta}(\overline{x}_{n})\mathcal{R}_{tb,\theta}(\overline{x}_{n})|^{2}+\sum_{n=1}^{N_{tb}}\omega_{n}^{tb,2}|\Delta\mathcal{R}_{tb,\theta}(\overline{x}_{n})|^{2}$
$\displaystyle\quad+\sum_{n=1}^{N_{sb}}\omega_{n}^{sb,0}|\mathcal{R}_{sb,\theta}(\widehat{x}_{n},\widehat{t}_{n})|^{2}+\sum_{n=1}^{N_{sb}}\omega_{n}^{sb,1}|\partial_{t}\mathcal{R}_{sb,\theta}(\widehat{x}_{n},\widehat{t}_{n})|^{2}\sum_{n=1}^{N_{sb}}\omega_{n}^{sb,2}|\partial_{t}^{2}\mathcal{R}_{sb,\theta}(\widehat{x}_{n},\widehat{t}_{n})|^{2},\quad\quad\quad$

where the coefficients

$\omega^{d,k}_{n},\ \omega^{int,k}_{n},\ \omega^{tb,j}_{n},\
\omega^{sb,j}_{n},\ k=0,1,\ j=0,1,2$

are the quadrature weights. It is easy to see that the error for the loss
function is

$\displaystyle|J_{\lambda}(\theta,\eta)-J_{\lambda}^{N}(\theta,\eta)|\leq
C\min\\{N_{d}^{-\alpha_{d,k}},N_{int}^{-\alpha_{int,k}},N_{tb}^{-\alpha_{tb,j}},N_{sb}^{-\alpha_{sb,j}}:k=0,1,\
j=0,1,2\\},\quad$ (4.5)

where $C$ depends on the continuous norm $\|\cdot\|_{C(\Omega)}$ of the
integrals, the rate $\alpha^{d,k}$, $\alpha^{int,k}$, $\alpha^{tb,j}$,
$\alpha^{sb,j}$ ($k=0,1,\ j=0,1,2$) are positive and depend on the regularity
of the underlying integrand i.e, on the space $C(\Omega)$. Therefore, the
underlying solutions and neural networks should be sufficiently regular such
that the residuals can be approximated to a high accuracy by the quadrature
rule.

Now, we define the generalization errors as

$\displaystyle\begin{cases}\mathcal{E}_{G,q}:=\left\|q-q^{*}\right\|_{L^{2}(\Omega)},\\\
\mathcal{E}_{G,u}:=\left\|u-u^{*}\right\|_{C\left([0,T];L^{2}(\Omega)\right)},\end{cases}$
(4.6)

where $u^{*}:=u_{\theta^{*}},\ q^{*}:=q_{\eta^{*}}$ with
$(\theta^{*},\eta^{*})$ is the minimizer of the functional (4.4). Also, we
estimate generalization errors in terms of the following training errors:

  * •

The measurement data training errors:
$\mathcal{E}_{T,d}:=\mathcal{E}_{T,d,0}+\mathcal{E}_{T,d,1}$, where

$\displaystyle\begin{cases}\mathcal{E}_{T,d,0}:=\left(\sum_{n=1}^{N_{d}}\omega_{j}^{d,0}\left|q_{\eta}\left(x_{n}\right)\mathcal{R}_{d,\theta^{*}}\left(x_{n}\right)\right|^{2}\right)^{\frac{1}{2}},\\\
\mathcal{E}_{T,d,1}:=\left(\sum_{n=1}^{N_{d}}\omega_{j}^{d,1}\left|\Delta\mathcal{R}_{d,\theta^{*}}\left(x_{n}\right)\right|^{2}\right)^{\frac{1}{2}}.\end{cases}$
(4.7)

  * •

The interior PDE training errors:
$\mathcal{E}_{T,int}:=\mathcal{E}_{T,int,0}+\mathcal{E}_{T,int,1}$, where

$\displaystyle\begin{cases}\mathcal{E}_{T,int,0}:=\left(\sum_{n=1}^{N_{int}}\omega_{n}^{int,0}|\mathcal{R}_{{int},\theta^{*},\eta^{*}}(\widetilde{x}_{n},\widetilde{t}_{n})|^{2}\right)^{\frac{1}{2}},\\\
\mathcal{E}_{T,int,1}:=\left(\sum_{n=1}^{N_{int}}\omega_{n}^{int,1}|\partial_{t}\mathcal{R}_{{int},\theta^{*},\eta^{*}}(\widetilde{x}_{n},\widetilde{t}_{n})|^{2}\right)^{\frac{1}{2}}.\end{cases}$
(4.8)

  * •

The initial condition training errors:
$\mathcal{E}_{T,tb}:=\mathcal{E}_{T,tb,0}+\mathcal{E}_{T,tb,1}+\mathcal{E}_{T,tb,2}$,
where

$\displaystyle\begin{cases}\mathcal{E}_{T,tb,0}:=\left(\sum_{n=1}^{N_{tb}}\omega_{n}^{tb,0}|\mathcal{R}_{tb,\theta^{*}}(\overline{x}_{n})|^{2}\right)^{\frac{1}{2}},\\\
\mathcal{E}_{T,tb,1}:=\left(\sum_{n=1}^{N_{tb}}\omega_{n}^{tb,1}|q_{\eta}(\overline{x}_{n})\mathcal{R}_{tb,\theta^{*}}(\overline{x}_{n})|^{2}\right)^{\frac{1}{2}},\\\
\mathcal{E}_{T,tb,2}:=\left(\sum_{n=1}^{N_{tb}}\omega_{n}^{tb,2}|\Delta\mathcal{R}_{tb,\theta^{*}}(\overline{x}_{n})|^{2}\right)^{\frac{1}{2}}.\end{cases}$
(4.9)

  * •

The spatial boundary condition training errors:
$\mathcal{E}_{T,sb}:=\mathcal{E}_{T,sb,0}+\mathcal{E}_{T,sb,1}+\mathcal{E}_{T,sb,2}$,
where

$\displaystyle\begin{cases}\mathcal{E}_{T,sb,0}:=\left(\sum_{n=1}^{N_{sb}}\omega_{n}^{sb,0}|\mathcal{R}_{sb,\theta^{*}}(\widehat{x}_{n},\widehat{t}_{n})|^{2}\right)^{\frac{1}{2}},\\\
\mathcal{E}_{T,sb,1}:=\left(\sum_{n=1}^{N_{sb}}\omega_{n}^{sb,1}|\partial_{t}\mathcal{R}_{sb,\theta^{*}}(\widehat{x}_{n},\widehat{t}_{n})|^{2}\right)^{\frac{1}{2}},\\\
\mathcal{E}_{T,sb,2}:=\left(\sum_{n=1}^{N_{sb}}\omega_{n}^{sb,2}|\partial_{t}^{2}\mathcal{R}_{sb,\theta^{*}}(\widehat{x}_{n},\widehat{t}_{n})|^{2}\right)^{\frac{1}{2}}.\end{cases}$
(4.10)

###  4.2 Proofs of the estimates.

Now we can state the theorem about the generalization error estimates.

######  Theorem 4.1.

Recall the errors defined in (4.6)-(4.10). Under Assumption 1, there exists a
unique solution to the inverse problem (1.1)-(1.3). Moreover, for the
approximate solution $(u^{*},q^{*})$ of the inverse problem with
$(\theta^{*},\eta^{*})$ being a global minimizer of the loss function
$J_{\lambda}^{N}(\theta,\eta)$, we have the following generalization error
estimates

$\displaystyle\mathcal{E}_{G,q}$ $\displaystyle\leq
C\Big{(}\mathcal{E}_{T,d}+\mathcal{E}_{T,int}+\mathcal{E}_{T,sb,1}+\mathcal{E}_{T,sb,2}+\mathcal{E}_{T,tb,1}+\mathcal{E}_{T,tb,2}+C_{q}^{\frac{1}{2}}N^{\frac{-\alpha}{2}}+O(\delta^{1/3})\Big{)},$
(4.11) $\displaystyle\mathcal{E}_{G,u}$ $\displaystyle\leq
C\Big{(}\mathcal{E}_{T,d}+\mathcal{E}_{T,int}+\mathcal{E}_{T,sb}+\mathcal{E}_{T,tb}+C_{q}^{\frac{1}{2}}N^{\frac{-\alpha}{2}}+O(\delta^{1/3})\Big{)},$

where

$\displaystyle N$
$\displaystyle=\min\left\\{N_{d},N_{int},N_{sb},N_{tb}\right\\},$
$\displaystyle\alpha$
$\displaystyle=\min\left\\{\alpha_{int,0},\alpha_{int,1},\alpha_{sb,0},\alpha_{sb,1},\alpha_{sb,2},\alpha_{tb,0},\alpha_{tb,1},\alpha_{d}\right\\},$

in (4.5), and

$C_{q}=\max\left\\{C_{q,0},C_{q,1},C_{qs,0},C_{qs,1},C_{qs,2},C_{qt,0},C_{qt,1},C_{qd}\right\\},$

with

$\displaystyle
C_{qd}=C_{qd}(\|\mathcal{L^{*}}\mathcal{R}_{d,\theta^{*}}\|_{C(\Omega)}),$
$\displaystyle
C_{q,0}=C_{q,0}(\|\mathcal{R}_{int,\theta^{*},\eta^{*}}\|_{C(\Omega_{T})}),$
$\displaystyle
C_{q,1}=C_{q,1}(\|\partial_{t}\mathcal{R}_{int,\theta^{*},\eta^{*}}\|_{C(\Omega_{T})}),$
$\displaystyle
C_{qs,0}=C_{qs,0}(\|\mathcal{R}_{sb,\theta^{*}}\|_{C(\partial\Omega_{T})}),$
$\displaystyle
C_{qs,1}=C_{qs,1}(\|\partial_{t}\mathcal{R}_{sb,\theta^{*}}\|_{C(\partial\Omega_{T})}),$
$\displaystyle
C_{qs,2}=C_{qs,2}(\|\partial_{t}^{2}\mathcal{R}_{sb,\theta^{*}}\|_{C(\partial\Omega_{T})}),$
$\displaystyle
C_{qt,0}=C_{qt,0}(\|\mathcal{R}_{tb,\theta^{*}}\|_{C(\Omega)}),$
$\displaystyle
C_{qt,1}=C_{qt,1}(\|\mathcal{L^{*}}\mathcal{R}_{tb,\theta^{*}}\|_{C(\Omega)}).$

The constant $C$ depends on $\|q^{*}\|_{L^{\infty}(\Omega)},$ $\Omega$ and
$T$.

###### Proof.

First, we introduce $\hat{u}:=u^{*}-u$ and realize that

$\begin{cases}\begin{aligned}
(\partial_{t}-\Delta+q^{*}(x))\hat{u}(x,t)&=\mathcal{R}_{{int},\theta^{*},\eta^{*}}(x,t)+(q-q^{*})u[q](x,t),&&(x,t)\in\Omega_{T},\\\
\hat{u}(x,t)&=\mathcal{R}_{sb,\theta^{*}}(x,t),&&(x,t)\in\partial\Omega_{T},\\\
\hat{u}(x,0)&=\mathcal{R}_{tb,\theta^{*}}(x),&&x\in\Omega,\end{aligned}\end{cases}$

with the final condition

$\displaystyle\hat{u}(x,T)=u[q^{*}](x,T)-u[q](x,T)=\mathcal{R}_{d,\theta^{*}}(x)-(\varphi-
G_{\epsilon}\varphi^{\delta}).$

We make the decomposition $\hat{u}:=\hat{u}_{1}+\hat{u}_{2}$, where
$\hat{u}_{1},\ \hat{u}_{2}$ satisfy

$\begin{cases}\begin{aligned}
(\partial_{t}-\Delta+q^{*}(x))\hat{u}_{1}(x,t)&=(q^{*}-q)(x)u[q](x,t),&&(x,t)\in\Omega_{T},\\\
\hat{u}_{1}(x,t)&=0,&&(x,t)\in\partial\Omega_{T},\\\
\hat{u}_{1}(x,0)&=0,&&x\in\Omega,\end{aligned}\end{cases}$ (4.12)

with

$\displaystyle\hat{u}_{1}(x,T)=\mathcal{R}_{d,\theta^{*}}(x)-(\varphi(x)-G_{\epsilon}\varphi^{\delta}(x))-\hat{u}_{2}(x,T),$
(4.13)

and

$\begin{cases}\begin{aligned}
(\partial_{t}-\Delta+q^{*}(x))\hat{u}_{2}(x,t)&=\mathcal{R}_{{int},\theta^{*},\eta^{*}}(x,t),&&(x,t)\in\Omega_{T},\\\
\hat{u}_{2}(x,t)&=\mathcal{R}_{sb,\theta^{*}}(x,t),&&(x,t)\in\partial\Omega_{T},\\\
\hat{u}_{2}(x,0)&=\mathcal{R}_{tb,\theta^{*}}(x),&&x\in\Omega,\end{aligned}\end{cases}$

respectively.

Define the operator $\mathcal{L}^{*}$ as
$\mathcal{L}^{*}\psi=(-\Delta+q^{*})\psi.$ With Assumption 1, we can apply
Lemma 3.1 to (4.12)-(4.13) and deduce that

$\displaystyle\|q^{*}-q\|_{L^{2}(\Omega)}$ (4.14) $\displaystyle\leq
C\|\mathcal{L}^{*}\hat{u}_{1}(\cdot,T)\|_{L^{2}(\Omega)}$
$\displaystyle=C\|\mathcal{L}^{*}\mathcal{R}_{d,\theta^{*}}-\mathcal{L}^{*}(\varphi-
G_{\epsilon}\varphi^{\delta})-\mathcal{L}^{*}\hat{u}_{2}(\cdot,T)\|_{L^{2}(\Omega)}$
$\displaystyle=C\|\mathcal{L}^{*}\mathcal{R}_{d,\theta^{*}}+\Delta\varphi-\Delta
G_{\epsilon}\varphi^{\delta}-q^{*}(\varphi-
G_{\epsilon}\varphi^{\delta})-\mathcal{L}^{*}\hat{u}_{2}(\cdot,T)\|_{L^{2}(\Omega)}$
$\displaystyle\leq
C\left(\|\mathcal{L}^{*}\mathcal{R}_{d,\theta^{*}}\|_{L^{2}(\Omega)}+\|\Delta\varphi-\Delta
G_{\epsilon}\varphi^{\delta}\|_{L^{2}(\Omega)}+\|q^{*}(\varphi-
G_{\epsilon}\varphi^{\delta})\|_{L^{2}(\Omega)}+\|\mathcal{L}^{*}\hat{u}_{2}(\cdot,T)\|_{L^{2}(\Omega)}\right),$

with $C=C(\|q^{*}\|_{L^{\infty}(\Omega)},\Omega,T)$. Using Lemma 2.1, we get

$\|\Delta\varphi-\Delta G_{\epsilon}\varphi^{\delta}\|_{L^{2}(\Omega)}\leq
C\epsilon(1+\delta\epsilon^{-3}).$ (4.15)

Also, we have

$\displaystyle|(\varphi-G_{\epsilon}\varphi^{\delta})|$
$\displaystyle\leq\int_{|x-y|\leq\epsilon}\rho_{\epsilon}(|x-y|)~{}|\varphi(x)-\varphi^{\delta}(y)|\;dy$
$\displaystyle\leq\int_{|x-y|\leq\epsilon}\rho_{\epsilon}(|x-y|)~{}|\varphi(x)-\varphi(y)|dy+\int_{|x-y|\leq\epsilon}\rho_{\epsilon}(|x-y|)~{}|\varphi(y)-\varphi^{\delta}(y)|\;dy$
$\displaystyle\leq C\epsilon+\delta.$

Thus, there holds

$\|q^{*}(\varphi-G_{\epsilon}\varphi^{\delta})\|_{L^{2}(\Omega)}\leq
C\epsilon+\delta.$ (4.16)

By straightforward computations, we have

$\displaystyle\|\mathcal{L}^{*}\hat{u}_{2}(\cdot,T)\|_{L^{2}(\Omega)}$
$\displaystyle\leq\|\partial_{t}\hat{u}_{2}(\cdot,T)\|_{L^{2}(\Omega)}+\|\mathcal{R}_{{int},\theta^{*},\eta^{*}}(\cdot,T)\|_{L^{2}(\Omega)}$
$\displaystyle\leq\|\partial_{t}\hat{u}_{2}\|_{L^{\infty}(0,T;L^{2}(\Omega))}+\|\mathcal{R}_{{int},\theta^{*},\eta^{*}}(\cdot,T)\|_{L^{2}(\Omega)}.$
(4.17)

Setting $w(x,t):=\partial_{t}\hat{u}_{2}(x,t)$, it satisfies

$\begin{cases}\begin{aligned}
(\partial_{t}-\Delta+q^{*})w(x,t)&=\partial_{t}\mathcal{R}_{{int},\theta^{*},\eta^{*}}(x,t),&&(x,t)\in\Omega_{T},\\\
w(x,t)&=\partial_{t}\mathcal{R}_{sb,\theta^{*}}(x,t),&&(x,t)\in\partial\Omega_{T},\\\
w(x,0)&=\mathcal{R}_{{int},\theta^{*},\eta^{*}}(x,0)-\mathcal{L}^{*}\mathcal{R}_{tb,\theta^{*}}(x),&&x\in\Omega.\end{aligned}\end{cases}$
(4.18)

Using the regularity theory for the direct problem (4.18), we obtain

$\displaystyle\|\partial_{t}\hat{u}_{2}\|_{L^{\infty}(0,T;L^{2}(\Omega))}$
$\displaystyle=\|w\|_{L^{\infty}(0,T;L^{2}(\Omega))}$ (4.19)
$\displaystyle\leq
C\big{(}\|\partial_{t}\mathcal{R}_{{int},\theta^{*},\eta^{*}}\|_{L^{2}(\Omega_{T})}+\|\partial_{t}\mathcal{R}_{{sb},\theta^{*}}\|_{H^{1}(0,T;L^{2}(\partial\Omega))}$
$\displaystyle\qquad\
+\|\mathcal{R}_{{int},\theta^{*},\eta^{*}}(\cdot,0)\|_{L^{2}(\Omega)}+\|\mathcal{L}^{*}\mathcal{R}_{tb,\theta^{*}}\|_{L^{2}(\Omega)}\big{)}.$

Combining (4.14)-(4.19) together and using the Sobolev embedding theorem, we
get

$\displaystyle\mathcal{E}_{G,q}$ $\displaystyle=\|q^{*}-q\|_{L^{2}(\Omega)}$
$\displaystyle\leq
C\big{(}\|\mathcal{R}_{{int},\theta^{*},\eta^{*}}\|_{H^{1}(0,T;L^{2}(\Omega))}+\|\mathcal{L}^{*}\mathcal{R}_{d,\theta^{*}}\|_{L^{2}(\Omega)}+\|\mathcal{L}^{*}\mathcal{R}_{{tb},\theta^{*}}\|_{L^{2}(\Omega)}$
$\displaystyle\
\qquad+\|\partial_{t}\mathcal{R}_{{sb},\theta^{*}}\|_{H^{1}(0,T;L^{2}(\partial\Omega))}+\epsilon(1+\delta\epsilon^{-3})+\epsilon+\delta\big{)}$
$\displaystyle\leq
C\big{(}\|\mathcal{R}_{{int},\theta^{*},\eta^{*}}\|_{H^{1}(0,T;L^{2}(\Omega))}+\|q^{*}\mathcal{R}_{d,\theta^{*}}\|_{L^{2}(\Omega)}+\|\Delta\mathcal{R}_{d,\theta^{*}}\|_{L^{2}(\Omega)}+\|q^{*}\mathcal{R}_{{tb},\theta^{*}}\|_{L^{2}(\Omega)}$
$\displaystyle\qquad\
+\|\Delta\mathcal{R}_{{tb},\theta^{*}}\|_{L^{2}(\Omega)}+\|\partial_{t}\mathcal{R}_{{sb},\theta^{*}}\|_{H^{1}(0,T;L^{2}(\partial\Omega))}+\epsilon(1+\delta\epsilon^{-3})+\epsilon+\delta\big{)},$

with $C=C(\|q^{*}\|_{L^{\infty}(\Omega)},\Omega,T)>0$. Picking
$\epsilon=O(\delta^{1/3})$, we achieve the estimate

$\displaystyle\mathcal{E}_{G,q}$ $\displaystyle=\|q^{*}-q\|_{L^{2}(\Omega)}$
(4.20) $\displaystyle\leq
C\big{(}\|\mathcal{R}_{{int},\theta^{*},\eta^{*}}\|_{H^{1}(0,T;L^{2}(\Omega))}+\|q^{*}\mathcal{R}_{d,\theta^{*}}\|_{L^{2}(\Omega)}+\|\Delta\mathcal{R}_{d,\theta^{*}}\|_{L^{2}(\Omega)}$
$\displaystyle\qquad\
+\|q^{*}\mathcal{R}_{{tb},\theta^{*}}\|_{L^{2}(\Omega)}+\|\Delta\mathcal{R}_{{tb},\theta^{*}}\|_{L^{2}(\Omega)}+\|\partial_{t}\mathcal{R}_{{sb},\theta^{*}}\|_{H^{1}(0,T;L^{2}(\partial\Omega))}+O(\delta^{1/3})\big{)}.$

Finally, we will evaluate the generalization error for the unknown $u(x,t)$
employing the obtained generalization error (4.20) of the potential function.
From the classical regularity theory for PDE, if $F(x,t)$ is sufficiently
smooth, then it holds that $u\in{L^{2}(0,T;L^{\infty}(\Omega))}$.
Consequently,

$\displaystyle\mathcal{E}_{G,u}$
$\displaystyle=\|\hat{u}\|_{C([0,T];L^{2}(\Omega))}$ $\displaystyle\leq
C\big{(}\|\mathcal{R}_{{int},\theta^{*},\eta^{*}}+(q^{*}-q)u[q]\|_{L^{2}(\Omega_{T})}+\|\mathcal{R}_{{sb},\theta^{*}}\|_{H^{1}(0,T;L^{2}(\partial\Omega))}+\|\mathcal{R}_{{tb},\theta^{*}}\|_{L^{2}(\Omega)}\big{)}$
$\displaystyle\leq
C\big{(}\|\mathcal{R}_{{int},\theta^{*},\eta^{*}}\|_{L^{2}(\Omega_{T})}+\|q^{*}-q\|_{L^{2}(\Omega)}\;\|u\|_{L^{2}(0,T;L^{\infty}(\Omega))}+\|\mathcal{R}_{{sb},\theta^{*}}\|_{H^{1}(0,T;L^{2}(\partial\Omega))}$
$\displaystyle\qquad\
+\|\mathcal{R}_{{tb},\theta^{*}}\|_{L^{2}(\Omega)}\big{)}$ $\displaystyle\leq
C\big{(}\|\mathcal{R}_{{int},\theta^{*},\eta^{*}}\|_{H^{1}(0,T;L^{2}(\Omega))}+\|q^{*}\mathcal{R}_{d,\theta^{*}}\|_{L^{2}(\Omega)}+\|\Delta\mathcal{R}_{d,\theta^{*}}\|_{L^{2}(\Omega)}+\|\mathcal{R}_{{tb},\theta^{*}}\|_{L^{2}(\Omega)}$
$\displaystyle\qquad\
+\|q^{*}\mathcal{R}_{{tb},\theta^{*}}\|_{L^{2}(\Omega)}+\|\Delta\mathcal{R}_{{tb},\theta^{*}}\|_{L^{2}(\Omega)}+\|\mathcal{R}_{{sb},\theta^{*}}\|_{H^{2}(0,T;L^{2}(\partial\Omega))}+O(\delta^{1/3})\big{)}.$

The proof is complete. ∎

The estimate (4.11) demonstrates that well-trained neural networks will
produce small generalization errors for the inverse problem. Specifically,
when all components of the training errors, including the interior PDE errors,
measurement errors as well as the initial and boundary value ones, are
sufficiently small, and the training sampling is large enough, the
generalization errors for inverse problem using neural networks can be limited
well. This differs from classical stability results that rely solely on the
knowledge of data. In this work, the generalization error estimates reflect
stability due to both the model itself and the reconstruction algorithm. From
Theorem 4.1, we see that we can limit the errors of both the inverse and
forward problems by controlling the residuals and the mollified parameter.
This provides important insights into the mathematical properties of our
approach and plays an important role on the construction of algorithms.

##  5 Numerical reconstructions.

###  5.1 Reconstruction algorithm.

The neural networks $u_{\theta}(x,t)$ and $q_{\eta}(x)$ depend on the
parameters $\theta$ and $\eta$ describing the networks information for
specific activation functions. Within the standard paradigm of deep learning,
one trains the networks by finding the optimal parameters
$(\theta^{*},\eta^{*})$ such that the loss function (4.4) is minimized. Our
target is the unknown solution of the inverse problem (1.1)-(1.3) and we wish
to find the trainable parameters $(\theta^{*},\eta^{*})$ such that the
corresponding neural networks $(u_{\theta^{*}},q_{\eta^{*}})$ approximate
$(u,q)$ well. More precisely, to solve (1.1)-(1.3) we first parameterize $u$
and $q$ by deep neural networks $u_{\theta}$ and $q_{\eta}$ with network
parameters $(\theta,\eta)$ respectively. Then, we design an appropriate loss
function, which is minimized to determine the parameters $(\theta,\eta)$.
Finally, a gradient-based method is applied to alternately update the network
parameters so that $(u_{\theta},q_{\eta})$ gradually approximates $(u,q)$ for
our inverse problem.

We provide a schematic of the neural networks in Figure 2. The left part
visualizes two unknowns as two standard neural networks parameterized by
$\theta$ and $\eta$, respectively. The right part applies the given physical
laws to the networks. B.C., I.C. and D are the boundary condition, initial
status and the measurement data obtained from random sample points in training
sets $\mathcal{S}_{sb}$, $\mathcal{S}_{tb}$ and $\mathcal{S}_{d}$
respectively. The training points in $\mathcal{S}_{int}$ are randomly sampled
as the PDE residuals points in interior spatio-temporal domain. The loss
function with some Sobolev norm is computed on the sample points, which can be
done efficiently through automatic differentiation (AD) in case of derivative
information. Minimizing the loss with respect to the parameters
$(\theta,\eta)$ alternately produces $u_{\theta^{\ast}}$ and
$q_{\eta^{\ast}}$, which serves as the approximation to the solution of the
inverse problem.

Figure 2: The schematic of the deep neural networks for solving the inverse
potential problem.

With the support of Theorem 4.1, we can construct the proposed Algorithm 1 for
solving the inverse problem (1.1)-(1.3).

Algorithm 1 Data-driven solution of the inverse problem (1.1)-(1.3).

0: Noisy data $\varphi^{\delta}$ for the inverse problem.

Initialize Network architectures $\left(q_{\eta},u_{\theta}\right)$ and
parameters $(\eta,\theta)$.

for $j=1,\cdots,K$ do

Sample $\mathcal{S}_{d},\ \mathcal{S}_{int},\ \mathcal{S}_{sb},\
\mathcal{S}_{tb}$.

$\eta\leftarrow\operatorname{Adam}\left(-\nabla_{\eta}J_{\lambda}^{N}(\eta,\theta),\tau_{\eta},\lambda_{\eta}\right)$,$\theta\leftarrow\operatorname{Adam}\left(-\nabla_{\theta}J_{\lambda}^{N}(\eta,\theta),\tau_{\theta},\lambda_{\theta}\right),$

end for

$\left(q_{\eta},u_{\theta}\right)$.

The above minimization problem is to search a minimizer of a possibly non-
convex function $J_{\lambda}^{N}(\theta,\eta)$ over
$\Theta\subset\mathbb{R}^{\mathcal{M}}$ for possibly very large $\mathcal{M}$.
The hyper-parameters $(\tau_{\eta},\tau_{\theta})$ are learning rates and
$(\lambda_{\eta},\lambda_{\theta})$ are balance hyper-parameters between PDE
and measurement data residuals. The robust analysis for hyper-parameters
$\lambda$ and the architectures of neural networks are studied in the next
subsection.

The optimizer in Algorithm 1 is Adam (Adaptive Moment Estimation), which is an
optimization algorithm commonly used in deep learning for training neural
networks. The key idea of Adam is to adaptively adjust the learning rate for
each parameter based on estimates of both the first-order moment (the mean)
and the second-order moment (the uncentered variance) of the gradients. This
adaptation helps Adam to perform well in different types of optimization
problems. The algorithm maintains an exponentially moving average of gradients
($m_{t}$) and squared gradients ($V_{t}$) for each parameter. At each
iteration, Adam updates the parameters using a combination of these moving
average estimates. It incorporates bias correction to account for the fact
that the estimates are biased towards zero at the beginning of training. We
set $g_{t}$ be gradients w.r.t. stochastic objective at timestep $t$,
$\beta_{1},\beta_{2}\in[0,1)$ be the exponential decay rates for the moment
estimates and $\tau$ be the initial learning rate. Good default settings for
the tested machine learning problems are $\tau=0.001$, $\beta_{1}=0.9$,
$\beta_{2}=0.999$ and fuzzy factor $\epsilon=10^{-8}$. The updates are
calculated as follows:

  * (1)

Initialize the first moment vector $m$ and the second moment vector $V$ with
zeros for each parameter:

$m_{0}=V_{0}=0.$

  * (2)

Update the first moment estimate $m$ using a weighted average of the current
gradient $g_{t}$ and the previous first moment estimate $m_{t-1}$:

$m_{t}=\beta_{1}m_{t-1}+(1-\beta_{1})g_{t}.$

  * (3)

Update the second moment estimate $V$ using a weighted average of the squared
gradients and the previous second moment estimate $V_{t-1}$:

$V_{t}=\beta_{2}V_{t-1}+(1-\beta_{2})g_{t}^{2}.$

  * (4)

Calculate the bias-corrected first and second moment estimate to correct for
their initialization bias:

${\hat{m}}_{t}={\frac{m_{t}}{1-(\beta_{1})^{t}}},\quad{\hat{V}}_{t}={\frac{V_{t}}{1-(\beta_{2})^{t}}}.$

  * (5)

Update the parameters $\xi$ by moving in the direction of the first moment
estimate, where the learning rate is $\tau$ divided by the square root of the
second moment estimate:

$\xi_{t}=\xi_{t-1}-\tau\frac{\hat{m}_{t}}{\sqrt{\hat{V}_{t}}+\epsilon}.$

The hyper-parameters in Adam include the learning rate and the decay rates for
the moving averages. These hyper-parameters need to be tuned based on the
specific problem and dataset to achieve optimal performance. Adam has several
advantages that make it popular in deep learning:

  * (a)

Adaptive learning rate: Adam automatically adapts the learning rate for each
parameter based on the estimated first and second moments. This adaptive
behavior helps in effectively navigating the optimization landscape and can
lead to faster convergence.

  * (b)

Efficiency: Adam uses the moving averages to maintain a history of gradients,
which eliminates the need to store and compute gradients for each iteration
separately. This makes Adam memory-efficient and allows for efficient
parallelization during training.

  * (c)

Robustness: Adam performs well across a wide range of optimization problems
and is less sensitive to hyper-parameter tuning compared to some other
optimizers. It can handle sparse gradients and noisy data effectively.

The proposed algorithm, which utilizes (4.2) as the loss function, exhibits
superior performance in recovering smooth solutions due to the high regularity
of the PDEs residual term. This regularity term promotes smoother solutions
and is an important factor in achieving higher accuracy. Furthermore, the use
of automatic differentiation (AD) implementations enables the efficient
calculation of the necessary derivatives. This feature is a significant
advantage of our approach as it allows for the accurate optimization of the
objective function, which is crucial for effective solution of inverse
problems. To validate the effectiveness of the proposed algorithm and to
substantiate our claims, we conduct a series of numerical experiments.

###  5.2 Numerical experiments.

In this subsection, we will present several numerical examples for the spatial
domain $\Omega\subset\mathbb{R}^{d}$ with $d=2,3$. We define the following
relative errors for exact solutions $(u,q)$ and numerical approximations
$(u^{*},q^{*})$ as

$\displaystyle{Re}_{u}:=\frac{\|u-u^{*}\|_{L^{2}(\Omega_{T})}}{\|u\|_{L^{2}(\Omega_{T})}},\quad{Re}_{q}:=\frac{\|q-q^{*}\|_{L^{2}(\Omega)}}{\|q\|_{L^{2}(\Omega)}},\quad{Re}_{\Delta
u}:=\frac{\|\Delta u-\Delta u^{*}\|_{L^{2}(\Omega_{T})}}{\|\Delta
u\|_{L^{2}(\Omega_{T})}}.$

#### Example 1 (two-dimensional experiment):

For equation (1.1), we set the exact solution $u$ and the domain $\Omega_{T}$
as

$u(x,y,t)=(x^{2}+y^{2}+1)\exp(t),\quad(t,x,y)\in\Omega_{T}=[0,1]^{3}.$

The exact potential $q$ is given as

$q(x,y)=\sin(\pi x)\sin(\pi y).$

The initial and boundary conditions can be calculated from the representation
of $u$ straightforwardly. The exact measurement will be

$\displaystyle u(x,y,1)=\varphi(x,y)=(x^{2}+y^{2}+1)\exp(1),$

and in our experiments the noisy data is set as

$\displaystyle\varphi^{\delta}(x,y):=\varphi(x,y)+\delta\cdot(2\;\text{rand}(\text{shape}(\varphi(x,y)))-1),$
(5.1)

where $\text{rand}(\text{shape}(\varphi))$ is a random variable generated by
uniform distribution in $[0,1]$.

For the implementation details, we use a fully connected neural network for
$u_{\theta}$ and $q_{\eta}$ with 3 hidden layers, each with a width of 20. We
take

$\text{N}=N_{int}+N_{sb}+N_{tb}+N_{d}=256+256\times 4+256+256=1792$

as the number of collocation points, which are randomly sampled in four
different domains, i.e., interior spatio-temporal domain, spatial and temporal
boundary domain, and additional measurement domain. The activation function is
$\tanh(x)$ and the hyper-parameter is $\lambda=0.01$. The number of training
epochs is set to be $5\times 10^{4}$, and the initial learning rates
$\tau_{\theta},\tau_{\eta}$ both start with $0.001$ and shrink $10$ times
every $2\times 10^{4}$ iteration. The test sets are chosen by a uniform mesh

$\displaystyle\mathcal{T}:=\\{(t_{k},x_{i},y_{j}):k,i,j=0,1,\cdots,49\\}\subset\Omega_{T}.$
(5.2)

Since the noisy level of the measurement data affects the reconstruction
accuracy, in this simulation, we test the training performance for various
noisy levels. Figure 3 records the training process, i.e., the training loss,
the relative error for the reconstruction of $q$ and the relative error for
the recovery of $u$ with respect to the iterations for different noise levels
$\delta=0.1\%,\;1\%,\;5\%,\;10\%$ by the proposed scheme. After training, we
test the reconstruction result on test sets $\mathcal{T}$. The distribution of
the temperature field $u(x,t)$ also depends on the time $t$, Figure 4 shows
the time-series relative error of the recovered $u$ with various noise levels
after logarithmic re-scaling. As shown in these figures, the training
performance deteriorates as the noise level of the measurement data
increasing.

Figure 3: The training loss (left), the relative error for $q$ (middle), the
relative error for $u$ (right) after logarithmic re-scaling.  Figure 4: The
time series relative error (test) for $u$ for different noisy data.

Figure 5 shows the exact solution of the ptential term. Figure 6 shows the
reconstruction results for $q(x)$ by optimizing the proposed loss function
(first line) and the corresponding absolute pointwise error for various noisy
level $\delta=0.1\%,5\%,10\%$ (second line). Meanwhile, Figure 7 presents the
reconstruction solution $u$ (first line) and corresponding absolute pointwise
error (second line) for various noisy level measurement data at $t=1/7$. We
can see that the reconstruction accuracy for $q$ deteriorates as the noise
level of the measurement data increasing, but the performance for $u$ is still
satisfactory.

Figure 5: The exact potential function $q$. Figure 6: The reconstruction of
potential function $q$ (upper) and corresponding absolute pointwise error
$|q-q^{*}|$ (bottom) for various noisy level measurement data. Figure 7: The
reconstruction of solution $u$ (upper) and corresponding absolute pointwise
error $|u-u^{*}|$ (bottom) for various noisy level measurement data at
$t=1/7$.

Table 1 presents the recovery results solved by two schemes: (I) the proposed
frameworks with the loss function (4.2), (II) the DGM frameworks with the loss
function (4.3). We record the generalization error of $q$, $u$ and $\Delta u$
in $L^{2}$-error from the noisy input data with $\delta=0.01$. Due to the
random sampling of the training data points, the inversion results have some
stochasticity. Thus, we perform Algorithm 1 with the loss function in the
formulation (4.2) and formulation (4.3) five times, respectively. The relative
errors (mean and standard deviation) for the recovery of $q$, $u$ and $\Delta
u$ are shown in Table 1. As observed, optimizing the loss function proposed in
this paper leads to more accurate recovery results, especially for the
reconstruction of $q$ compared with DGM frameworks. Moreover, although the
reconstruction accuracy of $u$ in $L^{2}$-error for both two frameworks are
relatively close, the accuracy of $\Delta u$ in $L^{2}$-error for proposed
scheme in this paper performs better. This suggests that the proposed
frameworks are better able to capture smooth solutions.

Table 1:  The inversion results with noisy measuremrnt data ($\delta=0.01$) solved by two schemes: (I) the proposed frameworks with the loss function (4.2), (II) the DGM frameworks with the loss function (4.3). | $Re_{q}$ | $Re_{u}$ | $Re_{\Delta u}$
---|---|---|---
I | $0.5011\%\pm 0.0102\%$ | $0.0215\%\pm 0.0049\%$ | $0.2550\%\pm 0.0169\%$
II | $12.4545\%\pm 4.8016\%$ | $0.1148\%\pm 0.0133\%$ | $2.3468\%\pm 0.8133\%$

#### Example 2 (two-dimensional experiment):

For equation (1.1), we set the exact solution $u$ and the domain $\Omega_{T}$
as

$u(x,y,t)=t\exp(x+y),\quad(t,x,y)\in\Omega_{T}=[0,1]\times[0,2]^{2}.$

The exact potential $q$ is given as

$\displaystyle q(r)$ $\displaystyle=\begin{cases}15\left(\cos
r-\sqrt{3}/2\right)+2,\quad 0\leq r\leq\pi/6,\\\ 2,\quad\text{ otherwise
},\end{cases}$ $\displaystyle r(x,y)$
$\displaystyle=\sqrt{(x-1)^{2}+(y-1)^{2}}.$

The exact measurement will be

$\displaystyle u(x,y,1)=\varphi(x,y)=\exp(x+y),$

and the noisy data $\varphi^{\delta}$ is generated by (5.1). The network
architectures and hyper-parameters such as activation function, balance hyper-
parameter $\lambda$ are all the same as Example 1. The number of training
epochs is set to be $1\times 10^{5}$, and the initial learning rates
$\tau_{\theta},\tau_{\eta}$ both start with $0.001$ and shrink $10$ times
every $2\times 10^{4}$ iteration. The test sets are chosen by a uniform mesh
as (5.2)

In this simulation, we evaluate the training performance under various levels
of measurement noise. The training process under Algorithm 1 is recorded in
Figure 8, which includes the training loss, and the relative errors for the
reconstruction of $q$ and $u$ during training process for different noise
levels ($\delta=0,\ 1\%,\ 5\%,\ 10\%$). Figure 9 displays the exact potential
function $q$, while the approximated $q^{*}$ under different noise level
measurements ($\delta=1\%,\ 5\%,\ 10\%$) are shown in Figure 10. We can see
that the numerical reconstructions still satisfy the theoretical results even
with zero initial condition and nonsmooth exact potential. This means that in
numerical reconstructions we may release the conditions of Assumption 1 to
some extent.

Figure 8: The training loss (left), the relative error for $q$ (middle), the
relative error for $u$ (right) after logarithmic re-scaling.  Figure 9: The
exact potential function $q$. Figure 10: The reconstruction of potential
function $q$ (upper) and corresponding absolute pointwise error $|q-q^{*}|$
(bottom) for various noisy level measurement data.

Now, we start to verify the convergence of the iteration in Theorem 4.1 with
different neural network architectures. In the experiments, for a fixed number
of per-layer neurons $\text{NN}=20$, we compute the reconstruction errors for
$q$ versus the noise level $\delta$ using logarithmic re-scaling with various
hidden layers $\text{NL}=3,4,6$. The results of these experiments are
presented in the left of Figure 11. The theoretical estimate $O(\delta^{1/3})$
is shown by the black line. Similarly, fixing hidden layer $\text{NL}=6$, the
reconstruction errors for $q$ under various per-layer neurons
$\text{NN}=10,15,20$ are given in the right of Figure 11. From this figure, we
see that the error could be bounded by the rate $\delta^{1/3}$ to some extent,
which supports the theoretical analysis in Theorem 4.1.

Figure 11: The reconstruction error of $q$ versus noise level $\delta$ under
various layer numbers for fixed per-layer neuron numbers $\text{NN}=20$
(left), and under various per-layer neuron for fixed layer numbers
$\text{NL}=6$ (right).

In order to evaluate the effectiveness of the proposed scheme in terms of
hyper-parameters and network structure, a series of experiments are conducted.
Specifically, we examine the impact of the balance hyper-parameter $\lambda$
in (4.2) and the network structure, including the number of hidden layers and
neurons. For a fixed number of hidden layers ($\text{NL}=3$) and a fixed
number of neurons per-layer ($\text{NN}=20$), we compute the reconstruction
errors (mean and standard deviation) for $q$ and $u$ using various values of
$\lambda$, such as $\lambda=10^{j},\ -4\leq j\leq 1$. The results of these
experiments are presented in Table 2, which indicates that the performance of
the inverse problem is highly dependent on the balance hyper-parameter
$\lambda$. Specifically, we find that the relative reconstruction errors are
optimized when $\lambda$ is set to $10^{-2}$. Furthermore, we observe that the
reconstruction errors increase significantly as $\lambda$ exceeds this optimal
value. These results suggest that the selection of the balance hyper-parameter
is critical to achieving good performance in this inverse problem.

Next we experiment with various combinations of hidden layers and neuron
numbers for the inverse problem using Algorithm 1. We set the dimension to
$d=2$ and try a total of $16$ combinations of hidden layers (NL) and per-layer
neuron numbers (NN), with

$\text{NL}=3,6,9,14,\quad\text{NN}=10,15,20,25.$

For each combination, we run Algorithm 1 for $1\times 10^{5}$ iterations and
record the relative errors (mean and standard deviation) for $(q,u)$ in Table
3. It indicates that deeper (larger NL) and/or wider (larger NN) neural
networks tend to yield lower reconstruction errors, although this causes
higher computational cost. However, we also observe that for fixed neuron
number NN, increasing the number of hidden layers NL, for example
$\text{NL}\geq 15$, causes the algorithm fail to converge as the number of
iterations increases. This suggests that increasing the number of layers
and/or neurons can enhance the representation capacity of neural networks. But
it may also introduce more parameters to train, and lead to longer training
times and potential overfitting of the representation.

Table 2: The relative errors for $q$ (first line) and $u$ (second line) using various hyper-parameter $\lambda$ from noisy measurement data with $\delta=0.01$, where the hidden layers $\text{NL}=3$, per-layer neurons $\text{NN}=20$.  Error | $\lambda=10^{-4}$ | $\lambda=10^{-3}$ | $\lambda=10^{-2}$
---|---|---|---
${Re}_{q}$ | $2.5113\%\pm 1.1339\%$ | $0.9942\%\pm 0.1121\%$ | $0.7124\%\pm 0.0639\%$
${Re}_{u}$ | $0.1374\%\pm 0.0652\%$ | $0.1216\%\pm 0.0320\%$ | $0.0560\%\pm 0.0204\%$
Error | $\lambda=10^{-1}$ | $\lambda=1$ | $\lambda=10$
${Re}_{q}$ | $0.8388\%\pm\underline{}0.2448\%$ |  $1.0809\%\pm 0.1741$% | $2.2891\%\pm 1.3302\%$
${Re}_{u}$ | $0.0331\%\pm 0.0094\%$ | $0.0478\%\pm 0.0122\%$ | $0.0874\%\pm 0.0110\%$
Table 3: The relative errors for $q$ (top) and $u$ (bottom) using various combinations of (NL,NN) for problem dimension $d=2$ from noisy measurement data with $\delta=0.01$, where the hyper-parameter $\lambda=10^{-2}$.  $NL$ | $NN=10$ | $NN=15$ | $NN=20$ | $NN=25$
---|---|---|---|---
$3$ | $1.7258\%\pm 0.6258\%$ | $1.1664\%\pm 0.3531\%$ | $0.7124\%\pm 0.0639\%$ | $0.5978\%\pm 0.1071\%$
$6$ | $1.5015\%\pm 0.6640\%$ | $0.7234\%\pm 0.2407\%$ | $0.4482\%\pm 0.1065\%$ | $0.3093\%\pm 0.0310\%$
$9$ | $0.7797\%\pm 0.1577\%$ | $0.4788\%\pm 0.1598\%$ | $0.4796\%\pm 0.0900\%$ | $0.4874\%\pm 0.2034\%$
$14$ | $0.5602\%\pm 0.2118\%$ | $0.4077\%\pm 0.1990\%$ | $0.3958\%\pm 0.0552\%$ | $0.4108\%\pm 0.0637\%$
$3$ | $0.0841\%\pm 0.0216\%$ | $0.0615\%\pm 0.0290\%$ | $0.0560\%\pm 0.0204\%$ | $0.0306\%\pm 0.0105\%$
$6$ | $0.0575\%\pm 0.0238\%$ | $0.0490\%\pm 0.0062\%$ | $0.0280\%\pm 0.0071\%$ | $0.0358\%\pm 0.0090\%$
$9$ | $0.0751\%\pm 0.0129\%$ | $0.0433\%\pm 0.0037\%$ | $0.0523\%\pm 0.0166\%$ | $0.0314\%\pm 0.0132\%$
$14$ | $0.0806\%\pm 0.0103\%$ | $0.0612\%\pm 0.0144\%$ | $0.0593\%\pm 0.0145\%$ | $0.0422\%\pm 0.0115\%$

#### Example 3 (three-dimensional experiment):

We also take the following 3-dimensional experiment. We set the exact solution
and the domain of equation (1.1) as

$\displaystyle u(x,y,t)=t\exp(x+y+z),\quad(t,x,y,z)\in\Omega_{T}=[0,1]^{4}.$

The exact potential $q$ is given as

$\displaystyle q(x,y,z)=x+y+z.$

We also employ fully-connected neural networks with $\text{NL}=4,\
\text{NN}=20$ for both $u_{\theta}$ and $q_{\eta}$. The number of training
points are

$\text{N}=N_{int}+N_{sb}+N_{tb}+N_{d}=256+256\times 6+256+256=2304,$

which are randomly sampled from four different training sets. The other
network architectures and hyper-parameters such as activation function,
balance hyper-parameter $\lambda$, number of training epochs, and initial
learning rate are all the same as Example 1. The test sets are chosen by a
uniform mesh

$\mathcal{T}:=\\{(t_{k},x_{i},y_{j},z_{l}):\
k,i,j,l=0,1,\cdots,49\\}\subset\Omega_{T}.$

Figure 12 shows the exact potential function $q$ and the relative errors
versus iterations during training process for different noise scales
$\delta=0,\ 1\%,\ 5\%,\ 10\%$. Figure 13 presents the potential functions
$q_{\eta}$ recovered from different noise levels and the corresponding point
by point absolute errors on test sets. The inversion results are satisfactory
and reasonable overall.

Figure 12: The exact potential function $q$ (left), the relative error for
recovered $q$ with different noise levels (right). Figure 13: The
reconstruction of potential function $q$ (first line), absolute pointwise
error $|q-q^{*}|$ (second line) for various noise levels.

Finally, we conduct the experiments to evaluate the robustness of the proposed
scheme in terms of network structure (number of hidden layers and per-layer
nurons). More specifically, we run Algorithm 1 with per-layer neuron numbers
$\text{NN}=5,10,20,25$ for fixed hidden layers $\text{NL}=6$ and with hidden
layers $\text{NL}=4,6,8,10$ for fixed per-layer neuron numbers $\text{NN}=25$,
respectively. The reconstruction errors are presented in Figure 14. It also
seems that larger NL and/or larger NN neural networks yield lower
reconstruction error. In this example, for fixed hidden layer $\text{NL}=6$,
we test the per-layer neuron numbers $\text{NN}\geq 30$ and find that there
will be bad reconstruction result, that is, the relative error for $q$ with
$\text{NL}=6,\text{NN}=30$ is larger than the case $\text{NN}=10,20,25$ as the
iterations increasing. Therefore, more layers and/or neurons with much more
parameters to train, yield longer training time, and may result in overfitting
of the reconstruction.

Figure 14: The reconstruction error of $q$ under various per-layer neuron
numbers with fixed layer numbers NL=6 (left), the reconstruction error of $q$
under various layer numbers with fixed per-layer neuron number NN=25 (right).

##  6 Concluding remarks.

In this work, a deep neural network-based reconstruction scheme has been
proposed to solve an inverse potential problem in the parabolic equation. The
proposed method has shown superior performance in high-dimensional space. We
prove the uniqueness of the inverse potential problem. A new loss function has
been introduced, which includes regularization terms that depend on the
derivatives of the residuals for both the partial differential equation and
the measurement data. These regularization terms aim to address the ill-
posedness of the inverse problem and enhance the regularity of the solution.
Additionally, the mollification method has been employed to improve the
regularity of the noisy data, where it can reduce the perturbation errors
caused by numerical differentiation on the noisy data. Generalization
estimates based on the conditional stability of linear inverse source problems
and the mollification error estimate on noisy data have been established,
which provide a measure of the stability and accuracy of the proposed method
in solving the inverse potential problem. Numerical experiments have been
conducted to evaluate the performance of proposed method, which indicate the
efficiency the approach in this work.

## Acknowledgments

Mengmeng Zhang is supported by Foundation of Hebei University of Technology
(Grant No.282022550) and Foundation of Tianjin Education Commission Research
Program(Grant No.2022KJ102). Zhidong Zhang is supported by National Natural
Science Foundation of China (Grant No. 12101627).

## References

  * [1] Gang Bao, Xiaojing Ye, Yaohua Zang, and Haomin Zhou.  Numerical solution of inverse problems by weak adversarial networks.  Inverse Problems, 36(11):115003, 31, 2020.  URL: https://doi.org/10.1088/1361-6420/abb447, doi:10.1088/1361-6420/abb447.
  * [2] K. Cao and D. Lesnic.  Reconstruction of the space-dependent perfusion coefficient from final time or time-average temperature measurements.  J. Comput. Appl. Math., 337:150–165, 2018.  URL: https://doi.org/10.1016/j.cam.2018.01.010, doi:10.1016/j.cam.2018.01.010.
  * [3] De-Han Chen, Daijun Jiang, and Jun Zou.  Convergence rates of Tikhonov regularizations for elliptic and parabolic inverse radiativity problems.  Inverse Problems, 36(7):075001, 21, 2020.  URL: https://doi.org/10.1088/1361-6420/ab8449, doi:10.1088/1361-6420/ab8449.
  * [4] Qun Chen and Jijun Liu.  Solving an inverse parabolic problem by optimization from final measurement data.  J. Comput. Appl. Math., 193(1):183–203, 2006.  URL: https://doi.org/10.1016/j.cam.2005.06.003, doi:10.1016/j.cam.2005.06.003.
  * [5] Yuyao Chen, Lu Lu, George Em Karniadakis, and Luca Dal Negro.  Physics-informed neural networks for inverse problems in nano-optics and metamaterials.  Opt. Express, 28(8):11618–11633, 2020.  URL: https://opg.optica.org/oe/abstract.cfm?URI=oe-28-8-11618, doi:10.1364/OE.384875.
  * [6] Mourad Choulli and Masahiro Yamamoto.  Generic well-posedness of an inverse parabolic problem—the Hölder-space approach.  Inverse Problems, 12(3):195–205, 1996.  URL: https://doi.org/10.1088/0266-5611/12/3/002, doi:10.1088/0266-5611/12/3/002.
  * [7] Mourad Choulli and Masahiro Yamamoto.  An inverse parabolic problem with non-zero initial condition.  Inverse Problems, 13(1):19–27, 1997.  URL: https://doi.org/10.1088/0266-5611/13/1/003, doi:10.1088/0266-5611/13/1/003.
  * [8] Zui-Cha Deng, Liu Yang, and Jian-Ning Yu.  Identifying the radiative coefficient of heat conduction equations from discrete measurement data.  Appl. Math. Lett., 22(4):495–500, 2009.  URL: https://doi.org/10.1016/j.aml.2008.06.023, doi:10.1016/j.aml.2008.06.023.
  * [9] Chenguang Duan, Yuling Jiao, Yanming Lai, Dingwei Li, Xiliang Lu, and Jerry Zhijian Yang.  Convergence rate analysis for deep Ritz method.  Commun. Comput. Phys., 31(4):1020–1048, 2022.  URL: https://doi.org/10.4208/cicp.oa-2021-0195, doi:10.4208/cicp.oa-2021-0195.
  * [10] Weinan E and Bing Yu.  The deep Ritz method: a deep learning-based numerical algorithm for solving variational problems.  Commun. Math. Stat., 6(1):1–12, 2018.  URL: https://doi.org/10.1007/s40304-018-0127-z, doi:10.1007/s40304-018-0127-z.
  * [11] Lawrence C. Evans.  Partial differential equations, volume 19 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, second edition, 2010.  URL: https://doi.org/10.1090/gsm/019, doi:10.1090/gsm/019.
  * [12] Ian Goodfellow, Yoshua Bengio, and Aaron Courville.  Deep Learning.  MIT Press, 2016.  http://www.deeplearningbook.org.
  * [13] Qingguo Hong, Jonathan W. Siegel, and Jinchao Xu.  Rademacher complexity and numerical quadrature analysis of stable neural networks with applications to numerical pdes.  arXiv preprint.  doi:arXiv:2104.02903.
  * [14] Victor Isakov.  Inverse problems for partial differential equations, volume 127 of Applied Mathematical Sciences.  Springer, Cham, third edition, 2017.  URL: https://doi.org/10.1007/978-3-319-51658-5, doi:10.1007/978-3-319-51658-5.
  * [15] Yuling Jiao, Yanming Lai, Dingwei Li, Xiliang Lu, Fengru Wang, Yang Wang, and Jerry Zhijian Yang.  A rate of convergence of physics informed neural networks for the linear second order elliptic PDEs.  Commun. Comput. Phys., 31(4):1272–1295, 2022.  URL: https://doi.org/10.4208/cicp.oa-2021-0186, doi:10.4208/cicp.oa-2021-0186.
  * [16] Bangti Jin, Xiliang Lu, Qimeng Quan, and Zhi Zhou.  Convergence rate analysis of Galerkin approximation of inverse potential problem.  Inverse Problems, 39(1):Paper No. 015008, 26, 2023.  URL: https://doi.org/10.1088/1361-6420/aca70e, doi:10.1088/1361-6420/aca70e.
  * [17] Barbara Kaltenbacher and William Rundell.  Recovery of multiple coefficients in a reaction-diffusion equation.  Journal of Mathematical Analysis and Applications, 481(1):123475, 2020.  URL: https://www.sciencedirect.com/science/article/pii/S0022247X19307437, doi:https://doi.org/10.1016/j.jmaa.2019.123475.
  * [18] V. L. Kamynin and A. B. Kostin.  Two inverse problems of the determination of a coefficient in a parabolic equation.  Differ. Uravn., 46(3):372–383, 2010.  URL: https://doi.org/10.1134/S0012266110030079, doi:10.1134/S0012266110030079.
  * [19] Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis.  Learning nonlinear operators via deeponet based on the universal approximation theorem of operators.  Nature machine intelligence, 3(3):218–229, 2021.
  * [20] Lu Lu, Xuhui Meng, Zhiping Mao, and George Em Karniadakis.  DeepXDE: a deep learning library for solving differential equations.  SIAM Rev., 63(1):208–228, 2021.  URL: https://doi.org/10.1137/19M1274067, doi:10.1137/19M1274067.
  * [21] Yulong Lu, Jianfeng Lu, and Min Wang.  A priori generalization analysis of the deep ritz method for solving high dimensional elliptic partial differential equations.  volume 134 of Proceedings of Machine Learning Research, pages 3196–3241. PMLR, 2021.  URL: https://proceedings.mlr.press/v134/lu21a.html.
  * [22] Tao Luo and Haizhao Yang.  Two-layer neural networks for partial differential equations: Optimization and generalization theory, 2020.  arXiv:2006.15733.
  * [23] Siddhartha Mishra and Roberto Molinaro.  Estimates on the generalization error of physics-informed neural networks for approximating a class of inverse problems for PDEs.  IMA J. Numer. Anal., 42(2):981–1022, 2022.  URL: https://doi.org/10.1093/imanum/drab032, doi:10.1093/imanum/drab032.
  * [24] Aleksei Ivanovich Prilepko and Andrew Borisovich Kostin.  On certain inverse problems for parabolic equations with final and integral observation.  Matematicheskii Sbornik, 183(4):49–68, 1992.  URL: https://dx.doi.org/10.1070/SM1993v075n02ABEH003394, doi:10.1070/SM1993v075n02ABEH003394.
  * [25] Aleksey I. Prilepko, Dmitry G. Orlovsky, and Igor A. Vasin.  Methods for solving inverse problems in mathematical physics, volume 231 of Monographs and Textbooks in Pure and Applied Mathematics.  Marcel Dekker, Inc., New York, 2000.
  * [26] M. Raissi, P. Perdikaris, and G. E. Karniadakis.  Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  J. Comput. Phys., 378:686–707, 2019.  URL: https://doi.org/10.1016/j.jcp.2018.10.045, doi:10.1016/j.jcp.2018.10.045.
  * [27] Maziar Raissi, Alireza Yazdani, and George Em Karniadakis.  Hidden fluid mechanics: A navier-stokes informed deep learning framework for assimilating flow visualization data, 2018.  arXiv:1808.04327.
  * [28] William Rundell.  The determination of a parabolic equation from initial and final data.  Proc. Amer. Math. Soc., 99(4):637–642, 1987.  URL: https://doi.org/10.2307/2046467, doi:10.2307/2046467.
  * [29] Tim De Ryck, Ameya D. Jagtap, and Siddhartha Mishra.  Error estimates for physics informed neural networks approximating the navier-stokes equations, 2023.  arXiv:2203.09346.
  * [30] Yeonjong Shin, Jérôme Darbon, and George Em Karniadakis.  On the convergence of physics informed neural networks for linear second-order elliptic and parabolic type PDEs.  Commun. Comput. Phys., 28(5):2042–2074, 2020.  URL: https://doi.org/10.4208/cicp.oa-2020-0193, doi:10.4208/cicp.oa-2020-0193.
  * [31] Yeonjong Shin, Zhongqiang Zhang, and George Em Karniadakis.  Error estimates of residual minimization using neural networks for linear pdes, 2020.  arXiv:2010.08019.
  * [32] Khemraj Shukla, Patricio Clark Di Leoni, James Blackshire, Daniel Sparkman, and George Em Karniadakis.  Physics-informed neural network for ultrasound nondestructive quantification of surface breaking cracks.  Journal of Nondestructive Evaluation, 39:1–20, 2020.
  * [33] Justin Sirignano and Konstantinos Spiliopoulos.  DGM: a deep learning algorithm for solving partial differential equations.  J. Comput. Phys., 375:1339–1364, 2018.  URL: https://doi.org/10.1016/j.jcp.2018.08.029, doi:10.1016/j.jcp.2018.08.029.
  * [34] Dumitru Trucu, Derek B. Ingham, and Daniel Lesnic.  Space-dependent perfusion coefficient identification in the transient bio-heat equation.  J. Engrg. Math., 67(4):307–315, 2010.  URL: https://doi.org/10.1007/s10665-009-9319-6, doi:10.1007/s10665-009-9319-6.
  * [35] Masahiro Yamamoto and Jun Zou.  Simultaneous reconstruction of the initial temperature and heat radiative coefficient.  volume 17, pages 1181–1202. 2001.  Special issue to celebrate Pierre Sabatier’s 65th birthday (Montpellier, 2000).  URL: https://doi.org/10.1088/0266-5611/17/4/340, doi:10.1088/0266-5611/17/4/340.
  * [36] Yaohua Zang, Gang Bao, Xiaojing Ye, and Haomin Zhou.  Weak adversarial networks for high-dimensional partial differential equations.  J. Comput. Phys., 411:109409, 14, 2020.  URL: https://doi.org/10.1016/j.jcp.2020.109409, doi:10.1016/j.jcp.2020.109409.
  * [37] Mengmeng Zhang, Qianxiao Li, and Jijun Liu.  On stability and regularization for data-driven solution of parabolic inverse source problems.  J. Comput. Phys., 474:Paper No. 111769, 20, 2023.  URL: https://doi.org/10.1016/j.jcp.2022.111769, doi:10.1016/j.jcp.2022.111769.
  * [38] Mengmeng Zhang and Jijun Liu.  Recovery of non-smooth radiative coefficient from nonlocal observation by diffusion system.  J. Inverse Ill-Posed Probl., 28(3):389–410, 2020.  URL: https://doi.org/10.1515/jiip-2019-0029, doi:10.1515/jiip-2019-0029.
  * [39] Zhengqi Zhang, Zhidong Zhang, and Zhi Zhou.  Identification of potential in diffusion equations from terminal observation: analysis and discrete approximation.  SIAM J. Numer. Anal., 60(5):2834–2865, 2022.  URL: https://doi.org/10.1137/21M1446708, doi:10.1137/21M1446708.