Error analysis for hybrid finite element/neural network discretizations
We describe and analyze a hybrid finite element/neural network method for
predicting solutions of partial differential equations. The methodology is
designed for obtaining fine scale fluctuations from neural networks in a local
manner. The network is capable of locally correcting a coarse finite element
solution towards a fine solution taking the source term and the coarse
approximation as input. Key observation is the dependency between quality of
predictions and the size of training set which consists of different source
terms and corresponding fine & coarse solutions. We provide the a priori error
analysis of the method together with the stability analysis of the neural
network. The numerical experiments confirm the capability of the network
predicting fine finite element solutions. We also illustrate the generalization
of the method to problems where test and training domains differ from each
other.
Uladzislau Kapustsin, Utku Kaya, Thomas Richter
# Error analysis for hybrid finite element/neural network discretizations

###### Abstract

We describe and analyze a hybrid finite element/neural network method for predicting solutions of partial differential equations. The methodology is designed for obtaining fine scale fluctuations from neural networks in a local manner. The network is capable of locally correcting a coarse finite element solution towards a fine solution taking the source term and the coarse approximation as input. Key observation is the dependency between quality of predictions and the size of training set which consists of different source terms and corresponding fine & coarse solutions. We provide the a priori error analysis of the method together with the stability analysis of the neural network. The numerical experiments confirm the capability of the network predicting fine finite element solutions. We also illustrate the generalization of the method to problems where test and training domains differ from each other.

Keywords:Neural networks PDE approximation Finite element method Msc: 65Y10 65N12 65N30 

## 1 Introduction and motivation

The use of neural networks to approximate solutions to partial differential equations has made considerable progress in recent years. In particular the class of Physics Inspired Neural Networks (PINN) [27], like the Deep Ritz method, has been investigated intensely [8]. Especially for high dimensional problems [18; 4; 9] or parameter dependent partial differential equations [1] these approaches offer a new point of view and promise a substantial increase in efficiency. There are further methods in the literature which exploit the weak formulations of the PDEs [14] or incorporate finite element (FE) spaces in the loss function [2; 22; 24]. See [28] for a review.

However, when standard problems such as three dimensional fluid mechanics are considered, neural network based approaches have to compete with highly sophisticated and refined classical discretization methods. Finally, the training of the
network often remains the crucial problem. Even though the theoretical approximation properties of neural networks may be superior, the efficient solution of the associated optimization problems is still an open problem. On the other hand, there are highly efficient Newton-Krylov space methods, possibly with optimal multigrid methods for preconditioning.

The main drawback of simple PINN's is the need to re-train when the problem parameters change. DeepONet [17] is one approach to this issue. Here, instead of training one net to represent a specific problem, two nets are used to represent possible solutions as well as the solution operator. But even this approach could not yet prove the efficiency and accuracy of established methods.

An alternative approach we follow is to combine established methods for coarse representation of a solution with a neural network to resolve fine scales whose representation is often not possible. This approach has many potential applications, for example in fluid mechanics, where finite element or finite volume methods can reproduce the coarse structure with great accuracy while respecting conservation principles, but simultaneous resolution of fine-scale turbulent processes is often not possible or would simply be too expensive. A super-resolution methodology in this direction was introduced in [11]. With the _Deep Neural Network Multigrid Solver_ (DNN-MG) [21; 19; 20] we have introduced a concept which embeds a neural network fluently into a multigrid hierarchy, solves the coarse grid levels directly, e.g. with a finite element multigrid method, and predicts the corrections on fine grid levels locally by a neural network. An application to stationary linear problems including a simplified error analysis of this hybrid approach has been studied in our preliminary work [12].

Here, the neural network intervenes only locally: the grid is decomposed into patches, e.g., in 2d a range of \(p\times p\) (\(p\in\mathbb{N}\) is small here, usually less than four) elements, and on each of these patches the correction to a finer solution is pulled from the network. This approach was shown to increase efficiency for standard flow problems compared to established methods. Furthermore, the local design, i.e. the application of a mesh to all patches, allows a very good generalizability. The DNN-MG solver can be regarded as a numerical solution method in the sense of a domain decomposition method rather than an approximation method. The network never sees the whole solution but always only small sections.

The mathematical analysis of PINNs is already well advanced. In particular, the aspect of approximation with neural networks is very well established, starting with the universal approximation theorem in the late 80s [3; 6]. Especially for the application to partial differential equations relevant results are available [10; 18; 25] and optimal estimates in \(W^{n,p}\)-spaces are known, see [7] for a review neural network's approximation properties. Of the PINNs, the Deep Ritz method [8] in particular is well studied. It is based on the direct approximation of the energy functional (in the case of symmetric differential operators, e.g., Laplace or Stokes) with neural networks and Monte Carlo integration. Here, quite comprehensive a priori [25] as well as a posteriori error estimation are available [23].

The goal of this work is to investigate hybrid approaches that enrich a finite element solution on coarse grids with fine scale fluctuations from a neural network in terms of the DNN-MG method. We restrict ourselves to the simple linear Poisson equation and give the complete a priori error analysis of the hybrid method.

In the next section we will briefly describe the finite element discretization and introduce some notation. Section 3 then introduces the hybrid approximation
method and describes the training of the neural networks. We start the analysis of the method with the version where global network updates are used and then refine the a priori analysis for the local neural network updates. Numerical demonstrations follow in Section 4.

## 2 Preliminaries

We start with describing the model problem and then introduce the notation and finite element discretization. Let \(\Omega\subset\mathbb{R}^{d},\ d\in\{2,3\}\) be a domain with polygonal boundary. For \(f\in H^{-1}(\Omega)\) let \(u\in H^{1}_{0}(\Omega)\) be the weak solution to the Poisson equation

\[-\Delta u=f,\quad u|_{\partial\Omega}=0. \tag{1}\]

### Finite element discretization

By \(\Omega_{h}\) we denote a finite element mesh of the domain \(\Omega\), where \(h\) stands for the diameter of the largest cell. For simplicity we assume that \(\Omega_{h}\) is a simplical mesh that satisfies the common assumptions of structural regularity and shape regularity [5]. We assume that there is a hierarchy of finite element meshes

\[\Omega_{H}:=\Omega_{0}\preccurlyeq\Omega_{1}\prec\cdots\preccurlyeq\Omega_{L}=: \Omega_{h}, \tag{2}\]

where we denote by \(\Omega_{l-1}\preccurlyeq\Omega_{l}\), that each element of the fine mesh \(T\in\Omega_{l}\) originates from the uniform refinement of a coarse element \(T^{\prime}\in\Omega_{l-1}\), for instance, uniform splitting of a quadrilateral or triangular element into four and of a hexahedral or tetrahedral element into eight smaller ones, respectively.

On \(\Omega_{h}\) let \(V_{h}\) be the space of piecewise polynomials of degree \(r\geq 1\) satisfying the homogeneous Dirichlet condition on the boundary \(\partial\Omega\)

\[V_{h}:=\left\{\phi\in C(\bar{\Omega}),\ \phi|_{T}\in P^{(r)}(T)\ \forall T\in \Omega_{h},\ \phi|_{\partial\Omega}=0\right\}\]

where \(P^{(r)}(T)\) is the space of polynomials of degree \(r\) on a cell \(T\in\Omega_{h}\). On the hierarchy of meshes (2) we hence define the hierarchy of spaces \(V_{h}^{(l)}\) on \(\Omega_{l}\) for \(l=0,\ldots,L\). These spaces are nested

\[V_{h}^{(l-1)}\subset V_{h}^{(l)},\quad l=1,\ldots,L.\]

If not necessary for understanding of the specific context, we will skip the index \(l\) referring to the mesh level. Then, \(u_{h}\in V_{h}\) is the finite element solution to

\[(\nabla u_{h},\nabla\phi_{h})=(f,\phi_{h})\quad\forall\phi_{h}\in V_{h}, \tag{3}\]

with the \(L^{2}\) inner product \((\cdot,\cdot)\). For two right hand sides \(f,g\in H^{-1}(\Omega)\) and the corresponding finite element solutions \(u_{h}^{f},u_{h}^{g}\in V_{h}\) it holds

\[\|\nabla(u_{h}^{f}-u_{h}^{g})\|_{L^{2}(\Omega)}=\|f-g\|_{H^{-1}(\Omega)} \tag{4}\]

Given sufficient regularity of the right hand side, namely \(f\in H^{r-1}(\Omega)\) (using for simplicity the notation \(H^{0}(\Omega):=L^{2}(\Omega)\)) and of the domain, either having a
\(C^{r+1}\)-boundary or, in the case \(r=1\) being convex, the standard a priori error estimate

\[\|u-u_{h}\|_{L^{2}(\Omega)}+h\|\nabla(u-u_{h})\|_{L^{2}(\Omega)}\leq ch^{r+1}\|f \|_{H^{r-1}(\Omega)} \tag{5}\]

holds.

### Notation

Over subdomains \(\omega\subseteq\Omega\) the \(L^{2}-\) and \(H^{s}-\) norms are denoted by \(\|\cdot\|_{\omega}\) and \(\|\cdot\|_{s,\omega}\), respectively. We omit the index \(\omega\), if the norm is considered on the whole domain \(\Omega\). With \(\|\cdot\|_{2}\) we denote the Euclidean norm for vectors and the spectral norm for matrices. With \(X_{\omega}^{h}\) and \(X_{\omega}^{H}\) we denote the nodes of the meshes \(\Omega_{h}\) and \(\Omega_{H}\) that lie in the subdomain \(\omega\), respectively. Moreover, for \(v\in C(\overline{\omega})\) we define

