\(^\ast \)This work was supported by the project "The Development of Advanced and Applicative Research Competencies in the Logic of STEAM + Health” /POCU/993/6/13/153310, project co-financed by the European Social Fund through The Romanian Operational Programme Human Capital 2014-2020.

\(^\S \)Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy, Cluj-Napoca, Romania, e-mail: mihai.nechita@ictp.acad.ro, orcid.org/0000-0002-6448-912X

\(^\bullet \)Department of Mathematics, Babeș-Bolyai University, Cluj-Napoca, Romania

1 Introduction

Physics-informed neural networks [ 16 ] have recently emerged as an effective way of approximating ill-posed/inverse problems for partial differential equations (PDEs), which are challenging to solve numerically due their inherent instability [ 11 ] . One of their main advantages is the versatility with which they can solve both forward and inverse problems. As reported in [ 9 ] for well-posed boundary value problems, PINNs are not able, however, to outperform finite elements methods for low dimensional forward problems (2d and 3d).

In this paper we explore PINNs for an important class of inverse problems, namely the unique continuation (data assimilation) problem in which measurements are available in a subset of the domain and no boundary conditions are given. The goal is to find a solution that extends the data subject to the PDE considered. Numerical methods for such problems typically consider regularization at the continuous level (Tikhonov or quasi-reversibility). More recently, other kinds of methods have been proposed that make use of regularization at the discrete level in the framework of stabilized finite element methods (FEM). In this latter approach, conditional stability estimates can be used to prove error bounds and convergence in terms of the degree of ill-posedness and the approximation order, see e.g. [ 4 , 5 , 15 ] . Conditional stability estimates have also been employed in [ 14 ] to obtain bounds for the generalization error of PINNs solving unique continuation problems for Poisson, Stokes, heat and wave equations.

We are interested in the ill-posed unique continuation problem for the Helmholtz equation, whose stability properties (in terms of the frequency) depend on the geometry of the data set relative to the geometry of the target set where the solution is reconstructed [ 6 ] . We consider PINNs for approximating this problem and we first give a bound in section 3 on the generalization error that takes into account the frequency dependence. We then explore numerically in section 4 how different frequencies and geometric configuration impact the effectiveness of PINNs. We see that PINNs give good approximations when the frequency is low and the target set is inside the convex hull of the data set (robust stability in terms of frequency as discussed in section 2.1). This includes results where data perturbations are present. When the target set is outside of the convex hull of the data set (stability highly sensitive to frequency) or when the frequency is large we observe poor results with vanilla PINNs. This indicates that more sophisticated techniques are needed to obtain PINNs that are robust in frequency for inverse problems subject to the Helmholtz equation.

2 Unique continuation for the Helmholtz equation

Let \(\Omega \subset \mathbb R^2\) be a bounded domain (open and connected), and let \(\omega \subset \Omega \) be an open and connected subset. We consider the ill-posed unique continuation problem for the Helmholtz equation: find \(u\in H^1(\Omega )\) such that

with wave number \(k{\gt}0\), source term \(f\in L^2(\Omega )\) and datum \(g\in H^1(\omega )\). In other words, partial measurements are given in a subset \(\omega \subset \Omega \) of the domain for the solution of the Helmholtz equation, while no boundary conditions are prescribed. If problem ?? has a solution, then its uniqueness is guaranteed by the unique continuation principle for elliptic operators. We will assume that the function \(g\) is the restriction on \(\omega \) of a solution to the Helmholtz equation with source term \(f\).

Problem ?? is ill-posed in the sense of Hadamard: there is no uniform stability with respect to the data \(f\) and \(g\), see e.g. [ 1 ] or [ 11 ] . A small perturbation in data can lead to a large change in the solution. Nonetheless, assuming an additional a priori bound, the solution can be bounded by the data in the following way.

2.1 Conditional stability

