Finite element interpolated neural networks for solving forward and
  inverse problems
We propose a general framework for solving forward and inverse problems
constrained by partial differential equations, where we interpolate neural
networks onto finite element spaces to represent the (partial) unknowns. The
framework overcomes the challenges related to the imposition of boundary
conditions, the choice of collocation points in physics-informed neural
networks, and the integration of variational physics-informed neural networks.
A numerical experiment set confirms the framework's capability of handling
various forward and inverse problems. In particular, the trained neural network
generalises well for smooth problems, beating finite element solutions by some
orders of magnitude. We finally propose an effective one-loop solver with an
initial data fitting step (to obtain a cheap initialisation) to solve inverse
problems.
Santiago Badia, Wei Li, Alberto F. MartÃ­n
# Finite element interpolated neural networks for solving forward and inverse problems

###### Abstract.

We propose a general framework for solving forward and inverse problems constrained by partial differential equations, where we interpolate neural networks onto finite element spaces to represent the (partial) unknowns. The framework overcomes the challenges related to the imposition of boundary conditions, the choice of collocation points in physics-informed neural networks, and the integration of variational physics-informed neural networks. A numerical experiment set confirms the framework's capability of handling various forward and inverse problems. In particular, the trained neural network generalises well for smooth problems, beating finite element solutions by some orders of magnitude. We finally propose an effective one-loop solver with an initial data fitting step (to obtain a cheap initialisation) to solve inverse problems.

Key words and phrases:neural networks, PINNs, finite elements, PDE approximation, inverse problems

## 1. Introduction

Many problems in science and engineering are modelled by _low-dimensional_ (e.g., 2 or 3 space dimensions plus time) partial differential equations (PDEs). Since PDEs can rarely be solved analytically, their solution is often approximated using numerical methods, among which the finite element method (FEM) has been proven to be effective and efficient for a broad range of problems. The FEM enjoys a very solid mathematical foundation [1]. For many decades, advanced discretisations have been proposed, e.g., preserving physical structure [2], and optimal (non)linear solvers that can efficiently exploit large-scale supercomputers have been designed [3, 4, 5].

Grid-based numerical discretisations can readily handle forward problems. In a forward problem, all the data required for the PDE model to be well-posed is provided (geometry, boundary conditions and physical parameters), and the goal is to determine the state of the model. In an inverse problem setting, however, the model parameters are not fully known, but one can obtain some observations, typically noisy and/or partial, of the model state. Inverse problem solvers combine the partially known model and the observations to infer the information which is missing to complete the model. Inverse problems can be modelled using PDE-constrained minimisation [6].

Traditional numerical approximations for low-dimensional PDEs, like FEM, are linear. Finite element (FE) spaces are finite-dimensional vector spaces in which one seeks for the best approximation in some specific measure. As a result, the method/grid is not adapted to local features (e.g., sharp gradients or discontinuities) and convergence can be slow for problems that exhibit multiple scales. Although adaptive FE methods can efficiently handle this complexity, they add an additional loop to the simulation workflow (the mark and refine loop) and problem-specific robust error estimates have to be designed [7].

When the FEM is used to solve PDE-constrained inverse problems, the unknown model parameters are usually described using FE-like spaces, even though neural network (NN) representations have recently been proposed [8, 9]. The loss function accounts for the misfit term between the observation and the state, which in turns depends on the unknown model parameters. The gradient of the loss function with respect to the unknown model parameters requires a chain rule that involves the solution of the forward problem. An efficient implementation of this gradient relies on the adjoint method [10]. Inverse problem solvers add an additional loop to the simulation workflow, which involves the solution of the full forward problem and the adjoint of its linearisation at each iteration. There is usually a burden of computation cost in the first stages of the adjoint method, when full forward problems are solved despite being far from the desired solution.
The tremendous success of NNs in data science has motivated many researchers to explore their application in PDE approximation. Physics-informed NNs (PINNs) have been proposed in [11] to solve forward and inverse problems. A NN approximates the PDE solution, while the loss function evaluates the strong PDE residual on a set of randomly selected collocation points. NNs can also be combined with a weak statement of the PDE (see, e.g., Deep Ritz Method [12] or variational PINNs (VPINNs) [13]). NNs have some very interesting properties that make them perfectly suited for the approximation of forward and inverse PDE problems. First, NNs are genuinely _nonlinear approximations_, similar to, e.g., free-knot B-splines [14]. The solution is sought in a nonlinear manifold in the parameter space, which automatically adapts to the specific problem along the training process.1 Unlike FE bases, NNs are also perfectly suited (and originally designed) for data fitting. The NN parameters usually have a _global_ effect on the overall solution. As a result, one can design solvers for PDE-constrained inverse problems in which both the unknown model parameters and state variables are learnt along the same training process [16]. The loss function includes the data misfit and a penalised PDE residual term. State and unknown model parameters are not explicitly linked by the forward problem, and thus no forward problems are involved in each iteration of the optimisation loop. As a result, one can use NNs to design adaptive forward and inverse problems with a one-loop solver.

Footnote 1: It is illustrative to observe how e.g. linear regions in NNs with ReLU activation functions adapt to the solution being approximated [15]. The decomposition of the physical domain into linear regions is a polytopal conforming mesh.

Despite all these efforts, PINNs and related methods have not been able to outperform traditional numerical schemes for low dimensional PDEs; see, e.g., the study in [17]. There are some (intertwined) reasons for this. First, nonlinear approximability comes at the price of non-convex optimisation at the training process. Currently, non-convex optimisation algorithms for NN approximation of PDEs are costly and unreliable, especially when the PDE solutions contain multi-scale features or shocks [18, 19]. As a result, despite the enhanced expressivity of NNs, this improvement is overshadowed by poor and costly training. Second, the integration of the PDE residual terms is not exact, and the error in the integration is either unbounded or not taken into account. Poor integration leads to poor convergence to the desired solution (due to a wrong cost functional) and one can find examples for which the optimal solution is spurious [20]. In [15], the authors propose adaptive quadratures for NN in low dimensions that are proven to be more accurate that standard Monte Carlo, especially for sharp features. Finally, the usual PDE residual norms being used in the loss function are ill-posed at the continuous level in general. It is well-known that such a variational crime has negative effects in the convergence of iterative solvers for FE discretisations [21], and it will also hinder the non-convex optimisation at the training process. These issues prevent a solid mathematical foundation of these methods, and strong assumptions are required to prove partial error estimates [22, 23].

Additionally, NNs have not been designed to strongly satisfy Dirichlet boundary conditions. Thus, the loss function must include penalty terms that account for the boundary conditions, which adds an additional constraint to the minimisation and has a very negative effect on the training [24]. Such imposition of boundary conditions is not consistent, and Nitsche's method comes with the risk of ending up with an ill-posed formulation.2 Recently, some authors have proposed to multiply the NN with a distance function that vanishes on the Dirichlet boundary [25]. However, this arguably complicates the geometrical discretisation step compared to grid-based methods. The computation of such distance functions is complex in general geometries and has only been used for quite simple cases in 2D. Furthermore, it is unclear how to use this approach for non-homogeneous boundary conditions, which require a lifting of the Dirichlet values inside the domain. In comparison, (unstructured) mesh generation is a mature field and many mesh generators are available [26]. Unfitted FEs have become robust and general schemes that can handle complex geometries on Cartesian meshes [27]. With a mesh, the definition of the lifting is trivial, e.g., one can use a FE offset function.

Footnote 2: The coefficient in Nitsche’s method must be _large enough_ for stability, which can be mathematically quantified in FEM using inverse inequalities. However, NNs nature do not enjoy inverse inequalities; gradients can be arbitrarily large, and can only be indirectly bounded via regularisation.

In [28], the authors propose to interpolate NNs onto a FE space to overcome the integration issues and design a well-posed PDE-residual loss functional. The imposition of the boundary conditions is handled via distance functions to overcome the training problems related to boundary conditions [29]. The authors compare the solution of the interpolated NN (a FE function) with different standard PINN formulations. Despite the fact that the solution belongs to a fixed FE space (and cannot exploit nonlinear approximation, compared to the other PINN strategies), the results are superior in general. This technique, coined interpolated VPINNs (IVPINNs), has been applied to forward coercive grad-conforming PDEs on rectangular domains.
Unlike other PINNs, a priori error bounds have been obtained [28], even though suboptimal compared to the FEM solution. One can argue what is the benefit of getting sub-optimal FE solutions (measured in the energy norm) using a far more expensive non-convex optimisation solver. However, IVPINNs shed light on the negative impact that integration, residual definition, imposition of boundary conditions, lack of well-posedness, and training have on a straightforward approximation of PDEs using NNs.

In this work, we build upon IVPINNs ideas. However, instead of enforcing boundary conditions at the NN, we propose to strongly impose the boundary conditions at the FE space level. This allows us to readily handle complex geometries without the need to define, e.g., distance functions. To distinguish the two approaches, we coin the proposed method FE interpolated NNs (FEINNs). Besides, we explore the benefits of considering the trained NN (instead of the FE interpolation) as the final solution of the problem, i.e., evaluate how the trained NN _generalises_. We also discuss different PDE residual norms and suggest to use Riesz preconditioning techniques to end up with a well-posed formulation in the continuous limit. We perform a numerical analysis of the method, and prove that the proposed formulation can recover (at least) the _optimal_ FE bounds. Next, we apply these techniques to inverse problems, using a one-loop algorithm, as it is customary in PINNs. We exploit the excellent properties of NNs to fit data. We propose a first step in which we get a state initial guess by data fitting. In a second step, we learn the unknown model parameters by PDE-residual minimisation for a fixed state. The previous steps provide an initialisation for a third fully coupled step with a mixed data-PDE residual cost function.

We carry out a comprehensive set of numerical experiments for forward problems. We check that expressive enough NNs can return FE solutions for different polynomial orders. For smooth problems, the generalisation results for the trained NNs are striking. The solution obtained with the non-interpolated FEIN solution can be orders of magnitude more accurate than the FE solution on the same mesh, while IVPINNs do not generalise that well. The definition of the residual norm (and its preconditioned version) can have a tremendous impact in the convergence of the minimisation algorithm. Finally, we test the proposed algorithm for inverse problems. Unlike standard inverse solvers for grid-based methods, we can solve inverse problems with effective and cheap initialisation and one-loop algorithms, even without any kind of regularisation terms.

The outline of the article is the following. Sec. 2 states the model elliptic problem that we tackle, its FE discretisation, the NN architecture, and the proposed loss functions in the FEINM discretisation. Sec. 3 proves that the interpolation of an expressive enough NN recovers the FE solution. In Sec. 4, the proposed discretisation is applied to inverse problems, by defining a suitable loss function that includes data misfit and a multi-step minimisation algorithm. Sec. 5 describes the implementation of the methods and Sec. 6 presents the numerical experiments on several forward and inverse problems. Finally, Sec. 7 draws conclusions and lists potential directions for further research.

## 2. Forward problem discretisation using neural networks

### Continuous problem

In this work, we aim to approximate elliptic PDEs using a weak (variational) setting. As a model problem, we consider a convection-diffusion-reaction equation, even though the proposed methodology can readily be applied to other coercive problems. The problem reads: find \(u\in H^{1}(\Omega)\) such that

\[-\boldsymbol{\nabla}\cdot(\boldsymbol{\kappa}\boldsymbol{u})+(\boldsymbol{ \beta}\cdot\boldsymbol{\nabla})u+\sigma u=f\quad\text{in }\Omega,\quad u=g\quad\text{on }\Gamma_{D}, \quad\kappa\boldsymbol{n}\cdot\boldsymbol{\nabla}u=\eta\quad\text{on }\Gamma_{N}, \tag{1}\]