\[\|v\|_{L^{2}(\omega)}:=\Big{(}\sum_{x\in X_{\omega}^{h}}v(x)^{2}\Big{)}^{\frac {1}{2}}. \tag{6}\]

Definition 1 (Patch): A patch \(\mathcal{P}\in\Omega_{h}\) is defined to be a subdomain that is geometrically identical to one certain cell \(M\) of the coarse mesh \(\Omega_{H}\). We exploit that a patch is not only identified with the degrees of freedoms of the element \(M\) but also by the cells assembling it, \(\mathcal{P}=\{T\in\Omega_{h}\ :T\subset M\}\).

By \(V_{\mathcal{P}}\) we denote the local finite element subspace

\[V_{\mathcal{P}}:=\text{span}\left\{\phi_{h}|_{\mathcal{P}},\ \phi_{h}\in V_{h} \right\}.\]

By \(R_{\mathcal{P}}:V_{h}\to V_{\mathcal{P}}\) we denote the restriction to the local patch space, defined via

\[R_{\mathcal{P}}(u_{h})(x)=u_{h}(x)\quad\forall x\in X_{\mathcal{P}}^{h}.\]

By \(P_{\mathcal{P}}:V_{\mathcal{P}}\to V_{h}\) we denote the prolongation defined by

\[P_{\mathcal{P}}(u_{\mathcal{P}})(x)=\begin{cases}\frac{1}{n(x)}u_{\mathcal{P} }(x)&x\in X^{h},\\ 0&\text{otherwise},\end{cases}\]

where \(n(x)\in\mathbb{N}\) is the number of patches that contain the degree of freedom \(x\).

## 3 Hybrid finite element neural network discretization

Consider two finite dimensional spaces \(V_{H}\) and \(V_{h}\) that are built on coarse and fine meshes \(\Omega_{h}\) and \(\Omega_{H}\), respectively. The idea of the paper is to determine an approximate solution on the coarse mesh \(\Omega_{H}\) with the finite element method and then to obtain the fine mesh fluctuations in forms of neural network updates. In other words, we seek hybrid solutions \(u_{\mathcal{N}}\) which are found by augmenting \(u_{H}\in V_{H}\) with a neural network update \(w_{\mathcal{N}}\in V_{h}\), i.e.

\[u_{\mathcal{N}}:=u_{H}+w_{\mathcal{N}}\in V_{h}.\]
The neural network predicts the finite element coefficients on the fine mesh \(\Omega_{h}\) and as input it receives all data that is available: the coarse solution \(u_{H}\in V_{H}\), the problem data, i.e. the source term \(f\). The output \(\mathcal{N}(u_{H},f)\) is a vector, which determines the coefficients of \(w_{\mathcal{N}}\in V_{h}\).

Moreover, these updates are to be obtained locally in such a way that the network is not acting on the aforementioned data on the whole domain \(\Omega\). Instead the network acts on the patches \(\mathcal{P}\) separately in order to obtain the values of an update \(w_{\mathcal{N}}|_{\mathcal{P}}\) by providing the coefficient vector \(\mathcal{N}(u_{H}|_{\mathcal{P}},f|_{\mathcal{P}})\) of \(V_{\mathcal{P}}\). An illustration of the local updates is given in Fig. 1.

Definition 2 (Hybrid solution): Given a subdivision of the domain \(\Omega\) into a set of patches, a neural network function \(\mathcal{N}(\cdot)\) acting on each patch and the coarse mesh solution \(u_{H}\in V_{H}\), we define the hybrid solution \(u_{\mathcal{N}}\in V_{h}\) as

\[u_{\mathcal{N}}:=u_{H}+\sum_{\mathcal{P}}P_{\mathcal{P}}w_{N}^{\mathcal{P}},\]

where \(w_{N}^{\mathcal{P}}=\sum_{i=1}^{N_{dof}}\mathcal{N}(R_{\mathcal{P}}u_{H},R_{ \mathcal{P}}f)_{i}\phi_{i}^{\mathcal{P}}\). Here, \(\mathcal{N}(R_{\mathcal{P}}u_{H},R_{\mathcal{P}}f)_{i}\) is the \(i\)-th output of the network and \(\phi_{i}^{\mathcal{P}},\;i\in\{1,\cdots,N_{dof}\}\) are the basis functions of \(V_{\mathcal{P}}\).

As networks we will only consider fully connected multilayer perceptrons:

Definition 3 (Multilayer perceptron): The function \(\mathcal{N}:\mathbb{R}^{N_{0}}\rightarrow\mathbb{R}^{N_{L}}\) defined via

\[\mathcal{N}=l_{L}\circ\sigma\circ l_{L-1}\circ\cdots\circ\sigma\circ l_{1} \tag{7}\]

is called a _multilayer perceptron (MLP)_ of _depth_\(L\in\mathbb{N}\) with an _activation function_\(\sigma:\mathbb{R}\rightarrow\mathbb{R}\). The \(l_{i}(x):\mathbb{R}^{N_{i-1}}\rightarrow\mathbb{R}^{N_{i}}\) are called _layers_ and defined as

\[l_{i}(x)=W_{i}x+b_{i},\;i=1,\ldots,L \tag{8}\]

Figure 1: Illustration of the hybrid solver. The finite element solution \(u_{H}\) is approximated on the coarse mesh (left). On each patch (one or multiple elements) this solution is locally extracted and interpolated into a refined mesh. Together with the fine mesh right hand side information \(f_{h}\) it is the input of a neural network. The output \(w_{\mathcal{N}}\) is the local correction towards an improved solution and is prolongated back onto the fine global mesh.
where \(W_{i}\in\mathbb{R}^{N_{i-1}\times N_{i}}\) are _weights_ and \(b_{i}\in\mathbb{R}^{N_{i}}\) are _biases_.

The input to the neural network (on each patch) is a vector \(x_{\mathcal{P}}\in\mathbb{R}^{N_{0}}\) and we decompose it into \(x_{\mathcal{P}}=(x_{\mathcal{P}}^{u_{H}},x_{\mathcal{P}}^{f})\), where

\[x_{\mathcal{P}}^{u_{H}}=\big{(}R_{\mathcal{P}}u_{H}(x)\big{)}_{x\in X_{ \mathcal{P}}^{H}}\;,\quad x_{\mathcal{P}}^{f}=\big{(}R_{\mathcal{P}}f(x)\big{)} _{x\in X_{\mathcal{P}}^{h}}\;. \tag{9}\]

### Training of the neural network

Now, let's explain how we pre-train a neural network \(\mathcal{N}\). First of all, it is necessary to select a set of training problems. We generate the training data by selecting a set \(\mathcal{F}\subset H^{-1}(\Omega)\) of right hand side functions \(f\). Then, we solve the Poisson equation for each \(f\in\mathcal{F}\) for both coarse and fine meshes \(\Omega_{H}\) and \(\Omega_{h}\), respectively. The input data is evaluated as described in (9). The output data \(z_{\mathcal{P}}\in\mathbb{R}^{N_{L}}\) is given as

\[z_{\mathcal{P}}=\big{(}u_{h}(x)-u_{H}(x)\big{)}_{x\in X_{\mathcal{P}}^{h}}^{T}\]

which is simply a difference between fine and coarse solution on each fine mesh node \(x\) belonging to the patch \(\mathcal{P}\). We would like to underline that the network update of the solution \(w_{\mathcal{N}}\) is given by the value of the network applied not globally to the whole domain \(\Omega\), but locally to each patch \(\mathcal{P}\) (see Definition 1).

Once the training data is available, we need to solve the following optimization problem

\[\min_{\begin{subarray}{c}W_{i}\in\mathbb{R}^{N_{i-1},N_{i}},b_{i}\in\mathbb{R }^{N_{i}},\\ i\in\{0,\ldots,L\}\end{subarray}}\quad\frac{1}{N_{T}N_{\mathcal{P}}}\sum_{ \mathcal{P}}\lVert z_{\mathcal{P}}-\mathcal{N}(y_{\mathcal{P}})\rVert_{2}^{2} \tag{10}\]

where \(N_{T}\) is the number of source terms in the training set \(\mathcal{F}\) and \(N_{\mathcal{P}}\) is the number of patches. In order to solve this problem we use one of the stochastic gradient descent based methods. After the network has been trained, we can finally apply it to other problems that were not in the training data. For this, firstly we again construct input data in the same way as described above and compute network

Figure 2: Multilayer perceptron
predictions. Then, we construct a complete network solution \(u_{\mathcal{N}}\) by summing up (and averaging) these predictions as described in Definition 2.

Considering everything we have developed so far from the finite element perspective we note that each coarse solution \(u_{H}\in V_{H}\) also belongs to \(V_{h}\) and there it takes the form

\[u_{H}=\sum_{i=1}^{N_{\mathit{dof}}}U_{h}^{i}\phi_{h}^{i}\]

where \(\{\phi_{h}^{i}\}_{i=1}^{N_{\mathit{dof}}}\) is the basis of the fine finite element space \(V_{h}\) and where \(U_{h}^{i}\in\mathbb{R}\) are corresponding coefficients. As a consequence of the fact that, we update the coarse solution \(u_{H}\) only on fine mesh nodes, we can consider this whole procedure as a simple update of fine mesh coefficients \(U_{h}^{i}\), i.e.

\[u_{\mathcal{N}}=\sum_{i=1}^{N_{\mathit{dof}}}(U_{h}^{i}+W_{\mathcal{N}}^{i}) \phi_{h}^{i}\in V_{h}.\]