We denote by \(B\subset \Omega \) a target set containing \(\omega \) such that \(B\setminus \bar{\omega } \subset \Omega \), i.e. \(B\setminus \omega \) does not touch the boundary of \(\Omega \). A standard result for quantitative unique continuation of elliptic operators, see e.g. [ 1 ] , states that there exist constants \(C_{st}(k){\gt}0\) and \(\alpha \in (0,1)\) such that the following Hölder stability estimate holds

for any \(u\in H^1(\Omega )\) satisfying ??, where the stability constant \(C_{st}(k)\) depends on the frequency \(k\). The exponent \(\alpha \in (0,1)\) encodes the degree of ill-posedness for the continuation problem: as \(\alpha {\lt}1\) decreases the Hölder stability deteriorates, while \(\alpha =1\) would give a well-posed problem with Lipschitz stability. Both \(C_{st}\) and \(\alpha \) depend on the geometric configuration in a nontrivial way, and we give details below regarding \(C_{st}\).

Frequency dependence for the stability constant. Data set \(\omega \) (dark grey) and target set \(B\) (light grey). \(\Omega \) is the whole square.

An important aspect for the stability of this ill-posed Helmholtz problem is the dependence of the stability constant \(C_{st}(k)\) on the frequency \(k\). For example, when there is a straight line that intersects \(B\) but not \(\bar{\omega }\), as in figure 1, it was proven in [ 6 , Example 4 ] that for any \(N\in \mathbb N\), \(C_{st}(k) \le k^N\) cannot hold uniformly in \(k\); this means that the stability constant in 2 grows superpolynomially in the frequency. Also, for three-ball inequalities (where \(\omega , B, \Omega \) are concentric balls) it was recently shown in [ 2 ] that in the maximum norm \(C_{st}(k)\) grows exponentially in \(k\) and this dependence is optimal.

Bounds that have a different behavior with respect to the frequency can be obtained under a convexity condition of the target domain \(B\) relative to the data set \(\omega \), essentially that \(B\) is included in the convex hull of \(\omega \), as for example in figure 2. Such a condition was first considered in [ 10 ] , where it was shown that the stability of the solution in the \(L^2\)-norm can actually improve in a certain sense as the frequency increases. In this vein, similar results that make use of a convexity condition were proven in [ 6 ] for a particular geometric setting prototypical for continuation inside the convex hull of \(\omega \). It was shown in [ 6 , Corollary 2 ] that there exist constants \(C{\gt}0\) and \(\alpha \in (0,1)\) such that

for any \(u\in H^1(\Omega )\) satisfying ??. Moreover, the norms in which data is measured can be weaken [ 6 , Corollary 3 and Lemma 2 ] : there exist constants \(C{\gt}0\) and \(\alpha \in (0,1]\) such that

for any \(u\in H^1(\Omega )\) satisfying ??. Note that the bound is robust in the \(L^2\)-norm, while the frequency dependence is linear for the \(H^1\)-seminorm. Apart from the good dependence on the frequency, this kind of estimate is particularly suitable for numerical analysis since it can be directly applied to the error equation.

3 Physics-informed neural networks (PINNs)

Let us first recall the setup of a feed-forward neural network, which is based on successively applying linear and nonlinear transformations to the inputs. Let \(\mathcal{N}^L(\mathbf{x}): \mathbb {R}^{d_{\text{in}}} \to \mathbb {R}^{d_{\text{out}}}\) be an \(L\)-layer neural network with \((L-1)\) hidden layers and \(N_\ell \) neurons in the \(\ell \)-th layer (\(N_0 = d_{\text{in}}\), \(N_L = d_{\text{out}}\)). Let us denote the weight matrix and bias vector in the \(\ell \)-th layer by \(\bm@general \boldmath \m@ne \mv@bold \bm@command {W}^\ell \in \mathbb {R}^{N_\ell \times N_{\ell -1}}\) and \(\mathbf{b}^\ell \in \mathbb {R}^{N_\ell }\), respectively. For a nonlinear activation function \(\sigma \) applied element wise, the feed-forward neural network is given by:

We denote the parameters of the network (weights and biases) by

and the parameter set by \(\Theta \), hence

The output of a feed-forward neural network,

