Solving the inverse potential problem in the parabolic equation by the
  deep neural networks method
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.
Mengmeng Zhang, Zhidong Zhang
# Solving the inverse potential problem in the parabolic equation by the deep neural networks method

###### 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.

### Mathematical model.

The following parabolic system is considered in this work:

\[\left\{\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}\right. \tag{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

\[u(x,T):=\varphi(x),\quad x\in\Omega. \tag{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

\[\|\varphi^{\delta}-\varphi\|_{L^{\infty}(\Omega)}\leq\delta. \tag{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.

### 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^{2,1}_{p}(\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 Holder 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].

### 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.

### 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

\[s_{\xi}(z):=W_{K}l_{K-1}\circ\cdots\circ l_{1}(z)+b_{K}, \tag{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

\[\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}\).

### 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.

Figure 1: The fully connected neural networks.
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

\[\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. \tag{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_

\[\|\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_

\[|\varphi(x)-\varphi(y)| \leq C_{\varphi}|y-x|,\] \[|\Delta\varphi(x)-\Delta\varphi(y)| \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_

\[\|\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,

\[I_{2} =\Delta\Big{[}\int_{R^{d}}\rho_{\epsilon}(|x-y|)\ (\varphi(y)- \varphi^{\delta}(y))\ dy\Big{]}\] \[=\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.

### 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

\[\left\{\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}\right. \tag{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

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

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

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

\[\|u\|_{W_{2}^{2,1}(\Omega_{T})}=\sqrt{\sum_{|\alpha|\leq 2}\|D^{\alpha}u\|_{L^{2 }(\Omega_{T})}^{2}+\|u_{t}\|_{L^{2}(\Omega_{T})}^{2}},\]

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_

\[\|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 \(\|\overline{q}\|_{L^{\infty}(\Omega)}\), \(\nu\), \(\Omega\) and \(T\)._

### 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

\[\left\{\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}\right. \tag{3.3}\]

with

\[w(x,T)=0,\quad x\in\Omega. \tag{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

\[\left\{\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}\right.\]

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.

### 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 \[\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 \[\mathcal{R}_{sb,\theta}(x,t):=u_{\theta}(x,t)-b(x,t),\quad(x,t)\in\partial \Omega_{T}.\]
* Initial status residual \[\mathcal{R}_{tb,\theta}(x):=u_{\theta}(x,0)-u_{0}(x),\quad x\in\Omega.\]
* Data residual \[\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:

\[\begin{split} J_{\lambda}(\theta,\eta)=&\|q_{ \eta}\mathcal{R}_{d,\theta}\|_{L^{2}(\Omega)}^{2}+\|\Delta\mathcal{R}_{d, \theta}\|_{L^{2}(\Omega)}^{2}+\lambda\|\mathcal{R}_{int,\theta,\eta}\|_{H^{1 }(0,T;L^{2}(\Omega))}^{2}\\ &+\|\mathcal{R}_{tb,\theta}\|_{L^{2}(\Omega)}^{2}+\|q_{\eta} \mathcal{R}_{tb,\theta}\|_{L^{2}(\Omega)}^{2}+\|\Delta\mathcal{R}_{tb,\theta} \|_{L^{2}(\Omega)}^{2}+\|\mathcal{R}_{sb,\theta}\|_{H^{2}(0,T;L^{2}(\partial \Omega))}^{2},\end{split} \tag{4.2}\]

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_

\[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} \tag{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:

\[\mathcal{S}_{d} :=\left\{(x_{n},T):x_{n}\in\Omega,\quad n=1,2,\cdots,N_{d}\right\},\] \[\mathcal{S}_{int} :=\left\{(\widetilde{x}_{n},\widetilde{t}_{n}):(\widetilde{x}_{n},\widetilde{t}_{n})\in\Omega_{T},\quad n=1,2,\cdots,N_{int}\right\},\] \[\mathcal{S}_{tb} :=\left\{(\widetilde{x}_{n},0):\widetilde{x}_{n}\in\Omega,\quad n =1,2,\cdots,N_{tb}\right\},\] \[\mathcal{S}_{sb} :=\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

\[J_{\lambda}^{N}(\theta,\eta) =\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}\] \[\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}\] \[\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}\] \[\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}, \tag{4.4}\]

where the coefficients

\[\omega_{n}^{d,k},\ \omega_{n}^{int,k},\ \omega_{n}^{tb,j},\ \omega_{n}^{sb,j}, \ k=0,1,\ j=0,1,2\]

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

\[|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\}, \tag{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

\[\begin{cases}\mathcal{E}_{G,q}:=\|q-q^{*}\|_{L^{2}(\Omega)}\,,\\ \mathcal{E}_{G,u}:=\|u-u^{*}\|_{C([0,T];L^{2}(\Omega))}\,,\end{cases} \tag{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 \[\begin{cases}\mathcal{E}_{T,d,0}:=\left(\sum_{n=1}^{N_{d}}\omega_{j}^{d,0}\,|q _{\eta}\left(x_{n}\right)\mathcal{R}_{d,\theta^{*}}\left(x_{n}\right)|^{2} \right)^{\frac{1}{2}},\\ \mathcal{E}_{T,d,1}:=\left(\sum_{n=1}^{N_{d}}\omega_{j}^{d,1}\,|\Delta\mathcal{ R}_{d,\theta^{*}}\left(x_{n}\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 \[\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 \[\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_{n}( \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 \[\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)

### 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_

\[\begin{split}\mathcal{E}_{G,q}&\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{)},\\ \mathcal{E}_{G,u}&\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{)},\end{split} \tag{4.11}\]

_where_

\[N =\min\left\{N_{d},N_{int},N_{sb},N_{tb}\right\},\] \[\alpha =\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_

\[C_{qd}=C_{qd}(\|\mathcal{L}^{*}\mathcal{R}_{d,\theta^{*}}\|_{C( \Omega)}), C_{q,0}=C_{q,0}(\|\mathcal{R}_{int,\theta^{*},\eta^{*}}\|_{C(\Omega_{ T})}),\] \[C_{q,1}=C_{q,1}(\|\partial_{t}\mathcal{R}_{int,\theta^{*},\eta^ {*}}\|_{C(\Omega_{T})}), C_{qs,0}=C_{qs,0}(\|\mathcal{R}_{sb,\theta^{*}}\|_{C( \partial\Omega_{T})}),\] \[C_{qs,1}=C_{qs,1}(\|\partial_{t}\mathcal{R}_{sb,\theta^{*}}\|_{ C(\partial\Omega_{T})}), C_{qs,2}=C_{qs,2}(\|\partial_{t}^{2}\mathcal{R}_{sb,\theta^{*}}\|_{C( \partial\Omega_{T})}),\] \[C_{qt,0}=C_{qt,0}(\|\mathcal{R}_{tb,\theta^{*}}\|_{C(\Omega)}), 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

\[\left\{\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}\right.\]

with the final condition

\[\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

\[\left\{\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}\right. \tag{4.12}\]

with

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

and

\[\left\{\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}\right.\]

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

\[\|q^{*}-q\|_{L^{2}(\Omega)} \tag{4.14}\] \[\leq C\|\mathcal{L}^{*}\hat{u}_{1}(\cdot,T)\|_{L^{2}(\Omega)}\] \[=C\|\mathcal{L}^{*}\mathcal{R}_{d,\theta^{*}}-\mathcal{L}^{*}( \varphi-G_{\epsilon}\varphi^{\delta})-\mathcal{L}^{*}\hat{u}_{2}(\cdot,T)\|_{ L^{2}(\Omega)}\] \[=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)}\] \[\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}). \tag{4.15}\]

Also, we have

\[|(\varphi-G_{\epsilon}\varphi^{\delta})| \leq\int_{|x-y|\leq\epsilon}\rho_{\epsilon}(|x-y|)\ |\varphi(x)-\varphi^{\delta}(y)|\ dy\] \[\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\] \[\leq C\epsilon+\delta.\]
Thus, there holds

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

By straightforward computations, we have

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

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

\[\left\{\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}\right. \tag{4.18}\]

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

\[\|\partial_{t}\hat{u}_{2}\|_{L^{\infty}(0,T;L^{2}(\Omega))} =\|w\|_{L^{\infty}(0,T;L^{2}(\Omega))} \tag{4.19}\] \[\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))}\] \[\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

\[\mathcal{E}_{G,q} =\|q^{*}-q\|_{L^{2}(\Omega)}\] \[\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)}\] \[\qquad+\|\partial_{t}\mathcal{R}_{sb,\theta^{*}}\|_{H^{1}(0,T;L^ {2}(\partial\Omega))}+\epsilon(1+\delta\epsilon^{-3})+\epsilon+\delta\big{)}\] \[\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)}\] \[\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

\[\mathcal{E}_{G,q} =\|q^{*}-q\|_{L^{2}(\Omega)} \tag{4.20}\] \[\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)}\] \[\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,

\[\mathcal{E}_{G,u} =\|\hat{u}\|_{C([0,T];L^{2}(\Omega))}\] \[\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{)}\] \[\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))}\] \[\qquad+\|\mathcal{R}_{tb,\theta^{*}}\|_{L^{2}(\Omega)}\big{)}\] \[\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)}\] \[\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.

### 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^{*}}\) and \(q_{\eta^{*}}\), which serves as the approximation to the solution of the inverse problem.

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

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}.\]

Figure 2: The schematic of the deep neural networks for solving the inverse potential problem.
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:

1. 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.
2. 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.
3. 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.

### 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

\[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

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

and in our experiments the noisy data is set as

\[\varphi^{\delta}(x,y):=\varphi(x,y)+\delta\cdot(2\;\text{rand}(\text{shape}( \varphi(x,y)))-1), \tag{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

\[\mathcal{T}:=\{(t_{k},x_{i},y_{j}):k,i,j=0,1,\cdots,49\}\subset\Omega_{T}. \tag{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 5 shows the exact solution of the pretential 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.
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

Figure 4: The time series relative error (test) for \(u\) for different noisy data.

Figure 5: The exact potential function \(q\).

Figure 3: The training loss (left), the relative error for \(q\) (middle), the relative error for \(u\) (right) after logarithmic re-scaling.
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.

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\).
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

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

The exact measurement will be

\[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.

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.

\begin{table}
\begin{tabular}{l c c c} \hline \hline  & \(Re_{q}\) & \(Re_{u}\) & \(Re_{\Delta u}\) \\ \hline I & \(0.5011\%\pm 0.0102\%\) & \(0.0215\%\pm 0.0049\%\) & \(0.2550\%\pm 0.0169\%\) \\ \hline II & \(12.4545\%\pm 4.8016\%\) & \(0.1148\%\pm 0.0133\%\) & \(2.3468\%\pm 0.8133\%\) \\ \hline \hline \end{tabular}
\end{table}
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).
Figure 8: The training loss (left), the relative error for \(q\) (middle), the relative error for \(u\) (right) after logarithmic re-scaling.

Figure 10: The reconstruction of potential function \(q\) (upper) and corresponding absolute pointwise error \(|q-q^{*}|\) (bottom) for various noisy level measurement data.

Figure 9: The exact potential function \(q\).
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.

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).
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

\[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

\[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

\begin{table}
\begin{tabular}{c c c c} \hline \hline Error & \(\lambda=10^{-4}\) & \(\lambda=10^{-3}\) & \(\lambda=10^{-2}\) \\ \hline \(Re_{q}\) & \(2.5113\%\pm 1.1339\%\) & \(0.9942\%\pm 0.1121\%\) & \(0.7124\%\pm 0.0639\%\) \\ \hline \(Re_{u}\) & \(0.1374\%\pm 0.0652\%\) & \(0.1216\%\pm 0.0320\%\) & \(0.0560\%\pm 0.0204\%\) \\ \hline Error & \(\lambda=10^{-1}\) & \(\lambda=1\) & \(\lambda=10\) \\ \hline \(Re_{q}\) & \(0.8388\%\pm 0.2448\%\) & \(1.0809\%\pm 0.1741\%\) & \(2.2891\%\pm 1.3302\%\) \\ \hline \(Re_{u}\) & \(0.0331\%\pm 0.0094\%\) & \(0.0478\%\pm 0.0122\%\) & \(0.0874\%\pm 0.0110\%\) \\ \hline \hline \end{tabular}
\end{table}
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\).

\begin{table}
\begin{tabular}{c c c c c} \hline \hline \(NL\) & \(NN=10\) & \(NN=15\) & \(NN=20\) & \(NN=25\) \\ \hline
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\%\) \\ \hline
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\%\) \\ \hline \hline \end{tabular}
\end{table}
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}\).
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.

Finally, we conduct the experiments to evaluate the robustness of the proposed scheme in terms of network structure (number of hidden layers and per-layerurons). More specifically, we run Algorithm 1 with per-layer neuron numbers \(\text{NN}=5,10,20,25\) for fixed hidden layers

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.
\(\mathrm{NL}=6\) and with hidden layers \(\mathrm{NL}=4,6,8,10\) for fixed per-layer neuron numbers \(\mathrm{NN}=25\), respectively. The reconstruction errors are presented in Figure 14. It also seems that larger \(\mathrm{NL}\) and/or larger \(\mathrm{NN}\) neural networks yield lower reconstruction error. In this example, for fixed hidden layer \(\mathrm{NL}=6\), we test the per-layer neuron numbers \(\mathrm{NN}\geq 30\) and find that there will be bad reconstruction result, that is, the relative error for \(q\) with \(\mathrm{NL}=6,\mathrm{NN}=30\) is larger than the case \(\mathrm{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.

## 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.

Figure 14: The reconstruction error of \(q\) under various per-layer neuron numbers with fixed layer numbers \(\mathrm{NL}\)=6 (left), the reconstruction error of \(q\) under various layer numbers with fixed per-layer neuron number \(\mathrm{NN}\)=25 (right).
## 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).