### Error estimate using a global network update

We will start with a simplified setting, where there is only one patch that covers the complete domain, \(\mathcal{P}=\Omega_{h}\), compare Definition 1. Here, the hybrid finite element solution is directly given by

\[u_{\mathcal{N}}=u_{H}+\mathcal{N}(R_{\Omega_{h}}u_{H},R_{\Omega_{h}}f).\]

For brevity of notation we simply write \(u_{\mathcal{N}}=u_{H}+\mathcal{N}(f)\), implicitly assuming that the network also receives values of the coarse mesh solution \(u_{H}\in V_{H}\) corresponding to \(f\) as an input.

Theorem 3.1 (A priori finite element error for the single-patch solution): _Let \(\mathcal{N}\) be the network trained on the training set \(\mathcal{F}=\{f_{1},\ldots,f_{N_{T}}\}\) such that the loss function (10) is reduced to the error order of \(\epsilon^{2}\). For \(f\in H^{r-1}(\Omega)\), let \(u\in H^{r+1}(\Omega)\cap H_{0}^{1}(\Omega)\) be the solution to the Poisson problem, \(u_{H}\in V_{H}\) be the coarse finite element approximation of polynomial degree \(r\geq 1\) and \(u_{\mathcal{N}}=u_{H}+\mathcal{N}(f)\in V_{h}\) be the hybrid solution. It holds_

\[\|\nabla(u-u_{\mathcal{N}})\|\leq C\left(h^{r}\|f\|_{r-1}+\min_{f_{i}\in \mathcal{F}}\left\{\|f-f_{i}\|_{-1}+\|\nabla(\mathcal{N}(f)-\mathcal{N}(f_{i} ))\|\right\}+\epsilon\right),\]

_where \(h<H\) is the mesh size of the fine space \(V_{h}\) and \(\epsilon\) is the network approximation and training error for the training data set._

Proof: For arbitrary \(f_{i}\in\mathcal{F}\) we split the error

\[\|\nabla(u-u_{\mathcal{N}})\|\leq\underbrace{\|\nabla(u-u_{h})\|}_{=(I)}+ \underbrace{\|\nabla(u_{h}-u_{h}^{f_{i}})\|}_{=(II)}+\underbrace{\|\nabla(u_{ h}^{f_{i}}-u_{\mathcal{N}}^{f_{i}})\|}_{=(III)}+\underbrace{\|\nabla(u_{\mathcal{N}}^{f_{i}}-u_{ \mathcal{N}})\|}_{=(IV)}\]

into the fine mesh finite element error

\[(I)=\|\nabla(u-u_{h})\|\leq Ch^{r}\|f\|_{r-1},\]
where \(u_{h}\in V_{h}\) is the finite element solution in the resolved space to the right hand side \(f\). Next, into the data approximation error

\[(II)=\|\nabla(u_{h}-u_{h}^{f_{i}})\|\leq\|f-f_{i}\|_{-1},\]

where \(u_{h}^{f_{i}}\in V_{h}\) is the fine finite element solution to an arbitrary \(f_{i}\in\mathcal{F}\) from the training data set and into the network approximation and optimization error

\[(III)=\|\nabla(u_{h}^{f_{i}}-u_{\mathcal{N}}^{f_{i}})\|\leq C\epsilon.\]

Finally, using the composition \(u_{\mathcal{N}}=u_{H}+\mathcal{N}(f)\), the generalization error of the network and a further error term depending on the richness of the data set remains

\[(IV)=\|\nabla(u_{\mathcal{N}}^{f_{i}}-u_{\mathcal{N}})\| \leq\|\nabla(u_{H}^{f_{i}}-u_{H})\|+\|\nabla(\mathcal{N}(f_{i})- \mathcal{N}(f))\|\] \[\leq\|f-f_{i}\|_{-1}+\|\nabla(\mathcal{N}(f_{i})-\mathcal{N}(f))\|.\]

Combining the above gives the result.

This lemma shows that the hybrid approach is able to reduce the error up to the accuracy of the fine finite element space \(V_{h}\), the tolerance of the neural network training, the richness of the training data set \(\mathcal{F}\) and the stability of the neural network that governs it ability to generalize beyond the training data set. This last term

\[\|\nabla(\mathcal{N}(f_{i})-\mathcal{N}(f))\|\]

will depend on the design of the neural network and it will be investigated in Section 3.4.

### Error estimates for local neural network updates

We now tend to the discussion of a local application of the neural network. The global analysis in Section 3.2 is still valid, if the domain \(\Omega^{tr}\), where the training data is generated and the domain \(\Omega\), where the actual simulation is run are the same and if the same meshes are used. The localized approach will however allow us for more flexibility when it comes to choosing an approximation of the right hand side to one of the training data \(f\approx f_{i}\), as this can be done individually on each patch. Furthermore, if we allow for a generalization of the domains, i.e. \(\Omega\neq\Omega^{tr}\) the analysis also must be refined.

To keep the notation simple, we will make the following assumptions on the meshes, see also Fig. 3.

Assumption (Uniform compatability of meshes and patches)Let \(\Omega\) be the domain of application with meshes \(\Omega_{h}\) and \(\Omega_{H}\). By \(\Omega^{tr}\) we denote the domain where the training data is generated and by \(\Omega_{h}^{tr}\) and \(\Omega_{H}^{tr}\) the corresponding meshes. We assume that the patches \(\{\mathcal{P}_{1},\ldots,\mathcal{P}_{N_{\mathcal{P}}}\}\) of \(\Omega_{h}\) and \(\Omega_{H}\) and the patches \(\{\mathcal{P}_{1}^{tr},\ldots,\mathcal{P}_{N_{\mathcal{P}}^{tr}}^{tr}\}\) are compatible in the following sense: For each \(\mathcal{P}\in\Omega_{h}\), there exists a training patch \(\mathcal{P}^{tr}\in\Omega_{h}^{tr}\) which is a translation and/or rotation of \(\mathcal{P}\).
This assumption is restrictive and only allows fully uniform meshes. We refer to Remark 1 for some hints on how the approach can be extended to more general settings.

Moreover, for each patch \(\mathcal{P}\in\Omega_{h}\) we define an enlarged domain \(\tilde{\mathcal{P}}\) with \(\mathcal{P}\subset\subset\tilde{\mathcal{P}}\) that could be the union of all patches \(\mathcal{P}^{\prime}\in\Omega_{h}\) that overlap with \(\mathcal{P}\), i.e.

\[\tilde{\mathcal{P}}:=\bigcup_{\begin{subarray}{c}\mathcal{P}^{\prime}\in \Omega_{h}\\ \tilde{\mathcal{P}}^{\prime}\cap\tilde{\mathcal{P}}\neq\emptyset\end{subarray}} \mathcal{P}^{\prime} \tag{11}\]

For the distance between \(\mathcal{P}\) and \(\partial\tilde{\mathcal{P}}\setminus\partial\Omega\) it holds \(d(\mathcal{P},\tilde{\mathcal{P}})=\mathrm{dist}(\mathcal{P},\partial\tilde{ \mathcal{P}}\setminus\partial\Omega)=\mathcal{O}(h_{\mathcal{P}})\) where \(h_{\mathcal{P}}\) is the diameter of a patch. This distance is relevant for local error estimates such as (13). Now we present the main result of this paper which we will prove later on.

Theorem 2.2 (A priori finite element error for the hybrid solution based on local patches): _For \(f\in H^{r-1}(\Omega)\), let \(u\in H^{r+1}(\Omega)\cap H^{1}_{0}(\Omega)\) be the solution to the Poisson problem, \(u_{H}\in V_{H}\) be the coarse finite element approximation of degree \(r\geq 1\) and \(u^{i}_{H}\in V_{H}\) be the coarse finite element approximation solution that corresponds to \(f_{i}\in L^{2}(\Omega)\), \(i=1,\ldots,N_{T}\). Let Assumption 3.3 be satisfied and \(\mathcal{N}\) be the network trained on the training set_

\[\mathcal{F}_{\mathcal{P}}=\{(f_{1},u^{1}_{H},u^{1}_{h})|_{\mathcal{P}},\ldots, (f_{N_{T}},u^{N_{T}}_{H},u^{N_{T}}_{h})|_{\mathcal{P}}:\mathcal{P}\in\Omega_{ h}\}\]

_such that the loss (10) is reduced below the error order of \(\epsilon^{2}\)._

_For the hybrid solution \(u_{\mathcal{N}}=u_{H}+\sum_{\mathcal{P}}P_{\mathcal{P}}\mathcal{N}(f)\in V_{h}\) it holds_

\[\|\nabla(u_{h}-u_{\mathcal{N}})\|\leq c\Big{[}\big{(}h^{r}+H^{2r} \big{)}\|f\|_{r-1}+H^{2r}\max_{1\leq i\leq N_{T}}\|f_{i}\|_{r-1,\Omega^{tr}}\] \[+\Big{(}\sum_{\mathcal{P}\in\Omega_{h}}\min_{1\leq i\leq N_{T}} \Big{\{}\|f\!-\!f_{i}\|^{2}_{-1,\tilde{\mathcal{P}}}\!+\!\|\nabla(u_{H}\!-\!u^ {i}_{H})\|_{\tilde{\mathcal{P}}}^{2}\!+\!\|\nabla(\mathcal{N}(f_{i})\!-\! \mathcal{N}(f))\|_{\mathcal{P}}^{2}\Big{\}}\Big{)}^{\frac{1}{2}}\Big{]}+\epsilon. \tag{12}\]