where \(\Omega\subset\mathbb{R}^{d}\) is a Lipschitz polyhedral domain, \(\Gamma_{D}\) and \(\Gamma_{N}\) are a partition of its boundary such that \(\text{meas}(\Gamma_{D})>0\), \(\kappa\), \(\sigma\in L^{\infty}(\Omega)\), \(\boldsymbol{\beta}\in W^{1,\infty}(\Omega)^{d}\) such that \(\sigma-\boldsymbol{\nabla}\cdot\boldsymbol{\beta}>0\) and \(\boldsymbol{\beta}\cdot\boldsymbol{n}_{|\Gamma_{N}}\geq 0\), \(f\in H^{-1}(\Omega)\), \(g\in H^{1/2}(\Gamma_{D})\), and \(\eta\in H^{-1/2}(\Gamma_{N})\).

Consider the space \(U\doteq H^{1}(\Omega)\), \(\tilde{U}\doteq H^{1}_{0,\Gamma_{D}}(\Omega)\doteq\left\{v\in U\ :\ v_{|\Gamma_{D}}=0\right\}\), a continuous lifting \(\bar{u}\in U\) of the Dirichlet boundary condition (i.e., \(\bar{u}=g\) on \(\Gamma_{D}\)), and the forms

\[a(u,v)=\int_{\Omega}\kappa\boldsymbol{\nabla}u\cdot\boldsymbol{\nabla}v+( \boldsymbol{\beta}\cdot\boldsymbol{\nabla})uv+\sigma uv,\quad\ell(v)=\int_{ \Omega}fv+\int_{\Gamma_{N}}\eta v.\]

(We use the symbol \(\overline{\cdot}\) to denote trial functions and spaces with zero traces.) The variational form of the problem reads: find \(u=\bar{u}+\bar{u}\) where

\[\bar{u}\in\tilde{U}\ :\ a(\bar{u},v)=\ell(v)-a(\bar{u},v),\quad\forall v\in\tilde{U}. \tag{2}\]

In this setting, the problem with a non-homogeneous Dirichlet boundary condition is transformed into a homogeneous one via the lifting and a modification of the right-hand side (RHS). The well-posedness of the
problem relies on the coercivity and continuity of the forms:

\[a(u,u)\geq\gamma\|u\|_{U}^{2},\quad a(u,v)\leq\xi\|u\|_{U}\|v\|_{U},\quad\ell(v) \leq\chi\|v\|_{U}.\]

Below, we will make use of the PDE residual

\[\mathcal{R}(\tilde{u})\,\doteq\,\ell(\cdot)-a(\tilde{u}+\tilde{u},\cdot)\in \tilde{U}^{\prime}. \tag{3}\]

### Finite element approximation

Next, we consider a family of conforming shape-regular partitions \(\{\mathcal{T}_{h}\}_{h>0}\) of \(\Omega\) such that their intersection with \(\Gamma_{N}\) and \(\Gamma_{D}\) is also a partition of these lower-dimensional manifolds; \(h\) represents a characteristic mesh size. On such partitions, we can define a trial FE space \(U_{h}\subset U\) of order \(k_{U}\) and the subspace \(\tilde{U}_{h}\doteq U_{h}\cap H^{1}_{0,\Gamma_{D}}(\Omega)\) of FE functions with zero trace.

We define a FE interpolant \(\pi_{h}:C^{0}\to U_{h}\) obtained by evaluation of the degrees of freedom (DoFs) of \(U_{h}\). In this work, we consider grad-conforming Lagrangian (nodal) spaces (and thus composed of piece-wise continuous polynomials), and DoFs are pointwise evaluations at the Lagrangian nodes. Analogously, we define the interpolant \(\tilde{\pi}_{h}\) onto \(\tilde{U}_{h}\). We can pick a FE lifting \(\tilde{u}_{h}\in U_{h}\) such that \(\tilde{u}_{h}=\pi_{h}(g)\) on \(\Gamma_{D}\). (The interpolant is restricted to \(\Gamma_{D}\) and could be, e.g., a Scott-Zhang interpolant if \(g\) is non-smooth.) Usually in FEM, \(\tilde{u}_{h}\) is extended by zero on the interior.

Using the Galerkin method, the test space is defined as \(V_{h}\doteq\tilde{U}_{h}\), and let \(k_{V}\) be its order. Following [28], we also explore Petrov-Galerkin discretisations. To this end, we consider \(k_{V}\) such that \(s=k_{U}/k_{V}\in\mathbb{N}\), and a family of partitions \(\mathcal{T}_{h/s}\) obtained after \(s\) levels of uniform refinement of \(\mathcal{T}_{h}\). In this case, we choose \(V_{h}\) to be the FE space of order \(k_{V}\) on \(\mathcal{T}_{h/s}\) with zero traces on \(\Gamma_{D}\). We note that the dimension of \(\tilde{U}_{h}\) and \(V_{h}\) are identical. In this work, we only consider \(k_{V}=1\), i.e., a _linearised_ test FE space. The well-posedness of the Petrov-Galerkin discretisation is determined by the discrete inf-sup condition:

\[\inf_{u_{h}\in\tilde{U}_{h}}\sup_{v_{h}\in V_{h}}\frac{a(u_{h},v_{h})}{\|u_{h} \|_{U}\|v_{h}\|_{U}}\geq\beta>0.\]

In both cases, the problem can be stated as: find \(u_{h}=\tilde{u}_{h}+\tilde{u}_{h}\) where

\[\tilde{u}_{h}\in\tilde{U}_{h}\ :\ a(\tilde{u}_{h},v_{h})=\ell(v_{h})-a( \tilde{u}_{h},v_{h}),\quad\forall v_{h}\in V_{h}. \tag{4}\]

We represent with \(\mathcal{R}_{h}\) the restriction \(\mathcal{R}|_{U_{h}\times V_{h}}\). Given \(u_{h}\in U_{h}\), \(\mathcal{R}(u_{h})\in V_{h}^{\prime}\). \(V_{h}^{\prime}\) is isomorphic to \(\mathbb{R}^{N}\), where \(N\) is the dimension of \(V_{h}\) (and \(\tilde{U}_{h}\)). This representation depends on the basis chosen to span \(V_{h}\).

### Neural networks

We consider a fully-connected, feed-forward NN, obtained by the composition of affine maps and nonlinear activation functions. The network architecture is represented by a tuple \((n_{0},\ldots n_{L})\in\mathbb{N}^{(L+1)}\), where \(L\) is the number of layers and \(n_{k}\) is the number of neurons on layer \(1\leq k\leq L\). We take \(n_{0}=d\) and, for scalar-valued PDEs, we have \(n_{L}=1\). In this work, we use \(n_{1}=n_{2}=...=n_{L-1}=n\), i.e. all the hidden layers have an equal number of neurons \(n\).

At each layer \(1\leq k\leq L\), we represent with \(\mathbf{\Theta}_{k}:\mathbb{R}^{n_{L-1}}\rightarrow\mathbb{R}^{n_{k}}\) the affine map at layer \(k\), defined by \(\mathbf{\Theta}_{k}\mathbf{x}=\mathbf{W}_{k}\mathbf{x}+\mathbf{b}_{k}\) for some weight matrix \(\mathbf{W}_{k}\in\mathbb{R}^{n_{k}\times n_{k-1}}\) and bias vector \(\mathbf{b}_{k}\in\mathbb{R}^{n_{k}}\). The activation function \(\rho:\mathbb{R}\rightarrow\mathbb{R}\) is applied element-wise after every affine map except for the last one. Given these definitions, the network is a parametrizable function \(\mathcal{N}(\mathbf{\theta}):\mathbb{R}^{d}\rightarrow\mathbb{R}\) defined as:

\[\mathcal{N}(\mathbf{\theta})=\mathbf{\Theta}_{L}\circ\rho\circ\mathbf{\Theta}_{L-1}\circ \ldots\circ\rho\circ\mathbf{\Theta}_{1}, \tag{5}\]

where \(\mathbf{\theta}\) stands for the collection of all the trainable parameters \(\mathbf{W}_{k}\) and \(\mathbf{b}_{k}\) of the network. Although the activation functions could be different at each layer or even trainable, we apply the same, fixed activation function everywhere. However, we note that the proposed methodology is not restricted to this specific NN architecture. In this work, we denote the NN architecture with \(\mathcal{N}\) and a realisation of the NN with \(\mathcal{N}(\mathbf{\theta})\).

### Finite element interpolated neural networks

In this work, we propose the following discretisation of (2), which combines the NN architecture in (5) and the FE problem in (4). Let us consider a norm \(\|\cdot\|_{Y}\) for the discrete residual (choices for this norm are discussed below). We aim to find

\[u_{\mathcal{N}}\in\arg\min_{w_{\mathcal{N}}\in\mathcal{N}}\mathcal{L}(w_{ \mathcal{N}}),\qquad\mathcal{L}(w_{\mathcal{N}})\,\doteq\,\|\mathcal{R}_{h}( \tilde{\pi}_{h}(w_{\mathcal{N}}))\|_{Y}. \tag{6}\]

The computation of \(u_{\mathcal{N}}\) involves a non-convex optimisation problem (due to the nonlinear dependence of \(u_{\mathcal{N}}\) on \(\mathbf{\theta}\)). We prove in the next section that the \(\tilde{\pi}_{h}(u_{\mathcal{N}})\) equal to the FE solution is a global minimum of this functional.

In this method, the NN is _free_ on \(\Gamma_{D}\), the imposition of the Dirichlet boundary conditions relies on a FE lifting \(\tilde{u}_{h}\) and the interpolation \(\tilde{\pi}_{h}\) onto \(\tilde{U}_{h}\) applied to the NN (thus vanishing on \(\Gamma_{D}\)). Conceptually, the
proposed method trains a NN _pinned_ on the DoFs of the FE space \(\tilde{U}_{h}\), with a loss function that measures the FE residual of the interpolated NN for a given norm. The motivation behind the proposed method is to eliminate the Dirichlet boundary condition penalty term in standard PINNs and related methods [11, 12], while avoiding enforcing the conditions at the NN level (see, e.g., [25] for PINNs and [29] for VPINNs). It also solves the issues related to Monte Carlo integration [20] and avoids the need to use adaptive quadratures [15]. Using standard element-wise integration rules, the integrals in \(\mathcal{R}_{h}\) can be exactly computed (or, at least, its error can be properly quantified for non-polynomial physical parameters and forcing terms). Moreover, in the current setting, we can consider different alternatives for the residual norm and better understand the deficiencies and variational crimes related to standard choices.

### Loss function

As discussed above, the loss function involves the norm of the FE residual. The residual is isomorphic to the vector \([\mathbf{r}_{h}(w_{h})]_{i}=\left\langle\mathcal{R}(w_{h}),\varphi^{i} \right\rangle\doteq\mathcal{R}(w_{h})(\varphi^{i})\), where \(\{\varphi^{i}\}_{i=1}^{N}\) are the FE shape functions that span the test space \(V_{h}\). As a result, we can consider the loss function:

\[\mathcal{L}(u_{\mathcal{N}})=\|\mathbf{r}_{h}(\tilde{\pi}_{h}(u_{\mathcal{N}} ))\|_{\ell^{2}}.\]

This is the standard choice (possibly squared) in the methods proposed so far in the literature that rely on variational formulations [13, 28, 30]. However, as it is well-known in the FE setting, this quantity is ill-posed in the limit \(h\downarrow 0\)[21]. At the continuous level, the norm of \(\mathcal{R}(u)\) is not defined.

If the problem is smooth enough and \(\mathcal{R}\) is well-defined on \(L^{2}(\Omega)\) functions, we can define its \(L^{2}(\Omega)\) projection onto \(V_{h}\) as follows:

\[\mathcal{M}_{h}^{-1}\mathcal{R}_{h}(w_{h})\in V_{h}\ :\ \int_{\Omega} \mathcal{M}_{h}^{-1}\mathcal{R}_{h}(w_{h})v_{h}=\mathcal{R}_{h}(w_{h})(v_{h}), \quad\forall v_{h}\in V_{h}.\]

Next, one can define the cost function

\[\mathcal{L}(u_{N})=\|\mathcal{M}_{h}^{-1}\mathcal{R}_{h}(\tilde{\pi}_{h}(u_{N }))\|_{L^{2}(\Omega)},\]

which, for quasi-uniform meshes, is equivalent (up to a constant) to the scaling of the Euclidean norm, i.e., \(h^{d}\|\mathbf{r}_{h}(\tilde{\pi}_{h}(u_{N}))\|_{\ell^{2}}\). However, for non-smooth solutions, the \(L^{2}\) norm of the residual still does not make sense at the continuous level, and thus, the convergence must deteriorate as \(h\downarrow 0\). One can instead define a discrete Riesz projector \(\mathcal{B}_{h}^{-1}:V_{h}^{\prime}\to V_{h}\) such that

\[\mathcal{B}_{h}^{-1}\mathcal{R}_{h}(w_{h})\in V_{h}\ :\ \left(\mathcal{B}_{h}^{-1} \mathcal{R}_{h}(w_{h}),v_{h}\right)_{U}=\mathcal{R}(w_{h})(v_{h}),\quad\forall v _{h}\in V_{h}.\]

For the model case proposed herein, \(\|\cdot\|_{U}\) is the \(H^{1}\) or \(H^{1}_{0,\Gamma_{D}}\)-norm and \(\mathcal{B}_{h}^{-1}\) is the inverse of the discrete Laplacian. Then, one can consider the cost function:

\[\mathcal{L}(u_{\mathcal{N}})=\|\mathcal{B}_{h}^{-1}\mathcal{R}_{h}(\tilde{\pi }_{h}(u_{\mathcal{N}}))\|_{L^{2}(\Omega)}, \tag{7}\]

or

\[\mathcal{L}(u_{\mathcal{N}})=\|\mathcal{B}_{h}^{-1}\mathcal{R}_{h}(\tilde{\pi }_{h}(u_{\mathcal{N}}))\|_{H^{1}(\Omega)}.\]

These cost functions are well-defined in the limit \(h\downarrow 0\). In practice, one can replace \(\mathcal{B}_{h}^{-1}\) by any spectrally equivalent approximation in order to reduce computational demands. For example, in the numerical experiments section, we consider several cycles of a geometric multigrid (GMG) preconditioner.

## 3. Analysis

In this section, we first show that the proposed loss functions are differentiable. Next, we show that the interpolation of the NN architecture can return any FE function in a given FE space. Combining these two results, we observe that there exists a global minimum of the FEIN problem in (6) such that its interpolation is the solution of the FE problem (4).

**Proposition 3.1**.: _The loss function is differentiable for \(\mathcal{C}^{0}\) activation functions._

Proof.: Using the chain rule, we observe that

\[\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}\boldsymbol{\theta}}=\frac{\mathrm{d} \mathcal{L}}{\mathrm{d}\mathbf{r}_{h}}\frac{\mathrm{d}\mathbf{r}_{h}}{\mathrm{ d}\mathbf{u}_{h}}\frac{\mathrm{d}\mathbf{u}_{h}}{\mathrm{d}\boldsymbol{\theta}},\]

for \(\mathbf{u}_{h}\) being the DoFs of the FE space \(U_{h}\). The first derivative in the RHS simply involves the squared root of a quadratic functional. The second derivative is the standard Jacobian of the FE problem. The third
derivative is the vector of derivatives of the NN at the nodes of \(U_{h}\), which is well-defined for \(C^{0}\) activation functions. As a result, \(\mathcal{L}\) is differentiable. 

Consequently, one can use gradient-based minimisation techniques. We note that this is not the case when the NN is evaluated without FE interpolation. For instance, refer to [15] for a simple example that shows ReLU activation functions cannot be used for PDE approximation using PINNs and related methods. In PINNs, one must compute \(\boldsymbol{\nabla}_{\boldsymbol{\theta}}\boldsymbol{\nabla}_{\boldsymbol{x} }N\), which poses additional smoothness requirements on the activation function. However, in the proposed methodology (as in [28]), the spatial derivatives are computed by the interpolated function, not the NN, and thus not affected by this constraint. For simplicity, we prove the result for the ReLU activation function.

**Proposition 3.2**.: _Let \(U_{h}\) be a FE space on a mesh \(\mathcal{T}_{h}\) in \(\mathbb{R}^{d}\) with DoFs equal to \(N\simeq h^{1/d}\). Let \(\mathcal{N}\) be a neural network architecture with 3 layers, \((3dN,dN,N)\) neurons per layer, and a ReLU activation function. For any \(u_{h}\in U_{h}\), there exists a choice of the NN parameters \(\boldsymbol{\theta}\) such that \(\pi_{h}(u_{\mathcal{N}})=u_{h}\)._

Proof.: At each node \(\boldsymbol{n}\in\mathcal{T}_{h}\), one can define a box \(\mathcal{B}_{\boldsymbol{n}}\) centred at \(\boldsymbol{n}\) that only contains this node of the mesh. Let us consider first the 1D case. For ReLU activation functions, one can readily define a hat function with support in \([0,1]\) as follows. First, we consider

\[f_{1}(x)=2x,\quad f_{2}(x)=4x-2,\quad f_{3}(x)=2x-2.\]

One can check that \(f=\rho(f_{1})-\rho(f_{2})+\rho(f_{3})\) is a hat function with value 1 at \(x=1/2\) and support in \([0,1]\). One can readily consider a scaling and translation to get \(\operatorname{supp}(f)\subset\mathcal{B}_{\boldsymbol{n}}\). This way, assuming one has \(3N\) neurons in the first layer and \(N\) neurons in the second layer, one can emulate the 1D FE basis in the second layer.

In 2D, one can create the 1D functions for both \(x\) and \(y\) directions. It requires \(6N\) neurons in the first layer and \(2N\) neurons in the second layer. Thus, for each node, we have two hat functions, namely \(b_{1}\) and \(b_{2}\), that depend on \(x\) and \(y\), respectively. Now, in a third layer with \(N\) neurons, we can compute \(\rho(b_{1}+b_{2}-1)\) at each node. We can generalise this construction to an arbitrary dimension \(d\). We need \(3dN\) neurons in the first layer to create the 1D functions in all directions. The hat functions are created in a second layer with \(dN\) neurons. The final functions are combined as \(\rho(\sum_{i=1}^{d}b_{i}-d+1)\). We note that, by the construction of \(b_{i}\), these functions have value one in the corresponding node and their support is contained in the corresponding box.

In the last layer, we end up with a set of functions \(\psi_{i}\) that are equal to 1 on one node and zero on the rest. Besides, the FE function can also be expressed as \(u_{h}=\sum_{i=1}^{N}u^{i}\boldsymbol{\varphi}^{i}(\boldsymbol{x})\) and \(\pi_{h}(u_{\mathcal{N}})=\sum_{i=1}^{N}u_{\mathcal{N}}(\boldsymbol{n}_{i}) \boldsymbol{\varphi}^{i}(\boldsymbol{x})\). Linearly combining the last layer functions with the DoF values \(\{u^{i}\}_{i=1}^{N}\) we construct a NN realisation that proves the proposition. 

**Remark 3.3**.: _For other activation functions like tanh or sigmoid, it is not possible to construct localised functions with compact support as in the proof above. However, one can consider a piecewise polynomial approximation of these activation functions (e.g., using B-splines) with this property [31]. Then, one can use a similar construction as in ReLU._

We note that this construction can be further optimised by exploiting the structure of the underlying FE mesh \(\mathcal{T}_{h}\). For instance, for a structured mesh of a square with \(n\) parts per direction (\(N=n^{d}\)), only \(3n\) neurons are needed. We can exploit the fact that many nodes share the same coordinates in some directions. For the same reason, only \(dn\) neurons are required in the second layer. On the other hand, for more than 3 layers, the computations can be arranged among neurons/layers in different ways. For simplicity, in the proposition, we consider a worst-case scenario situation (no nodes share coordinate components and we only consider the arrangement in the proposition statement).

**Proposition 3.4**.: _Let us assume that the FE problem (4) is well-posed and admits a unique solution \(\tilde{u}_{h}\). The FEIN problem (6) admits a global minimiser \(u_{\mathcal{N}}\) such that \(\tilde{\pi}_{h}(u_{\mathcal{N}})=\tilde{u}_{h}\)._

Proof.: First, we note that the loss function differentiable (by Prop. 3.1) and positive. Besides, from the statement of the problem and Prop. 3.2, one can readily check that there exists a \(u_{\mathcal{N}}\) such that \(\mathcal{R}_{h}(\tilde{\pi}_{h}(u_{\mathcal{N}}))=\mathcal{R}_{h}(\tilde{u}_{h })=0\) and thus \(\mathcal{L}(u_{\mathcal{N}})=0\), i.e. \(u_{\mathcal{N}}\) is a global minimum of the cost function. 

As a result, the FEINN method can exhibit optimal convergence rates (the ones of FEM), provided the NN is expressive enough compared to the FE space. In Sec. 6.1, we experimentally observe this behaviour. This
analysis is different from the one in [28], which, using a completely different approach, proves sub-optimal results in a different setting. The numerical experiments in [28] and in Sec. 6.1 show that IVPINNs can also recover optimal convergence rates. In fact, the results above can straightforwardly be extended to IVPINNs. The sub-optimality in [28] is related to the choice of the residual norm, the \(\ell^{2}\) norm of the residual vector. Sharper estimates could likely be obtained with the new residual norms suggested in Sec. 2.5.

## 4. Inverse problem discretisation using neural networks

In this section, we consider a PDE-constrained inverse problem that combines observations of the state variable \(u\) and a partially known model (1). Let us represent with \(\mathbf{\Lambda}\) the collection of unknown model parameters. It can include the physical coefficients, forcing terms and Dirichlet and Neumann boundary values. We parametrise \(\mathbf{\Lambda}\) with one or several NNs, e.g., as the ones proposed for the state variable in Sec. 2.4, which will be represented with \(\mathbf{\Lambda}_{N}\). Again, \(n_{0}=d\), while \(n_{L}\) depends on whether the unknown model parameter of the specific problem is a scalar-valued (\(n_{L}=1\)), vector-valued (\(n_{L}=d\)) or tensor-valued (\(n_{L}=d^{2}\)) field.

Let us denote with \(\mathcal{R}(\mathbf{\Lambda},u)\) the PDE residual in (3), where we make explicit its dependence with respect to the unknown model parameters (idem for \(\mathcal{R}_{h}\)). For integration purposes, we consider the interpolation of the model parameters onto FE spaces, which we represent with \(\boldsymbol{\pi}_{h}(\mathbf{\Lambda}_{N})\). The discrete model parameter FE spaces can in general be different to \(U_{h}\) (just as their infinite-dimensional counterpart spaces might be different to \(U\)) and do not require imposition of boundary conditions. Besides, the interpolation can be restricted to different boundary regions for Dirichlet and Neumann values. If we consider a discontinuous nodal FE space with nodes on the quadrature points of the Gaussian quadrature being used for integration (as in the numerical experiments), the interpolated and non-interpolated methods are equivalent. Thus, the interpolant simply accounts for the integration error being committed when integrating the NNs for the unknown model parameters.