depends on the tuning parameter \(\theta \). Training the network means using an optimization algorithm to find the parameters \(\theta \in \Theta \) that minimize a certain loss function \(\mathcal{L}_\theta \).

Let us now briefly describe physics-informed neural networks (PINNs) as introduced in [ 16 ] , based on feed-forward neural networks. The idea is to consider a loss function with two components: the residual of the differential equation and the data residual. For typical boundary value problems, the data residual contains boundary conditions, initial conditions, and any other available measurements of the solution. In our case, since no boundary conditions are given, the data residual will only involve the interior partial measurements.

To solve the ill-posed unique continuation problem ?? with PINNs we consider the residual of the Helmholtz equation

and the data residual

and aim to minimize the two residuals over the admissible set of tuning parameters \(\Theta \) by considering the loss function

To approximate the integrals in the setting of a neural network, we introduce the following collocation points: \(\mathcal{T}_\Omega := \{ \textbf{x}_i^\Omega \} _{i=1}^{|\mathcal{T}_\Omega |}\) for the whole domain \(\Omega \) and \(\mathcal{T}_\omega := \{ \textbf{x}_i^\omega \} _{i=1}^{|\mathcal{T}_\omega |}\) for the data region \(\omega \subset \Omega \) in which observations of the solution are given as \(\{ g_i(\textbf{x}^\omega )\} _{i=1}^{|\mathcal{T}_\omega |}\). The loss function for the PINN is hence defined as

with the PDE residual

and the data residual

By training the neural network we aim to

We will denote by

the PINN solution obtained by this algorithm, see figure 3 for a sketch.

As the loss function is highly non-linear and non-convex, we will use the standard approach of minimizing by gradient-based methods such as ADAM [ 12 ] or L-BFGS [ 7 ] , more details are given in section 4.

3.1 Estimating the generalization error

Such PINNs for unique continuation problems have been considered in [ 14 ] where the authors prove for the first time estimates for the generalization error by employing conditional stability estimates; details of the analysis and numerical examples were given for the Poisson, Stokes, heat and wave equations. In brief, one uses conditional stability estimates to bound the error in a target domain (generalization error) by the residuals, which are then bounded in terms of the training error and the number of training samples. We now adapt the analysis in [ 14 , Theorem 2.4 ] to the case of the Helmholtz equation by using the stability results discussed in section 2.1.

in a weak sense. Applying the conditional stability estimate ?? we obtain

The conclusion follows by assuming quadrature rules with approximation bounds as in [ 14 , Eqs (2.9) and (2.11) ] , where for some \(\tau , \tau _d {\gt} 0\) one has that

with constants \(C_q = C_q\left(\| \mathring {\mathcal{R}}_{\theta ^{\ast }}\| _{C^{n-2}(\Omega )}\right)\) and \(C_{qd} = C_{qd}\left(\| \mathring {\mathcal{R}}_{\theta ^{\ast },d}\| _{C^{n}(\omega )}\right)\).

4 Numerical experiments

We present numerical experiments for the Helmholtz unique continuation problem ?? solved with physics-informed neural networks (PINNs) described in section 3 and implemented using the open-source library DeepXDE [ 13 ] . The test case we will focus on has been considered as a benchmark test for this problem solved with primal-dual stabilized finite element methods: with a conforming discretization in [ 6 ] and with a high-order hybridized discontinuous Galerkin method in [ 3 ] . It represents the Helmholtz version of the classical Hadamard example for ill-posed elliptic equations.

Let \(n \in \mathbb {N}\) and consider the Cauchy problem

whose solution for \(n{\gt}k\) is given by

for \(n=k\) by \(u(x,y) = \sin (kx) y\), and for \(n{\lt}k\) by

For such Hadamard-type solutions, we consider the interior datum \(g=u|_\omega \) and study two geometric configurations of \(\omega \) and \(B\): one with frequency-robust stability bounds and one exponentially sensitive to the frequency, as discussed in section 2.1, namely

similar to figure 2, and

sketched in figure 1.