Figure 3: Training and application domains \(\Omega^{tr}\) and \(\Omega\) as well as the corresponding meshes \(\Omega^{tr}_{h}\) and \(\Omega_{h}\) can differ. Both however must be split into the same kind of patches. A patch \(\mathcal{P}\) and the extended patch \(\tilde{\mathcal{P}}\) is marked in both domains in orange and blue, respectively.
#### 3.3.1 Local finite element error estimates

To prepare the proof to Theorem 2.2 we state several auxiliary estimates. We start by citing an a priori error estimate for the local finite element error. Applied to each patch it holds

\[\|\nabla(u-u_{h})\|_{\mathcal{P}}\leq c\Big{(}h^{r}\|u\|_{r+1,\tilde{\mathcal{P} }}+d(\mathcal{P},\tilde{\mathcal{P}})^{-1}\|u-u_{h}\|_{\tilde{\mathcal{P}}} \Big{)}. \tag{13}\]

See (26, Theorem 5.1) for details and the proof. The second term in (13) is a local \(L^{2}\) error term which can be further analyzed as follows. For the Poisson problem we consider the adjoint solution \(z\in H^{1}_{0}(\Omega)\)

\[-\Delta z=\frac{u-u_{h}}{\|u-u_{h}\|}\chi_{\tilde{\mathcal{P}}}\text{ in }\Omega,\;z=0\text{ on }\partial\Omega,\]

where \(\chi_{\tilde{\mathcal{P}}}\) is the characteristic function of \(\tilde{\mathcal{P}}\). Hereby we get

\[\|u-u_{h}\|_{\tilde{\mathcal{P}}}\leq ch^{r}\|u\|_{r+1}\cdot h\|\chi_{\tilde{ \mathcal{P}}}\|\leq ch^{r+1}h_{\tilde{\mathcal{P}}}^{\frac{d}{2}}\|u\|_{r+1},\]

as \(|\tilde{\mathcal{P}}|^{\frac{1}{2}}=\mathcal{O}(h_{\tilde{\mathcal{P}}}^{ \frac{d}{2}})\). Combination with (13) and noting that \(d(\mathcal{P},\tilde{\mathcal{P}})=\mathcal{O}(h_{\tilde{\mathcal{P}}})\) gives

\[\|\nabla(u-u_{h})\|_{\mathcal{P}}\leq ch^{r}\big{(}\|u\|_{r+1,\tilde{\mathcal{ P}}}+h\cdot h_{\mathcal{P}}^{\frac{d}{2}-1}\|u\|_{r+1}\big{)}. \tag{14}\]

In addition, we will require local finite element error estimates in Sobolev norms with negative index.

Lemma 1 (Local finite element estimates in negative norms): _It holds_

\[\|u-u_{h}\|_{1-r,\mathcal{P}}\leq ch^{2r}\Big{(}\|u\|_{r+1,\tilde{\mathcal{P}} }+h\cdot h_{\mathcal{P}}^{\frac{d}{2}-1}\|u\|_{r+1}\Big{)}. \tag{15}\]

Proof: Let \(\psi\in H^{r-1}(\mathcal{P})\) and \(E\psi\in H^{r-1}(\tilde{\mathcal{P}})\) be its natural extension to a smooth domain \(\hat{\mathcal{P}}\) satisfying \(\mathcal{P}\subset\subset\hat{\mathcal{P}}\subset\subset\tilde{\mathcal{P}}\). Further, let \(z\in H^{r+1}(\tilde{\mathcal{P}})\cap H^{1}_{0}(\tilde{\mathcal{P}})\) be the solution to the adjoint problem

\[(\phi,E\psi)=(\nabla\phi,\nabla z)_{\tilde{\mathcal{P}}}\quad\forall\phi\in H ^{1}_{0}(\tilde{\mathcal{P}})\]

and therefore it satisfies \(\|z\|_{r+1}\leq c\|E\psi\|_{r-1,\tilde{\mathcal{P}}}\leq c\|\psi\|_{r-1, \mathcal{P}}\). Then, we have

\[(u-u_{h},\psi)_{\mathcal{P}}\leq\|\nabla(u-u_{h})\|_{\mathcal{P}}ch^{r}\|\psi \|_{r-1,\mathcal{P}},\]

which, together with (14) gives (15).
#### 3.3.2 Localized finite element solutions

The comparison of the solution on the application mesh \(\Omega_{h}\) with training data obtained on \(\Omega_{h}^{tr}\) will be by means of local problems that are defined only on the surrounding of a single patch.

Definition 4 (Local problems): Let \(\mathcal{P}\subset\subset\tilde{\mathcal{P}}\) be a patch and a slightly enlarged domain matching the mesh \(\Omega_{h}\). By \(V_{h}(\tilde{\mathcal{P}})\), \(V_{H}(\tilde{\mathcal{P}})\) and \(V_{h,0}(\tilde{\mathcal{P}})\) we denote local finite element spaces, \(V_{h,0}\) having zero boundary data on \(\partial\tilde{\mathcal{P}}\). For \(u_{H}\in V_{H}(\tilde{\mathcal{P}})\) we define \(v_{h}\in V_{h,0}(\tilde{\mathcal{P}})\) via

\[(\nabla(u_{H}+v_{h}),\nabla\phi_{h})_{\tilde{\mathcal{P}}}=(f,\phi_{h})\quad \forall\phi_{h}\in V_{h,0}(\tilde{\mathcal{P}}). \tag{16}\]

Lemma 2 (Local problems): _Let \(u_{h}\in V_{h}\) and \(u_{H}\in V_{H}\) be solutions to_

\[(\nabla u_{h},\phi_{h})=(f,\phi_{h})\quad\forall\phi_{h}\in V_{h},\quad(\nabla u _{H},\phi_{H})=(f,\phi_{H})\quad\forall\phi_{h}\in V_{H}. \tag{17}\]

_For the solution to the local problem \(v_{h}\in V_{h,0}(\tilde{\mathcal{P}})\) with finite elements of degree \(r\in\mathbb{N}\) it holds_

\[\|\nabla v_{h}\|_{\tilde{\mathcal{P}}} \leq\|\nabla(u_{h}-u_{H})\|_{\tilde{\mathcal{P}}} \tag{18}\] \[\|v_{h}\|_{1-r,\tilde{\mathcal{P}}} \leq c\big{(}h^{r}+H^{r})\|\nabla(u_{h}-u_{H})\|_{\tilde{\mathcal{P }}}. \tag{19}\]

Proof: Testing (16) with \(\phi_{h}=v_{h}\) and inserting \(\pm(\nabla u_{h},\nabla v_{h})_{\tilde{\mathcal{P}}}\) gives

\[\|\nabla v_{h}\|_{\tilde{\mathcal{P}}}^{2}=(f,v_{h})_{\tilde{\mathcal{P}}}-( \nabla u_{h},\nabla v_{h})_{\tilde{\mathcal{P}}}+(\nabla(u_{H}-u_{h}),\nabla v _{h})_{\tilde{\mathcal{P}}}.\]

The first term is zero using (17), as \(v_{h}\) can be extended to \(\hat{v}_{h}\in V_{h}\) by zero outside of \(\tilde{\mathcal{P}}\). Estimating with Cauchy-Schwarz and dividing by \(\|\nabla v_{h}\|_{\tilde{\mathcal{P}}}\) gives the energy norm estimate (18).

Next, for \(\psi\in H^{r-1}(\tilde{\mathcal{P}})\) we define the adjoint solution \(z\in H^{1}_{0}(\tilde{\mathcal{P}})\)

\[(\phi,\psi)=(\nabla\phi,\nabla z)\quad\forall\phi\in H^{1}_{0}(\tilde{ \mathcal{P}}) \tag{20}\]

For \(\phi:=v_{h}\in V_{h,0}(\tilde{\mathcal{P}})\subset H^{1}_{0}(\tilde{\mathcal{P}})\) it holds by means of (16) for all \(z_{h}\in V_{h,0}(\tilde{\mathcal{P}})\)

\[(v_{h},\psi) =(\nabla v_{h},\nabla z)_{\tilde{\mathcal{P}}}-\Big{(}(\nabla v_ {h},\nabla z_{h})_{\tilde{\mathcal{P}}}+(\nabla u_{H},\nabla z_{h})_{\tilde{ \mathcal{P}}}-(f,z_{h})_{\tilde{\mathcal{P}}}\Big{)}\] \[=\big{(}\nabla v_{h},\nabla(z-z_{h})\big{)}_{\tilde{\mathcal{P}}} -\Big{(}(f,z_{h})_{\tilde{\mathcal{P}}}-(\nabla u_{H},\nabla z_{h})_{\tilde{ \mathcal{P}}}\Big{)}\] \[=\big{(}\nabla v_{h},\nabla(z-z_{h})\big{)}_{\tilde{\mathcal{P}}} -\Big{(}(\nabla u_{h},\nabla z_{h})_{\tilde{\mathcal{P}}}-(\nabla u_{H},\nabla z _{h})_{\tilde{\mathcal{P}}}\Big{)}\] \[=\big{(}\nabla v_{h},\nabla(z-z_{h})\big{)}_{\tilde{\mathcal{P}}} -(\nabla(u_{h}-u_{H}),\nabla(z_{h}-z_{H}))_{\tilde{\mathcal{P}}},\]