Let us consider a measurement operator \(\mathcal{D}:U\to\mathbb{R}^{M}\) and the corresponding vector of observations \(\mathbf{d}\in\mathbb{R}^{M}\). The loss function for the inverse problem must contain the standard data misfit term and a term that accounts for the PDE residual. The method is understood as a (PDE-)constrained minimisation problem. As a result, the PDE residual is weighted by a (dynamically adapted) penalty coefficient. We consider the loss functional:

\[\mathcal{L}(\mathbf{\Lambda},u)\doteq\|\mathbf{d}-\mathcal{D}(u)\|_{\ell^{2} }+\alpha\|\mathcal{R}_{h}(\boldsymbol{\pi}_{h}(\mathbf{\Lambda}),\tilde{\pi} _{h}(u))\|_{Y}, \tag{8}\]

for any of the choices of the residual norm discussed above and \(\alpha\in\mathbb{R}^{+}\) is a penalty coefficient for the weak imposition of the PDE constraint. The inverse problem reads:

\[u_{N},\mathbf{\Lambda}_{N}\in\operatorname*{arg\,min}_{w_{N},\Xi_{N}\in \mathbf{\Lambda}_{N}\times\mathbf{\Lambda}_{N}}\mathcal{L}(\boldsymbol{\Xi}_ {N},w_{N}). \tag{9}\]

We refer to [32] for an application of penalty methods to inverse problems. However, their approach is more akin to the adjoint method, where they eliminate the state. We note that our approach is a _one-loop_ minimisation algorithm, i.e., one can minimise for both the state and unknown model parameters at the same time. This differs from adjoint methods, in which the loss function and the minimisation is in terms of \(\mathbf{\Lambda}\) only, but the state \(u(\mathbf{\Lambda})\) is constrained to be the solution of the (discrete) PDE at each iterate of \(\mathbf{\Lambda}\).

To alleviate the challenges associated with the training of the loss function described in (9) and enhance the robustness of our method, we propose the following algorithm. The motivation behind its design is to exploit the excellent properties of NNs for data fitting. First, we train the state NN with the observations. Next, we train the unknown model parameters NNs with the PDE residual, but freeze the state variable to the value obtained in the previous step. These steps are computationally lightweight because they do not involve differential operators in the training processes. These two initial steps are finally used as initialisation for the one-loop minimisation in (9). We summarise the algorithm below:

* Step 1 (Data fitting): Train the state neural network to fit the observed data, using standard NN initialisation: \[u_{N}^{0}=\operatorname*{arg\,min}_{w_{N}\in\mathbf{\Lambda}_{N}}\|\mathbf{d} -\mathcal{D}(w_{N})\|_{\ell^{2}}.\]
* Step 2 (Unknown model parameters initialisation): Train the model parameter NNs with the PDE residual for the fixed state \(u_{N}^{0}\) computed in Step 1, using standard NN initialisation: \[\mathbf{\Lambda}_{N}^{0}=\operatorname*{arg\,min}_{\Xi_{N}\in\mathbf{\Lambda}_ {N}}\|\mathcal{R}_{h}(\boldsymbol{\pi}_{h}(\mathbf{\Xi}_{N}),\tilde{\pi}_{h}( u_{N}^{0}))\|_{Y}.\]
* Step 3 (Fully coupled minimisation): Train both the state and model parameter NNs the full loss function (8), starting from \(u_{\mathcal{N}}^{0}\) and \(\mathbf{\Lambda}_{\mathcal{N}}^{0}\).

It is important to point out that the three-step training process is facilitated by the incorporation of NNs. Our attempts to train a FE function in the data step have not been successful, especially when the number of observations is much smaller than the DoFs of the FE space. This is attributed to the local support of FE functions, which limits the adjustment of the values of the free nodes that are directly influenced by the observations. In contrast, NNs with their global support, allow for parameter tuning across the entire domain in the data step.

## 5. Implementation

We rewrite (8) in the following algebraic form

\[\mathcal{L}(\mathbf{\theta}_{\lambda},\mathbf{\theta}_{u})=\left\|\mathbf{e}(\mathbf{ u}_{h}(\mathbf{\theta}_{u}))\right\|_{\ell^{2}}+\alpha\left\|\mathbf{r}_{h}( \mathbf{u}_{h}(\mathbf{\theta}_{u}),\mathbf{\lambda}_{h}(\mathbf{\theta}_{\lambda}))\right\| _{\ell^{1}}, \tag{10}\]

where \(\mathbf{e}\doteq\mathbf{d}-\mathcal{D}_{h}\mathbf{u}_{h}\) is the data misfit error, \(\mathbf{r}_{h}\) is the variational residual vector, \(\mathbf{u}_{h}\), \(\mathbf{\lambda}_{h}\) are the vectors of DoFs of \(\tilde{\mathbf{\pi}}_{h}(u_{\mathcal{N}}(\mathbf{\theta}_{u}))\) and \(\mathbf{\pi}_{h}(\mathbf{\Lambda}_{\mathcal{N}}(\mathbf{\theta}_{\lambda}))\) of the NN realisations \(u_{\mathcal{N}}(\mathbf{\theta}_{u})\) and \(\mathbf{\Lambda}_{\mathcal{N}}(\mathbf{\theta}_{\lambda})\) for the arrays of parameters \(\mathbf{\theta}_{u}\) and \(\mathbf{\theta}_{\lambda}\), respectively. We have chosen the \(\ell^{1}\) residual norm in (10) because it is the one we have used in the numerical tests for inverse problems in Sec. 6. However, the proposed implementation is general and can be easily adapted to other choices of residual norms proposed above.

We describe below an implementation of FEINNs using Julia packages, even though the proposed implementation is general. In Julia, we rely on the existing packages Flux.jl[33, 34] for the neural network part and Gridap.jl[35, 36] for the FEM part. We employ ChainRules.jl[37] to automatically propagate user-defined rules across the code.

To minimise the loss function (10) with gradient-based training algorithms, these gradients are required:

\[\frac{\partial\mathcal{L}}{\partial\mathbf{\theta}_{u}}=\left(\frac{\partial \mathcal{L}}{\partial\mathbf{r}_{h}}\frac{\partial\mathbf{r}_{h}}{\partial \mathbf{u}_{h}}+\frac{\partial\mathcal{L}}{\partial\mathbf{e}}\frac{\partial \mathbf{e}}{\partial\mathbf{u}_{h}}\right)\frac{\partial\mathbf{u}_{h}}{ \partial\mathbf{\theta}_{u}},\qquad\frac{\partial\mathcal{L}}{\partial\mathbf{\theta}_ {\lambda}}=\frac{\partial\mathcal{L}}{\partial\mathbf{r}_{h}}\frac{\partial \mathbf{r}_{h}}{\partial\mathbf{\lambda}_{h}}\frac{\partial\mathbf{\lambda}_{h}}{ \partial\mathbf{\theta}_{\lambda}}.\]

Existing chain rules in ChainRules.jl can readily handle \(\partial\mathcal{L}/\partial\mathbf{r}_{h}\) and \(\partial\mathcal{L}/\partial\mathbf{e}\). We need to define specific rules for the automatic differentiation of the following tasks:

* The interpolation of a NN onto a FE space in \(\partial\mathbf{u}_{h}/\partial\mathbf{\theta}_{u}\) and \(\partial\mathbf{\lambda}_{h}/\partial\mathbf{\theta}_{\lambda}\);
* The computation of the FE residual in \(\partial\mathbf{r}_{h}/\partial\mathbf{u}_{h}\) and \(\partial\mathbf{r}_{h}/\partial\mathbf{\lambda}_{h}\);
* The measurement operator \(\mathcal{D}\) on the FE state in \(\partial\mathbf{e}/\partial\mathbf{u}_{h}\).

It is important to highlight that we never explicitly construct the global Jacobian matrices in our implementation. To evaluate the gradient \(\partial\mathcal{L}/\partial\mathbf{\theta}_{\lambda}\), we utilise Gridap.jl to compute the Jacobian \(\partial\mathbf{r}_{h}/\partial\mathbf{\lambda}_{h}\) cell-wise (i.e., at each cell of \(\mathcal{T}_{h}\) separately), and restrict the vector \(\partial\mathcal{L}/\partial\mathbf{r}_{h}\) to each cell. By performing the vector Jacobian product (VJP) within each cell for \(\partial\mathcal{L}/\partial\mathbf{r}_{h}\) and \(\partial\mathbf{r}_{h}/\partial\mathbf{\lambda}_{h}\), we obtain the cell-wise vectors that can be assembled to form \(\partial\mathcal{L}/\partial\mathbf{\lambda}_{h}\). With the help of Flux.jl, we can calculate the gradient \(\partial\mathcal{L}/\partial\mathbf{\theta}_{\lambda}\) by performing the VJP for \(\partial\mathcal{L}/\partial\mathbf{\lambda}_{h}\) and \(\partial\mathbf{\lambda}_{h}/\partial\mathbf{\theta}_{\lambda}\), without explicitly constructing the Jacobian \(\partial\mathbf{\lambda}_{h}/\partial\mathbf{\theta}_{\lambda}\). This cell-wise approach recasts most of the floating point operations required to compute the gradients in terms of dense matrix-vector products. This results in a reduction of the computational times and memory requirements.

The gradient \(\partial\mathcal{L}/\partial\mathbf{\theta}_{u}\) has two contributions, corresponding to the FE residual and data misfit terms. The same process described above is applied to compute the former contribution. The contribution of the data misfit term involves the computation of \(\partial\mathcal{L}/\partial\mathbf{e}\)\(\partial\mathbf{e}/\partial\mathbf{u}_{h}\), which has not been discussed so far. In our implementation, it also follows an efficient cell-wise approach. In particular, we identify those cells with at least one observation point and, for these cells, we evaluate the cell shape functions at the observation points. This is nothing but the restriction of \(\partial\mathbf{e}/\partial\mathbf{u}_{h}\) to the observation points and DoFs of the cell. We then restrict the vector \(\partial\mathcal{L}/\partial\mathbf{e}\) to these cells, and compute the VJP among these vector and Jacobian restrictions. Finally, we assemble the resulting cell-wise vector contributions to obtain the data misfit global contribution vector to the vector \(\partial\mathcal{L}/\partial\mathbf{u}_{h}\).

Once all the rules for Jacobian computations are appropriately defined, ChainRules.jl seamlessly combine them, enabling smooth gradient computation during the training process.

Let us finish the section with a discussion about computational cost. In FEINNs and IVPINNs, one computes the spatial derivatives in the residual on FE functions in \(\partial\mathbf{r}_{h}/\partial\mathbf{u}_{h}\) and the derivatives of pointwise evaluations of the NN with respect to parameters in \(\partial\mathbf{u}_{h}/\partial\mathbf{\theta}\) separately. The expression of the polynomial derivatives is straightforward and the parameter differentiation is the one required in standard data fitting (and
thus, highly optimised in machine learning frameworks). On the contrary, in standard PINNs the residual is not evaluated with the projection \(\mathbf{u}_{h}\) but the NN itself. One must compute \(\partial\mathbf{r}_{h}/\partial\boldsymbol{\theta}\) directly. It involves _nested_ derivatives (in terms of parameters and input features) that are more expensive (and less common in data science).

## 6. Numerical experiments

### Forward problems

We use the standard \(L^{2}\) and \(H^{1}\) error norms to evaluate the precision of the approximation \(u^{id}\) for forward problems:

\[e_{L^{2}(\Omega)}(u^{id})=\left\|u-u^{id}\right\|_{L^{2}(\Omega)},\qquad e_{H^ {1}(\Omega)}(u^{id})=\left\|u-u^{id}\right\|_{H^{1}(\Omega)},\]