To assess the effect of increasing the frequency, we will take exact solutions ?? with \(n=5,k=1\) and \(n=7,k=5\), both having \(\sqrt{n^2-k^2} = \sqrt{24}\) .

Hyper-parameters. The PINNs described in section 3 need to be configured with the following hyper-parameters: number of hidden layers \(L-1\) (depth), number of neurons in each hidden layer \(\ell \) (width), learning rate \(\lambda \), activation function \(\sigma \). In order to find good configurations for these, we draw on previous numerical experiments for forward and inverse problems presented in the DeepXDE tutorials [ 13 ] and the numerical experiments for unique continuation subject to the Poisson equation in [ 14 ] . We also validate the choices by comparing them with the results given by the hyper-parameter optimization in [ 8 ] , which uses Gaussian processes-based Bayesian optimization. For this we consider the search space \(L-1\in [4,50],N_\ell \in [20,150],\lambda \in [10^{-4}, 10^{-2}],\sigma \in \{ \sin ,\tanh \} \). In these ways we obtain some optimal or near-optimal configurations to which the results presented below correspond. From different numerical experiments, it turns out that a good choice for the activation function is \(\sigma = \sin \) and for the learning rate \(\lambda = 10^{-3}\), which will be fixed from this point onwards.

Training. We sample \(N=N_\Omega =N_\omega \) points on Cartesian grids from the PDE domain \(\Omega \) and the measurement domain \(\omega \). We train the model for 50000 iterations with the ADAM optimizer [ 12 ] and then we train again with L-BFGS [ 7 ] , following [ 13 ] . Since PINNs approximations are obtained by non-convex minimization which in general does not have a unique solution, the method might converge to different solutions depending on the network’s initial values. We use the standard strategy of training the PINNs with different random starting values for the optimizer and average the errors over 30 retrainings. We note that the smallest training error is considerably better.

Results. We first show in figure 4 the exact solution \(u\) for the Hadamard-type function ?? with \(n=7,k=5\), and the PINN solution \(u^\ast \) when considering a network generated with \(N=400\) training points, \(L-1=4\) hidden layers and \(N_\ell =24\) neurons in each layer. As discussed above, the activation function is \(\sigma = \sin \) and the learning rate \(\lambda = 10^{-3}\).

We study the efficienty of PINNs by considering the \(L^2\) and \(H^1\) relative percentage errors for the number of training points \(N=20^2,40^2,80^2\) in the geometric configuration ?? with robust bounds in theorem 2. table 1 shows the errors in the target domain \(B\), while table 2 shows the global errors. In both tables we observe similar behaviours for the error. Moreover, we notice that even for very few training points the generalization errors are small (around 3%). Note that the reported errors correspond to averages over 30 retrainings. The generalization errors corresponding to the network with the smallest training error are typically smaller by a factor of \(3\) to \(5\). The training error in these examples is in between \(10^{-3}\) and \(10^{-2}\). The training of such a network is very fast, on a single Quadro GPU it takes about 1-2 minutes.

We also test the performance of the method when the measurements are perturbed with random values sampled from the normal distribution with zero mean and standard deviation 0.01 (representing approximately 2.5% of the average of the exact solution ?? with \(n=5,k=1\) considered in that example). The results given in table 3 show that the numerical approximation is robust to such perturbations in data.

In table 4 we keep the same (convex) geometric configuration but increase the frequency to \(k=5\), with \(n=7\) such that \(\sqrt{n^2-k^2} = \sqrt{24}\) as before. We observe that the approximation improves as the frequency increases: the errors decrease compared to table 2. This is probably an artefact of the particular conditions of the experiment, but we note that such a surprising phenomenon for inverse Helmholtz problems has been previously noticed and theoretically discussed in [ 10 ] . However, when increasing the wave number \(k{\gt}10\) we report that the PINNs no longer provide a good approximation to the solution.

We now consider continuation outside the convex hull of the data set in ??. We observe in table 5 that even though the wave number is small \(k=1\), the PINNs fail to find a good approximation of the solution. This could be caused by a bad search space for hyper-parameter optimization or it could be an inherent limitation of vanilla PINNs.