where we used that \(u_{h}\) is solution to (17) and also Galerkin orthogonality with respect to \(V_{H}\) and \(V_{h}\). Hence, taking \(z_{h}=I_{h}z\in V_{h,0}(\tilde{\mathcal{P}})\) and \(z_{H}=I_{H}z\in V_{H,0}(\tilde{\mathcal{P}})\) as the interpolations

\[\big{|}(v_{h},\psi)_{\tilde{\mathcal{P}}}\big{|}\leq\|\nabla v_{h}\|_{\tilde{ \mathcal{P}}}ch^{r}\|\psi\|_{r-1,\tilde{\mathcal{P}}}+\|\nabla(u_{h}-u_{H})\|_{ \tilde{\mathcal{P}}}c_{i}\big{(}h^{r}+H^{r}\big{)}\|\psi\|_{r-1,\tilde{\mathcal{ P}}}.\]

Using the energy norm estimate (18) and taking the supremum over \(\psi\in H^{r-1}(\tilde{\mathcal{P}})\) gives the estimate.

After these preparations we conclude with the proof to the main theorem.
Proof (Proof of Theorem 2.2): As in the single-patch case we first introduce the fine mesh solution \(u_{h}\in V_{h}\) to the right hand side \(f\in L^{2}(\Omega)\)

\[\|\nabla(u-u_{\mathcal{N}})\|\leq Ch^{r}\|f\|_{r-1}+\|\nabla(u_{h}-u_{\mathcal{N }})\|. \tag{21}\]

As \(u_{\mathcal{N}}=u_{H}+\sum P_{\mathcal{P}}\mathcal{N}(f)\) is composed of local updates we, from here on, consider the local contributions \(\|\nabla(u_{h}-u_{\mathcal{N}})\|_{\mathcal{P}}\). On each patch \(\mathcal{P}\) we introduce the following local solutions to local problems given by Definition 4

\[\begin{split}& v_{h}\in V_{h,0}(\tilde{\mathcal{P}})\qquad(\nabla(u_{ H}+v_{h}),\nabla\phi_{h})_{\tilde{\mathcal{P}}}=(f,\phi_{h})_{\tilde{ \mathcal{P}}}\qquad\forall\phi_{h}\in V_{h,0}(\tilde{\mathcal{P}})\\ & v_{h}^{i}\in V_{h,0}(\tilde{\mathcal{P}})\qquad(\nabla(u_{H}^{i }+v_{h}^{i}),\nabla\phi_{h})_{\tilde{\mathcal{P}}}=(f_{i},\phi_{h})_{\tilde{ \mathcal{P}}}\qquad\forall\phi_{h}\in V_{h,0}(\tilde{\mathcal{P}})\end{split} \tag{22}\]

The index \(i\) will refer to an element of the training data. The first local problem in (22) is to estimate the local error between \(u_{h}\) and \(u_{H}^{h}:=u_{H}+v_{h}\) on the application domain \(\Omega\), whereas the second local problem estimates the local error between the training data \(u_{h}^{i}\) and \(u_{H}^{h,i}:=u_{H}^{i}+v_{h}^{i}\) on the training mesh. We split the error as

\[\|\nabla(u_{h}-u_{\mathcal{N}})\|_{\mathcal{P}}\leq\|\nabla(u_{h} -u_{H}^{h})\|_{\mathcal{P}}+\|\nabla(u_{H}^{h}-u_{H}^{h,i})\|_{\mathcal{P}}\\ +\|\nabla(u_{h}^{i}-u_{H}^{h,i})\|_{\mathcal{P}}+\|\nabla(u_{h}^{i }-u_{\mathcal{N}})\|_{\mathcal{P}} \tag{23}\]

First and third term can be estimated by (29, Theorem 5.2), where local finite element algorithms are analyzed that decompose the solution \(u_{h}\) into a global coarse solution \(u_{H}\in V_{H}\subset V_{h}\) and into a local fine mesh solution. However, we refine their proof and basically only use (29, Lemma 1), introduce \(\pm u_{H}\) and use Lemma 2 above to get

\[\begin{split}\|\nabla(u_{h}-u_{H}^{h})\|_{\mathcal{P}}& \leq c\|u_{h}-u_{H}^{h}\|_{1-r,\tilde{\mathcal{P}}}\\ &\leq c\Big{(}\|u_{h}-u_{H}\|_{1-r,\tilde{\mathcal{P}}}+H^{r}\| \nabla(u_{h}-u_{H})\|_{\tilde{\mathcal{P}}}\Big{)}.\end{split} \tag{24}\]

We insert \(\pm u\) and using the local error estimates (14) (on each \(\mathcal{P}\) making up \(\tilde{\mathcal{P}}\)) and the local negative norm estimate (15) to bound

\[\|\nabla(u_{h}-u_{H}^{h})\|_{\mathcal{P}}\leq cH^{2r}\big{(}\|u\|_{r+1,\tilde{ \mathcal{P}}}+H\cdot h_{\mathcal{P}}^{\frac{d}{2}-1}\|u\|_{r+1}\big{)}. \tag{25}\]

Likewise, for \(u_{h}^{i}-u_{H}^{h,i}\) in (23) we get

\[\|\nabla(u_{h}^{i}-u_{H}^{h,i})\|_{\mathcal{P}}\leq cH^{2r}\big{(}\|u^{i}\|_{r+ 1,\tilde{\mathcal{P}}}+H\cdot h_{\mathcal{P}}^{\frac{d}{2}-1}\|u^{i}\|_{r+1, \Omega^{tr}}\big{)}. \tag{26}\]

Next, considering the definition of (22), the second term in (23) can be estimated as an error that measures the richness of the training data

\[\|\nabla(u_{H}^{h}-u_{H}^{h,i})\|_{\mathcal{P}}\leq\|\nabla(u_{H}-u_{H}^{i}) \|_{\tilde{\mathcal{P}}}+\|f-f_{i}\|_{-1,\tilde{\mathcal{P}}} \tag{27}\]

The last term in (23) is split into

\[\|\nabla(u_{h}^{i}-u_{\mathcal{N}})\|_{\mathcal{P}}\leq\|\nabla(u_{h}^{i}-u_ {\mathcal{N}}^{i})\|_{\mathcal{P}}+\|\nabla(u_{\mathcal{N}}^{i}-u_{\mathcal{N }})\|_{\mathcal{P}}, \tag{28}\]

where \(u_{\mathcal{N}}^{i}:=u_{H}^{i}+\sum_{\mathcal{P}}P_{\mathcal{P}}\mathcal{N}(f_{ i})\) is the local network approximation for the training data element \(u_{H}^{i}\). The first term

\[\epsilon_{net}^{i}:=\|\nabla(u_{h}^{i}-u_{\mathcal{N}}^{i})\|_{\mathcal{P}} \tag{29}\]
depends on the expressitivity of the neural network and the optimization error. The second term of (28) is estimated as

\[\|\nabla(u_{\mathcal{N}}^{i}-u_{\mathcal{N}})\|_{\mathcal{P}}\leq\|\nabla(u_{H}^{ i}-u_{H})\|_{\mathcal{P}}+c\|\nabla(\mathcal{N}(f_{i})-\mathcal{N}(f))\|_{ \mathcal{P}} \tag{30}\]

and again consists of the data error \(\|\nabla(u_{H}-u_{H}^{i})\|_{\mathcal{P}}\) and, finally, the local network generalization error.

We combine (23)-(30) and sum over all patches to get

\[\|\nabla(u_{h}-u_{\mathcal{N}})\|^{2}\leq c\Big{[}\big{(}h^{2r}+H^ {4r}\big{)}\|f\|_{r-1}^{2}\\ +\sum_{\mathcal{P}\in\Omega_{h}}\min_{i}\Big{\{}\|f-f_{i}\|_{-1, \mathcal{\bar{P}}}^{2}+\|\nabla(u_{H}-u_{H}^{i})\|_{\mathcal{\bar{P}}}^{2}+( \epsilon_{net}^{i})^{2}+H^{4r}\|u^{i}\|_{r+1,\mathcal{\bar{P}}}^{2}\\ +H^{4r+2}\cdot h_{\mathcal{P}}^{d-2}\big{(}\|u\|_{r+1,\Omega}^{2 }+\|u^{i}\|_{r+1,\Omega^{tr}}^{2}\big{)}+\|\nabla(\mathcal{N}(f_{i})-\mathcal{ N}(f))\|_{\mathcal{P}}^{2}\Big{\}}\Big{]}\]

Taking the square root and pulling the \(\|u^{i}\|_{H^{r+1}(\mathcal{\bar{P}})}\leq\|f^{i}\|_{r-1,\Omega^{tr}}\) term out of the minimum we get

\[\|\nabla(u_{h}-u_{\mathcal{N}})\|\] \[+\max_{i}\Big{\{}\big{(}1+H\cdot h_{\mathcal{P}}^{\frac{d}{2}-1} \sqrt{N_{\mathcal{P}}}\big{)}H^{2r}\|f_{i}\|_{r-1,\Omega^{tr}}\Big{\}}\] \[+\!\left(\sum_{\mathcal{P}\in\Omega_{h}}\min_{i}\Big{\{}\|f-f_{i} \|_{\mathcal{P}}^{2}+\|\nabla(u_{H}-u_{H}^{i})\|_{\mathcal{P}}^{2}+\|\nabla( \mathcal{N}(f_{i})-\mathcal{N}(f))\|_{\mathcal{P}}^{2}\Big{\}}\right)^{\frac{ 1}{2}}\Big{]}\!+\!\epsilon_{net},\]