where \(u\) is the true state, \(\left\|\cdot\right\|_{L^{2}(\Omega)}=\sqrt{\int_{\Omega}\left|\cdot\right|^{2}}\), and \(\left\|\cdot\right\|_{H^{1}(\Omega)}=\sqrt{\int_{\Omega}\left|\cdot\right|^{2} +\left|\boldsymbol{\nabla}(\cdot)\right|^{2}}\). The integrals in these terms are evaluated with Gauss quadrature rule, and a sufficient number of quadrature points are used to guarantee accuracy. Note that, in the forward problem experiments, \(u^{id}\) can either be a NN or its interpolation onto a suitable FE space. We will specify which representation is used explicitly when necessary.

As for the experiments, we first compare FEINNs with IVPINNs by solving the forward convection-diffusion-reaction problem (1). Next, we shift to the Poisson equation, i.e., problem (1) with \(\boldsymbol{\beta}=\boldsymbol{0}\) and \(\sigma=0\), and analyse the impact of preconditioning on accelerating convergence during the training process. Finally, we showcase the effectiveness of FEINNs in complex geometries by solving a Poisson problem in a domain characterised by irregular shapes. It is worth noting that a comprehensive comparison in terms of computational cost and accuracy between IVPINNs, PINNs, and VIPINNs has already been conducted in [28]. In these experiments, the accuracy of IVPINNs is similar or better than the other PINNs being analysed for a given number of NN evaluations. The computational cost of IVPINNs is reported to be lower than standard PINN approaches, which is explained by the different cost of differentiation in each case, as explained in Sec. 5. As a result, we restrict ourselves to the comparison between FEINNs and IVPINNs and refer the reader to [28] for the relative merit of FEINNs over PINNs.

In all the experiments in this section, we adopt the NN architecture in [28], namely \(L=5\) layers, \(n=50\) neurons for each hidden layer, and \(\rho=\tanh\) as activation function, so that we can readily compare these results with the ones in [28] for standard PINNs. Besides, this choice strikes a good balance between the finest FE resolution being used and the NN expressivity. Indeed, we have experimentally observed that increasing the expressiveness of the NN (additional number of layers and/or neurons per layer) for the finest FE mesh being used in our experiments does not noticeably improve the results.

In addition, we employ Petrov-Galerkin discretisations, i.e., we use a linearised test space \(V_{h}\) as defined in Sec. 2.2. Unless otherwise specified, we adopt the \(\ell^{2}\) norm in the loss function (6). In all the experiments in this section and Sec. 6.2, we use the Glorot uniform method [38] for NN parameter initialisation and the BFGS optimiser in Optim.jl[39].3

Footnote 3: We have experimentally observed that L-BFGS is not as effective as BFGS for the problems considered in this paper.

#### 6.1.1. Convection-diffusion-reaction equation with a smooth solution