where \(N_{\mathcal{P}}\) is the number of patches. Even if each coarse mesh element is chosen as single patch, it holds \(N_{\mathcal{P}}=\mathcal{O}(H^{-d})\) and \(h_{\mathcal{P}}=H\) such that \(H^{\frac{d}{2}}\sqrt{N_{\mathcal{P}}}=\mathcal{O}(1)\) and the a priori estimate follows.

Remark 1 (Extension to more general settings): Assumption 3.3 gives little freedom in the generalization of the domain \(\Omega^{tr}\). Basically, all domains have to be put together by blocks of patches that are found in the training data. The path for a further generalization would be by means of introducing a parametric setup: data is still kept local on the patches \(\mathcal{P}\), but the training of the network and the application of the network is always by means of transformation to a reference patch \(\mathcal{P}_{r}\). In setting up the neural network approach, the training data then must include a sufficient variety of different patch geometries and sizes that are also found in the application domain \(\Omega\). This approach would also be a first step towards an application to locally refined finite element methods.

### Stability of the neural network

We continue with a stability result which will be used to further analyze the error estimates from Theorems 1 and 2.

Lemma 3 (Network stability): _Let the activation function \(\sigma:\mathbb{R}\rightarrow\mathbb{R}\) of the MLP (see Def. 3) satisfy_

\[|\sigma(y)-\sigma(y_{i})|\leq c_{0}|y-y_{i}|\]
_with \(c_{0}>0\). Then, for each patch \(\mathcal{P}\), for the inputs \(y\) and \(y^{f_{i}}\) and the corresponding finite element functions \(w_{\mathcal{N}}\) and \(w_{\mathcal{N}}^{f_{i}}\), obtained from the network updates the following inequality holds_

\[\|\nabla(w_{\mathcal{N}}-w_{\mathcal{N}}^{f_{i}})\|_{\mathcal{P}}\leq c\cdot c_ {W}\cdot c_{0}^{L}\cdot h^{\frac{d}{2}-1}\|y-y^{f_{i}}\|_{2} \tag{31}\]

_where_

\[c_{W}:=\prod_{j=1}^{L}\|W^{j}\|_{2}. \tag{32}\]

Proof: There exist constants \(c_{1},c_{2}>0\) independent of \(h_{\mathcal{P}}\) and \(n_{\mathcal{P}}\), which is the number of coefficients in each patch, such that

\[c_{1}\|v\|_{\mathcal{P}}^{2}\leq\frac{h_{\mathcal{P}}^{d}}{n_{\mathcal{P}}}\|v \|_{l^{2}(\mathcal{P})}^{2}\leq c_{2}\|v\|_{\mathcal{P}}^{2}\quad\forall v\in V _{\mathcal{P}}.\]

The number of fine mesh coefficients in each patch \(n_{\mathcal{P}}\) scales like \(n_{\mathcal{P}}=\mathcal{O}(h^{-d}\cdot h_{\mathcal{P}}^{d})\), hence

\[c_{1}\|v\|_{\mathcal{P}}^{2}\leq h^{d}\|v\|_{l^{2}(\mathcal{P})}^{2}\leq c_{2 }\|v\|_{\mathcal{P}}^{2}\quad\forall v\in V_{\mathcal{P}}. \tag{33}\]

Together with the inverse inequality and (33), for each \(v\in V_{h}\) it holds

\[\|\nabla v\|_{\mathcal{P}}^{2}\leq c_{inv}^{2}h^{-2}\|v\|_{\mathcal{P}}^{2} \leq ch^{d-2}\|v\|_{\ell^{2}(\mathcal{P})}^{2}.\]

Using the definition of the network gives

\[\|w_{\mathcal{N}}-w_{\mathcal{N}}^{f_{i}}\|_{\ell^{2}(\mathcal{P})}=\|z_{L}(y )-z_{L}(y^{f_{i}})\|_{2} \tag{34}\]

where \(z_{i}(y)=(l_{i}\circ\sigma\circ l_{i-1}\circ\cdots\circ\sigma\circ l_{1})(y)\) and \(l_{i}\) are as defined in (8). By using the definition of \(z_{j}\) and from the assumption we obtain for an arbitrary layer \(j\)

\[\|z_{j}(y)-z_{j}(y^{f_{i}})\|_{2} =\|W^{j}((\sigma\circ z_{j-1})(y)-(\sigma\circ z_{j-1})(y^{f_{i}} ))\|_{2} \tag{35}\] \[\leq\|W^{j}\|_{2}\|(\sigma\circ z_{j-1})(y)-(\sigma\circ z_{j-1}) (y^{f_{i}})\|_{2}\] \[\leq c_{0}\|W^{j}\|_{2}\cdot\|z_{j-1}(y)-z_{j-1}(y^{f_{i}})\|_{2}.\]

Then, by applying (35) recursively to the last layer we obtain

\[\|z_{L}(y)-z_{L}(y^{f_{i}})\|_{2}\leq c_{0}^{L}\prod_{i=1}^{L}\|W_{j}\|_{2} \cdot\|y-y^{f_{i}}\|_{2}.\]

Hence, by applying it to (34) we arrive to

\[\|\nabla(w_{\mathcal{N}}-w_{\mathcal{N}}^{f_{i}})\|_{\mathcal{P}}\leq c_{0}^{ L}h^{\frac{d}{2}-1}\prod_{j=1}^{L}\|W^{j}\|_{2}\cdot\|y-y^{f_{i}}\|_{2}.\]

Note that more superior bounds can be obtained by using methods based on a relaxation to a polynomial optimization problem [16]. As next we show how the difference between inputs are related to the differences of corresponding source terms and the coarse finite element solutions.
Lemma 4: _Given input vectors \(y\) and \(y^{f_{i}}\) in form of (9) together with \(f,f_{i}\in C(\Omega)\) and \(u_{H},u_{H}^{f_{i}}\in V_{H}\). It holds_

\[\left\|y-y^{f_{i}}\right\|_{2}\leq c(h^{-\frac{d}{2}}\|u_{H}-u_{H}^{f_{i}}\|_{ \mathcal{P}}+\|f-f_{i}\|_{\ell^{2}(\mathcal{P})})\]

Proof: Per definition of \(y\) and \(y^{f_{i}}\) we have

\[\left\|y-y^{f_{i}}\right\|_{2}^{2}=\left\|u_{H}-u_{H}^{f_{i}}\right\|_{\ell^{2 }(\mathcal{P})}^{2}+\|f-f_{i}\|_{\ell^{2}(\mathcal{P})}^{2}.\]

The claim follows by applying inequality (33) to the first term above.

By combining Lemmas 3 and 4 we obtain the following Corollary.

Corollary 1: _Under the same assumptions of Lemma 3 it holds_

\[\|\nabla(w_{\mathcal{N}}-w_{\mathcal{N}}^{f_{i}})\|_{\mathcal{P}}\leq c_{ \mathcal{N}}(h^{-1}\|u_{H}-u_{H}^{f_{i}}\|_{\mathcal{P}}+h^{\frac{d}{2}-1}\|f- f_{i}\|_{\ell^{2}(\mathcal{P})})\]

_where \(c_{\mathcal{N}}=\mathcal{O}(c_{W}\cdot c_{0}^{L})\)._

As next, an immediate consequence of previous results is presented.

Corollary 2 (A priori error estimate for the hybrid solution (single-patch case)): _Under the assumptions of Theorem 1 and Corollary 1 it holds_

\[\|\nabla(u-u_{\mathcal{N}})\| \leq c\bigg{[}h^{r}\|f\|_{r-1}+\epsilon\] \[+\min_{f_{i}\in\mathcal{F}}\Bigl{\{}(1+c_{\mathcal{N}}c_{P}h^{-1} )\|f-f_{i}\|_{-1}+c_{\mathcal{N}}h^{\frac{d}{2}-1}\|f-f_{i}\|_{\ell^{2}(\Omega )}\Bigr{\}}\bigg{]}\]

_where \(c_{\mathcal{N}}=\mathcal{O}(c_{W}\cdot c_{0}^{L})\) and \(c_{P}\) is the Poincare constant._

Proof: The claim follows from Theorem 1 and Corollary 1, where we have used the Poincare inequality and stability of the weak Poisson problem

\[h^{-1}\|u_{H}-u_{H}^{f_{i}}\|\leq c_{P}h^{-1}\|\nabla(u_{H}-u_{H}^{f_{i}})\| \leq c_{P}h^{-1}\|f-f_{i}\|_{-1}\]

Likewise, for the multi-patch case the following result is obtained, which we state without proof.

Corollary 3 (A priori error estimate for the hybrid solution (multi-patch case)): _Under the assumptions of Theorem 2 and Corollary 1 it holds_