We replicate most of the experiment settings in [28, Convergence test #1], allowing the interested reader to check how other PINNs perform in similar experiments by looking at this reference. Specifically, the problem is defined on a square domain \(\Omega=[0,1]^{2}\), \(\Gamma_{D}\) are spanned by the left and right sides, and \(\Gamma_{N}\) by the the top and bottom ones. We choose the following analytical functions for the model parameters:

\[\kappa(x,y)=2+\sin(x+2y),\qquad\boldsymbol{\beta}(x,y)=\left[\sqrt{x-y^{2}+5},\ \sqrt{y-x^{2}+5}\right]^{\mathrm{T}},\qquad\sigma(x,y)=\mathrm{e}^{\frac{k}{3}- \frac{\gamma}{3}}+2,\]

and pick \(f\), \(g\) and \(\eta\) such that the exact solution is:

\[u(x,y)=\sin(3.2x(x-y))\cos(x+4.3y)+\sin(4.6(x+2y))\cos(2.6(y-2x)).\]

We discretise the domain using uniform meshes of quadrilateral elements of equal size.

It is crucial to emphasize that IVPINNs and FEINNs share a fundamental idea at their core: the interpolation of NNs (or their product with a function for IVPINNs) onto a corresponding FE space. The primary distinction lies in the approach used to impose the Dirichlet boundary condition. FEINNs rely on the trial FE space to enforce the boundary condition, using an interpolation that enforces zero trace on \(\Gamma_{D}\). IVPINNs, however, rely on an offset function \(\tilde{u}\) and a distance function \(\Phi\), where \(\tilde{u}\) is as smooth as \(u\) and satisfies the Dirichlet boundary condition and \(\Phi\in\tilde{U}\). The authors propose in [28] to train an auxiliary neural
network or use data transfinite interpolation to compute the lifting \(\bar{u}\) of the Dirichlet data \(g\). So, the true state can be expressed as \(\Phi\circ u_{\mathcal{N}}+\bar{u}\). IVPINNs interpolate this expression onto the FE space. The interpolated NN composition now belongs to \(\tilde{U}_{h}\) due to the property of \(\Phi\), and the full expression (approximately) satisfies the Dirichlet boundary condition because of the existence of \(\bar{u}\). The loss function of the method reads:

\[u_{\mathcal{N}}\in\arg\min_{w_{N}}\|\hat{\mathbf{r}}_{h}(\pi_{h}(\Phi\circ w_{ \mathcal{N}}+\bar{u}))\|,\quad[\hat{\mathbf{r}}_{h}]_{i}\doteq\ell(\varphi^{i} )-a(\pi_{h}(\Phi\circ w_{\mathcal{N}}+\bar{u}),\varphi^{i}),\]

where \(\{\varphi^{i}\}_{i=1}^{N}\) are the shape functions that span \(V_{h}\). In our numerical experiments, we have considered the training of an auxiliary NN to approximate \(\bar{u}\), but the results were not satisfactory. (Probably, because we are computing a function in \(\Omega\) with data on \(\Gamma_{D}\) only.) We have considered instead a discrete harmonic extension (i.e., a FE approximation of the Poisson problem with \(g\) on \(\Gamma_{D}\)) to approximate \(\bar{u}\). Since we consider a trivial square domain, the distance function \(\Phi\) can readily be defined as the product of the linear polynomials, i.e., \(\Phi(x,y)=x(1-x)\); note that \(\Gamma_{D}\) only includes the left and right sides of the squared domain.

In addition to evaluating the performance of FEINNs and IVPINNs, we also examine how the NNs generalise. For IVPINNs, we compute the error of the NN composition \(\Phi\circ u_{\mathcal{N}}+\bar{u}\), while, in the case of FEINNs, we compute the error of \(u_{\mathcal{N}}\) directly. We emphasise that this setting aligns with the principles of NN training: we train the NN with data in a set of points (the nodes of the mesh), and if the training is effective, we expect the NN to yield low error on the whole domain \(\bar{\Omega}\).

For the first experiment, we investigate the impact of mesh refinement on the approximation error. Keeping \(k_{U}=6\) fixed, we discretise the domain using a uniform mesh of quadrilaterals with different levels of refinement. To account for the impact of NN initialisation on both FEINNs and IVPINNs, we run 10 experiments with different initialisations for each mesh resolution. Fig. 0(a) and 0(b) illustrate the \(L^{2}\) errors and \(H^{1}\) errors, respectively, for the different methods versus mesh size. The curves labelled as "FEM" refer to the errors associated to the FEM solution, those labelled as "FEINN" and "IVPINN" to the errors of the interpolated NNs resulting from either method, and, finally, the label tag "(NN only)" is used to refer to the (generalisation) error associated to the NN itself (i.e., not to its FE interpolation). Due to the negligible variance in the errors of the interpolated NNs for both FEINNs and IVPINNs, we present the average error among those obtained for the 10 experiments. We also provide the slopes of the FEM convergence curves in Fig. 0(a) and 0(b). The computed slopes closely match the expected theoretical values, validating the FEM solution and the accuracy of the error computation. Based on the observations from Fig. 1, FEINNs not only generalise better compared to IVPINNs, they also have the potential to outperform FEM. This capability of FEINNs is not coincidental, as all errors associated to the NNs resulting from FEINNs, consistently remain below the FEM convergence curve. Additionally, we observe that as the mesh becomes finer, IVPINNs starts to struggle. While more training iterations may reduce the errors of IVPINNs, it is worth noting that the number of training iterations reaches the prescribed limit of 30,000 for the three finest mesh resolutions. It is also interesting to compare the distribution of errors among the NNs resulting from IVPINNs and FEINNs. We observe a high sensitivity of the errors to NN initialisation for IVPINNs, whereas the errors tend to cluster for FEINNs. We also observe that the \(L^{2}\) error of FEM gets closer to that of the non-interpolated NN resulting from FEINN as the mesh is refined. This behaviour is expected, since the FE mesh is being refined while the NN architecture is fixed. There is a point in which the NN is not expressive enough to represent the optimal FE solution and thus, Prop. 3.4 does not hold any more.

Since \(u\in\mathcal{C}^{\infty}(\bar{\Omega})\) in this problem, similar to FEM, we can also explore at which rate the error decays as we increase the polynomial order of the trial space (i.e., the NN interpolation space). We maintain a fixed mesh consisting of \(15\times 15\) quadrilaterals, and increase \(k_{U}\) from 1 up to 6. We perform 10 experiments for each order, with a different NN initialisation for each experiment. Fig. 1(a) and 1(b) depict the \(L^{2}\) and \(H^{1}\) errors, respectively, against \(k_{U}\). Once again, we observe that FEINNs have comparable performance to FEM, and more importantly, the non-interpolated NNs resulting from FEINNs have lower errors than FEM. In some cases, these can outperform FEM by more than two orders of magnitude. On the same mesh, the NN obtained with FEINNs is comparable to the FE solution obtained using between one and two orders more. Overall, IVPINNs demonstrate a comparable level of performance to FEM, with the exception occurring at \(k_{U}=6\). After 30,000 training iterations, it fails to achieve the same performance as FEM. Notably, the non-interpolated NN compositions from IVPINNs yield satisfactory results at lower orders, but as the order increases, they fail to reach the accuracy of FEM. The same comment about the expressivity limit of the NN architecture applies here. As we increase the order, the improvement of FEINN becomes less pronounced, since we are keeping fix the NN architecture.
In the sequel, we investigate how the computational cost and convergence rates of FEINNs and IVPINNs compare. To this end, we solve the same problem so far in this section by training FEINNs and IVPINNs for a fixed number of iterations, and then visualise at which rate the \(L^{2}\) and \(H^{1}\) errors decay with time. This is reported in Fig. 3. We used \(k_{U}=6\), a mesh consisting of \(15\times 15\) quadrilaterals (resulting in a problem with 8,099 DoFs), and a fixed number of 30,000 BFGS training iterations. Besides, we consider two different choices of the offset function to study its impact on the performance of IVPINNs. In particular, the curve labelled "(smooth offset)" in Fig. 3 denotes the same offset function used so far in this section (namely, a discrete harmonic extension, as per suggested in [28]), while the one labelled "(standard offset)" denotes the offset function that one naturally uses in FEM (and we also use here with FEINNs). We observe that both FEINNs and IVPINNs have roughly the same computational cost per iteration, around 0.023 seconds per iteration on a GeForce RTX 3090 GPU. Moreover, as shown in Fig. 1(a) and 1(b), IVPINNs converge consistently slower than FEINNs. Fig. 3 also illustrates that the choice of offset function greatly influences the convergence rate of IVPINNs. Indeed, IVPINN with the smooth offset function converges much faster

Figure 1. Convergence of errors with respect to the mesh size of the trial space for the forward convection-diffusion-reaction problem with a smooth solution.

Figure 2. Convergence of errors with respect to the order of trial bases for the forward convection-diffusion-reaction problem with a smooth solution.
than with the standard offset function, while for FEINNs, we readily obtain a faster convergence rate without the need for a special offset function. It is also worth noting that the authors in [28] observe that IVPINNs are less computationally expensive than PINNs and VPINNs, and thus for this problem, which indicates that FEINNs are also more efficient than the latter two methods.

#### 6.1.2. Convection-diffusion-reaction equation with a singular solution

The second problem we solve is still (1), but with a singular solution. We adopt most of the settings in [28, Convergence test #2]. The domain and boundaries are the same as those in Sec. 6.1.1. The coefficients are \(\kappa=1\), \(\boldsymbol{\beta}=[2,3]^{\mathsf{T}}\), and \(\sigma=4\). We pick \(f\), \(\eta\), and \(g\) such that the true state is, in polar coordinates,

\[u(r,\theta)=r^{\frac{2}{3}}\sin(\frac{2}{3}(\theta+\frac{\pi}{2})).\]

Since \(u\in H^{5/3-\epsilon}(\Omega)\) for any \(\epsilon>0\), the expected \(H^{1}\) error decay rate is around \(2/3\). Consequently, increasing \(k_{U}\) is unlikely to effectively reduce the error. Therefore, we keep \(k_{U}=2\), and focus our study on the impact of mesh refinement on error reduction. Fig. 4 depicts how \(L^{2}\) and \(H^{1}\) errors decay as we increase the mesh size. The errors of the non-interpolated NNs are not displayed in the plots, because their performance is relatively poor. This observation is consistent with previous findings in [28], which

Figure 4. Convergence of errors with respect to the mesh size of the trial space for the forward convection-diffusion-reaction problem with a singular solution.

Figure 3. Comparison among FEINNs and IVPINNs in terms of \(L^{2}\) and \(H^{1}\) errors versus computational time. Training was performed in both methods for a fixed number of 30,000 iterations.
highlight the inferior performance of PINNs and VPINNs compared to IVPINNs in this singular solution scenario. Fig 4 show how the \(L^{2}\) and \(H^{1}\) errors change as the mesh size changes. We obtain the expected error decay rate in Fig. 3(b). We conclude that both FEINNs and IVPINNs perform well in addressing this singular problem, and they successfully overcome the limitations in NNs in this particular situation.

#### 6.1.3. The effect of preconditioning on Poisson equation with a singular solution

In this experiment, we investigate whether preconditioning can effectively accelerate the training process, and examine the potential of leveraging widely used GMG preconditioners from FEM to aid in the training of FEINNs.

We only consider the \(L^{2}\)-norm of the preconditioned loss, i.e., (7). At a purely algebraic level, we can rewrite (7) as:

\[\mathcal{L}(\mathbf{\theta}_{u})=\left\|\mathbf{B}^{-1}\mathbf{A}\mathbf{u}_{h}( \mathbf{\theta}_{u})-\mathbf{B}^{-1}\mathbf{f}\right\|_{\ell^{2}}, \tag{11}\]

where \(\mathbf{B}\) is the preconditioner, \(\mathbf{A}\) is the coefficient matrix resulting from discretisation, \(\mathbf{\theta}_{u}\) are the parameters for \(u_{N}\), \(\mathbf{u}_{h}\) is the vector of DoFs of \(\bar{U}_{h}\), and \(\mathbf{f}\) is the RHS vector. (We note that, since the mesh being used is (quasi-)uniform, we can replace the \(L^{2}\)-norm by the Euclidean \(\ell^{2}\)-norm; they differ by a scaling.)

We consider three types of preconditioners. The first one is \(\mathbf{B}_{\text{inv}}=\mathbf{A}\). Plugged into (11), the loss becomes \(\left\|\mathbf{u}_{h}(\mathbf{\theta}_{u})-\mathbf{A}^{-1}\mathbf{f}\right\|_{ \ell^{2}}\). This loss resembles the loss in data fitting tasks, and it should theoretically be easier for NNs to minimise. We also consider another preconditioner \(\mathbf{B}_{\text{inv\_lin}}\), which is defined as the matrix resulting from discretisation with \(V_{h}\) as both trial and test FE spaces. Note that \(V_{h}\) is built out of a mesh resulting from the application of \(k_{U}\) levels of uniform refinement to the mesh associated to \(U_{h}\), since \(k_{V}=1\). The preconditioner \(\mathbf{B}_{\text{inv\_lin}}\) is computationally cheaper to invert than \(\mathbf{B}_{\text{inv}}\), since it is symmetric positive definite (SPD) and involves linear FEM bases only. The last (and cheapest to invert) one is, as mentioned before, a GMG preconditioner \(\mathbf{B}_{\text{GMG}}\) of \(\mathbf{B}_{\text{inv\_lin}}\).

We now change to the Poisson equation. The problem is defined on \(\Omega=[0,1]^{2}\), with \(\Gamma_{\text{D}}=\partial\Omega\) and \(\kappa=1\). Choose \(f\) and \(g\) such that the true state is the same as the singular \(u\) in Sec. 6.1.2. We divide the domain uniformly into \(64\times 64\) quadrilaterals, and then employ \(k_{U}=2\) or \(k_{U}=4\) for the NN interpolation space. To evaluate the effectiveness of the aforementioned preconditioners, we perform four experiments for each order. Three of them employ the \(\mathbf{B}_{\text{inv}}\), \(\mathbf{B}_{\text{inv\_lin}}\), and \(\mathbf{B}_{\text{GMG}}\) preconditioners, respectively, while the fourth experiment serves as a baseline without any preconditioning, denoted as \(\mathbf{B}_{\text{none}}\). In all experiments, we use the same initial parameters for the NNs to ensure a fair comparison.

Fig. 5 shows the \(L^{2}\) error history of FEINNs using different preconditioners during training for the first 1,000 iterations. We can extract several findings from the figure. Firstly, as \(k_{U}\) increases, the training for the unpreconditioned loss becomes more challenging. This is evident from the flatter error curve for \(\mathbf{B}_{\text{none}}\) in Fig. 4(b) compared to Fig. 4(a). Then, the preconditioners contribute to faster convergence as the error curves of the preconditioned FEINNs are much steeper compared to the one without any preconditioner. Next, the cheaper \(\mathbf{B}_{\text{inv\_lin}}\) preconditioner is surprisingly as effective as the \(\mathbf{B}_{\text{inv}}\) preconditioner. Lastly, \(\mathbf{B}_{\text{GMG}}\) leads to a substantial acceleration of \(L^{2}\) convergence for both order 2 and order 4. Specifically, in Fig. 4(a), for \(k_{U}=2\), the standard unpreconditioned loss function requires more than 800 iterations to reduce the \(L^{2}\)-error below \(10^{-3}\), while \(\mathbf{B}_{\text{inv\_lin}}\) and \(\mathbf{B}_{\text{inv}}\) loss functions attain the same error in around 100 iterations and

Figure 5. \(L^{2}\) error history during training of FEINNs for the forward Poisson problem with a singular solution using different preconditioners.
\(\mathbf{B}_{\text{GMG}}\) requires around 300 iterations. Besides, the GMG-preconditioned FEINN achieves an error level that closely matches the FEINN preconditioned by the other preconditioners after around 900 iterations. Overall, the difference between the errors of these preconditioned FEINNs and the error of the unpreconditioned FEINN exceeds one order of magnitude after enough iterations, and reach two others of magnitude in many cases. Similarly, for \(k_{U}=4\) as shown in Fig (b)b, the unpreconditioned case requires around 1,000 iterations to reduce the \(L^{2}\)-error below \(10^{-3}\), \(\mathbf{B}_{\text{GMG}}\) preconditioned FEINN needs around 300 iterations, and \(\mathbf{B}_{\text{inv\_lin}}\) and \(\mathbf{B}_{\text{inv}}\) preconditioned FEINNs only require around 100 iterations. In this second case, the difference between preconditioned and unpreconditioned training exceeds two orders of magnitude. Although the GMG preconditioner may not be as effective as the other preconditioners, it still exhibits a remarkable reduction in the \(L^{2}\) error compared to the unpreconditioned FEINN, reaching approximately two orders of magnitude after 400 iterations, while being a very cheap preconditioner.

#### 6.1.4. Poisson equation on a complex geometry

In this section, we demonstrate the capabilities of FEINNs in solving forward Poisson problems defined on general domains. We focus on a slightly modified version of [40, Example 4]. As shown in Fig. (a)a, the computational domain \(\Omega\) features a bone-shaped region, which is parameterised by \((x(\theta),y(\theta))\). The parametric equations are defined as \(x(\theta)=0.6\cos(\theta)-0.3\cos(3\theta)\) and \(y(\theta)=0.7\sin(\theta)-0.07\sin(3\theta)+0.2\sin(7\theta)\) with \(\theta\in[0,2\pi]\). Finding appropriate \(\Phi\) and \(\bar{u}\) for IVPINNs is very challenging for this irregular domain, so we only examine the performance of FEINNs in our experiments. We consider the Poisson problem with \(\kappa(x,y)=2+\sin(xy)\), \(\Gamma_{\text{D}}=\partial\Omega\). We choose \(f\) such that the solution is \(u(x,y)=\mathrm{e}^{x}(x^{2}\sin(y)+y^{2})\).

In this study, our focus is on examining the impact of mesh refinement on FEINNs. Consequently, we fix \(k_{U}=2\) and discretise \(\Omega\) using unstructured triangular meshes with an increasing number of cells. In the loss function, we employ the \(\ell^{1}\)-norm for the residual vector. Although the \(\ell^{2}\)-norm is equally effective, we aim to showcase the flexibility in choosing the norm and to provide evidence supporting the suitability of the \(\ell^{1}\) norm for the PDE loss. This is particularly relevant, as we consistently use the \(\ell^{1}\) norm for the PDE part in the subsequent experiments for inverse problems. Similar to the previous sections, we conduct 10 experiments for each mesh resolution, each with distinct initialisations of the NNs.

Fig. 6 illustrates the changes in \(L^{2}\) and \(H^{1}\) errors as the DoFs in the FE interpolation space increase. Overall, FEINNs demonstrate almost identical performance to FEM. Importantly, similar to the findings for the forward convection-diffusion-reaction problem with a smooth solution, the non-interpolated NNs consistently outperform FEM. Notably, when the mesh is "fine enough", there is a remarkable two-order-of-magnitude difference in \(H^{1}\) errors between the NNs and FEM, as illustrated in Fig. (b)b.

To further confirm the superior performance of the NNs in terms of \(H^{1}\) error, we present the point-wise gradient error magnitudes for the FEINN solution and the NN solution in Fig. (b)b and (c)c, respectively. These

Figure 6. Convergence of errors with respect to DoFs of the trial space for the forward Poisson problem on a bone-shaped geometry.
figures correspond to one of our experiments conducted on the finest mesh. Notably, we observe a significant two-order-of-magnitude reduction in error magnitude for the NN solution compared to the FEINN solution across most regions of the domain. Additionally, the lack of smoothness of the gradient of the interpolated solution in FEINN on a low order \(\mathcal{C}^{0}\) space is evident in Fig. (b)b. In contrast, one can observe the smoothness of the error of the FEINN trained NN in Fig. (c)c.

### Inverse problems

In the experiments for inverse problems, we introduce the following relative \(L^{2}\) and \(H^{1}\) errors to measure the accuracy of an identified solution \(z^{id}\):

\[\varepsilon_{L^{2}(\Omega)}(z^{id})=\frac{\left\|z^{id}-z\right\|_{L^{2}( \Omega)}}{\left\|z\right\|_{L^{2}(\Omega)}},\qquad\varepsilon_{H^{1}(\Omega)} (z^{id})=\frac{\left\|z^{id}-z\right\|_{H^{1}(\Omega)}}{\left\|z\right\|_{H^{1 }(\Omega)}},\]

where \(z\) is the ground truth.

The optimisation involving the penalty term (see (8)) occurs at Step 3, requiring the selection of the norm for \(\mathcal{R}_{h}\) and the corresponding coefficient \(\alpha\). In [41, Ch. 17], the authors provide insights into the distinction between utilising \(\ell^{1}\) and \(\ell^{2}\) norms. According to [41, Theorem 17.1], when employing the \(\ell^{2}\)-norm for \(\mathcal{R}_{h}\), the minimiser of (9) becomes a global solution to the inverse problem as \(\alpha\) approaches to infinity. Furthermore, [41, Theorem 17.3] states that there exists an \(\alpha^{*}\) such that the minimiser in (9) for the \(\ell^{1}\)-norm of \(\mathcal{R}_{h}\) is compelled to coincide with the solution of the inverse problem for any \(\alpha\geq\alpha^{*}\). To avoid choosing an arbitrarily large \(\alpha\), we opt to use the \(\ell^{1}\) norm. Moreover, the authors propose [41, Framework 17.2] for adjusting the coefficient \(\alpha\). Following this, we partition Step 3 into several sub-steps. We use a sequence of \(\{\alpha_{k}\}\) for these sub-steps, where \(\alpha_{k}>\alpha_{k-1}\) for \(k>1\).4

Footnote 4: Against common experience in the inverse problem community [32], penalty coefficients for the PDE residual term are usually considered fixed in PINNs and related methods (see [11, 13]). Similarly, for forward problems, the Dirichlet penalty term (which is also a constraint in the minimisation) is usually kept fixed in these formulations.

As mentioned before, we split the training process into three steps. Although we have extensively tested training only (10), the three-step strategy consistently yielded superior results. As a result, all the experiments in this section will follow this training process. We introduce the notation \([n_{1},n_{2},k\times n_{3}]\) to represent the number of iterations for each step: \(n_{1}\) iterations for the data fitting step, followed by \(n_{2}\) iterations for the model parameter initialisations step, and \(k\) sub-steps in the coupled step, with each sub-step consisting of \(n_{3}\) iterations. The sub-steps simply represent a new value of the penalty coefficient. We use the notation \(\alpha=[\alpha_{1},\alpha_{2},...,\alpha_{k}]\), where \(\alpha_{1}\), \(\alpha_{2}\),..., and \(\alpha_{k}\) are the penalty coefficients at each sub-step.

We employ the softplus activation function for FEINNs in our inverse problem experiments, even though the tanh activation generally performs comparably or even better. We aim to explore alternative activation functions for NNs in the context of solving PDE-constrained problems using FEINNs. We use linear FE interpolation space for FEINNs.5 In the remaining experiments, unless otherwise specified, we consider \(z^{id}\) to be the FE interpolation of the NN \(z_{N}\).

Figure 7. True state and gradient error magnitude in FEINN and NN solutions for the forward Poisson problem on a bone-shaped domain.
#### 6.2.1. Poisson equation with partial observations

We begin our inverse problem experiments with a Poisson equation involving partial observations. Following the experiment presented in [9, Sec. 3.1.3], we consider the computational domain \([0,1]^{2}\) with Dirichlet boundary conditions on the left, bottom, and top sides, and a Neumann boundary condition on the right side. The unknown state and diffusion coefficient (Fig. 7(d)) are:

\[u(x,y)=\sin(\pi x)\sin(\pi y),\qquad\kappa(x,y)=1+0.5\sin(2\pi x)\sin(2\pi y).\]

Fig. 7(a) illustrates the true state, and our observations are limited to every DoF inside the white box located at the center of the figure. The objectives of this experiment are to reconstruct the partially known state and to recover the unknown diffusion coefficient.

We discretise the domain by \(50\times 50\) quadrilaterals. Both NNs, \(u_{\mathcal{N}}\) and \(\kappa_{\mathcal{N}}\), have the same structure with \(L=2\) layers and each hidden layer has \(n=20\) neurons. To ensure the positivity of the diffusion coefficient, we apply a rectification function \(r(x)=|x|+0.01\) as the activation for the output layer of \(\kappa_{\mathcal{N}}\). Although \(r(x)=x^{2}+0.01\) also produces satisfactory results, it is more common to use an output layer with linear features. The training iterations are \([400,400,3\times 400]\), and the penalty coefficients are \(\alpha=[0.1,0.3,0.9]\).

In Fig. 8, we display the FEINNs solutions along with their corresponding errors in comparison to the true solutions. The identified state \(u^{id}\) in Fig. 7(b) closely resembles the true state \(u\) in Fig. 7(a), accompanied by very small point-wise errors in Fig. 7(c). These observations highlight the effectiveness of FEINNs at completing the partial observations. Fig. 7(f) displays the small point-wise error of \(\kappa^{id}\), further confirming the accuracy of our approach on discovering the unknown diffusion coefficient.

In order to also consider the relative merits of FEINNs compared to other approaches proposed in the literature, Fig. 9 reports the relative error history for the state and coefficient throughout the training process for FEINNs and our Julia implementation of the adjoint-based NN method (adjoint NN). While adjoint NN approximates the unknown diffusion coefficient with a neural network, it still approximates the state using a FE space, and uses the adjoint solver to compute the gradient of the data misfit with respect to the

Figure 8. Comparison of the true solutions (first column), the FEINN solutions (second column), and corresponding point-wise errors (third column) for the inverse Poisson problem with partial observations. The presented results are from a specific experiment. The first row depicts the state, while the second row represents the coefficient. The observations of \(u\) are limited to the white box.
NN parameters, resulting in a two-loop optimisation process. Adjoint NN was first introduced in [8], and then further explored in [9]. In the experiment, both methods ran their optimisations for 2,000 iterations with identical \(\kappa_{N}\) structures and initialisations. The gaps in the state error curves in Fig. (a)a and (b)b for FEINNs correspond to the second model parameter initialisation step, where \(u_{N}\) is not trained. Similarly, the coefficient error curve in Fig. (c)c for FEINNs starts at iteration 401 as \(\kappa_{N}\) is not trained in the initial data fitting step. The shapes of the curves for FEINNs align with the motivation behind the three-step training process, where the first and second steps aim to lead \(u_{N}\) and \(\kappa_{N}\) to a good initialisation, while the third step focuses on further improving the accuracy. The experiments were performed on a single core of an AMD Ryzen Threadripper 3960X CPU, and we also report the computational cost for both methods: for FEINNs, the average training time is 0.028 seconds per iteration, whereas adjoint NN requires 0.040 seconds per iteration.6 Notably, benefiting from our three-step training strategy and an additional network for state approximation, Fig. 9 reveals that FEINNs have the potential to yield superior accuracy compared to adjoint NN, as all FEINN curves remain below the error curves of adjoint NN after approximately 800 iterations. Moreover, we also plot the error curves for the non-interpolated NNs in Fig. 9. Similar to the findings in Sec. 6.1.1, the smoothness of NN contributes to improved \(H^{1}\) accuracy of a smooth \(u\).

Footnote 6: It is important to note that the computational cost of these two methods is not easily comparable, since they have different computational requirements. The adjoint method involves (non)linear solvers per external iteration, while FEINNs must compute the differentiation of the NN with respect to parameters not only for the physical coefficients but also the state variable. Thus, the relative cost of these methods will be influenced by various factors, including the structure of NNs, the implementation of the (non)linear solver, the specific problem being addressed, etc.

Figure 10. Comparison among FEINNs and adjoint NN in terms of relative errors (depicted using box plots) with different initialisations for the inverse Poisson problem with partial observations.

Figure 9. Comparison among FEINNs and adjoint NN in terms of relative errors during training for the inverse Poisson problem with partial observations. The optimisation loop was run for 2,000 iterations in both cases.
In the sequel, we also compare the robustness with respect to NN initialisation of FEINNs and adjoint NN. We solve the inverse problem 100 times with different NN initialisations. The same NN structure and parameter initialisation of \(\kappa_{N}\) were used for both methods in order to have a fair comparison. In Fig. 10, we depict with box plots the relative errors for the state and the diffusion coefficient from these 100 experiments. Whiskers in the box plot represent the minimum and maximum values within 1.5 times the interquartile range.

Let us first comment on the results obtained with FEINNs. Most of the errors for the state \(\varepsilon_{L^{2}(\Omega)}(u^{id})\) and \(\varepsilon_{H^{1}(\Omega)}(u^{id})\) are very small, with the largest \(\varepsilon_{L^{2}(\Omega)}(u^{id})\) below 0.8%. Besides, the majority of the relative coefficient errors \(\varepsilon_{L^{2}(\Omega)}(\kappa^{id})\) are below 1%. Consequently, we conclude that FEINNs are robust with respect to initialisation in solving this inverse problem with partial observations. Again, the label tag "(NN only)" of FEINNs denotes the errors of the NNs themselves. We observe that the \(L^{2}\) errors for \(\kappa_{N}\) are nearly equivalent to their interpolated counterparts. However, consistent with the findings in Fig. 9, since \(u\) is smooth, \(u_{N}\) surpasses their interpolations in \(H^{1}\) accuracy, with potential for improved \(L^{2}\) accuracy.

The results corresponding to adjoint NN are presented in Fig. 10 as box plots labelled "AdjointNN". During training, we observe that a good initialisation for \(\kappa_{N}\) is imperative, otherwise the optimisation quits prematurely as the gradient norm drops below \(10^{-10}\). This occurrence results in considerably adverse outcomes, at times with \(\varepsilon_{L^{2}(\Omega)}(\kappa^{id})\) exceeding \(10^{8}\). To enhance visual clarity, when constructing the box plots, errors surpassing 1 are standardised to 1. As Fig. 10 indicates, we also explored the activation function tanh as proposed in [9] (tagged as "(\(\rho_{\kappa}=\tanh\))"). However, neither of these configurations produce results outperforming those achieved by FEINNs. Therefore, the adjoint NN method clearly shows less robustness than FEINNs in this partial observations situation.

#### 6.2.2. Poisson equation with noisy observations

In this experiment, we explore the effectiveness of FEINNs in solving an inverse Poisson problem with noisy data. Following the settings in [9, Sec. 3.1.2], we consider the true state (Fig. 10(a)) and diffusion coefficient (Fig. 10(d)) as:

\[u(x,y)=\sin(\pi x)\,\sin(\pi y),\qquad\kappa(x,y)=\frac{1}{1+x^{2}+y^{2}+(x-1) ^{2}+(y-1)^{2}}.\]

The domain \(\Omega\), its discretisation and boundary conditions remain the same as in Sec. 6.2.1. The state at each DoF is known but contaminated with Gaussian noise \(\epsilon\sim N(0,0.05^{2})\). The objectives of this experiment are to reconstruct the state from the noisy data and to estimate the unknown diffusion coefficient.

The structures for \(u_{N}\) and \(\kappa_{N}\) are the same as the ones in Sec. 6.2.1. We again apply \(r(x)=|x|+0.01\) to the output layer of \(\kappa_{N}\) to ensure a positive diffusion coefficient. The training iterations are \([300,100,2\times 300]\), with a total of 1,000, matching the setup in [9, Sec. 3.1.2]. The penalty coefficients are \(\alpha=[1.0,3.0]\).

In Fig. 10, the last two columns display the outcomes from one of our experiments. The identified state \(u^{id}\) in Fig. 10(b) and its low point-wise error in Fig. 10(c) validate the FEINNs capability of recovering the state despite the presence of noise in the data. Fig. 10(e) shows the identified diffusion coefficient \(\kappa^{id}\), which, although visually slightly different from \(\kappa\) in Fig. 10(d), still captures its pattern very well. The point-wise error in Fig. 10(f) further confirms that FEINNs effectively predict the values of the diffusion coefficient.

The error history plots for FEINNs and adjoint NN when applied to the inverse Poisson problem with noisy observations are shown in Fig. 12. The optimisation loop was run in both cases up to 1,000 iterations, and we used the same \(\kappa_{N}\) architecture and parameter initialisation. Consistent with the findings in [9], the loss function in the adjoint method requires explicit regularisation. This is evident as the error curves corresponding to adjoint NN with no regularisation (label tag "(no reg)") start increasing very shortly after the optimisation begins, while the results are much improved by using the regularisation proposed in [9].7 Notably, even _without any regularisation_, the FEINN errors are very stably decreasing. FEINNs could possibly benefit from effective regularisation, but we have not explored this option to keep the method simple and less tuning-dependent. In terms of computational cost, FEINNs demand 0.025 seconds per iteration, while adjoint NN takes 0.043 seconds per iteration. Additionally, in terms of \(u\), the errors of the (non-interpolated) NNs are frequently below their interpolation counterparts during training. This indicates that \(u_{N}\) possesses the capacity to improve \(u\) accuracy despite the noisy observations.

Footnote 7: We use \(\ell^{1}\) regularisation on \(\kappa_{N}\). After testing various regularisation coefficients, we have concluded that the best results are obtained for \(10^{-3}\).

Let us assess the robustness of FEINNs with respect to NN initialisation. We generate the Gaussian noise with the same random seed and repeat the experiment 100 times with differently initialised NNs. The
resulting box plots are shown in Fig. (a)a, where label "FEINN" is for the interpolated NNs and "FEINN (NN only)" is for the non-interpolated ones. We observe that the smoothness of \(u_{N}\) contributes to enhanced accuracy, as both boxes of \(\varepsilon_{L^{2}(\Omega)}(u^{id})\) and \(\varepsilon_{H^{1}(\Omega)}(u^{id})\) for \(u_{N}\) are positioned lower than their interpolation counterparts. Besides, FEINNs generally produce very good results, with \(\varepsilon_{L^{2}(\Omega)}(u^{id})\) mostly below 0.6%, and \(\varepsilon_{L^{2}(\Omega)}(\kappa^{id})\) mostly under 2%.

In this experiment, To compare the performance of FEINNs against adjoint NN with regularisation (as described above), we provide the results for the latter method in the same figure (labelled as "AdjointNN"). We observe that the boxes of \(\varepsilon_{H^{1}(\Omega)}(u^{id})\) and \(\varepsilon_{L^{2}(\Omega)}(\kappa^{id})\) of adjoint NN are positioned higher than that of FEINNs, suggesting that FEINNs generally achieve better accuracy in terms of these two relative errors.

Figure 11. Comparison of the true solutions (first column), the FEINN solutions (second column), and corresponding point-wise errors (third column) for the inverse Poisson problem with noisy observations. The presented results are from a specific experiment. The first row depicts the state, while the second row represents the coefficient.

Figure 12. Comparison among FEINNs and adjoint NN in terms of relative errors during training for the inverse Poisson problem with noisy observations. The optimisation loop was run for 1,000 iterations in both cases.
Furthermore, the state NN in FEINNs generalises well and is far more accurate than the the FE interpolation. In contrast, adjoint NN relies on a FE function for state approximation, lacking such capability of FEINNs.

In this example, we are also interested in exploring how the variability of noise affects FEINNs accuracy. We fix the NN initialisation and the distribution of the Gaussian noise (\(\epsilon\sim N(0,0.05^{2})\)), and repeat the experiment 100 times with different random noise seeds. The resulting box plots are shown in Fig. (b)b with label "FEINN (var noise)". We observe that the noise randomness impacts the accuracy of FEINNs more than NN initialisation randomness, with broader error boxes. Nonetheless, FEINNs are still robust in this scenario, since most \(\varepsilon_{L^{2}(\Omega)}(u^{id})\) are below \(0.8\%\) and most \(\varepsilon_{L^{2}(\Omega)}(\kappa^{id})\) are less than \(3\%\).

#### 6.2.3. Inverse heat conduction problem

In our final experiment for this paper, we attack an inverse heat conduction problem (HCP). In many heat transfer applications, the boundary values are either unavailable or difficult to measure over the entire surface. The goal of HCPs is to estimate the surface temperature (Dirichlet boundary value), and/or heat flux (Neumann boundary value), based on temperature data measured at certain points within the domain [42]. Our example combines the challenges in [42] and [43], where we consider a two-layered half-tube cross-section as the computational domain \(\Omega\), as shown in Fig. (a)a. The domain \(\Omega\) can be described in polar coordinates as \(\theta\in[0,\pi]\) and \(r\in[0.05,0.11]\). The tube is composed of two layers of media, with a diffusion coefficient of \(\kappa_{1}=1\) for \(r\in[0.05,0.08]\), and \(\kappa_{2}=100\) for \(r\in[0.08,0.11]\). The unknown boundary values are, in polar coordinates,

\[g(\theta,r)=200,\ r=0.05,\qquad\eta(\theta,r)=-100-50\sin(\theta),\ r=0.11.\]

The horizontal section of the tube is also a Neumann boundary, with known \(\eta=0\).

We discretise the domain with \(2\times 50\times 50\) triangles and solve the forward problem using FEM with the aforementioned boundary conditions. Fig. (a)a shows the FEM solution of the temperature, and we use the temperature at the yellow dots as our observations. Since the temperature has different patterns in the two layers due to the discontinuity in the diffusion coefficient, we use a deeper NN with 6 layers and 20 neurons for each hidden layer (\(L_{u}=6\), \(n_{u}=20\)) as \(u_{N}\). The Dirichlet boundary value \(g\) is just a part of \(u\), so in this problem, \(u_{N}\) is defined over the whole domain, including the Dirichlet boundary. We evaluate \(u^{id}\) on the Dirichlet boundary to obtain \(g^{id}\). We train another NN \(\eta_{N}\) with \(L_{\eta}=3\) and \(n_{\eta}=20\) on the Neumann boundary to predict the Neumann boundary value. We set the number of training iterations to \([700,200,3\times 700]\), and use the penalty coefficients \(\alpha=[0.001,0.003,0.009]\).

Fig. (b)b shows the relative point-wise errors of FEINNs solutions for the boundary values \(\eta\) and \(g\), obtained from one of our experiments. The errors at most of the Neumann boundary points are below \(1\%\), indicating accurate recovery. Besides, the identified Dirichlet value \(g^{id}\) is even more accurate, with a maximum error of approximately \(3\times 10^{-6}\). Overall, FEINNs excel at accurately reconstructing the boundary values. Fig. (a)a depicts the history of relative errors during training from the same experiment. We observe

Figure 13. Left: comparison among FEINNs and adjoint NN in terms of relative errors (depicted using box plots) with different NN initialisations for the inverse Poisson problem with noisy observations. Right: box plots depicting the relative errors of FEINNs trained with Gaussian noise generated by different random seeds.
that the data step and model parameter initialisation step reduce corresponding errors as expected, and the errors steadily decrease after a few hundred iterations of adjustment in the coupled step.

To study FEINNs' reliability in solving IHCPs, we repeat the experiment 100 times with different NN initialisations. The resulting errors are presented in Fig. (b)b as box plots along with the original data points. Even though the number of observations (100) is much smaller that the DoFs (2,601) of the trial space, FEINNs recover the temperature distribution accurately, with most \(\varepsilon_{L^{2}(\Omega)}(\mu^{id})\) errors below \(10^{-5}\), and only a few outliers with higher errors. The proposed formulation demonstrates robustness despite a significant discontinuity in the diffusion coefficient, a limited number of observations, and no regularisation. Furthermore, the majority of experiments (at least 90%) yield remarkably low errors.

## 7. Conclusions

In this paper, we propose a general framework, called FEINNs, to approximate forward and inverse problems governed by _low-dimensional_ PDEs, by combining NNs and FEs to overcome some of the limitations (numerical integration error, treatment of Dirichlet boundary conditions, lack of solid mathematical foundations) of existing approaches proposed in the literature to approximate PDEs with NNs, such as, e.g., PINNs. For forward problems, we interpolate the NN onto the FE space with zero traces (non-homogeneous Dirichlet boundary conditions are enforced via a standard offset FE function), and evaluate the FE residual for the resulting FE function. The loss function is the norm of the FE residual. We propose different

Figure 14. True temperature distribution and the relative point-wise errors of FEINNs on the Dirichlet and Neumann boundaries for the IHCP. Yellow dots on the temperature figure indicate observation locations.

Figure 15. The error history during training from a specific experiment and the box plots of the relative errors of FEINNs with different initialisations for the IHCP.
norms, and suggest the use of standard FE preconditioners (e.g., a fixed number of GMG cycles) to end up with a well-posed loss function in the limit \(h\downarrow 0\). For inverse problems, the unknown model parameters are parametrised via NNs, which can also be interpolated onto FE spaces. The loss function in this case combines the data misfit term with a penalty term for the PDE residual. We propose a three step algorithm to speed up the training process of the resulting formulation, where we perform two cheap data fitting steps (no differential operators involved) to provide a good initialisation for a fully coupled minimisation step.

We have conducted numerous numerical experiments to assess the computational performance and accuracy of FEINNs. We use forward convection-diffusion-reaction problems to compare FEINNs against IVPINNs, a recently proposed related method which mainly differs in the treatment of Dirichlet boundary conditions and has been proven to be superior to other PINN formulations in certain situations [28]. The computational cost per iteration of IVPINNs and FEINNs is virtually the same. However, IVPINNs struggle to keep the convergence of FEINNs (and reach the FEM error) as we increase mesh resolution or polynomial order. Additionally, the (non-interpolated) NNs trained with FEINNs exhibits excellent generalisation, with superior performance compared to the FE solution and the non-interpolated NN composition of IVPINNs. For singular solutions, both FEINNs and IVPINNs have comparable performance to FEM. We evaluate the effect of the residual norm and show how preconditioned norms accelerate the training. Moreover, experiments performed on a non-trivial geometry highlights the capability FEINNs handling complex geometries and Dirichlet boundary condition effortlessly, which is not the case of IVPINNs or standard PINNs.

In the experiments for the inverse problems, we show that FEINNs are capable of estimating unknown diffusion coefficient from partial or noisy observations of the state and recovering the unknown boundary values from discrete observations. We additionally compare the performance of FEINNs against the adjoint-based NN method [8, 9]. The numerical results demonstrate that FEINNs exhibit greater robustness for partial observations and are comparable for noisy observations. However, adjoint methods require the tunning of the regularisation term to be effective, while FEINNs are robust without regularisation. The conducted experiments also prove that the three-step training process employed by FEINNs is a sound strategy.

This work can be extended in many directions. First, one could consider transient and/or nonlinear PDEs, in which NNs and non-convex optimisation have additional benefits compared to standard linearisation and iterative linear solvers in FEM. Besides, while this work concentrates on problems in \(H^{1}\), the framework can be extended to problems in \(H(\text{curl})\) and \(H(\text{div})\) spaces, combined with compatible FEM [2, 44]. To target large scale problems, one could design domain decomposition [45, 46] and partition of unity methods [47] to end up with suitable algorithms for massively parallel distributed-memory platforms and exploit existing parallel FE frameworks GridapDistributed.jl[48]. Lastly, we want to explore in the future the usage of adaptive meshes to exploit the nonlinear approximability of NNs within the same training loop.

## 8. Acknowledgments

This research was partially funded by the Australian Government through the Australian Research Council (project numbers DP210103092 and DP220103160). This work was also supported by computational resources provided by the Australian Government through NCI under the NCMAS and ANU Merit Allocation Schemes. W. Li acknowledges the support from the Laboratory for Turbulence Research in Aerospace and Combustion (LTRAC) at Monash University through the use of their HPC Clusters.