\[\|\nabla(u_{h}-u_{\mathcal{N}})\| \leq c\Big{[}(h^{r}+H^{2r})\|f\|_{r-1}+H^{2r}\max_{i}\|f_{i}\|_{r- 1,\Omega^{tr}}\] \[\qquad+\Big{(}\sum_{\mathcal{P}\in\Omega_{h}}\min_{i}\left\{\|f- f_{i}\|_{-1,\tilde{\mathcal{P}}}^{2}+c_{\mathcal{N}}^{2}h^{d-2}\|f-f_{i}\|_{\ell^{2} (\mathcal{P})}^{2}\right.\] \[\qquad+\left\|\nabla(u_{H}-u_{H}^{i})\right\|_{\tilde{\mathcal{P }}}^{2}+c_{\mathcal{N}}^{2}h^{-2}\|u_{H}-u_{H}^{i}\|_{\mathcal{P}}^{2}\Big{\}} \Big{]}+\epsilon\]

_where \(c_{\mathcal{N}}=\mathcal{O}(c_{W}\cdot c_{0}^{L})\)._
The error estimates from Theorems 1 and 2, as well as the two corollaries, show a balanced error. In particular, they prove the practical usefulness of the method: the total error is given as a balance of richness of the training data, i.e., number and distribution of the training elements, and quality of the training data, i.e., the resolution of the training data. Both can be controlled. In addition, there is the network and optimization error, which in principle can be influenced by the architecture, depth and width of the network. The practical efficiency of the method will always depend on the specific example and especially on how much the effort of offline phase and online phases differ. For typical problems in fluid mechanics, the 3d flow around obstacles, we observed substantial increases in efficiency for relevant generalizations of the training data [20]. Numerically we could not yet identify the role of the coarse mesh contribution \(O(H^{2r})\) that suggests the relation \(h=H^{2}\) between the two grid levels.

## 4 Numerical experiments

In the following paragraphs we will document different numerical simulations to explore the performance of the hybrid finite element neural network approach. We will start by an in-depth analysis of the algorithm and the sensitivity of the approximation properties on the various aspects like network size, variety of training data or optimization and regularization procedure. Then, we add a second test case to explore the generalization capacity of the approach.

### Configuration of the test case

We start by describing our experimental setup. Let \(\Omega=(0,1)^{2}\) be the unit square. We consider the two-dimensional Poisson equation with homogeneous Dirichlet boundary conditions

\[-\Delta u=f\ \text{in}\ \Omega,\quad u=0\ \text{on}\ \partial\Omega. \tag{36}\]

We consider right hand sides from the following set of functions

\[\mathcal{F}\ :=\ \Big{\{}f(x_{1},x_{2})\ =\ \sin(2\pi\,\cdot\,C_{1}(x_{1}\,+ \,C_{2}))\,\cdot\,\sin(2\pi\,\cdot\,C_{3}(x_{2}\,+\,C_{4}))\Big{\}} \tag{37}\]

and training data is declared by picking random and uniformly distributed samples from this set. For each such \(f\in\mathcal{F}\) we will then compute finite element solutions \(u_{h}\) and \(u_{H}\) as training data.

### Neural network setup and optimization

The neural network is a multilayer perceptron as described previously in Section 3, see Definition 3. We consider networks with \(L-1=4\) hidden layers, each having a width of \(512\), if not stated otherwise. The total number of trainable parameters depends on the size of the input and output vectors, see Table 1. In order to perform \(k\) refinements, the network receives 4 nodal values of the coarse solution
and \((2^{k}+1)^{2}\) nodal values of the source term. The output vector contains \((2^{k}+1)^{2}\) nodal values which predicts the difference between fine and coarse finite element solutions.

#### 4.2.1 Generation of training data

For creating training and test data sets we randomly generate finite sets of functions taken from \(\mathcal{F}\) as given in (37). To generate such a function it suffices to generate four real numbers, namely \(C_{1},C_{2},C_{3}\) and \(C_{4}\). Since \(C_{2}\) and \(C_{4}\) are the phase shifts of respective terms in \(f\) divided by \(2\pi\) we sample them uniformly at random from [0.0,1.0]. \(C_{1}\) and \(C_{3}\) are frequencies of respective terms, so restrict them to be only from [1.0,1.5].

The computations are performed for a coarse step size \(H\) and for the fine mesh sizes \(h=H/2\), \(h=H/4\) and \(h=H/8\). Hereby, we generate sets of training and test triples \((f_{i},u_{H}^{i},u_{h}^{i})\in\mathcal{F}\times V_{H}\times V_{h}\) for \(i=1,\ldots,N_{T}\). In our experiments, we set \(H=2^{-3}\) which corresponds to 64 cells.

Due to the local neural network setup, each triple \((f_{i},u_{H}^{i},u_{h}^{i})\) generates a larger number of training patches \(N_{\mathcal{P}}\), one for each patch of the triangulation. The number of patches are always the same as the number of elements in the coarse mesh. Hence, \(H=h_{\mathcal{P}}\) and we have \(N_{\mathcal{P}}=64\).

A second test data set of the same size is additionally chosen from \(\mathcal{F}\).

#### 4.2.2 Optimization

For the training of the network we use the Adam optimizer [15]. The loss function is the mean square error divided by the number of patches and the number of training data as defined in (10). As Lemma 3 shows, stability of the network is tightly related with the multiplications of spectral norms of the weights. Therefore, we also perform experiments with a modified loss function by including a regularization term which we report in Section 4.4.

### Accuracy of the hybrid finite element neural network solver

In this first test case, we analyze the accuracy of the hybrid solver. For this purpose, we consider the mean error over all test data. Starting from the coarse solution with grid size \(H=2^{-3}\), the hybrid method is enriched with neural networks. In doing so, we predict 1, 2 or 3 grid levels, thus trying to achieve accuracies \(h=H/2\), \(h=H/4\) and \(h=H/8\).

\begin{table}
\begin{tabular}{c|c|c|c|c} \hline H/h & \(N_{0}\) & \(N_{1,\ldots,4}\) & \(N_{5}\) & \# parameters \\ \hline \(2^{2}\) & 4 + 25 & 512 & 25 & 816 153 \\ \(2^{3}\) & 4 + 81 & 512 & 81 & 873 553 \\ \(2^{4}\) & 4 + 289 & 512 & 289 & 1 086 753 \\ \hline \end{tabular}
\end{table}
Table 1: The number of neurons \(N_{l}\) at each layer \(l\in\{0,\ldots,5\}\) and the total number of trainable parameters (weights and biases) \(\sum_{l=1}^{L=5}N_{l-1}\cdot N_{l}+N_{l}\).
In this first test, we only want to investigate the approximation ability of the neural network and we use an excessive number of training data to keep the error term \(\min_{i}\left\{\|f-f_{i}\|+\|\nabla(u_{H}-u_{H}^{i})\|\right\}\) small. Thus, from (12), considering linear finite elements \(r=1\), it remains

\[\|\nabla(u-u_{\mathcal{N}})\|=\mathcal{O}\big{(}h+H^{2}\big{)}\]

Since we do not know the analytical solution we use a reference solution instead, which is calculated on a mesh with element size \(h/4\). The numerical results shown in Table 2 confirm the estimate and the hybrid neural network solution recovers the accuracy of the fine mesh solution.

#### 4.3.1 Dependency on the training data set

In Fig. 4 we show the approximation error of the hybrid simulation for the three prediction levels \(h=H/2\), \(h=H/4\) and \(h=H/8\) depending on the size of the training data set \(N_{T}\). The results indicate that the corresponding fine mesh

\begin{table}
\begin{tabular}{c c c} \hline \hline \(h\) & \(\|u-u_{h}\|\) & \(\|u-u_{\mathcal{N}}\|\) \\ \hline \(H/2\) & 1.25e-04 & 1.25e-04 \\ \(H/4\) & 3.21e-05 & 3.29e-05 \\ \(H/8\) & 8.13e-06 & 1.11e-05 \\ \hline \hline \end{tabular} 
\begin{tabular}{c c} \hline \hline - \(\blacktriangle\)- & Reference \(\|u-u_{h}\|\) & Hybrid \(\|u-u_{\mathcal{N}}\|\) \\ \hline \hline \end{tabular}
\end{table}
Table 2: Performance of the hybrid solution for three different refinement levels. We show the average error on the test data. The training data set \(\mathcal{F}\) is sufficiently big such that the data error is negligible.

Figure 4: Dependency of the prediction quality \(\|u-u_{\mathcal{N}}\|\) (for the test data) on the size of the training data set.
accuracy is reached in all three cases, however, only for an increase amount of training data. This is in accordance to the error estimate (12), which shows that the error term \(H^{2r}+h^{r}\) must be balanced with the data error \(\min_{i}\|f-f_{i}\|\) and \(\min_{i}\|\nabla(u_{H}-u_{H}^{i})\|\). Furthermore, Fig. 4 suggests the scaling \(\min_{i}\|f-f_{i}\|=\mathcal{O}(N_{T}^{-\frac{1}{2}})\).

Figure 5 indicates the constant \(c_{W}\) (32), see Lemma 3, and its dependency on the target step size \(h\) as well as on the amount of training data \(N_{T}\). Growth of this constant is slower for finer target step sizes. All in all, a linear dependency of \(c_{W}\) on \(N_{T}\) appears eminent. This result is not desirable as \(c_{W}\) enters the error estimate of Theorem 3.1 in the same way as the data error \(\min\|f-f_{i}\|\). This data error must decrease faster than the growth of the constant.

#### 4.3.2 Error depending on the network complexity

Next we study the dependency of the error on the complexity of the network. Again we consider the cases \(h=H/2,H/4\) and \(H/8\). We used \(N_{T}=2^{12}\) training problems for this experiment. First we choose the number of hidden layers of the perceptron between 1, 2, 4 or 8 with 512 neurons in each layer.

The results can be seen in Figure 6. We observe that most of the time using four layers seems to be optimal for cases \(h=H/2\) and \(H/4\). For the case \(h=H/8\) we observe that after three layers the network is more prone to overfitting, and while the train error keeps reducing, the test error is not getting any better.

Just like above, let's also consider how the constant \(c_{W}\) (32) from the bound changes as we vary the number of layers. This can be seen in Figure 7. Here we observe that the value of \(c_{W}\) increases with number of layers.

Next, we fix the number of layers to be 4 and vary the number of neurons per layer from 8 to 512. The results can be seen in Figure 8. Here we observe that error behaves similarly to the experiment where we varied the amount of training data - the error decreases as the number of neurons increases. So, we confirm that 512 is most likely indeed the optimal choice in the current setting. Results for the the bound can be seen in Figure 9. Same as above we observe that the value of \(c_{W}\) increases with number of neurons.

Figure 5: Bound vs number of training problems
### Impact of the regularization

Since experiments in the previous sections showed that the bound \(c_{W}\) from (32) grows fast with respect to the training data set size (see Fig. 5) and the network complexity (see Figs. 7 and 9), we aim to penalize it. To this end, we modify our loss function such that the optimization problem (10) is replaced with

\[\min_{\begin{subarray}{c}W_{i},b_{i},\\ i\in\{0,\dots,L\}\end{subarray}}\quad\frac{1}{N_{T}N_{\mathcal{P}}}\sum_{ \mathcal{P}}\lVert z_{\mathcal{P}}-\mathcal{N}(y_{\mathcal{P}})\rVert_{2}^{2} \;+\;\frac{\alpha}{N_{T}N^{param}}\prod_{i=1}^{N_{L}}\lVert W_{i}\rVert_{F} \tag{38}\]

where \(W_{i}\in\mathbb{R}^{N_{i-1}}\times\mathbb{R}^{N_{i}}\) stands for the weights and \(N^{param}\) for the number of parameters. We use the Frobenius norm \(\lVert\cdot\rVert_{F}\) since minimizing the singular norm would be extremely computationally expensive and the Frobenius norm bounds the singular norm.

Figure 6: Dependency of the prediction quality \(\lVert u-u_{\mathcal{N}}\rVert\) (for the test data) on the number of network layers. The size of the training data is chosen as \(N_{T}=2^{12}\) for all cases.

Figure 7: Corresponding to Fig. 6 we show the estimates for the bound for \(c_{W}\) depending on the number of layers.
In Fig. 10 we compare the resulting constant for fixed \(h=H/2\) and varying \(\alpha\). For the sake of comparison we also consider \(\alpha=0\), which means using the loss function without regularization (10). In Fig. 11 we show the effect on the quality of the hybrid approximation. As the choice \(\alpha=10^{-4}\) leads to a smaller error and a smaller bound \(c_{W}\) among other choices, we present further results on the dependency to the training data set size, but now with \(\alpha=10^{-4}\) in the modified loss function (38). These results can be seen in Figure 12 and Figure 13 respectively.

### Role of data preprocessing

Finally, we justify our specific choice of data preprocessing, so-called standardization, which we employed in every experiment presented in this paper. The data

Figure 8: Dependency of the prediction quality \(\|u-u_{\mathcal{N}}\|\) (for the test data) on the number of neurons in hidden layers. The size of the training data is chosen as \(N_{T}=2^{12}\) for all cases.

Figure 9: Corresponding to Fig. 8 we show the estimates for the bound for \(c_{W}\) depending on the number of layers.
preprocessing has a large effect on the optimization of the network and its approximation quality. In order to illustrate this, we consider fixed \(h=H/2\) and we vary the amount of training data from \(2^{12}\) to \(2^{16}\). As for preprocessing methods, we considered no preprocessing, min-max scaling to \([0,1]\) and standardization. These methods were applied to each input and output feature separately.

Min-max scaling works by uniformly rescaling a given feature to the interval \([0,1]\). In other words for a given set of samples \(X=(X_{1},X_{2},\ldots,X_{n})\in\mathbb{R}^{n}\) with minimum values \(X_{\min}\) and maximum value \(X_{\max}\), the min-max normalization of it can be defined as following:

\[\tilde{X}_{min-max}:=\frac{X-X_{\min}}{X_{\max}-X_{\min}}. \tag{39}\]

Figure 11: Approximation errors \(\|u-u_{\mathcal{N}}\|\) of the solutions obtained with regularization factors \(\alpha\in\{0,10^{-4},10^{-2},1\}\).

Figure 10: The bound \(c_{W}\) corresponding to different regularization factors \(\alpha\in\{0,10^{-4},10^{-2},1\}\)
Standardization is similar but instead of rescaling the data to \([0,1]\), it transforms the data to have mean value of \(0\) and variance of \(1\), i.e.

\[\tilde{X}_{std}=\frac{X-\text{Mean}(X)}{\text{Var}(X)}. \tag{40}\]

Figure 14 shows the results where \(y\)-axis corresponds to the average values of the validation metric. As it was mentioned above, validation metric is defined as the \(\ell_{2}\)-norm (6) of difference between given and a reference solution. \(x\)-axis here corresponds to the number of training data.

While the topmost line depicts the validation metric of coarse (input) solutions on the test data set, the bottommost line depicts the validation metric for the fine (target) solutions on both the test data sets. The rest of the lines depicts the validation metric for solutions obtained by using the proposed method with aforementioned preprocessing methods on train and test data sets.

Figure 12: Regularization with \(\alpha=10^{-4}:\) Dependency of the prediction quality \(\|u-u_{\mathcal{N}}\|\) (for the test data) on the size of the training data set.

Figure 13: Estimate for the bound for \(c_{W}\) obtained with regulariza- tion factor \(\alpha=10^{-4}\).
We observe that no preprocessing gives the worst results compared to other methods and the error of the resulting solution is even getting bigger with more training data. The second best preprocessing method turned out to be min-max scaling. It performs pretty well even with small amount of the training data and the error improves with more data.

The best preprocessing method is standardization. In this case, the error of the proposed method reaches the error of the target one slightly faster and can give the comparable error to min-max scaling while using less data for the training. Since standardization performed the best, we have used it in all our experiments. Due to the fact that standardization consists of just a few arithmetic operations performed on both input and output features, the computational overhead of it is negligible.

### Generalizations to different domains

In this section we discuss the application of the network to a domain different from the one, that was used for generation of the training data. In the previous sections we have used \(\Omega_{\rm train}=(0,1)^{2}\) for both training and testing. Now, we leave the experimental settings mostly untouched and after the training was performed on \(\Omega_{\rm train}\) we apply the network to the solutions of Equation 1 on \(\Omega_{\rm test}=(0,2)\times(0,1)\). In this section we consider a problems of the aforementioned form for our further tests. The equation we approximate solutions of is again the Poisson problem

\[-\Delta u=\sin(2\pi\cdot C_{1}(x_{1}+C_{2}))\cdot\sin(2\pi\cdot C_{3}(x_{2}+C_ {4}))\ \mbox{in}\ \Omega,\quad u=0\ \mbox{on}\ \partial\Omega. \tag{41}\]

where \(C_{1}=1.2,C_{2}=0.2,C_{3}=1.4,C_{4}=0.4\). Table 3 shows the results for this test-case. We observe optimal convergence for predicting one \(H/2\) or two mesh levels \(H/4\), but the accuracy is slightly non-optimal if three mesh levels are to be predicted, i.e. in the case \(h=H/8\). Figure 15 shows the solution and the error \(u-u_{\mathcal{N}}\) for the extended domain problem.

Figure 14: Importance of preprocessing: dependency of the prediction quality \(\|u-u_{\mathcal{N}}\|\) (for the test data) on the size of the training data set.
\begin{table}
\begin{tabular}{c c c c} \hline \hline \(h\) & \(\|u-u_{H}\|\) & \(\|u-u_{h}\|\) & \(\|u-u_{\mathcal{N}}\|\) \\ \hline \(H/2\) & 3.83e-04 & 1.00e-04 & 1.00e-04 \\ \(H/4\) & 3.82e-04 & 2.53e-05 & 2.68e-05 \\ \(H/8\) & 3.82e-04 & 6.37e-06 & 1.19e-05 \\ \hline \hline \end{tabular} 
\begin{tabular}{c c} \hline \hline \(\cdot\
## 5 Conclusion

In this paper we have introduce introduced and analyzed a local hybrid finite element neural network method and its application to the Poisson problem. We have employed neural networks to reduce the complexity of solving Poisson equation for a family of source terms. The network operates locally, which enables it to be domain agnostic, meaning that it can be applied to problems on different domains which are not included in the training data.

The most important component of this paper is the theoretical investigation of the proposed method. In particular we performed an a priori error analysis of the local hybrid finite element neural network method including the stability analysis of the network being used. Theoretical findings are accompanied with numerical results, which have shown that the error of the hybrid solution can be controlled by the amount of the training data and by tuning the hyperparameters of the network.

We have shown first generalization results to a slightly modified domain and found good performance. However, generalization of the method will fail if the character of the solution changes. For instance, generalization from the very regular square-domain setting to an L-shaped domain with reentrant corners will not give satisfactory results, as the singular behavior close to the corner is never seen during training of the network. Since neural networks essentially perform an interpolation of the data and processes available in the training, an appropriate enrichment of the training data will be necessary in this case.

Further steps will focus on the extension of the analysis to the time-dependent case, where the hybrid DNN-MG solver has shown very high accuracy in relevant problems [21; 19; 20]. The time-dependent case is appealing as it could enable one to achieve the accuracy of the fine temporal discretizations by calculation solutions from coarse discretizations only and receiving the temporal fluctuations from the network. This will be the focus of an upcoming work.

## Data availability

The Python scripts for reproducing the numerical test cases are published on Zenodo [13].

## Acknowledgements

The authors acknowledge the support of the GRK 2297 MathCoRe, funded by the Deutsche Forschungsgemeinschaft, Grant Number 314838170.

## Conflict of interest

The authors declare that they have no conflict of interest.