Adaptive quadratures for nonlinear approximation of low-dimensional PDEs
  using smooth neural networks
Physics-informed neural networks (PINNs) and their variants have recently
emerged as alternatives to traditional partial differential equation (PDE)
solvers, but little literature has focused on devising accurate numerical
integration methods for neural networks (NNs), which is essential for getting
accurate solutions. In this work, we propose adaptive quadratures for the
accurate integration of neural networks and apply them to loss functions
appearing in low-dimensional PDE discretisations. We show that at opposite ends
of the spectrum, continuous piecewise linear (CPWL) activation functions enable
one to bound the integration error, while smooth activations ease the
convergence of the optimisation problem. We strike a balance by considering a
CPWL approximation of a smooth activation function. The CPWL activation is used
to obtain an adaptive decomposition of the domain into regions where the
network is almost linear, and we derive an adaptive global quadrature from this
mesh. The loss function is then obtained by evaluating the smooth network
(together with other quantities, e.g., the forcing term) at the quadrature
points. We propose a method to approximate a class of smooth activations by
CPWL functions and show that it has a quadratic convergence rate. We then
derive an upper bound for the overall integration error of our proposed
adaptive quadrature. The benefits of our quadrature are evaluated on a strong
and weak formulation of the Poisson equation in dimensions one and two. Our
numerical experiments suggest that compared to Monte-Carlo integration, our
adaptive quadrature makes the convergence of NNs quicker and more robust to
parameter initialisation while needing significantly fewer integration points
and keeping similar training times.
Alexandre Magueresse, Santiago Badia
Adaptive Quadratures for Nonlinear Approximation of Low-Dimensional PDEs using Smooth Neural Networks

###### Abstract.

Physics-informed neural networks (PINNs) and their variants have recently emerged as alternatives to traditional partial differential equation (PDE) solvers, but little literature has focused on devising accurate numerical integration methods for neural networks (NNs), which is essential for getting accurate solutions. In this work, we propose adaptive quadratures for the accurate integration of neural networks and apply them to loss functions appearing in low-dimensional PDE discretisations. We show that at opposite ends of the spectrum, continuous piecewise linear (CPWL) activation functions enable one to bound the integration error, while smooth activations ease the convergence of the optimisation problem. We strike a balance by considering a CPWL approximation of a smooth activation function. The CPWL activation is used to obtain an adaptive decomposition of the domain into regions where the network is almost linear, and we derive an adaptive global quadrature from this mesh. The loss function is then obtained by evaluating the smooth network (together with other quantities, e.g., the forcing term) at the quadrature points. We propose a method to approximate a class of smooth activations by CPWL functions and show that it has a quadratic convergence rate. We then derive an upper bound for the overall integration error of our proposed adaptive quadrature. The benefits of our quadrature are evaluated on a strong and weak formulation of the Poisson equation in dimensions one and two. Our numerical experiments suggest that compared to Monte-Carlo integration, our adaptive quadrature makes the convergence of NNs quicker and more robust to parameter initialisation while needing significantly fewer integration points and keeping similar training times.

\({}^{\dagger}\)School of mathematics, Monash university, Clayton, Victoria 3800, Australia

\({}^{\sharp}\) Centre Internacional de Metodes Numerics a l'Enginveria, Campus Nord, 08034 Barcelona, Spain

_E-mail addresses: alexandre.magueresse@monash.edu, santiago.badia@monash.edu._

_Date:_ March 23, 2023.

_Key words and phrases._ adaptive quadrature, integration, PDE, neural networks.

\({}^{*}\)Corresponding author.

## 1. Introduction

In the last years, the use of Neural Networks (NNs) to approximate partial differential equations (PDEs) has been widely explored as an alternative to standard numerical methods such as the Finite Element Method (FEM). The idea behind these methods is to represent the numerical solution of a PDE as the realisation of a NN. The parameters of the NN are then optimised to minimise a form of the residual of the PDE or an energy functional. The original Physics-Informed Neural Network (PINN) method [1] considers the minimisation of the PDE in the strong form. Other frameworks recast the PDE as a variational problem and minimise the weak energy functional, like in Variational PINNs (VPINNs, [2]), the Deep Ritz Method (DRM, [3]) and the Deep Nitsche Method (DNM, [4]).

**Comparison between the FEM and NNs.** NNs activated by the \(\mathrm{ReLU}\) function (Rectified Linear Unit) constitute a valuable special case, as they can emulate linear finite elements spaces on a mesh that depends on the parameters of the network [5]. This enables one to draw close links between the well-established theory of the FEM and the recently introduced neural frameworks [6].

There are additional benefits to using the ReLU activation function, such as addressing the issue of exploding or vanishing gradients and being computationally cheap to evaluate. However, the numerical experiments in [7] suggest that the theoretical approximation bounds for \(\mathrm{ReLU}\) networks are not met in practice: even though there exist network configurations that realise approximations of a given accuracy, the loss landscape is rough because of the low regularity of \(\mathrm{ReLU}\). As a result, gradient-based optimisation algorithms are largely dependent on the initial parameters of the network and may fail to converge to a satisfying local minimum. It was also shown in [8] that training a network with a smooth activation function helps the network converge.

Compared to the FEM, NNs can be differentiated thanks to automatic differentiation rather than relying on a numerical scheme (and ultimately a discretisation of the domain). Another advantage of NNs is that they are intrinsically adaptive, in the sense that contrary to traditional PDE solvers that rely on a fixed mesh and fixed shape functions, NNs can perform nonlinear approximation. Finally, the machine learning approach proposes a unified setting to handle forward and inverse problems that can involve highly nonlinear differential operators.
However, properly enforcing initial and boundary conditions remains one of the major challenges for NNs [9], especially when the domain has a complex geometry.

**Numerical integration of NNs.** Current efforts in the literature are mainly focused on developing novel model architectures or solving increasingly complex PDEs, while there are comparatively few results on the convergence, consistency and stability of PINNs and related methods. These are the three main ingredients that provide theoretical guarantees on the approximation error, and they remain to be established for NNs.

In particular, the numerical quadratures that are used to evaluate the loss function have mostly been investigated experimentally, and the link between the integration error of the loss function, often referred to as the generalisation gap, and the convergence of the network is still not fully understood. The overwhelming majority of implementations of NNs solvers rely on Monte Carlo (MC) integration because it is relatively easy to implement. However, it is known that its convergence rate with respect to the number of integration of points is far from optimal in low-dimensional settings. For example, the integration error decays as \(n^{-2}\) for the trapezoidal rule in dimension one, whereas it is of the order of \(n^{-1/2}\) for MC regardless of the dimension. Furthermore, local behaviours of the function to integrate are generally poorly handled by MC, because the points are sampled uniformly across the domain.

A few recent works have considered more advanced integration methods for NN solvers, including Gaussian quadratures or quasi-MC integration. Some studies contemplate resampling the MC points during the training or adding more points where the strong, pointwise residual is larger. We refer the reader to [10] for a comprehensive review of sampling methods and pointwise residual-based resampling strategies applied to NNs. They show that the convergence of the model is significantly influenced by the choice of the numerical quadrature. A theoretical decay rate of the generalisation gap when the training points are carefully sampled was shown in [11], while [12] proved the consistency of PINNs and related methods in the sample limit. Finally, the generalisation error was bounded with the sum of the training error and the numerical integration error in [13]. However, this result requires strong assumptions on the solution and the numerical quadratures, and the constants involved in the bound are not explicit and may not be under control.

It was shown in [14] that "quadrature errors can destroy the quality of the approximated solution when solving PDEs using deep learning methods". The authors explored the use of a Continuous Piecewise Linear (CPWL) interpolation of the output of the network and observed better convergence compared to MC integration. Their experiments were performed in dimension one, and the networks were interpolated on fixed, uniform meshes.

**Motivations and contributions.** Our work aligns with current research in the theory and practice of NNs applied to approximating PDEs. We seek to equip NNs with some of the tools of the FEM in terms of numerical integration, for low-dimensional (\(d\leq 3\)) problems arising from computational physics or engineering.

The space of CPWL activation functions is of great interest for NNs owing to its structure of vector space and its stability by composition. This means that whenever a network is activated by a CPWL function, its realisation is CPWL too. Our idea is that provided that we can decompose the domain at hand into regions where the network is almost linear in each cell, we can integrate with high accuracy any integral that involves the network and its partial derivatives by relying on Gauss-like quadratures of a sufficiently high order. The linear and bilinear forms are decomposed on this mesh in a similar way as in the FEM, except that the underlying mesh is adaptive instead of fixed. We provide an algorithm to obtain such a decomposition of the domain based on the expression of the activation function and the parameters of the network.

Nevertheless, as will be shown in Sect. 3, CPWL activations are not smooth enough for the optimisation problem of a NN approximating a PDE to be well-behaved. If the activation function is CPWL, it needs to be regularised so that the network can be properly trained. In this work, all the networks that we train are activated by smooth functions. However, the quadrature points and the mesh are obtained from the projection of this smooth network onto a CPWL space, when the activation function is replaced by a CPWL approximation of the smooth activation. We design a procedure to obtain a CPWL approximation of a smooth function on \(\mathbb{R}\) as a whole.

Our extensive numerical experiments demonstrate that our proposed integration method enables the NN to converge faster and smoother, to lower generalisation errors compared to MC integration, and to be more robust to the initialisation of its parameters while keeping similar training times.

**Outline.** We recall the general setting of PINN and its variants and introduce some notations related to the variational problem we consider in Sect. 2. In Sect. 3, after explaining why NNs activated by CPWL functions are not suitable to approximate PDEs, we suggest several smooth approximations of \(\mathrm{ReLU}\). In Sect. 4, we introduce a class of activation functions that behave linearly at infinity, propose a method to approximate them
by CPWL functions on \(\mathbb{R}\) as a whole, and prove its quadratic convergence rate. Then, in Sect. 5, we provide an algorithm to decompose the physical domain into convex regions where the NN is almost linear. We also describe how to create Gaussian quadratures for convex polygons and analyse the integration error of our method. We present our numerical experiments in Sect. 6. We solve the Poisson equation in dimensions one and two on the reference segment and square, but also on a more complex domain, to show that our adaptive quadrature can adapt to arbitrary polygonal domains. We discuss our findings and conclude our work in Sect. 7. In order to keep the article more readable, we decided to write the proofs of our propositions in Apdx. A and our algorithms in Apdx. B.

## 2. Preliminaries

Consider the boundary condition problem

\[\left\{\begin{array}{ll}\mathcal{D}u=f&\text{in }\Omega\\ \mathcal{B}u=g&\text{on }\Gamma.\end{array}\right., \tag{1}\]

where \(\Omega\) is a domain in \(\mathbb{R}^{d}\) with boundary \(\Gamma\), \(\mathcal{D}\) is a domain differential operator, \(f:\Omega\to\mathbb{R}\) is a forcing term, \(\mathcal{B}\) is a trace/flux boundary operator and \(g:\Gamma\to\mathbb{R}\) the corresponding prescribed value. We approximate the solution \(u\) as the realisation of a NN of a given architecture. We introduce some notations and recall the key principles of the PINN framework in the rest of this section.

### Neural networks

Our solution ansatz is a fully-connected, feed-forward NN, which is obtained as the composition of affine maps and nonlinear activation functions applied element-wise. In this work, we separate the structure of a NN from its realisation when the activation function and the values of the weights and biases are given. We describe the architecture of a network by a tuple \((n_{0},\ldots n_{L})\in\mathbb{N}^{(L+1)}\), where \(L\) is the number of affine maps it is composed of, and for \(1\leq k\leq L\) the number of neurons on layer \(k\) is \(n_{k}\). We take \(n_{0}=d\) and since we only consider scalar-valued PDEs, we have \(n_{L}=1\).

For each layer number \(1\leq k\leq L\), we write \(\boldsymbol{\Theta}_{k}:\mathbb{R}^{n_{k-1}}\to\mathbb{R}^{n_{k}}\) the affine map at layer \(k\), defined by \(\boldsymbol{\Theta}_{k}\boldsymbol{x}=\boldsymbol{W}_{k}\boldsymbol{x}+ \boldsymbol{b}_{k}\) for some weight matrix \(\boldsymbol{W}_{k}\in\mathbb{R}^{n_{k}\times n_{k-1}}\) and bias vector \(\boldsymbol{b}_{k}\in\mathbb{R}^{n_{k}}\). We apply the activation function \(\rho:\mathbb{R}\to\mathbb{R}\) after every affine map except for the last one so as not to constrain the range of the model. The expression of the network is then

\[\mathcal{N}(\boldsymbol{\Theta},\rho)=\boldsymbol{\Theta}_{L}\circ\rho\circ \boldsymbol{\Theta}_{L-1}\circ\ldots\circ\rho\circ\boldsymbol{\Theta}_{1},\]

where \(\boldsymbol{\Theta}\) stands for the collection of trainable parameters of the network \(\boldsymbol{W}_{k}\) and \(\boldsymbol{b}_{k}\). Although the activation functions could be different at each layer or even trainable, we apply the same, fixed activation function everywhere.

### Abstract setting for the loss function

The loss function serves as a distance metric between the solution and the model. We write the loss in terms of the network itself to simplify notations, although it is to be understood as a function of its trainable parameters. Throughout this paper, we write \(\left\langle u,v\right\rangle_{\Omega}=\int_{\Omega}uv\,\mathrm{d}\Omega\) and \(\left\langle u,v\right\rangle_{\Gamma}=\int_{\Gamma}uv\,\mathrm{d}\Gamma\) the classical inner products on \(L^{2}(\Omega)\) and \(L^{2}(\Gamma)\). We write \(\left\|\cdot\right\|_{\Omega}\) and \(\left\|\cdot\right\|_{\Gamma}\) the corresponding \(L^{2}\) norms.

#### 2.2.1. Strong formulation

According to the original PINN formulation [1], the network is trained to minimise the strong residual of Eq. (1). The minimisation is performed in \(L^{2}\), so we suppose that \(\mathcal{D}u\) and \(f\) are in \(L^{2}(\Omega)\), and \(\mathcal{B}u\) and \(g\) are in \(L^{2}(\Gamma)\), This ensures that the residual

\[\mathcal{R}(u)\doteq\frac{1}{2}\left\|\mathcal{D}u-f\right\|_{\Omega}^{2}+ \frac{\beta}{2}\left\|\mathcal{B}u-g\right\|_{\Gamma}^{2},\]

is well-posed, where \(\beta>0\) is a weighting term for the boundary conditions.

#### 2.2.2. Weak formulation

The FEM derives a weak form of Eq. (1) that we can write in the following terms:

\[\left\{\begin{array}{ll}\text{Find }u\in U\text{ such that }\\ \forall v\in V,&a(u,v)=\ell(v).\end{array}\right.\]

In this formulation, \(a:U\times V\to\mathbb{R}\) is a bilinear form, \(\ell:V\to\mathbb{R}\) is a linear form and \(U\) and \(V\) are two functional spaces on \(\Omega\) such that \(a(u,v)\) and \(\ell(v)\) are defined for all \(u\in U\) and \(v\in V\). When the bilinear form \(a\) is symmetric, coercive and continuous in \(U\times V\), and \(U=V\), this weak setting can be recast as the minimisation of the energy functional

\[\mathcal{J}(u)\doteq\frac{1}{2}a(u,u)-\ell(u)\]
for \(u\in U\). In our case, \(U\) is the manifold of NNs with a given architecture defined on \(\Omega\). One difference with the FEM setting is that it is difficult to strongly enforce essential boundary conditions on NN spaces. Among other approaches, it is common to multiply the network by a function that vanishes on the boundary and add a lifting of the boundary condition [15]. Instead, we follow the penalty method and modify the residual: the boundary terms that come from integration by parts do not cancel and we suppose that it is possible to alter the bilinear form \(a\) into \(a_{\beta}\) and linear form \(\ell\) into \(\ell_{\beta}\) such that \(a_{\beta}\) is still symmetric and continuous. To ensure the coercivity of the bilinear form and weakly enforce Dirichlet boundary conditions, we add the term \(\beta\left\langle u,v\right\rangle_{\Omega}\) in the bilinear form and \(\beta\left\langle g,v\right\rangle_{\Gamma}\) in the linear form. Provided that the penalty coefficient \(\beta\) is large enough, \(a_{\beta}\) can be made coercive [16, p. 196] and we consider the energy minimisation problem

\[\text{Find }u\in U\text{ such that for all }v\in U,\mathcal{J}_{\beta}(u) \doteq\frac{1}{2}a_{\beta}(u,u)-\ell_{\beta}(u)\leq\mathcal{J}_{\beta}(v). \tag{2}\]

We see that the strong problem is a special case of the weak formulation. Indeed, by expanding the inner products and rearranging the terms, we obtain

\[\mathcal{R}(u)=\frac{1}{2}\left(\left\langle\mathcal{D}u,\mathcal{D}u \right\rangle_{\Omega}+\beta\left\langle\mathcal{B}u,\mathcal{B}u\right\rangle _{\Gamma}\right)-\left(\left\langle\mathcal{D}u,f\right\rangle_{\Omega}+\beta \left\langle\mathcal{B}u,g\right\rangle_{\Gamma}\right)+C,\]

where \(C\) is a constant with respect to \(u\) that we can disregard in the context of optimisation. If \(\mathcal{B}\) is a Dirichlet operator, this expression corresponds to the penalty method. After considering the functional derivative of \(\mathcal{R}\) and enforcing it to be zero, we recognise the notations of the weak form of problem Eq. (1), where the bilinear and linear form are defined as \(a_{\beta}(u,v)\doteq\left\langle\mathcal{D}u,\mathcal{D}v\right\rangle_{\Omega }+\beta\left\langle\mathcal{B}u,\mathcal{B}v\right\rangle_{\Gamma}\) and \(\ell_{\beta}(v)\doteq\left\langle\mathcal{D}v,f\right\rangle_{\Omega}+\beta \left\langle\mathcal{B}v,g\right\rangle_{\Gamma}\). In the remainder of this work, we keep the notations of the weak form as in Eq. (2) and simply write \(\mathcal{J}(u)\) the energy to be minimised.

### Evaluation of the loss function and optimisation

The bilinear and linear forms are sums of integrals over the domain \(\Omega\) and its boundary \(\Gamma\). In the general case, these integrals have to be approximated because no analytical expressions can be obtained, or they are numerically intractable. This is why most PINN implementations and their variants approximate the loss function with MC integration. The integration points can be resampled at a given frequency during training to prevent overfitting or to place more points where the pointwise residual is larger [10]. In this work, we are precisely designing a new integration method to replace MC.

The trainable parameters of the network are randomly initialised and a gradient descent algorithm can be used to minimise \(\mathcal{J}\) with respect to \(\boldsymbol{\Theta}\). The partial derivatives of \(u\) with respect to the space variable \(\boldsymbol{x}\), and the gradient of the loss function with respect to the parameters \(\boldsymbol{\Theta}\) are made available through automatic differentiation.

## 3. Regularisation of neural networks

It was shown in [5] that NNs activated by the \(\mathrm{ReLU}\) activation function can emulate first-order Lagrange elements, thus advocating the use of \(\mathrm{ReLU}\) networks to approximate PDEs. In this section, we show that the convergence of CPWL networks is slow and noisy whenever the energy functional \(\mathcal{J}\) involves a spatial derivative of \(u\). We then provide a reason for preferring smoother activations and suggest a few alternatives for the \(\mathrm{ReLU}\) function.

### Limitations of CPWL activation functions

When the activation is CPWL, the NN \(u\) is continuous and differentiable, and its partial derivatives are piecewise constant and discontinuous. Even though derivatives of order two and higher will be evaluated as zero, from a theoretical point of view however, this means that only the first-order derivatives of the network make pointwise sense, while derivatives of higher order do not.

Let us consider an illustrative example to clarify the kind of problems that arise when the energy functional involves a spatial derivative of \(u\). We consider the Poisson problem \(u^{\prime\prime}=f^{\prime\prime}\) in \(\Omega\) with \(u=f\) on \(\Gamma\), which we reframe into the following variational problem using the Nitsche method

\[a(u,v) =\left\langle u^{\prime},v^{\prime}\right\rangle_{\Omega}-\left\langle u ^{\prime}\cdot n,v\right\rangle_{\Gamma}-\left\langle u,v^{\prime}\cdot n \right\rangle_{\Gamma}+\beta\left\langle u,v\right\rangle_{\Gamma}\] \[\ell(v) =\left\langle-f^{\prime\prime},v\right\rangle_{\Omega}-\left\langle f,v^{\prime}\cdot n\right\rangle_{\Gamma}+\beta\left\langle f,v\right\rangle_{ \Gamma}.\]

Here \(n\) stands for the outward-pointing unit normal to \(\Gamma\). We consider the manufactured solution \(f=\mathrm{ReLU}\) on \(\Omega=[-1,+1]\). In this case \(f^{\prime\prime}=\delta_{0}\) is the Dirac delta centered at zero, which makes the first integral in \(\ell(v)\) equal to \(-v(0)\). We consider a one-layer network \(u(x)=f(x+c)\) to approximate the solution to this problem, and we wish to recover \(c=0\). The energy functional has the following expression:

\[\mathcal{J}(c)=\frac{1}{2}\left\{\begin{array}{ll}0&\text{if }c<-1\\ \beta c^{2}+|c|-(\beta-1)&\text{if }-1\leq c\leq+1\\ 2\beta c^{2}-2(\beta-1)c&\text{if }c>+1\end{array}\right.\]
We verify that \(\mathcal{J}\) shows a discontinuity at \(-1\). Besides, \(\mathcal{J}\) is decreasing on \([-1,0]\) and increasing on \([0,+1]\). We compute \(\mathcal{J}(0)=-\frac{1}{2}(\beta-1)\), which is negative whenever \(\beta>1\). This proves that \(\mathcal{J}\) has a global minimum at \(c=0\). However, the gradient of \(\mathcal{J}\) is discontinuous at \(0\), being equal to \(-1/2\) on the left and \(+1/2\) on the right. Because \(\mathrm{ReLU}\) is not differentiable at zero, any gradient-based optimiser will suffer from oscillations, implying slow and noisy convergence.

This example only involved one unknown, but one would already have to choose a low learning rate and apply a decay on the learning rate to make the network converge to an acceptable solution. One can imagine that when the optimisation problem involves a whole network, these oscillations can interfere with and severely hinder convergence. Ensuring that the activation function is at least \(\mathcal{C}^{1}\) enables one to interchange the derivative and integral signs, and if it is \(\mathcal{C}^{2}\) then the gradient of the loss function is continuous.

### Regularisation of the activation function

The regularity of a NN is entirely dictated by that of its activation function. In order to make the energy functional smoother while keeping the motivation of CPWL, our idea is to approximate the CPWL function by a smoother equivalent and replace every occurrence of the CPWL activation by its smoother counterpart in the network.

Our approach is similar to the "canonical smoothings" of \(\mathrm{ReLU}\) networks introduced in [17]. We take \(\mathrm{ReLU}\) as an example and point out several families of smooth equivalents. In each case, we introduce the normalised variable \(\bar{x}=x/\gamma\), where \(\gamma\) is a constant that only depends on the choice of the family.

* If we rewrite \(\mathrm{ReLU}(x)=\frac{1}{2}(x+|x|)\), we see that the lack of regularity comes from the absolute value being non-differentiable at zero. We can thus find a smooth equivalent of the absolute value and replace it in this expression. We set \(\gamma=2\varepsilon\) and define \(\rho_{\varepsilon}\) as \[\rho_{\varepsilon}(x)=\frac{1}{2}\left(x+\gamma\sqrt{1+\bar{x}^{2}}\right).\]
* The first derivative of \(\mathrm{ReLU}\) is discontinuous at zero. We can replace it with any continuous sigmoid function and obtain a smoothing of \(\mathrm{ReLU}\) by integrating it. We set \(\gamma=\frac{2}{\ln 2}\varepsilon\) and define \(\rho_{\varepsilon}\) as \[\rho_{\varepsilon}(x)=\frac{1}{2}\left(x+\gamma\ln(2\cosh(\bar{x}))\right).\]
* The second derivative of \(\mathrm{ReLU}\) is the Dirac delta at zero. We can draw from the well-established theory of mollifiers to obtain a smooth Dirac and integrate it twice to build a regularised version of \(\mathrm{ReLU}\). We set \(\gamma=2\sqrt{\pi}\varepsilon\) and define \(\rho_{\varepsilon}\) as \[\rho_{\varepsilon}(x)=\frac{1}{2}\left(x+x\,\mathrm{erf}(\bar{x})+\frac{ \gamma}{\sqrt{\pi}}\exp(-\bar{x}^{2})\right).\]

The following properties hold for the three classes of smooth equivalents of \(\mathrm{ReLU}\) introduced above and all values of \(\varepsilon>0\). The function \(\rho_{\varepsilon}\) is \(\mathcal{C}^{\infty}(\mathbb{R})\), convex and monotonically increasing in \(\mathbb{R}\). The graph of \(\rho_{\varepsilon}\) lies above that of \(\mathrm{ReLU}\). Similarly to \(\mathrm{ReLU}\), \(\rho_{\varepsilon}\) satisfies the property \(\rho_{\varepsilon}(x)-\rho_{\varepsilon}(-x)=x\). Finally, the parameter \(\varepsilon\) controls the pointwise distance between \(\rho_{\varepsilon}\) and \(\mathrm{ReLU}\), in the sense that \(\|\rho_{\varepsilon}-\rho\|_{L^{\infty}(\mathbb{R})}=\rho_{\varepsilon}(0)=\varepsilon\).

In particular, Lemma 2.5 in [17] applies to all the candidates above, and we refer to Proposition 2.6 in the same article for a bound of \(\|\mathcal{N}(\boldsymbol{\Theta},\mathrm{ReLU})-\mathcal{N}(\boldsymbol{ \Theta},\rho_{\varepsilon})\|_{L^{\infty}(\Omega)}\) in terms of \(\|\mathrm{ReLU}-\rho_{\varepsilon}\|_{L^{\infty}(\mathbb{R})}\).

The computational cost of these functions must also be taken into account and for this reason we choose to use the family obtained by the first method. We write it \(\mathrm{ReLU}_{\varepsilon}\) in the rest of this work.

Figure 1. Examples of regularising families for the \(\mathrm{ReLU}\) function and corresponding first and second derivatives.
## 4. Global approximation of an activation function by a CPWL function

In this section, we suppose that we are given a smooth activation \(\rho\) and we consider the NN \(u=\mathcal{N}(\boldsymbol{\Theta},\rho)\). We draw inspiration from the FEM and we would like to construct a mesh of the domain such that the behaviour of the NN is close to linear on each cell of the mesh. We design a method that approximates the smooth activation function \(\rho\) by a CPWL function \(\pi[\rho]\) and bound the distance between the two functions.

### Assumptions on the activation function

Since we have no prior bound on the parameters of the network, we need the approximation of \(\rho\) by \(\pi[\rho]\) to hold on \(\mathbb{R}\) as a whole. We impose that \(\pi[\rho]-\rho\) is square integrable on \(\mathbb{R}\), which implies that \(\rho\) has a linear behaviour at infinity. If \(\rho\) is twice differentiable on \(\mathbb{R}\), we know that \(\pi[\rho]-\rho\) behaves as \(x^{2}f^{\prime\prime}(x)\) at infinity. As a consequence, one way to enforce the integrability condition is to make sure that \(x^{2}f^{\prime\prime}(x)\) is of the order of \(|x|^{-\alpha}\) for a given \(\alpha>1/2\), which is equivalent to \(|x|^{\beta}f^{\prime\prime}(x)\) being bounded on \(\mathbb{R}\) for a given \(\beta>5/2\).

As will be made clear in the proof of Lem. 1, we also need \(f\) to be analytical in the neighbourhood of the zeros of \(f^{\prime\prime}\). Altogether, we restrict our approach to functions that belong to the functional space \(\mathcal{A}\) defined as

\[\mathcal{A}=\left\{\begin{array}{l}f\in\mathcal{C}^{2}(\mathbb{R}),|x|^{ \alpha}\,f^{\prime\prime}(x)\text{ is bounded for a given }\alpha>5/2,\\ f^{\prime\prime}(x)=0\implies f\text{ is analytical in a neighbourhood of }x\end{array}\right\}.\]

We point out that for all \(f\in\mathcal{A}\), \(\pi[f]-f\in L^{p}(\mathbb{R})\) for all \(p\geq 2\) and \(f^{\prime\prime}\in L^{p}(\mathbb{R})\) for all \(p\geq 2/5\).

### Approximation space

In our experiments, we observed that the way we approximate the activation function by a CPWL counterpart plays an important role in the performance of our adaptive quadrature. There are existing algorithms and heuristics to build \(\pi[\rho]\)[18] but we found that \(\pi[\rho]\) needs to satisfy extra properties for our method to be most efficient in practice. In particular, we noticed that it is preferable that the CPWL approximation of \(\rho\) lie below \(\rho\) in the regions where \(\rho\) is convex, and above \(\rho\) where \(\rho\) is concave. To that aim, we wish to decompose \(\mathbb{R}\) into intervals where \(\rho\) is either strictly convex, linear or strictly concave. We can then approximate \(\rho\) in each interval by connecting the tangents to \(\rho\) at free points. More formally, we define the partition induced by \(f\) as follows.

**Definition 1** (Partition induced by a function).: _Let \(f\in\mathcal{A}\). We define the set \(Z_{0}(f)=\{I\subset\mathbb{R},f^{\prime\prime}_{|I}=0\}\) that contains the intervals (possibly isolated points) where \(f^{\prime\prime}\) is zero and \(Z(f)=\partial Z_{0}(f)\cup\{-\infty,+\infty\}\) the set of points where \(f^{\prime\prime}\) is zero, excluding intervals but including \(\pm\infty\). We write \(z(f)=|Z(f)|\) and \(\xi_{1}<\ldots<\xi_{z(f)}\) the ordered elements of \(Z(f)\). The partition induced by \(f\) is defined as \(P(f)=\{[\xi_{i},\xi_{i+1}],1\leq i\leq z(f)-1\}\), that is the collection of intervals whose ends are two consecutive elements of \(Z(f)\). The elements of \(P(f)\) are called compatible intervals for \(f\). We finally introduce \(1\leq\kappa(f)\leq z(f)-1\) the number of segments in \(P(f)\) where \(f\) is not linear._

By way of example, the second derivative of \(\mathrm{ReLU}_{\varepsilon}\) is always positive so \(Z(\mathrm{ReLU}_{\varepsilon})=\varnothing\) and \(P(\mathrm{ReLU}_{\varepsilon})=\{\mathbb{R}\}\). However, \(Z(\tanh)=\{0\}\) and \(P(\tanh)=\{\mathbb{R}_{-},\mathbb{R}_{+}\}\).

For \(x\in\mathbb{R}\), we write \(T[f,x]\) the tangent to \(f\) at \(x\) and we extend this notation to \(x=\pm\infty\) by considering the asymptotic tangents to \(f\). From the construction of \(P(f)\), it is easy to see that for all \(I\in P(f)\) and \(x<y\in I\), the two tangents \(T[f,x]\) and \(T[f,y]\) intersect at a unique \(z\in|x,y[\). We write \(T[f,x,y]\) the CPWL function defined on \([x,y]\) as \(T[f,x]\) on \([x,z]\) and \(T[f,y]\) on \([z,y]\).

We are now ready to define our approximation space. For \(n\geq z(f)\), we define the space \(\mathcal{A}_{n}(f)\) of CPWL functions that connect the tangents to \(f\) at \(Z(f)\) as well as at \(n-z(f)\) other free points that do not belong to \(Z(f)\). As an example, the function plotted in Fig. (a)a belongs to \(\mathcal{A}_{7}(\tanh)\). There are four free points because \(z(\tanh)=3\). The square points are fixed and they are located at \(\pm\infty\) and zero.

### Convergence rate

For \(f\in\mathcal{A}\), \(n\geq z(f)\) and \(p\geq 2\), we write \(\pi_{n,p}[f]\) the best approximation of \(f\) in \(\mathcal{A}_{n}(f)\) in the \(L^{p}\) norm. In this paragraph, we bound the distance between \(f\) and \(\pi_{n,p}[f]\), and determine the convergence rate of this approximation. We start by finding an upper bound for the distance between \(f\) and \(T[f,x,y]\) for \(x<y\in I\), where \(I\) is a compatible interval for \(f\).

**Lemma 1**.: _Let \(f\in\mathcal{A}\) and \(p\in[2,\infty[\). There exist constants \(A_{p}(f)\), \(A_{\infty}(f)\), \(B_{p}(f)\) and \(B_{\infty}(f)\) such that for all interval \(I\) compatible for \(f\) and \(J=[x,y]\subset I\), it holds_

\[\left\|f-T[f,x,y]\right\|_{L^{p}(J)}\leq A_{p}(f)\left\|f^{\prime \prime}\right\|_{L^{p/(2p+1)}(J)}, \left\|f-T[f,x,y]\right\|_{L^{\infty}(J)}\leq A_{\infty}(f)\left\|f^{\prime \prime}\right\|_{L^{1/2}(J)}, \tag{4}\] \[\left\|f^{\prime}-T[f,x,y]^{\prime}\right\|_{L^{p}(J)}\leq B_{p}(f) \left\|f^{\prime\prime}\right\|_{L^{p/(p+1)}(J)}, \left\|f^{\prime}-T[f,x,y]^{\prime}\right\|_{L^{\infty}(J)}\leq B_{\infty}(f) \left\|f^{\prime\prime}\right\|_{L^{1}(J)}. \tag{3}\]
Thanks to this lemma, we can bound the distance between \(f\) and \(\pi_{n,p}[f]\) and prove a quadratic convergence for all \(p\).

**Proposition 1** (Best approximation on \(\mathcal{A}_{n}\)).: _Let \(f\in\mathcal{A}\), \(n\geq z(f)\) and \(p\in[2,\infty[\). Using the constants from Lem. 1, it holds_

\[\left\|f-\pi_{n,p}[f]\right\|_{L^{p}(\mathbb{R})}\leq\frac{A_{p} (f)\kappa(f)^{(2p+1)/p}\left\|f^{\prime\prime}\right\|_{L^{p/(2\pi+1)}(\mathbb{ R})}}{n^{2}}, \left\|f-\pi_{n,p}[f]\right\|_{L^{\infty}(\mathbb{R})}\leq\frac{A_{\infty}(f) \kappa(f)^{2}\left\|f^{\prime\prime}\right\|_{L^{1/2}(\mathbb{R})}}{n^{2}},\] \[\left\|f^{\prime}-\pi_{n,p}[f]^{\prime}\right\|_{L^{p}(\mathbb{R })}\leq\frac{B_{p}(f)\kappa(f)^{(p+1)/p}\left\|f^{\prime\prime}\right\|_{L^{p/ (p+1)}(\mathbb{R})}}{n}, \left\|f^{\prime}-\pi_{n,p}[f]^{\prime}\right\|_{L^{\infty}( \mathbb{R})}\leq\frac{B_{\infty}(f)\kappa(f)\left\|f^{\prime\prime}\right\|_{ L^{1}(\mathbb{R})}}{n}.\]

### Numerical implementation

In practice, we choose \(p=2\) and implement the least-square approximation of \(\rho\) on \(\mathcal{A}_{n}(\rho)\). We write \(\pi_{n}\) instead of \(\pi_{n,2}\). Since we are searching for the best approximant within \(\mathcal{A}_{n}[\rho]\), the unknowns to be determined are the locations \((\xi_{i})_{1\leq i\leq n}\) of the tangents. We write the total \(L^{2}\) norm between \(\rho\) and \(\rho_{n}\) in terms of \(\xi_{i}\) and seek to minimise this quantity via a gradient descent algorithm.

We are interested in two activation functions: \(\mathrm{ReLU}_{\varepsilon}\) and \(\tanh\). Both functions are symmetric around the origin, which allows us to restrict the optimisation domain. We also verify that both functions belong to \(\mathcal{A}\).

* The relationship \(\mathrm{ReLU}_{\varepsilon}(x)-\mathrm{ReLU}_{\varepsilon}(-x)=x\) holds for all \(x\in\mathbb{R}\). We impose the same condition on \(\pi_{n}[\mathrm{ReLU}_{\varepsilon}]\) and verify that \(\mathrm{ReLU}_{\varepsilon}-\pi_{n}[\mathrm{ReLU}_{\varepsilon}]\) is even. Since \(\mathrm{ReLU}_{\varepsilon}\) takes bounded values on \(\mathbb{R}_{-}\), we approximate \(\mathrm{ReLU}_{\varepsilon}\) by \(\pi_{n}[\mathrm{ReLU}_{\varepsilon}]\) on \(\mathbb{R}_{-}\) only and recover \(\pi_{n}[\mathrm{ReLU}_{\varepsilon}]\) as a whole by symmetry. The second derivative of \(\mathrm{ReLU}_{\varepsilon}\) is \(x\mapsto\frac{1}{2\gamma}(1+\overline{x}^{2})^{-3/2}\). It follows a polynomial decal rate with \(\alpha=3>5/2\). We verify that the second derivative of \(\mathrm{ReLU}_{\varepsilon}\) is never zero. The asymptotic tangents to \(\mathrm{ReLU}_{\varepsilon}\) intersect, so we can approximate \(\mathrm{ReLU}_{\varepsilon}\) with \(2\) pieces and this approximation is \(\mathrm{ReLU}\).
* The \(\tanh\) function is odd, so imposing that \(\pi_{n}[\tanh]\) is odd enforces that \(\tanh-\pi_{n}[\tanh]\) is even. We decide to perform the minimisation on \(\mathbb{R}_{+}\) and obtain an approximation of \(\tanh\) on \(\mathbb{R}\) as a whole by symmetry. The second derivative of \(\tanh\) is \(x\mapsto-2\tanh(x)(1-\tanh(x)^{2})\). This is asymptotically equivalent to \(-8\exp(-2x)\). In particular \(x\mapsto x^{3}\tanh^{\prime\prime}(x)\) is bounded on \(\mathbb{R}\). Moreover, \(\tanh^{\prime\prime}(x)=0\) if and only if \(x=0\), and \(\tanh\) is analytic in a neighbourhood of zero (its Taylor series converges within a radius of \(\pi/2\)). Since the asymptotic tangents to \(\tanh\) do not intersect, the minimum number of pieces to approximate \(\tanh\) on \(\mathbb{R}\) is three and the corresponding CPWL function is known as the hard \(\tanh\) activation.

As an example, we plot \(\pi_{7}[\tanh]\) in Fig. 1(a). The convergence plot in \(L^{2}\) norm of our method applied to both functions is shown in Fig. 1(b), together with the theoretical convergence rate of \(-2\). We verify that the convergence is asymptotically quadratic. We notice that the approximation of \(\mathrm{ReLU}_{\varepsilon}\) shows a preasymptotic behaviour in which the effective convergence rate is closer to \(1.5\).

## 5. Adaptive meshes and quadratures for neural networks

Let \(u=\mathcal{N}(\mathbf{\Theta},\rho)\) and \(\pi[u]=\mathcal{N}(\mathbf{\Theta},\pi[\rho])\) be two NNs sharing the same weights, activated by a smooth function and a CPWL approximation of it, obtained with the method explained above. In practice \(\pi[\rho]\) will be

Figure 2. Visualisation of \(\pi_{7}(\tanh)\) and convergence plot of our proposed method in \(L^{2}\) norm for \(\mathrm{ReLU}_{\varepsilon}\) and \(\tanh\). On (a), the square pink points are fixed, while the circle green points are chosen so as to minimise the \(L^{2}\) norm of \(\tanh-\pi_{7}[\tanh]\).
fixed, so we drop the subscript \(n\) in this section. We are interested in recovering a maximal decomposition of \(\Omega\) (resp. \(\Gamma\)) such that \(\pi[u]\) is linear in these regions. We say that this decomposition is a mesh of \(\Omega\) (resp. \(\Gamma\)) adapted to \(\pi[u]\). Since \(\pi[\rho]\) is fixed, this mesh depends only on \(\boldsymbol{\Theta}\). We introduce the notation \(\tau_{\Omega}(\boldsymbol{\Theta})\) and \(\tau_{\Gamma}(\boldsymbol{\Theta})\) to refer to these meshes, and \(\tau(\boldsymbol{\Theta})\) to refer to either of the two meshes. We equip each cell with a Gaussian quadrature to evaluate the loss function. Bounds for the integration error of our Adaptive Quadrature (AQ) are provided at the end of this section.

### Convexity of the mesh

We observe that each neuron corresponds to a scalar-valued linear map followed by the composition by \(\pi[\rho]\). The breakpoints of \(\pi[\rho]\) define hyperplanes in the input space of each neuron. Indeed, the composition of a linear map \(\boldsymbol{x}\to\boldsymbol{W}\boldsymbol{x}+\boldsymbol{b}\) by \(\pi[u]\) is CPWL and the boundaries of the regions where the composition is linear correspond to the equations \(\boldsymbol{W}\boldsymbol{x}+\boldsymbol{b}=\xi\), where \(\xi\) are the breakpoints of \(\pi[\rho]\).

Intuitively, we can obtain the adapted meshes \(\tau(\boldsymbol{\Theta})\) by considering the hyperplanes attached to every neuron of the network. The cells of the mesh are the largest sets of points that lie within exactly one piece of \(\pi[\rho]\) at each neuron. This is why these regions are also called the "activation patterns" of a NN, as we can label them by the pieces they belong to at each neuron [19].

The implementation of an algorithm that, given the parameters \(\boldsymbol{\Theta}\) of a NN, the linearisation \(\pi[\rho]\) and the input space \(\Omega\), outputs \(\tau_{\Omega}(\boldsymbol{\Theta})\) has to be robust to a large variety of corner cases, as in practice the cells can be of extremely small measures, or take skewed shapes. Besides, we are concerned with the complexity of this algorithm as it is meant to be run at every iteration (or at a given frequency) during the training phase. Fortunately enough, the only cells in the mesh that may not be convex must be in contact with the boundary.

**Lemma 2** (Convexity of the mesh).: _Let \(u:\mathbb{R}^{d}\to\mathbb{R}\) be a NN with weights \(\boldsymbol{\Theta}\) activated by a CPWL function and \(\Omega\subset\mathbb{R}^{d}\). If a cell in \(\tau_{\Omega}(\boldsymbol{\Theta})\) is not convex, it has to intersect with the boundary of \(\Omega\). In particular, if \(\Omega\) is convex, all the cells of \(\tau_{\Omega}(\boldsymbol{\Theta})\) are convex._

If \(\Omega\) is not convex, we decompose \(\Omega\) into a set of convex polytopes and construct an adapted mesh for each polytope independently. By Lem. 2 these submeshes are convex. The mesh composed of all the cells of the submeshes is adapted to \(u\) on \(\Omega\) and contains convex cells only. Thus without loss of generality, we suppose that \(\Omega\) is convex. Clipping a line by a polygon is made easier when the polygon is convex. Our method to build \(\tau_{\Omega}(\boldsymbol{\Theta})\) takes great advantage of the convexity of the mesh. Our algorithm is described in Alg. 1.

### Representation of a linear region

Let \(u_{k}\) denote the composition of the first \(k\) layers of the network \(u\), and \(\tau^{k}(\boldsymbol{\Theta})\) be the mesh adapted to \(u_{k}\). We represent a linear region of \(\tau^{k}(\boldsymbol{\Theta})\) and its corresponding local expression by \((R,\boldsymbol{W}_{R},\boldsymbol{b}_{R})\), where \(R\in\tau^{k}(\boldsymbol{\Theta})\) is a cell of the mesh, \(\boldsymbol{W}_{R}\in\mathbb{R}^{n_{k}\times d}\) and \(\boldsymbol{b}_{R}\in\mathbb{R}^{n_{k}}\) are the coefficients of the restriction of \(\pi[u_{k}]\) to \(R\): \(\pi[u_{k}]_{|R}(\boldsymbol{x})=\boldsymbol{W}_{R}\boldsymbol{x}+\boldsymbol {b}_{R}\).

### Initialisation of the algorithm

We initialise the algorithm with the region \((\Omega,\boldsymbol{I}_{d},\boldsymbol{0}_{d})\), where \(\boldsymbol{I}_{d}\) is the identity matrix of size \(d\) and \(\boldsymbol{0}_{d}\) is the zero vector of size \(d\). The mesh \(\{(\Omega,\boldsymbol{I}_{d},\boldsymbol{0}_{d})\}\) is adapted to \(u_{0}\).

### Composition by a linear map

Suppose that \((R,\boldsymbol{W}_{R},\boldsymbol{b}_{R})\) is adapted to \(u_{k}\). The composition of \(\pi[u_{k}]_{|R}\) by the linear map \(\boldsymbol{\Theta}_{k+1}\) remains linear on \(R\), only the coefficients of the map are changed. They become \(\boldsymbol{W}_{k+1}\boldsymbol{W}_{R}\) and \(\boldsymbol{W}_{k+1}\boldsymbol{b}_{R}+\boldsymbol{b}_{k+1}\).

### Composition by an activation function

Suppose that \((R,\boldsymbol{W}_{R},\boldsymbol{b}_{R})\) is adapted to \(u_{k}\). We compose \(\pi[u_{k}]_{|R}\) by \(\pi[\rho]\) componentwise. Let \(\boldsymbol{w}_{i}\) be the \(i\)-th row of \(\boldsymbol{W}_{R}\) and \(b_{i}\) be the \(i\)-th coordinate of \(\boldsymbol{b}_{R}\). We need to find the pre-images of the breakpoints of \(\pi[\rho]\), which amounts to consider the hyperplanes \(\boldsymbol{w}_{i}\boldsymbol{x}+b_{i}=\xi_{j}\), where \(\xi_{j}\) are the breakpoints of \(\pi[\rho]\). We explain our method in detail in Apdx. B.1.

The process of mesh extraction in dimension two is illustrated in Fig. 3. In Fig. 3, the dashed line represents the orientation of the activation function. The plain lines depict two hyperplanes corresponding to two different breakpoints of \(\pi[\rho]\). Fig. 3 pictures the hyperplanes resulting from several neurons, thus oriented in different directions. In this example, there would be at least four neurons on the layer at hand. We notice that since the ranges of the neurons are different, all neurons may not activate all the pieces of \(\pi[\rho]\), and thus they can give rise to different numbers of hyperplanes. The clipping operation is shown in Fig. 3, and the intersections of the hyperplanes against themselves are displayed in Fig. 3. In this example, the initial cell would be cut into \(13\) subcells.
### Gaussian quadratures for convex polygons

We decompose the integrals in the linear and bilinear forms on the cells of \(\tau_{\Omega}(\mathbf{\Theta})\) and \(\tau_{\Gamma}(\mathbf{\Theta})\). In these cells, the terms that only depend on the network \(\pi[u]\) and its spatial derivatives are polynomials. As a consequence, the linear and bilinear forms involving \(\pi[u]\) can be computed exactly using Gaussian quadratures on segments in dimension one, and polygons in dimension two.

In dimension one, the cells of the mesh are segments. Gaussian quadrature rules are known and made available through plenty of libraries.

However, in dimension two, the cells of the mesh can be arbitrary convex polygons. To handle the general case, one approach could consist in splitting each convex cell into triangles and then applying a known Gaussian quadrature for triangles. This approach is the least economical in terms of the number of quadrature points. At the opposite end of the spectrum, we could rely on Gaussian quadratures of reference \(n\)-gons. Still, the order of the mapping from the reference \(n\)-gon to the \(n\)-gon at hand is \(n-2\), which makes the Jacobian of this mapping costly to evaluate. In this work, we strike a tradeoff and decompose each convex cell into a collection of triangles and convex quadrangles. We use a recursive algorithm to obtain this decomposition, as explained in Apdx. B.2. We refer to [20] for numerically accurate symmetric quadrature rules for triangles and quadrangles.

### Alternatives

There exist computationally cheaper numerical integration methods based on evaluations of the integrands at the vertices of the mesh. Indeed, using Stokes theorem one can transform surface integrals on polygons into line integrals on their edges, and in turn, into pointwise evaluations at their vertices [21]. This approach requires knowing the local expression of \(\pi[u]\), that is the coefficients of the linear interpolation of \(u\) on each cell.

We have conducted preliminary experiments using this approach but we have observed that in addition to being numerically unstable, the overall cost including the interpolation step is not lower. Indeed, finding the best-fitting plane that passes through given points involves the inversion of a \(3\times 3\) system on each cell. The coefficients of these matrices are the integrals over each cell of the polynomials \(x^{2}\), \(xy\), \(y^{2}\), \(x\), \(y\) and \(1\). In many cases, the cells take skewed shapes so these matrices can be extremely ill-conditioned.

### Analysis of the integration error

Our proposed adaptive quadrature consists in approximating \(\mathcal{J}(u)\) by a quadrature on the cells of the mesh adapted to \(\pi[u]\). Conceptually, this is equivalent to choosing a suitable piecewise polynomial function defined on the mesh adapted to \(\pi[u]\), and approximating \(\mathcal{J}(u)\) by \(\mathcal{J}_{\text{num}}(v)\). Here, \(\mathcal{J}_{\text{num}}\) is a discrete version of \(\mathcal{J}\) in which the integrals have been replaced by numerical quadratures. We break down the numerical integration error into two parts as follows

\[\left|\mathcal{J}(u)-\mathcal{J}_{\text{num}}(v)\right|\leq\left|\mathcal{J}( u)-\mathcal{J}(v)\right|+\left|\mathcal{J}(v)-\mathcal{J}_{\text{num}}(v) \right|.\]

The second term is the numerical quadrature error and is commonplace in the field of the FEM. In particular, one can consider piecewise polynomial approximations of physical parameters, forcing term, and boundary conditions up to a given order and use a Gaussian quadrature that cancels the numerical quadrature error of the energy functional. The first term is the error incurred by the piecewise polynomial approximation of \(u\) by \(v\). Since we know that \(\mathcal{J}\) is a variational energy functional that comes from bounded linear and bilinear forms, we can bound the linearisation error in the following way

\[\left|\mathcal{J}(u)-\mathcal{J}(v)\right|\leq\frac{1}{2}\left|a(u-v,u+v) \right|+\left|\ell(u-v)\right|\leq\frac{1}{2}C_{a}\left\|u-v\right\|_{a}\left\| u+v\right\|_{a}+C_{\ell}\left\|u-v\right\|_{\ell},\]

Figure 3. Example of a mesh extraction. (a) Parallel hyperplanes associated with different breakpoints. (b) Hyperplanes linked to a whole layer. (c) Clipping hyperplanes with region boundary. (d) Intersection of hyperplanes.
where \(C_{a}\) and \(C_{\ell}\) are the constants that bound \(a\) and \(\ell\), and \(\left\lVert\cdot\right\rVert_{a}\) and \(\left\lVert\cdot\right\rVert_{\ell}\) are integral norms. For instance, when we consider a strong formulation of a weak Poisson problem, \(\left\lVert\cdot\right\rVert_{a}\) is an \(H^{2}\) norm and \(\left\lVert\cdot\right\rVert_{\ell}\) is an \(L^{2}\) norm. In the case of a weak formulation, \(\left\lVert\cdot\right\rVert_{a}\) and \(\left\lVert\cdot\right\rVert_{\ell}\) are \(H^{1}\) norms.

#### 5.8.1. Weak formulations

In Prop. 1, we showed that the \(L^{\infty}\) norm of \(\rho-\pi_{n}[\rho]\) decays as \(n^{-2}\) and that of \(\rho^{\prime}-\pi_{n}[\rho]^{\prime}\) decays as \(n^{-1}\). We consider the two networks \(u=\mathcal{N}(\boldsymbol{\Theta},\rho)\) and \(\pi_{n}[u]=\mathcal{N}(\boldsymbol{\Theta},\pi_{n}[\rho])\), that only differ by their activation function. We show that we can bound \(u-\pi_{n}[u]\) in terms of \(\rho-\pi_{n}[\rho]\) in the following proposition. We define the \(L^{\infty}\) norm of a vector-valued function as the maximum \(L^{\infty}\) norm of each of its components.

**Proposition 2** (Bounds for the difference of neural networks).: _Let \(\Omega\subset\mathbb{R}^{d}\) be a bounded domain, \(f\), \(g\) two continuous and almost everywhere differentiable functions such that \(f-g\), \(f^{\prime}\) and \(g^{\prime}\) are bounded on \(\mathbb{R}\). We further assume that \(f\) and \(f^{\prime}\) are Lipschitz continuous. Let \(u=\mathcal{N}(\boldsymbol{\Theta},f)\) and \(v=\mathcal{N}(\boldsymbol{\Theta},g)\) two NNs that only differ by their activation functions. Then there exists three constants \(C_{1}\), \(C_{2}\), \(C_{3}\) such that_

\[\left\lVert u-v\right\rVert_{L^{\infty}(\Omega)} \leq C_{1}\left\lVert f-g\right\rVert_{L^{\infty}(\mathbb{R})}, \tag{6}\] \[\left\lVert\nabla u-\nabla v\right\rVert_{L^{\infty}(\Omega)} \leq C_{2}\left\lVert f-g\right\rVert_{L^{\infty}(\mathbb{R})}+C _{3}\left\lVert f^{\prime}-g^{\prime}\right\rVert_{L^{\infty}(\mathbb{R})}. \tag{5}\]

_The three constants depend on \(\boldsymbol{\Theta}\) and the Lipschitz constant of \(f\). Additionally \(C_{2}\) and \(C_{3}\) depend on \(\left\lVert f^{\prime}\right\rVert_{L^{\infty}(\mathbb{R})}\), \(\left\lVert g^{\prime}\right\rVert_{L^{\infty}(\mathbb{R})}\) and the Lipschitz constant of \(f^{\prime}\)._

The proposition above can be used to bound \(\left\lVert u-\pi[u]\right\rVert_{a}\) and \(\left\lVert u-\pi[u]\right\rVert_{\ell}\) whenever their integrands can be expressed as bivariate polynomials involving \(u\) and \(\nabla u\). Indeed, for all \(n\geq 1\), \(\rho\) and \(\pi_{n}[\rho]\) are continuous and almost everywhere differentiable and their derivatives are bounded. Furthermore Prop. 2 showed that \(\rho-\pi_{n}[\rho]\) and \(\rho^{\prime}-\pi_{n}[\rho]^{\prime}\) are bounded. Since \(\rho\in\mathcal{A}\), we infer that \(\rho^{\prime}\) and \(\rho^{\prime\prime}\) are bounded, which implies that \(\rho\) and \(\rho^{\prime}\) are Lipschitz continuous. Thus Prop. 2 can be applied and we have

\[\left\lVert u-\pi_{n}[u]\right\rVert_{L^{p}(\Omega)} \leq\frac{C\left\lvert\Omega\right\rvert^{1/p}}{n^{2}}, \left\lVert\nabla u-\nabla\pi_{n}[u]\right\rVert_{L^{p}(\Omega)} \leq\frac{D\left\lvert\Omega\right\rvert^{1/p}}{n},\]

where \(C=C_{1}A_{\infty}(\rho)\kappa(\rho)^{2}\left\lVert\rho^{\prime\prime}\right\rVert _{L^{1/2}(\mathbb{R})}\) and \(D=C_{2}A_{\infty}(\rho)\kappa(\rho)\left\lVert\rho^{\prime\prime}\right\rVert _{L^{1/2}(\mathbb{R})}+C_{3}B_{\infty}(\rho)\kappa(\rho)\left\lVert\rho^{ \prime\prime}\right\rVert_{L^{1}(\mathbb{R})}\). We used the fact that \(n\geq\kappa(\rho)\) to bound \(\kappa(\rho)/n\) in \(D\). We conclude that \(\left\lVert u-\pi_{n}[u]\right\rVert_{H^{1}(\Omega)}\) decays as \(n^{-1}\).

#### 5.8.2. Strong formulations

Since \(\pi_{n}[u]\) is only CPWL, the \(H^{2}\) norm of \(u-\pi_{n}[u]\) is not defined, which means that we cannot use the same approach as above to bound the integration error. Instead one could apply a similar method to that used in the FEM, where \(u\) is approximated by a piecewise polynomial function that can be integrated exactly. The integration error would involve the mesh size.

## 6. Numerical experiments

We carry out thorough numerical experiments to validate our integration method. Choosing the best hyperparameters to train a NN usually requires trying several combinations and picking the one that leads to the lowest generalisation error. In addition to the usual learning rate and the number of epochs, we also need to tune the design of the loss function by specifying the penalty coefficient \(\beta\), and the integration method. In the case of MC integration, we select the number of points in the domain and on the boundary, as well as the resampling frequency. Our adaptive quadrature introduces a new set of hyperparameters, namely the number of pieces in the CPWL approximation of the activation function, the order of the quadrature in each cell as well as the remeshing frequency.

In this work, we are also concerned with reducing the computational budget involved with the training of a NN to solve a PDE. For this reason, we consider a two-hidden layer feed-forward NN with \(10\) neurons on both layers. The number of trainable parameters is \(141\) and \(151\) in dimensions one and two respectively. The optimiser is chosen as ADAM and all the networks are trained for \(5000\) epochs. To compensate for the low number of training iterations, we use a slightly higher learning rate compared to common practice in the literature. We set it to \(10^{-2}\) or \(10^{-3}\), whereas it seems to usually lie between \(10^{-4}\) and \(10^{-3}\). We specify the learning rate for each experiment in the sections below.

In our experiments, we compare models that have been trained with the same average number of integration points. In the case of a network trained with MC, the number of integration points is fixed during training, but this is not the case for the method we propose here. In practice, we observe that the number of points starts at an intermediate value, then slightly decreases, and finally rises above the initial value. We still take the average
number of integration points across the whole training as a reference because it is the average number of times the model is evaluated at each epoch and is a good indicator of the cost of a numerical method.

We compare AQ against MC integration to solve the Poisson equation in dimensions one and two, both in the strong and weak form using the Nitsche formulation. We report the relative \(L^{2}\) norm of the pointwise difference, defined as

\[E(u,\hat{u})=\frac{\|\hat{u}-u\|_{L^{2}(\Omega)}}{\|\hat{u}\|_{L^{2}(\Omega)}}.\]

The loss function corresponding to each problem is introduced in Subsect. 6.1 and the manufactured solutions we consider are presented in Subsect. 6.2. We give an outline of our experiments and short conclusions for each of them.

* We conduct a set of exploratory experiments in Subsect. 6.3 to identify the best learning rate and sampling frequency, as well as to evaluate and compare the robustness of AQ and MC to the number of integration points. We find that regardless of the integration method, the strong Poisson problem brings about higher levels of noise during the training phase compared to the weak Poisson formulation. In general, AQ can reach lower errors than MC and reduce this noise to a certain extent while using fewer integration points. This is especially true for the weak Poisson problem, where the convergence is noticeably faster and smoother with AQ.
* In Subsect. 6.4, we show that models trained with our proposed integration method are more robust to parameter initialisation compared to MC. Furthermore, AQ consistently leads to lower errors than MC even when it relies on fewer integration points.
* The next round of experiments in Subsect. 6.5 shows that it is possible to reduce the number of integration points by merging small regions with their neighbours.
* Finally in Subsect. 6.6 we solve a Poisson equation on a slightly more complex domain. We find that for a similar number of integration points, AQ reduces the error of MC by \(70\%\).

We provide a summary of our experiments in Subsect. 6.7 and discuss limitations and possible extensions of our method in Subsect. 6.8.

In all the following tables, \(N_{\Omega}\) and \(N_{\Gamma}\) stand for the number of integration points in the domain and on the boundary, respectively. The letters \(P\) and \(O\) denote the number of pieces in the approximation of the activation function by a CPWL function, and the order of the quadrature in each cell of the mesh adapted to the network.

### Study cases

In the following experiments, we consider two simple domains: the segment \(\Omega_{1}=[-1,+1]\) and the square \(\Omega_{2}=[-1,+1]^{2}\). We weakly enforce Dirichlet boundary conditions on the whole boundary of the domain. We introduce the bilinear and linear forms corresponding to each problem below.

**Strong formulation:**: Let \(f:\Omega\to\mathbb{R}\) be a function of class \(\mathcal{C}^{2}(\Omega)\). We want to solve the equation \(\Delta u=\Delta f\) with the Dirichlet boundary condition \(u_{|\Gamma}=f_{|\Gamma}\). The original PINN formulation corresponds to the strong form of the PDE and is associated with the loss

\[J(u)=\frac{1}{2}\|\Delta u-\Delta f\|_{\Omega}^{2}+\frac{\beta}{2}\|u-f\|_{ \Gamma}^{2},\]

where \(\beta>0\) is a weight factor for the boundary condition. We transform the squared \(L^{2}\) norms as inner products and obtain the following bilinear and linear forms

\[a(u,v) =\left\langle\Delta u,\Delta v\right\rangle_{\Omega}+\beta\left\langle u,v\right\rangle_{\Gamma},\] \[\ell(v) =\left\langle\Delta f,\Delta v\right\rangle_{\Omega}+\beta\left\langle f,v\right\rangle_{\Gamma}.\]
**Weak formulation:**: We also solve a weak form of the Poisson equation and use the Nitsche method to obtain a symmetric bilinear form. In this setting, we only require \(f\) to be of class \(\mathcal{C}^{1}(\Omega)\). The bilinear and linear forms are as follows

\[a(u,v) =\left\langle\nabla u,\nabla v\right\rangle_{\Omega}-\left\langle \nabla u\cdot n,v\right\rangle_{\Gamma}-\left\langle u,\nabla v\cdot n \right\rangle_{\Gamma}+\beta\left\langle u,v\right\rangle_{\Gamma},\] \[\ell(v) =\left\langle-\Delta f,v\right\rangle_{\Omega}-\left\langle f, \nabla v\cdot n\right\rangle_{\Gamma}+\beta\left\langle f,v\right\rangle_{ \Gamma}.\]

Here, \(n\) is the outward-pointing unit normal vector to \(\Gamma\) and \(\beta>0\) is a coefficient large enough so that the bilinear form is coercive.

In all our experiments, the relative \(L^{2}\) norm is evaluated with a quadrature of order \(10\) on a \(100\)-cell uniform mesh in 1D and on a \(100\times 100\) uniform grid in 2D.
### Activation functions and forcing terms

We solve the two problems we have just described with different forcing terms \(f\). Throughout our experiments, all the problems that involve a given forcing term are solved using the same activation function. We report two groups of forcing terms and their corresponding activation in Tab. 1. The notation \(\mathrm{sinc}\) stands for the cardinal sine function \(x\mapsto\sin(x)/x\), defined on \(\mathbb{R}\backslash\{0\}\) and continuously extended at zero by setting \(\mathrm{sinc}(0)=1\).

### Relationship between integration method, learning rate and sampling frequency

In this initial experiment, we focus on solving the four one-dimensional problems (\(\mathrm{ReLU}_{\varepsilon}\)/\(\tanh\), weak/strong) so that exploring a large space of hyperparameters remains computationally affordable. All the networks are trained from the same seed, which means that they are initialised from the same set of weights and biases, and the MC integration points are sampled in the same way across all experiments.

We study the connection between the learning rate, the frequency at which we resample the integration points and the integration method. We set the learning rate to \(10^{-2}\) or \(10^{-3}\) and resample the integration points every \(1\), \(10\), \(100\), \(500\), \(1000\) or \(5000\) epochs. Since the networks are trained for \(5000\) epochs, the last frequency corresponds to fixed integration points. We compare MC with \(50\) or \(100\) points against AQ with several choices of quadrature order and number of pieces in \(\pi[\rho]\). When the domain is one-dimensional, the boundary integrals reduce to evaluating the integrand at \(\Gamma_{1}=\{-1,+1\}\).

We also solve the weak Poisson problem in dimension two to further assess the effect of the integration method when the sampling frequency is fixed.

#### 6.3.1. Strong Poisson in 1D

The relative \(L^{2}\) norm after training the network for the two one-dimensional strong Poisson problems is shown in Tab. 2 and Tab. 3. We first observe that the sampling frequency and the number of integration points have a large influence on the final error. It is unclear whether increasing the number of integration points helps improve performance, as both trends appear for the two activation functions and integration methods.

Generally speaking, we note that AQ works best with high refresh rates (every \(100\) or fewer epochs), especially when the number of pieces and the order of the quadrature are low. However, it is almost always the case that AQ reaches similar or better accuracies than MC with fewer points. For instance, when the resampling rate is set to \(10\) epochs and the learning rate to \(10^{-2}\) for the \(\mathrm{ReLU}_{\varepsilon}\) activation function, a quadrature of order \(2\) with \(3\) pieces involves \(57\) points on average and leads to a final error of \(6.25\times 10^{-3}\) while all of the MC settings bring to higher errors. We notice that using a higher-order quadrature or more pieces in the CPWL approximation of the activation function reduces the error more consistently than increasing the number of MC points, even though it is not systematic in either case.

To illustrate the learning phase, we plot the learning curve and final pointwise error in two different settings in Fig. 4 and Fig. 5. We remark that the levels of noise in the two scenarios are very different. In the first case, the convergence with AQ is less perturbed by noise than that with MC, and the error seems to decay faster with AQ than with MC. In the second case, the training suffers from high levels of noise for all integration settings, and this is a representative example of most of the networks solving the strong Poisson problem.

To summarise, the strong Poisson problem brings about high levels of noise during the training phase for both integration methods. Still, AQ can reach lower errors than MC and reduce this noise to a certain extent while using fewer integration points.
\begin{table}
\begin{tabular}{c c c c c c c c c} \hline \hline \multirow{2}{*}{\(\eta\)} & \multirow{2}{*}{\(N_{\mathrm{\Pi}}\left(P,\,O\right)\)} & \multicolumn{6}{c}{\(\nu\)} \\ \cline{3-8}  & & 1 & 10 & 100 & 500 & 1000 & 5000 \\ \hline \multirow{8}{*}{\(10^{-2}\)} & \multirow{3}{*}{MC} & 50 & \(2.65\times 10^{-2}\) & \(3.50\times 10^{-2}\) & \(6.81\times 10^{-2}\) & \(4.01\times 10^{-1}\) & \(1.76\times 10^{-1}\) & \(5.79\times 10^{-2}\) \\  & & 100 & \(2.97\times 10^{-2}\) & \(7.99\times 10^{-2}\) & \(3.18\times 10^{-2}\) & \(1.48\times 10^{-2}\) & \(4.24\times 10^{-2}\) & \(2.51\times 10^{-1}\) \\ \cline{2-8}  & & 25 (2, 2) & \(1.11\times 10^{-1}\) & \(1.05\times 10^{-1}\) & \(3.00\times 10^{-2}\) & \(2.39\times 10^{-2}\) & \(4.47\times 10^{-2}\) & \(1.97\times 10^{-0}\) \\  & & 42 (2, 5) & \(5.13\times 10^{-2}\) & \(3.18\times 10^{-2}\) & \(9.43\times 10^{-2}\) & \(6.39\times 10^{-2}\) & \(1.52\times 10^{-2}\) & \(1.62\times 10^{-0}\) \\  & & 75 (2, 10) & \(1.84\times 10^{-3}\) & \(2.26\times 10^{-2}\) & \(1.22\times 10^{-2}\) & \(3.83\times 10^{-2}\) & \(7.25\times 10^{-3}\) & \(4.32\times 10^{-1}\) \\  & & 57 (3, 2) & \(2.55\times 10^{-2}\) & \(6.25\times 10^{-3}\) & \(1.70\times 10^{-1}\) & \(2.60\times 10^{-1}\) & \(2.46\times 10^{-1}\) & \(3.04\times 10^{-2}\) \\  & & 84 (3, 5) & \(1.08\times 10^{-2}\) & \(3.81\times 10^{-2}\) & \(1.32\times 10^{-2}\) & \(1.76\times 10^{-3}\) & \(3.76\times 10^{-2}\) & \(2.65\times 10^{-2}\) \\  & & 86 (5, 2) & \(7.56\times 10^{-2}\) & \(6.18\times 10^{-3}\) & \(6.66\times 10^{-3}\) & \(3.58\times 10^{-2}\) & \(1.04\times 10^{-2}\) & \(1.09\times 10^{-0}\) \\ \hline \multirow{8}{*}{\(10^{-3}\)} & \multirow{3}{*}{MC} & 50 & \(4.67\times 10^{-2}\) & \(1.15\times 10^{-1}\) & \(1.10\times 10^{-1}\) & \(1.06\times 10^{-1}\) & \(4.40\times 10^{-1}\) & \(5.77\times 10^{-1}\) \\  & & 100 & \(4.29\times 10^{-2}\) & \(1.17\times 10^{-1}\) & \(4.03\times 10^{-2}\) & \(3.18\times 10^{-2}\) & \(6.50\times 10^{-3}\) & \(1.44\times 10^{-1}\) \\ \cline{2-8}  & & 28 (2, 2) & \(1.91\times 10^{-2}\) & \(4.64\times 10^{-2}\) & \(6.49\times 10^{-2}\) & \(4.35\times 10^{-2}\) & \(7.09\times 10^{-2}\) & \(2.41\times 10^{-0}\) \\  & & 44 (2, 5) & \(1.80\times 10^{-2}\) & \(6.67\times 10^{-3}\) & \(1.88\times 10^{-2}\) & \(2.49\times 10^{-2}\) & \(2.81\times 10^{-2}\) & \(1.85\times 10^{-0}\) \\ AQ & 78 (2, 10) & \(1.01\times 10^{-2}\) & \(6.75\times 10^{-3}\) & \(2.23\times 10^{-2}\) & \(5.03\times 10^{-3}\) & \(5.74\times 10^{-3}\) & \(1.26\times 10^{-0}\) \\  & & 59 (3, 2) & \(3.87\times 10^{-2}\) & \(3.84\times 10^{-2}\) & \(1.32\times 10^{-2}\) & \(7.94\times 10^{-2}\) & \(1.71\times 10^{-1}\) & \(1.33\times 10^{-2}\) \\  & & 86 (3, 5) & \(1.40\times 10^{-2}\) & \(3.87\times 10^{-2}\) & \(1.61\times 10^{-2}\) & \(9.66\times 10^{-3}\) & \(6.76\times 10^{-3}\) & \(7.50\times 10^{-3}\) \\  & & 89 (5, 2) & \(6.32\times 10^{-3}\) & \(1.63\times 10^{-2}\) & \(3.50\times 10^{-2}\) & \(6.45\times 10^{-3}\) & \(3.03\times 10^{-2}\) & \(1.38\times 10^{-2}\) \\ \hline \hline \end{tabular}
\end{table}
Table 2. Comparison of the relative \(L^{2}\) norm for the strong Poisson problem with the \(\mathrm{ReLU}_{\varepsilon}\) activation, depending on learning rate, resampling frequency and integration method.

Figure 4. Learning curve (a) and pointwise error (b) for the strong Poisson problem with the \(\tanh\) activation, depending on learning rate, resampling frequency and integration method.

\begin{table}
\begin{tabular}{c c c c c c c c} \hline \hline \multirow{2}{*}{\(\eta\)} & \multirow{2}{*}{\(N_{\mathrm{\Pi}}\left(P,\,O\right)\)} & \multicolumn{6}{c}{\(\nu\)} \\ \cline{3-8}  & & 1 & 10 & 100 & 500 & 1000 & 5000 \\ \hline \multirow{8}{*}{\(10^{-2}\)} & \multirow{3}{*}{MC} & 50 & \(2.62\times 10^{-2}\) & \(2.81\times 10^{-2}\) & \(7.34\times 10^{-3}\) & \(5.55\times 10^{-2}\) & \(2.75\times 10^{-2}\) & \(4.61\times 10^{-1}\) \\  & & 100 & \(2.96\times 10^{-2}\) & \(1.17\times 10^{-2}\) & \(2.81\times 10^{-2}\) & \(1.72\times 10^{-3}\) & \(9.02\times 10^{-3}\) & \(4.71\times 10^{-2}\) \\ \cline{2-8}  & & 38 (3, 2) & \(2.86\times 10^{-2}\) & \(6.26\times 10^{-3}\) & \(5.95\times 10^{-2}\) & \(1.83\times 10^{-2}\) & \(4.11\times 10^{-2}\) & \(7.62\times 10^{-0}\) \\  & & 57 (3, 5) & \(3.13\times 10^{-3}\) & \(3.84\times 10^{-3}\) & \(2.65\times 10^{-3}\) & \(2.19\times 10^{-3}\) & \(1.37\times 10^{-3}\) & \(1.54\times 10^{-0}\) \\  & AQ & 118 (3, 10) & \(1.73\times 10^{-3}\) & \(1.54\times 10^{-4}\) & \(1.74\times 10^{-3}\) & \(1.41\times 10^{-2}\) & \(6.11\times 10^{-2}\) & \(2.92\times 10^{-0}\) \\  & & 76 (5, 2) & \(8.55\times 10^{-2}\) & \(2.06\times 10^{-3}\
sampling frequency when the network is trained with AQ. When relying on MC, the variance of the final error is much higher and the dependence between these the variance and the number of integration points is unclear.

When the sampling frequency is high enough (every \(100\) or fewer epochs), we observe that the number of integration points can very often be reduced while keeping the same levels of error. For instance, in the problem with the \(\tanh\) function with a learning rate of \(10^{-2}\) and a resampling frequency of \(10\), shifting from \(7\) pieces to \(5\) or \(3\) pieces while keeping order \(2\) does not deteriorate the performance of AQ. However, reducing the number of MC integration points increases the error in most cases.

To further compare the training phase with AQ and MC, we plot the learning curve and pointwise error for the \(\tanh\) activation in Fig. 6 and for the \(\operatorname{ReLU}_{\varepsilon}\) activation in Fig. 7. In both cases, we remark that the convergence of the network is significantly smoother with AQ than with MC, even with few integration points. Moreover, it is clear that the \(L^{2}\) norm decreases in the case of AQ, whereas it seems to plateau with MC. Finally, it is interesting to note that the learning curves follow the same trend at the beginning of the training. After a few hundred epochs, the learning curves corresponding to AQ keep on decreasing while the ones corresponding to MC break away from this trend and start to suffer from high levels of noise.

In conclusion, solving the weak Poisson problem with AQ is computationally cheaper and shows a considerably faster and smoother convergence than with MC.

#### 6.3.3. Weak Poisson in 2D

Based on the experiments in dimension one, we conclude that the integration points should be resampled around every \(10\) epochs and the learning rate should be set to \(10^{-2}\) to obtain optimal performances with both integration methods. We keep these hyperparameters in the rest of our experiments.

\begin{table}
\begin{tabular}{c c c c c c c c} \hline \hline \multirow{2}{*}{\(\eta\)} & \multirow{2}{*}{\(N_{\Omega}\left(P,O\right)\)} & \multicolumn{6}{c}{\(\nu\)} \\ \cline{3-8}  & & & \(1\) & \(10\) & \(100\) & \(500\) & \(1000\) & \(5000\) \\ \hline \multirow{4}{*}{MC} & \(50\) & \(4.82\times 10^{-2}\) & \(1.69\times 10^{-1}\) & \(1.41\times 10^{-0}\) & \(5.08\times 10^{-0}\) & \(8.96\times 10^{-0}\) & \(1.82\times 10^{+3}\) \\  & \(100\) & \(4.29\times 10^{-1}\) & \(6.39\times 10^{-1}\) & \(1.68\times 10^{-0}\) & \(1.40\times 10^{-0}\) & \(2.11\times 10^{+1}\) & \(7.26\times 10^{+3}\) \\ \cline{2-8}  & \(24\) (\(2,\,2\)) & \(2.87\times 10^{-2}\) & \(2.05\times 10^{-1}\) & \(8.80\times 10^{-2}\) & \(1.14\times 10^{-0}\) & \(3.54\times 10^{+2}\) & \(4.43\times 10^{+5}\) \\  & \(43\) (\(2,\,5\)) & \(1.04\times 10^{-2}\) & \(1.00\times 10^{-2}\) & \(4.90\times 10^{-2}\) & \(1.81\times 10^{-0}\) & \(2.04\times 10^{+2}\) & \(5.28\times 10^{+5}\) \\  & AQ & \(88\) (\(2,\,10\)) & \(8.51\times 10^{-3}\) & \(8.24\times 10^{-3}\) & \(6.89\times 10^{-2}\) & \(2.05\times 10^{-1}\) & \(1.10\times 10^{-0}\) & \(3.11\times 10^{+3}\) \\ \(10^{-2}\) & \(61\) (\(3,\,2\)) & \(8.12\times 10^{-3}\) & \(1.37\times 10^{-2}\) & \(1.59\times 10^{-2}\) & \(5.13\times 10^{-2}\) & \(4.83\times 10^{-0}\) & \(5.49\times 10^{+3}\) \\  & \(90\) (\(3,\,5\)) & \(8.74\times 10^{-3}\) & \(1.03\times 10^{-2}\) & \(1.02\times 10^{-2}\) & \(8.71\times 10^{-3}\) & \(7.79\times 10^{-3}\) & \(3.47\times 10^{+2}\) \\  & \(91\) (\(5,\,2\)) & \(9.67\times 10^{-3}\) & \(1.30\times 10^{-2}\) & \(6.35\times 10^{-3}\) & \(7.78\times 10^{-3}\) & \(1.49\times 10^{-2}\) & \(9.35\times 10^{+1}\) \\ \hline \multirow{4}{*}{MC} & \(50\) & \(1.91\times 10^{-1}\) & \(1.70\times 10^{-1}\) & \(9.88\times 10^{-1}\) & \(4.93\times 10^{-0}\) & \(4.97\times 10^{-0}\) & \(1.00\times 10^{+1}\) \\  & MC & \(100\) & \(2.73\times 10^{-3}\) & \(2.32\times 10^{-1}\) & \(1.37\times 10^{-0}\) & \(8.35\times 10^{-1}\) & \(6.46\times 10^{-1}\) & \(6.55\times 10^{+1}\) \\ \cline{2-8}  & \(26\) (\(2,\,2\)) & \(5.67\times 10^{-2}\) & \(3.50\times 10^{-2}\) & \(5.15\times 10^{-2}\) & \(1.26\times 10^{-0}\) & \(3.02\times 10^{-0}\) & \(2.32\times 10^{+3}\) \\ \cline{2-8}  & \(46\) (\(2,\,5\)) & \(8.66\times 10^{-3}\) & \(7.24\times 10^{-3}\) & \(1.98\times 10^{-2}\) & \(4.13\times 10^{-3}\) & \(1.70\times 10^{-1}\) & \(1.78\times 10^{+3}\) \\ \cline{2-8}  & \(93\) (\(2,\,10\)) & \(5.87\times 10^{-3}\) & \(1.04\times 10^{-2}\) & \(5.71\times 10^{-3}\) & \(5.57\times 10^{-3}\) & \(5.28\times 10^{+2}\) \\ \(10^{-3}\) & \(66\) (\(3,\,2\)) & \(2.30\times 10^{-2}\) & \(7.99\times 10^{-3}\) & \(1.42\times 10^{-2}\) & \(1.19\times 10^{-2}\) & \(3.26\times 10^{-2}\) & \(3.90\times 10^{-0}\) \\  & \(96\) (\(3,\,5\)) & \(1.97\times 10^{-2}\) & \(1.79\times 10^{-0}\) & \(3.78\times 10^{-2}\) & \(5.46\times 10^{-3}\) & \(5.53\times 10^{-3}\) & \(6.42\times 10^{-3}\) \\  & \(90\) (\(5,\,2\)) & \(1.78\times 10^{-2}\) & \(6.80\times 10^{-3}\) & \(6.61\times 10^{-3}\) & \(1.16\times 10^{-2}\) & \(9.72\times 10^{-3}\) & \(1.83\times 10^{-0}\) \\ \hline \hline \end{tabular}
\end{table}
Table 4. Comparison of the relative \(L^{2}\) norm for the weak Poisson problem with the \(\operatorname{ReLU}_{\varepsilon}\) activation, depending on learning rate, resampling frequency and integration method.

Figure 5. Learning curve (a) and pointwise error (b) for the strong Poisson problem with the \(\operatorname{ReLU}_{\varepsilon}\) activation. The learning rate is \(\eta=10^{-2}\) and the points are resampled every \(10\) epochs.
\begin{table}
\begin{tabular}{c c c c c c c c} \hline \hline \multirow{2}{*}{\(\eta\)} & \multirow{2}{*}{\(N_{\Omega}\left(P,O\right)\)} & \multicolumn{6}{c}{\(\nu\)} \\ \cline{3-8}  & & \multicolumn{2}{c}{1} & \multicolumn{2}{c}{10} & \multicolumn{2}{c}{100} & \multicolumn{2}{c}{500} & \multicolumn{2}{c}{1000} & \multicolumn{2}{c}{5000} \\ \hline \multirow{4}{*}{MC} & \(50\) & \(7.67\times 10^{-2}\) & \(4.99\times 10^{-1}\) & \(4.64\times 10^{-1}\) & \(2.86\times 10^{-0}\) & \(8.05\times 10^{-0}\) & \(3.45\times 10^{+2}\) \\  & \(100\) & \(1.75\times 10^{-1}\) & \(4.92\times 10^{-1}\) & \(9.31\times 10^{-1}\) & \(1.55\times 10^{-0}\) & \(1.64\times 10^{-0}\) & \(4.78\times 10^{-0}\) \\ \cline{3-8}  & \(40\,(3,2)\) & \(3.87\times 10^{-2}\) & \(1.49\times 10^{-2}\) & \(1.14\times 10^{-1}\) & \(1.34\times 10^{-0}\) & \(1.19\times 10^{+1}\) & \(4.21\times 10^{+2}\) \\  & \(46\,(3,5)\) & \(8.40\times 10^{-1}\) & \(3.06\times 10^{-2}\) & \(3.47\times 10^{-2}\) & \(2.19\times 10^{-2}\) & \(1.77\times 10^{-0}\) & \(6.64\times 10^{+1}\) \\  & AQ & \(110\,(3,10)\) & \(7.89\times 10^{-3}\) & \(4.56\times 10^{-3}\) & \(8.17\times 10^{-3}\) & \(5.09\times 10^{-2}\) & \(1.45\times 10^{-0}\) & \(3.92\times 10^{+2}\) \\ \(10^{-2}\) & \(84\,(5,2)\) & \(7.41\times 10^{-3}\) & \(1.36\times 10^{-2}\) & \(3.41\times 10^{-3}\) & \(4.18\times 10^{-2}\) & \(4.27\times 10^{-1}\) & \(4.34\times 10^{+2}\) \\  & \(108\,(5,5)\) & \(3.73\times 10^{-3}\) & \(4.77\times 10^{-3}\) & \(5.53\times 10^{-3}\) & \(9.31\times 10^{-3}\) & \(5.17\times 10^{-1}\) & \(3.05\times 10^{+2}\) \\  & \(98\,(7,2)\) & \(3.75\times 10^{-3}\) & \(9.94\times 10^{-3}\) & \(3.71\times 10^{-3}\) & \(5.95\times 10^{-2}\) & \(1.67\times 10^{-1}\) & \(1.66\times 10^{+2}\) \\ \hline \multirow{4}{*}{MC} & \(50\) & \(2.52\times 10^{-1}\) & \(2.54\times 10^{-1}\) & \(9.53\times 10^{-1}\) & \(2.76\times 10^{-0}\) & \(3.42\times 10^{-0}\) & \(1.47\times 10^{+1}\) \\  & MC & \(100\) & \(3.16\times 10^{-1}\) & \(3.77\times 10^{-1}\) & \(1.34\times 10^{-0}\) & \(2.56\times 10^{-0}\) & \(1.77\times 10^{-0}\) & \(1.87\times 10^{-0}\) \\ \cline{3-8}  & \(21\,(3,2)\) & \(1.59\times 10^{-1}\) & \(1.59\times 10^{-0}\) & \(1.21\times 10^{-1}\) & \(5.16\times 10^{-1}\) & \(5.55\times 10^{-0}\) & \(6.10\times 10^{+1}\) \\  & \(26\,(3,5)\) & \(5.44\times 10^{-1}\) & \(1.66\times 10^{-0}\) & \(3.47\times 10^{-1}\) & \(3.43\times 10^{-2}\) & \(1.54\times 10^{-1}\) & \(1.86\times 10^{-0}\) \\ AQ & \(57\,(3,10)\) & \(6.52\times 10^{-1}\) & \(6.20\times 10^{-1}\) & \(5.57\times 10^{-2}\) & \(4.92\times 10^{-2}\) & \(5.12\times 10^{-2}\) & \(3.86\times 10^{+1}\) \\  & \(50\,(5,2)\) & \(5.41\times 10^{-2}\) & \(3.90\times 10^{-2}\) & \(7.87\times 10^{-2}\) & \(5.77\times 10^{-2}\) & \(2.47\times 10^{-2}\) & \(5.02\times 10^{+1}\) \\  & \(71\,(5,5)\) & \(2.32\times 10^{-2}\) & \(2.31\times 10^{-2}\) & \(4.37\times 10^{-2}\) & \(6.62\times 10^{-2}\) & \(5.95\times 10^{-2}\) & \(2.18\times 10^{+1}\) \\  & \(41\,(7,2)\) & \(2.35\times 10^{-1}\) & \(1.81\times 10^{-1}\) & \(2.24\times 10^{-1}\) & \(1.08\times 10^{-1}\) & \(2.78\times 10^{-1}\) & \(1.50\times 10^{-0}\) \\ \hline \hline \end{tabular}
\end{table}
Table 5. Comparison of the relative \(L^{2}\) norm for the weak Poisson problem with the \(\tanh\) activation, depending on learning rate, resampling frequency and integration method.

Figure 6. Learning curve (a) and pointwise error (b) for the weak Poisson problem with the \(\tanh\) activation. The learning rate is \(\eta=10^{-2}\) and the points are resampled every \(10\) epochs.

Figure 7. Learning curve (a) and pointwise error (b) for the weak Poisson problem with the \(\operatorname{ReLU}_{\varepsilon}\) activation. The learning rate is set to \(10^{-2}\) and the points are resampled every \(10\) epochs.
We find that AQ can reach similar or lower errors with less than half the number of points MC needs. With the \(\mathrm{ReLU}_{\varepsilon}\) activation, AQ with fewer than \(1000\) points and MC with \(10000\) points reach comparable performances. In the case of the \(\tanh\) activation, AQ with around \(1660\) points wins over MC with \(10000\) points. However, we notice that increasing the order of the quadrature or the number of pieces does not significantly reduce the relative \(L^{2}\) norm.

Fig. 8 and Fig. 9 display the pointwise error between the ground truth and the solutions obtained with MC and AQ. In the case of \(\tanh\), we remark that MC struggles to model the transition between the high and low regions. Similarly for \(\mathrm{ReLU}_{\varepsilon}\), MC does not properly handle the region around the origin, where the solution shows larger variations. In both cases, the pointwise error for the model trained with AQ is much more homogeneous across the whole domain.

### Robustness to initialisation

In this second round of experiments, we study the robustness of our integration method to the initialisation of the network. We solve the weak Poisson problems in dimensions one and two. The learning rate is set to \(10^{-2}\) and we resample the integration points every \(10\) epochs. We report the minimum, maximum, average and standard deviation of the final \(L^{2}\) norm on \(10\) random initialisations, as well as the average training time in Tab. 6(a) and Tab. 6(b) for the \(\mathrm{ReLU}_{\varepsilon}\) and \(\tanh\) activations respectively. We use the same seeds for MC and AQ in both experiments. The number of integration points for MC is chosen so that it is larger than the maximum number of points for the same experiments run with AQ. The exact choices of quadrature order, number of pieces and integration points are reported in the corresponding tables.

We find that AQ has a consistently lower average \(L^{2}\) norm than MC. By comparing the standard deviation of the error, we infer that the solutions obtained by AQ are less dependent on the network initialisation compared to MC. Besides, we observe that the distributions of the errors obtained with AQ and MC do not overlap: the

\begin{table}

\end{table}
Table 6. Relative \(L^{2}\) norm for the weak Poisson problems in 2D with the \(\mathrm{ReLU}_{\varepsilon}\) (a) and \(\tanh\) (b) activations for various integration hyperparameters.

Figure 8. Comparison of the pointwise error for the weak Poisson problem with the \(\tanh\) activation. The learning rate is \(\eta=10^{-2}\) and the points are resampled every \(10\) epochs. MC is trained with \(5000\) points and AQ with \(5\) pieces and order \(3\) (\(2964\) points).
maximum error with AQ is lower than the minimum error with MC. More importantly, we find that AQ enables a reduction of the number of integration points while keeping similar or shorter training times.

We conclude that our proposed integration method is more robust to initialisation than MC and that it obtains higher accuracies than MC with fewer integration points.

### Reduction of the number of integration points

In this set of experiments, we want to show that it is possible to further reduce the number of integration points generated by AQ by merging small regions with their neighbours, without sacrificing the accuracy of our method. We keep the same hyperparameters as in the previous paragraph and focus on the two-dimensional problems.

Every time the mesh is updated, we compute the median of the size of the cells. We then identify the cells that are smaller than a given fraction of this median and merge them one by one with their neighbours. We always merge a cell with the largest of its neighbours. In dimension two, after all the cells are merged in this way, the mesh needs an extra post-processing step in order to account for non-convex cells. Indeed, the union of two convex polygons may be non-convex. We split non-convex aggregated cells with the ear clipping algorithm [22], and apply our splitting strategy to all resulting cells to obtain triangles and convex quadrangles.

We experiment with several merging thresholds and report our findings in Tab. 8. In dimension one, we find that merging regions up to \(50\%\) of the median of the regions size does not significantly affect the error while

\begin{table}

\end{table}
Table 7. Distribution of the relative \(L^{2}\) norm for the weak Poisson problem with the \(\mathrm{ReLU}_{\varepsilon}\) (a) and \(\tanh\) (b) activations on \(10\) random initialisations. The learning rate is set to \(10^{-2}\) and the points are resampled every \(10\) epochs.

Figure 9. Comparison of the pointwise error for the weak Poisson problem with the \(\mathrm{ReLU}_{\varepsilon}\) activation. The learning rate is set to \(10^{-2}\) and the points are resampled every \(10\) epochs. MC is trained with \(5000\) points and AQ with \(2\) pieces and order \(5\) (\(2681\) points).
allowing for a sizeable reduction of the number of integration points. In the two-dimensional case, raising the merging threshold above \(25\%\) harms the performance of the \(\mathrm{ReLU}_{\varepsilon}\) network, but not that of the \(\tanh\) network, as the number of integration points does not decrease.

It appears that increasing the merging threshold does not always reduce the number of integration points. Indeed, the aggregated regions may become less and less convex, so they will need to be split into several convex regions.

We illustrate the kind of meshes that a network generates and how they are merged in Fig. 10. These meshes are obtained from a trained \(\tanh\) network solving the weak Poisson problem. We found that the meshes from the previous round of experiments were easier to interpret, so we selected one of the \(10\) corresponding models. The boundary of the circular region can be easily identified. The lines that connect two sides of the square domain come from the first layer whereas the other lines come from the second layer. We observe that non-convex regions appear when small regions are merged.

### More complex domains

We complete our numerical experiments by solving a weak Poisson equation on a more complex two-dimensional domain, defined as the union of two rhombi. We solve the following Poisson

Figure 10. Final mesh, number of cells and number of integration points at several region merging thresholds for the weak Poisson problem and the \(\tanh\) activation. The \(\tanh\) activation is cut into \(5\) pieces and the quadrature is of order \(2\). The merging thresholds are \(0\%\) (a) and \(25\%\) (b).

\begin{table}
\begin{tabular}{c c c c} \hline \hline Dim. & & (\(N_{\Omega}\), \(N_{\Gamma}\)) (Threshold) & \(L^{2}\) norm \\ \hline \multirow{4}{*}{1} & \multirow{2}{*}{MC} & (200, 2) & \(6.39\times 10^{-1}\) \\  & & (100, 2) & \(4.55\times 10^{-1}\) \\ \cline{3-4}  & & (290, 2) (\(0\%\)) & \(8.25\times 10^{-3}\) \\  & & (198, 2) (\(50\%\)) & \(7.82\times 10^{-3}\) \\  & & (144, 2) (\(100\%\)) & \(4.41\times 10^{-2}\) \\ \hline \multirow{4}{*}{2} & \multirow{2}{*}{MC} & (5000, 500) & \(1.43\times 10^{-1}\) \\  & & (1000, 100) & \(2.35\times 10^{-1}\) \\ \cline{3-4}  & & (4041, 158) (\(0\%\)) & \(6.40\times 10^{-2}\) \\ \cline{3-4}  & & (3562, 148) (\(10\%\)) & \(7.27\times 10^{-2}\) \\ \cline{3-4}  & & (3217, 132) (\(25\%\)) & \(2.78\times 10^{-2}\) \\ \cline{3-4}  & & (2809, 115) (\(50\%\)) & \(1.41\times 10^{-1}\) \\ \hline \hline \end{tabular} \begin{tabular}{c c c} \hline Dim. & (\(N_{\Omega}\), \(N_{\Gamma}\)) (Threshold) & \(L^{2}\) norm \\ \hline \multirow{4}{*}{1} & \multirow{2}{*}{MC} & (200, 2) & \(2.34\times 10^{-1}\) \\  & & (100, 2) & \(4.92\times 10^{-1}\) \\ \cline{3-4}  & & (138, 2) (\(0\%\)) & \(3.08\times 10^{-3}\) \\  & AQ & (94, 2) (\(50\%\)) & \(4.18\times 10^{-3}\) \\  & & (85, 2) (\(100\%\)) & \(2.72\times 10^{-2}\) \\ \hline \multirow{4}{*}{2} & \multirow{2}{*}{MC} & (10000, 1000) & \(6.61\times 10^{-2}\) \\  & & (5000, 500) & \(8.51\times 10^{-2}\) \\ \cline{3-4}  & & (7798, 259) (\(0\%\)) & \(5.02\times 10^{-2}\) \\ \cline{3-4}  & AQ & (6895, 218) (\(10\%\)) & \(3.96\times 10^{-2}\) \\ \cline{3-4}  & & (8336, 244) (\(25\%\)) & \(5.26\times 10^{-2}\) \\ \cline{3-4}  & & (9088, 204) (\(50\%\)) & \(7.19\times 10^{-2}\) \\ \hline \end{tabular} \begin{tabular}{c c} \hline Dim. & (\(N_{\Omega}\), \(N_{\Gamma}\)) (Threshold) & \(L^{2}\) norm \\ \hline \multirow{4}{*}{1} & \multirow{2}{*}{MC} & (200, 2) & \(2.34\times 10^{-1}\) \\  & & (100, 2) & \(4.92\times 10^{-1}\) \\ \cline{3-4}  & AQ & (138, 2) (\(0\%\)) & \(3.08\times 10^{-3}\) \\ \cline{3-4}  & AQ & (94, 2) (\(50\%\)) & \(4.18\times 10^{-3}\) \\ \cline{3-4}  & & (85, 2) (\(100\%\)) & \(2.72\times 10^{-2}\) \\ \hline \multirow{4}{*}{2} & \multirow{2}{*}{MC} & (10000, 1000) & \(6.61\times 10^{-2}\) \\  & & (5000, 500) & \(8.51\times 10^{-2}\) \\ \cline{3-4}  & & (7798, 259) (\(0\%\)) & \(5.02\times 10^{-2}\) \\ \cline{3-4}  & AQ & (6895, 218) (\(10\%\)) & \(3.96\times 10^{-2}\) \\ \cline{3-4}  & & (8336, 244) (\(25\%\)) & \(5.26\times 10^{-2}\) \\ \cline{3-4}  & & (9088, 204) (\(50\%\)) & \(7.19\times 10^{-2}\) \\ \hline \end{tabular} \begin{tabular}{c c} \hline Dim. & (\(N_{\Omega}\), \(N_{\Gamma}\)) (Threshold) & \(L^{2}\) norm \\ \hline \multirow{4}{*}{1} & \multirow{2}{*}{MC} & (200, 2) & \(2.34\times 10^{-1}\) \\  & & (100, 2) & \(4.92\times 10^{-1}\) \\ \cline{3-4}  & AQ & (138, 2) (\(0\%\)) & \(3.08\times 10^{-3}\) \\ \cline{3-4}  & AQ & (94, 2) (\(50\%\)) & \(4.18\times 10^{-3}\) \\ \cline{3-4}  & & (85, 2) (\(100\%\)) & \(2.72\times 10^{-2}\) \\ \hline \multirow{4}{*}{2} & \multirow{2}{*}{MC} & (10000, 1000) & \(6.61\times 10^{-2}\) \\  & & (5000, 500) & \(8.51\times 10^{-2}\) \\ \cline{3-4}  & AQ & (7798, 259) (\(0\%\)) & \(5.02\times 10^{-2}\) \\ \cline{3-4}  & AQ & (6895, 218) (\(10\%\)) & \(3.96\times 10^{-2}\) \\ \cline{3-4}  & & (8336, 244) (\(25\%\)) & \(5.26\times 10^{-2}\) \\ \cline{3-4}  & & (9088, 204) (\(50\%\)) & \(7.19\times 10^{-2}\) \\ \hline \end{tabular} 
\begin{tabular}{c c} \hline Dim. & (\(N_{\Omega}\), \(N_{\Gamma}\)) (Threshold) & \(L^{2}\) norm \\ \hline \multirow{4}{*}{1} & \multirow{2}{*}{MC} & (200, 2) & \(2.34\times 10^{-1}\) \\  & & (100, 2) & \(4.92\times 10^{-1}\) \\ \cline{3-4}  & AQ & (138, 2) (\(0\%\)) & \(3.08\times 10^{-3}\) \\ \cline{3-4}  & AQ & (94, 2) (\(50\%\)) & \(4.18\times 10^{-3}\) \\ \cline{3-4}  & AQ & (85, 2) (\(100\%\)) & \(2.72\times 10^{-2}\) \\ \hline \multirow{4}{*}{2} & \multirow{2}{*}{MC} & (10000, 1000) & \(6.61\times 10^{-2}\) \\  & & (5000, 500) & \(8.51\times 10^{-2}\) \\ \cline{3-4}  & AQ & (779
equation

\[\left\{\begin{array}{rl}-\Delta u=x+y&\text{in }\Omega\\ u=0&\text{on }\Gamma\end{array}\right.. \tag{7}\]

We obtain an approximation of the solution with the FEM. We rely on the following weak formulation for the FEM solution:

\[\left\{\begin{array}{rl}\text{Find }u\in H_{0}^{1}(\Omega)\text{ such that for all }v\in H_{0}^{1}(\Omega),\\ \int_{\Omega}\Delta u\cdot\Delta v\,\mathrm{d}\Omega=\int_{\Omega}(x+y)v\, \mathrm{d}\Omega\end{array}\right.,\]

where \(H_{0}^{1}(\Omega)\) is the subset of functions in \(H^{1}(\Omega)\) that vanish on \(\Gamma\). We then obtain a mesh composed of \(3894\) cells from GMSH and solve this FEM problem with the Gridap software [23]. We rely on first-order Lagrange finite elements, and the total number of degrees of freedom is \(1852\). We use the FEM solution as the reference to compute the \(L^{2}\) norm for the solution obtained by a NN. The relative \(L^{2}\) norm is computed with a quadrature of order \(10\) on the cells of the mesh of the FEM solution.

We then approximate this equation with NNs using the same penalised weak form that we considered in the previous experiments. We increase the network architecture to \((2,20,20,1)\) and the activation function is \(\mathrm{ReLU}_{\varepsilon}\). We train our network for \(10000\) iterations using the ADAM optimiser with a learning rate of \(10^{-2}\). We cut the \(\mathrm{ReLU}_{\varepsilon}\) activation into \(3\) pieces, use a quadrature of order \(2\) in each cell, and do not merge small regions. The integration points are resampled every \(10\) epochs. On average, AQ involved \(5894\) integration points in the domain and \(292\) on the boundary. We match the computational budget of AQ for MC by sampling \(6000\) in-domain points and \(300\) points on the boundary for MC. To obtain a uniform distribution of integration points for MC, we first perform a triangulation of the domain. Then for as many points as we need, we draw a triangle at random by weighting them by their measure, and sample one point in the corresponding triangle.

We plot the solutions obtained by FEM, MC and AQ in Fig. 10(a). The pointwise difference between MC and FEM, and between AQ and FEM is shown in Fig. 10(b). The relative \(L^{2}\) norm is \(7.99\times 10^{-2}\) for MC and \(4.68\times 10^{-2}\) for AQ. We mention that the training times for the last experiment are \(700\) seconds with AQ and \(625\) seconds with MC. We believe that this \(12\%\) extension in time is well worth the \(70\%\) decrease in \(L^{2}\) norm. It is interesting to mention that the pointwise error is located around the same regions in the domain. However when using AQ, the magnitude of the local maxima of the pointwise error is much lower.

### Summary

Our numerical experiments confirmed the effectiveness of our proposed integration method in a wide range of scenarios. AQ enables to reach lower generalisation errors compared to MC while using fewer integration points. We have shown that this new integration method is more robust to the initial parameters of the network. Moreover, the convergence of the network is less noisy and the generalisation error decays quicker, as AQ can sustain higher learning rates without introducing extra noise in the learning phase. We are convinced that our proposed integration method can help NNs become a more reliable and predictable tool when it comes to approximating PDEs, so that they can be employed to solve real-world problems.

The most significant differences between MC and AQ were observed for the weak Poisson problem. One reason could be that the loss of the strong Poisson problems can be understood in a pointwise sense. In this case, MC minimises the pointwise distance between \(\Delta u\) and \(\Delta f\) at the location of the integration points. On the contrary, the loss of the weak Poisson problem is to be understood and minimised in a global sense. Our experiments also suggest that the benefits of our adaptive quadrature are more striking at the location of strong gradients.

### Discussion

Our approach suffers from the curse of dimensionality and is intended to be used for low-dimensional PDEs as a nonlinear alternative to the FEM. We noticed that the number of integration points in dimension one was of the order of a few tens to a few hundreds, when the usual number of integration points used for MC in the literature is rather a few hundreds to a few thousands. In dimension two, our method needs a few thousand integration points, whereas the models in the literature rely on a few tens of thousands of points. Even though our approach is still cheaper than MC in the two-dimensional case, the reduction factor is much smaller compared to the one-dimensional scenario. We found that merging linear regions can help reduce the number of linear regions and thus integration points. This phenomenon should amplify as the dimension increases, but we are only ever concerned with problems arising from physics, which involve space and time only. We would like to extend our method to the three-dimensional case in the future.

We have observed that the time spent on extracting the mesh is offset by the reduction in the number of integration points, as the overall training very often takes less time with AQ than with MC.

Our idea to regularise the activation function was motivated by the relationship between FEM and \(\mathrm{ReLU}\) networks, as well as the need for a smooth activation so that the energy functional itself is smooth. We believe
that it is more economical to use a smoothed CPWL activation function than to directly use a smooth activation such as \(\tanh\) because intuitively fewer pieces are necessary to reach a given tolerance in the former case. Whereas most of the current implementations of PINN and its variants can only handle tensor-product domains or balls, our proposed linearisation algorithm can naturally handle any polygonal domain. This is an important step towards the use of NNs solvers for practical issues.

One direction for future work is tackling the variational setting in the general case when it is not possible to recast the problem as an energy minimisation. We point out that our domain decomposition strategy and adaptive quadrature can be extended to higher-order approximations (e.g. splines) of the activation function, which would result in curved boundaries for the regions.

## 7. Conclusion

In this work, we introduce a procedure to define adaptive quadratures for smooth NNs suitable for low-dimensional PDE discretisation. It relies on a decomposition of the domain into regions where the network is almost linear. We also show the importance of smoothness in the training process and propose a regularisation of the standard ReLU activation function.

We carry out the numerical analysis of the proposed framework to obtain upper bounds of the integration error. Numerical experimentation shows that our integration method helps make the convergence less noisy and quicker compared to Monte Carlo integration. We observe that the number of integration points needed to

Figure 11. Approximation (a) and pointwise error (b) for the solution of Eq. (7).
compute the loss function can be significantly reduced with our method and that our integration strategy is more robust to the initialisation of the network. We illustrated the benefit of our adaptive integration on the Poisson equation, but it can be extended to any PDE involving a symmetric and coercive bilinear form.

While NN solvers are applied to increasingly intricate problems, they still lack reliability and theoretical guarantees before being used on real-world problems. Besides, most PINN frameworks can only handle tensor-product domains, while the FEM can be applied to much more complex domains. The adaptive quadratures we propose, the error bounds we prove, and the seamless possibility to handle polygonal domains of our method represent a few steps further in bridging these gaps.

For the sake of reproducibility, our implementation of the proposed adaptive quadrature is openly available at [24].

## Acknowledgements

This research was partially funded by the Australian Government through the Australian Research Council (project number DP220103160).

## Appendix A Proofs

### Proof of Lemma 1

Proof.: Let \(f\in\mathcal{A}\) and \(I\in P(f)\) be a compatible interval for the function \(f\). If \(f\) is linear in \(I\), then for all \(x<y\in I\), \(T[f,x]\) coincides with \(f\) on \([x,y]\), and so do \(T[f,y]\) and \(T[f,x,y]\). Besides \(f^{\prime\prime}=0\) so inequalities Eq. (3) and Eq. (4) hold trivially. We thus suppose that \(f\) is either strictly convex or strictly concave in \(I\), which means that \(f^{\prime\prime}\) is zero at the two ends of \(I\) and nowhere else in \(I\).

For \(\delta\geq 0\), we define the set \(E_{\delta}=\{(x,y)\in I^{2},y-x\geq\delta\}\). We also write \(E=E_{0}\). Let now \(\delta>0\) be fixed. For all \((x,y)\in E_{\delta}\), the quantity \(\|f^{\prime\prime}\|_{L^{p/(2p+1)}(]x,y]}\) is not zero since \(y>x\) and the integrand only vanishes on a set of zero measure. This enables us to consider the function

\[r:E_{\delta}\ni(x,y)\mapsto\frac{\|f-T[f,x,y]\|_{L^{p}(]x,y]}}{\|f^{\prime \prime}\|_{L^{p/(2p+1)}(]x,y]}}.\]

We want to show that \(r\) can be defined on \(E\) as a whole, and that \(r\) is bounded on \(E\).

_(i) Continuity and differentiability of \(c\) on \(E_{\delta}\)._ By definition of \(I\), for all \(x<y\in I\), the tangents to \(f\) at \(x\) and \(y\) intersect at an abscissa that we write \(c(x,y)\), with \(x<c(x,y)<y\). More precisely, for finite \(x<y\in I\), the function \(c\) has the following expression

\[c(x,y)=-\frac{[f(y)-yf^{\prime}(y)]-[f(x)-xf^{\prime}(x)]}{f^{\prime}(y)-f^{ \prime}(x)},\]

in which the denominator is not zero because \(f^{\prime}\) is injective in \(I\). We extend \(c\) to possible infinite values of \(x\) or \(y\) by setting \(c(x,+\infty)=\lim_{y\to+\infty}c(x,y)\) and \(c(-\infty,y)=\lim_{x\to-\infty}c(x,y)\). Since \(f\) and \(f^{\prime}\) are continuous on \(I\), and since \(f^{\prime}\) is injective on \(I\), we conclude that \(c\) is continuous on \(E_{\delta}\). By using the differentiation rule for quotients, similar arguments show that \(c\) is differentiable on \(E_{\delta}\).

_(ii) Continuity of \(r\) on \(E_{\delta}\)._ Since \(f\) and \(f^{\prime}\) are continuous on \(\mathbb{R}\), the map \(x\mapsto T[f,x]\) is continuous. Besides, \(c\) is continuous so the map \((x,y)\to T[f,x,y]\) is also continuous. Finally, any map \((x,y)\to\int_{x}^{y}\varphi(t,x,y)\,\text{d}t\) is continuous whenever \(\varphi\) is continuous. Here, \((x,y)\mapsto f-T[f,x,y]\) and \(f^{\prime\prime}\) are continuous, so we conclude that the numerator and the denominator of \(r\) are continuous. Moreover, since \(y-x\geq\delta>0\), the denominator of \(r\) is not zero, which makes \(r\) continuous on \(E_{\delta}\).

The fact that \(f\) belongs to \(\mathcal{A}\) ensures that the numerator of \(r\) is bounded. Besides, the denominator is a strictly increasing function of \(y-x\), which is bounded from below by \(\delta\,\|f^{\prime\prime}\|_{L^{\infty}(]x,y[]}\), and from above because \(f^{\prime\prime}\) is in \(L^{p/(2p+1)}(\mathbb{R})\). This proves that \(r(x,y)\) is bounded on \(E_{\delta}\). Thus the only case where \(r\) could be unbounded is when \(x\) and \(y\) get arbitrarily close (\(\delta\) goes to zero).

_(iii) Continuity and differentiability of \(c\) on \(E\)._ We now show that \(c\) can be continuously extended to \(E\) and that this extension is differentiable on \(E\). To do so, we prove that for all element \(x\in I\), \(\lim_{y\to x}c(x,y)=x\). One can prove that for all \(y\in I\), \(\lim_{x\to y}c(x,y)=y\) using a similar argument. Let \(\varepsilon\leq y-x\), so that \(x+\varepsilon\in I\). We start by rearranging the terms in the expression of \(c(x,x+\varepsilon)\) and obtain the following

\[c(x,x+\varepsilon)=x+\frac{f^{\prime}(x+\varepsilon)-\varepsilon^{-1}(f(x+ \varepsilon)-f(x))}{f^{\prime}(x+\varepsilon)-f^{\prime}(x)}\varepsilon.\]
Now, let \(k\geq 2\) be the smallest integer such that the \(k\)-th derivative of \(f\) at \(x\) is not zero. By definition of a compatible interval, \(k\) will be equal to \(2\) for all interior points of \(I\). Since \(f\) is analytical in the neighbourhood of the zeros of \(f^{\prime\prime}\), there also exists such a \(k>2\) for the ends of \(I\). We can always take \(\varepsilon\) as small as needed so that \(f\) is analytical around \(x\) in a neighbourhood of radius \(\varepsilon\). We obtain a Taylor expansion of the fraction above by expanding its numerator and denominator to order \(k-1\) around \(x\), expressed at \(x+\varepsilon\), and arrive at

\[c(x,x+\varepsilon)=x+\frac{k-1}{k}\varepsilon+O(\varepsilon^{2}).\]

This shows that the function \(\hat{c}\) defined by \(\hat{c}(x,y)=c(x,y)\) if \(x<y\) and \(\hat{c}(x,x)=x\) is continuous and differentiable on \(E\).

_(iv) Continuity of \(r\) on \(E\)._ Let \(x\in I\). We show that \(\varepsilon\mapsto r(x,x+\varepsilon)\) has a finite limit when \(\varepsilon\) goes to zero by finding an equivalent of the numerator and the denominator of \(r\). We write the numerator as \((I_{1}(\varepsilon)^{p}+I_{2}(\varepsilon)^{p})^{1/p}\), where \(I_{1}\) and \(I_{2}\) are the integrals on \([x,c(x,x+\varepsilon)]\) and \([c(x,x+\varepsilon),x+\varepsilon]\) respectively. The denominator of \(r\) is written \(I_{3}(\varepsilon)\). Here again, let \(k\geq 2\) be the smallest integer such that the \(k\)-th derivative of \(f\) at \(x\) is not zero.

We use the mean-value form of the remainder in the \(k\)-th-order Taylor expansion of \(f\) to obtain that for all \(t\in\left]x,c(x,x+\varepsilon)\right[\), there exists \(\eta_{t}\in[x,t]\) such that \(f(t)-T[f,x](t)=\frac{1}{k!}f^{(k)}(\eta_{t})(t-x)^{k}\). Then, owing to the mean-value theorem for integrals since \(f^{(k)}\) is continuous, there exists \(\eta\in\left]x,c(x,x+\varepsilon)\right[\) such that

\[I_{1}(\varepsilon)^{p}=\int_{x}^{c(x,x+\varepsilon)}\frac{1}{k!}|f^{(k)}(\eta _{t})(t-x)^{k}|^{p}\,\mathrm{d}t=\frac{1}{k!^{p}(kp+1)}|f^{(k)}(\eta)|^{p}(c(x,x+\varepsilon)-x)^{kp+1}.\]

We can now use the previous Taylor expansion of \(c\) around \((x,x)\) and also expand \(f^{(k)}(\eta)\) around \(x\) since \(\eta\) goes to \(x\) as \(\varepsilon\) approaches zero. We finally obtain the following expansion for \(I_{1}(\varepsilon)^{p}\)

\[I_{1}(\varepsilon)^{p}=\frac{1}{k!^{p}(kp+1)}\left(\frac{k-1}{k}\right)^{kp+1 }|f^{(k)}(x)|^{p}\varepsilon^{kp+1}+o(\varepsilon^{kp+1}).\]

We obtain a \(k\)-th order Taylor expansion of \(f-T[f,x+\varepsilon]\) as follows: there exist \(\eta_{1},\eta_{2}\in\left]x,x+\varepsilon\right[\), such that for all \(t\in\left]c(x,x+\varepsilon),x+\varepsilon\right[\), there exists \(\eta_{t}\) between \(x\) and \(t\) such that

\[f(t)-T[f,x+\varepsilon](t)=\frac{1}{k!}\left[f^{(k)}(\eta_{t})(t-x)^{k}-kf^{( k)}(\eta_{2})\varepsilon^{k-1}(t-x)+kf^{(k)}(\eta_{2})\varepsilon^{k}-f^{(k)}( \eta_{1})\varepsilon^{k}\right]+o(\varepsilon^{k}).\]

We observe that \(\eta_{1}\), \(\eta_{2}\) and \(\eta_{t}\) go to \(x\) as \(\varepsilon\) approaches zero. This means that the Taylor expansion of \(f^{(k)}(\eta_{t})\) is \(f^{(k)}(x)+O(\varepsilon)\), and similarly at \(\eta_{1}\) and \(\eta_{2}\). Besides, using the expansion of \(c(x,x+\varepsilon)\), we show that \(x+\frac{k-1}{k}\varepsilon+o(\varepsilon)\leq t\leq x+\varepsilon\), which means that \(t-x\) is of the order of \(\varepsilon\). After keeping the terms of order \(\varepsilon^{k}\), factoring by \(f^{(k)}(x)\varepsilon^{k}\) and introducing the change of variable \(u=\frac{t-x}{\varepsilon}\), we obtain

\[d_{2}(t)=\frac{f^{(k)}(x)}{k!}\left[u^{k}-ku+k-1\right]\varepsilon^{k}+o( \varepsilon^{k}).\]

The bounds of the integral become \(\frac{k-1}{k}+o(1)\) and \(1\), and we decompose the integration domain into \(\left[\frac{k-1}{k},1\right]\) and \(\left[\frac{k-1}{k}+o(1),\frac{k-1}{k}\right]\). The second integral is going to be negligible against the first because its integrand is bounded and its integration domain has a length of \(o(1)\). Putting everything together, we obtain the following expansion for \(I_{2}(\varepsilon)^{p}\)

\[I_{2}(\varepsilon)^{p}=\frac{1}{k!^{p}}\left(\int_{\frac{k-1}{k}}^{1}\left|u^{ k}-1-k(u-1)\right|^{p}\,\mathrm{d}u\right)|f^{(k)}(x)|^{p}\varepsilon^{kp+1}+o( \varepsilon^{kp+1}).\]

Using a similar argument, the Taylor expansion of the integrand of the denominator of \(r\) is \(\frac{1}{(k-2)!}f^{(k)}(x)(t-x)^{k-2}+o((t-x)^{k-2})\) and taking the \(L^{p/(2p+1)}\) norm of this relationship, we find that the expansion of \(I_{3}(\varepsilon)^{p}\) is

\[I_{3}(\varepsilon)^{p}=\frac{1}{(k-2)!^{p}}\left(\frac{2p+1}{kp+1}\right)^{2p+1 }|f^{(k)}(x)|^{p}\varepsilon^{kp+1}+o(\varepsilon^{kp+1}).\]

Altogether, in the expression of \(r(x,x+\varepsilon)^{p}\), we can factor the numerator and the denominator by \(|f^{(k)}(x)|^{p}\varepsilon^{kp+1}\) and taking the limit as \(\varepsilon\) goes to zero, we obtain

\[\lim_{\varepsilon\to 0}r(x,x+\varepsilon)^{p}=\frac{1}{(k(k-1))^{p}}\left(\frac{kp+1}{2p +1}\right)^{2p+1}\left(\frac{1}{kp+1}\left(\frac{k-1}{k}\right)^{kp+1}+\int_{ \frac{k-1}{k}}^{1}\left|u^{k}-1-k(u-1)\right|^{p}\,\mathrm{d}u\right).\]

This shows that the function \(\hat{r}\) defined by \(\hat{r}(x,y)=r(x,y)\) if \(x<y\) and \(\hat{r}(x,x)\) as the quantity above is continuous on \(E\) as a whole. Besides \(\hat{r}\) has finite limits on \(\partial E\), so we conclude that it is bounded in \(E\).
_(v) Proof of the other inequalities_. The same proof can be adapted for the bound on the derivatives, and the bound in \(L^{\infty}\) norms. The arguments are the same everywhere, except for the boundedness around the points of the type \((x,x)\), where other equivalents are found. It is still possible to define the corresponding function \(r\) on \(E\) as a whole and to show that it is bounded in \(E\). 

### Proof of Proposition 1

Proof.: This proof is inspired by the heuristics given in [18]. We keep the notations of the definition in Subsect. 4.2. Let \(I\) be a compatible interval of \(f\). If \(f\) is linear on \(I\), then \(f\) coincides with its tangent at either end of \(I\) so it is enough to define the restriction of \(\pi_{n,p}[f]\) to \(I\) as this tangent. We thus decompose the total squared \(L^{p}\) norm on the \(\kappa(f)\) compatible intervals of \(f\) where \(f\) is not linear. We decide to place \(n^{\prime}=\lfloor n/\kappa(f)\rfloor\) points in each interval.

Let \(I=[a,b]\) be a compatible interval for \(f\) where \(f\) is not linear. We write \(a<\xi_{1}<\ldots<\xi_{n^{\prime}}<b\) the free points in \(I\). For convenience we introduce \(\xi_{0}=a\) and \(\xi_{n^{\prime}+1}=b\). We split the integral on \(I\) on the intervals \([\xi_{i},\xi_{i+1}]\) and using Eq. (3) from Lem. 1, we readily have

\[\|f-\pi_{n,p}[f]\|_{L^{p}(I)}^{p}=\sum_{i=0}^{n^{\prime}}\|f-T[f,\xi_{i},\xi_{i +1}]\|_{L^{p}([\xi_{i},\xi_{i+1}])}^{p}\leq A_{p}(f)^{p}\sum_{i=0}^{n^{\prime} }\left\|f^{\prime\prime}\right\|_{L^{p/(2p+1)}([\xi_{i},\xi_{i+1}])}^{2p+1}.\]

We now define the \((\xi_{i})_{1\leq i\leq n^{\prime}}\) such that \(\left\|f^{\prime\prime}\right\|_{L^{p/(2p+1)}([\xi_{i},\xi_{i+1}])}^{p/(2p+1)}\) are all equal to the same constant that we write \(C_{I}\). By the linearity of the integral, we find that \((n^{\prime}+1)C_{I}=\left\|f^{\prime\prime}\right\|_{L^{p/(2p+1)}(I)}^{p/(2p+1 )}\) and \(\xi_{i}\) is such that

\[\int_{\xi_{0}}^{\xi_{i}}\left|f^{\prime\prime}(t)\right|^{p/(2p+1)}\,\mathrm{d }t=iC_{I}=\frac{i}{n^{\prime}+1}\int_{I}\left|f^{\prime\prime}(t)\right|^{p/(2p +1)}\,\mathrm{d}t.\]

The sum above evaluates to \((n^{\prime}+1)C_{I}^{2p+1}=(n^{\prime}+1)^{-2p}\left\|f^{\prime\prime}\right\| _{L^{p/(2p+1)}(I)}^{p}\). We apply the same method in all the compatible intervals for \(f\) that we write \(I_{k}\) for \(1\leq k\leq\kappa(f)\). Since \(2/5\leq p/(2p+1)\leq 1/2\) and \(f^{\prime\prime}\in L^{q}\) for all \(q\geq 2/5\), we bound the \(L^{p/(2p+1)}\) norm of \(f^{\prime\prime}\) on each \(I_{k}\) by that on \(\mathbb{R}\) as a whole. We finally obtain

\[\|f-\pi_{n,p}[f]\|_{L^{p}(\mathbb{R})}^{p} =\sum_{k=1}^{\kappa(f)}\|f-\pi_{n,p}[f]\|_{L^{p}(I_{k})}^{p}\leq \sum_{k=1}^{\kappa(f)}A_{p}(f)^{p}(n^{\prime}+1)C_{I_{k}}^{2p+1}\] \[\leq\kappa(f)\frac{A_{p}(f)^{p}\left\|f^{\prime\prime}\right\|_{L ^{p/(2p+1)}(\mathbb{R})}^{p}}{(n^{\prime}+1)^{2p}}\leq\frac{A_{p}(f)^{p} \kappa(f)^{2p+1}\left\|f^{\prime\prime}\right\|_{L^{p/(2p+1)}(\mathbb{R})}^{p }}{n^{2p}},\]

where we obtained the last inequality by multiplying the numerator and denominator by \(\kappa(f)^{2p}\) and used the fact that for all \(m,n>0\), \(m\left\lfloor n/m\right\rfloor\geq n-m\). We apply the same method to \(f^{\prime}\) using Eq. (3) from Lem. 1. In each \(I_{k}\), we sample the \((\xi_{i})_{1\leq i\leq n^{\prime}}\) such that \(\int_{\xi_{0}}^{\xi_{i}}\left|f^{\prime\prime}(t)\right|^{p/(p+1)}\,\mathrm{d}t =iD_{I_{k}}\) where \((n^{\prime}+1)D_{I_{k}}=\left\|f^{\prime\prime}\right\|_{L^{p/(p+1)}(I_{k})}^{p/ (p+1)}\). Then it holds

\[\|f-\pi_{n,p}[f]\|_{L^{p}(\mathbb{R})}^{p}\leq\sum_{k=1}^{\kappa(f)}B_{p}(f)^{ p}(n^{\prime}+1)D_{I_{k}}^{p+1}\leq\frac{B_{p}(f)^{p}\kappa(f)^{p+1}\left\|f^{ \prime\prime}\right\|_{L^{p/(p+1)}(\mathbb{R})}^{p}}{n^{p}}.\]

The same bounds are obtained for the \(L^{\infty}\) norms, by considering the maximum \(L^{\infty}\) norm over the \(n^{\prime}+1\) subsegments and sampling the \((\xi_{i})\) according to \(\int_{\xi_{0}}^{\xi_{i}}\left|f^{\prime\prime}\right|^{1/2}\) for the first bound and \(\int_{\xi_{0}}^{\xi_{i}}\left|f^{\prime\prime}\right|\) for the second bound. 

### Proof of Lemma 2

Proof.: We know that the cells in \(\tau_{\mathbb{R}^{d}}(\boldsymbol{\Theta})\) are all convex since they are defined as the interior of intersections of half-planes. The final mesh \(\tau_{\Omega}(\boldsymbol{\Theta})\) is the intersection of \(\tau_{\mathbb{R}^{d}}(\boldsymbol{\Theta})\) against \(\Omega\). Let \(K\in\tau_{\Omega}(\boldsymbol{\Theta})\). There exists \(K^{\prime}\in\tau_{\mathbb{R}^{d}}(\boldsymbol{\Theta})\) such that \(K=K^{\prime}\cap\Omega\). If \(K^{\prime}\subseteq\Omega\), then \(K=K^{\prime}\) is convex. By the contraposition principle, if \(K\) is not convex, then \(K^{\prime}\nsubseteq\Omega\). This shows that if \(K\) is a non-convex cell in \(\tau_{\Omega}(\boldsymbol{\Theta})\), then \(K\) intersects with both \(\Omega\) and \(\mathbb{R}^{d}\backslash\Omega\), which means that \(K\) intersects with the boundary of \(\Omega\). In the case where \(\Omega\) is convex, the intersection of two convex sets being convex implies that \(K=K^{\prime}\cap\Omega\) is always convex.
### Proof of Proposition 2

Proof.: We show the two bounds Eq. (5) and Eq. (6) on a single-layer network, and extend them to arbitrary NNs by induction. We also prove that \(\left\|\nabla u\right\|_{L^{\infty}(\Omega)}\leq C_{u}\left\|f^{\prime}\right\|_ {L^{\infty}(\mathbb{R})}^{\alpha}\) and \(\left\|\nabla v\right\|_{L^{\infty}(\Omega)}\leq C_{v}\left\|g^{\prime}\right\| _{L^{\infty}(\mathbb{R})}^{\alpha}\) for some constants \(C_{u}\) and \(C_{v}\) that only depend on \(\boldsymbol{\Theta}\), and some integer \(\alpha\).

To simplify notations, we write \(\left\|\cdot\right\|_{\Omega}\) instead of \(\left\|\cdot\right\|_{L^{\infty}(\Omega)}\) and similarly \(\left\|\cdot\right\|_{\mathbb{R}}\) instead of \(\left\|\cdot\right\|_{L^{\infty}(\mathbb{R})}\). For an affine map \(\boldsymbol{\Theta}\) associated with a matrix \(\boldsymbol{W}\), we write \(\left\|\boldsymbol{\Theta}\right\|_{\infty}=\max_{i,j}\left|\boldsymbol{W}_{i,j}\right|\). In this proof, we write \(C_{1}\), \(C_{2}\), \(C_{3}\), \(C_{u}\), \(C_{v}\) and \(\alpha\) the constants in the hypothesis of the induction and write the new constants after the induction with a prime symbol.

_(i) Initialisation._ We take an affine map \(\boldsymbol{\Theta}\) from \(\mathbb{R}^{d}\) to \(\mathbb{R}^{m}\) associated with a weight matrix \(\boldsymbol{W}\in\mathbb{R}^{m\times d}\) and we consider the functions \(u=f\circ\boldsymbol{\Theta}\) and \(v=g\circ\boldsymbol{\Theta}\). The partial derivatives of \(u\) are given by \(\partial_{j}u_{i}=\boldsymbol{W}_{i,j}f^{\prime}\circ\boldsymbol{\Theta}_{i}\). We readily obtain that \(\left\|u_{i}-v_{i}\right\|_{\Omega}\leq\left\|f-g\right\|_{\mathbb{R}}\) and \(\left\|\partial_{j}u_{i}-\partial_{j}v_{i}\right\|_{\Omega}\leq\left|\boldsymbol {W}_{i,j}\right|\left\|f^{\prime}-g^{\prime}\right\|_{\mathbb{R}}\). We also notice that \(\partial_{j}u_{i}\) is bounded by \(\left|\boldsymbol{W}_{i,j}\right|\left\|f^{\prime}\right\|_{\mathbb{R}}\). Similarly, \(\partial_{j}v_{i}\) is bounded by an analogous constant. We conclude by taking the maximum on \(1\leq i\leq m\) and \(1\leq j\leq d\). We can take \(C_{1}=1\), \(C_{2}\geq 0\), \(C_{3}=C_{u}=C_{v}=\left\|\boldsymbol{\Theta}\right\|_{\infty}\) and \(\alpha=1\).

_(ii) Induction with affine map._ Let \(u\) and \(v\) be two functions from \(\mathbb{R}^{d}\) to \(\mathbb{R}^{n}\) such that the two bounds hold and \(\nabla u\) and \(\nabla v\) are bounded in \(\Omega\). Let \(\boldsymbol{\Theta}\) be an affine map from \(\mathbb{R}^{n}\) to \(\mathbb{R}^{m}\). We want to show that the two bounds still hold for \(\boldsymbol{\Theta}\circ u-\boldsymbol{\Theta}\circ v\). We use the triangle inequality and the induction hypothesis to obtain that

\[\left\|(\boldsymbol{\Theta}\circ u)_{i}-(\boldsymbol{\Theta}\circ v)_{i} \right\|_{\Omega}\leq\sum_{k=1}^{n}\left|\boldsymbol{W}_{i,k}\right|\left\|u_ {k}-v_{k}\right\|_{\Omega}\leq nC_{1}\left\|\boldsymbol{\Theta}\right\|_{ \infty}\left\|f-g\right\|_{\mathbb{R}}.\]

The partial derivatives of \(\boldsymbol{\Theta}\circ u\) are given by \(\partial_{j}(\boldsymbol{\Theta}\circ u)_{i}=\sum_{k=1}^{n}\boldsymbol{W}_{i,k}\partial_{j}u_{k}\). Using the triangle inequality and the induction hypothesis, we arrive at

\[\left\|\partial_{j}(\boldsymbol{\Theta}\circ u)_{i}-\partial_{j}(\boldsymbol {\Theta}\circ v)_{i}\right\|_{\Omega}\leq\sum_{k=1}^{n}\left|\boldsymbol{W}_{i,k}\right|\left\|\partial_{j}u_{k}-\partial_{j}v_{k}\right\|_{\Omega}\leq n \left\|\boldsymbol{\Theta}\right\|_{\infty}(C_{1}\left\|f-g\right\|_{ \mathbb{R}}+C_{2}\left\|f^{\prime}-g^{\prime}\right\|_{\mathbb{R}}).\]

We verify that \(\left\|\partial_{j}(\boldsymbol{\Theta}\circ u)_{i}\right\|_{\Omega}\leq C_{u }\sum_{k=1}^{n}\left|\boldsymbol{W}_{i,k}\right|\left\|f^{\prime}\right\|_{ \mathbb{R}}^{\alpha}\) and the same argument shows that \(\partial_{j}v_{i}\) is bounded by a similar constant. We conclude by taking the maximum on \(1\leq i\leq m\) and \(1\leq j\leq d\). We can take \(C_{1}^{\prime}=C_{2}^{\prime}=nC_{1}\left\|\boldsymbol{\Theta}\right\|_{ \infty}\), \(C_{3}^{\prime}=nC_{2}\left\|\boldsymbol{\Theta}\right\|_{\infty}\), \(C_{u}^{\prime}=nC_{u}\left\|\boldsymbol{\Theta}\right\|_{\infty}\), \(C_{v}^{\prime}=nC_{v}\left\|\boldsymbol{\Theta}\right\|_{\infty}\) and \(\alpha^{\prime}=\alpha\).

_(iii) Induction with activation._ Let \(u\) and \(v\) be two functions from \(\mathbb{R}^{d}\) to \(\mathbb{R}^{n}\) such that the two bounds hold and \(\nabla u\) and \(\nabla v\) are bounded in \(\Omega\). We want to show that the two bounds still hold for \(f\circ u-g\circ v\). Using the triangle inequality, the fact that \(f\) is Lipschitz and the induction hypothesis, we obtain

\[\left\|(f\circ u)_{i}-(g\circ v)_{i}\right\|_{\Omega} \leq\left\|(f\circ u)_{i}-(f\circ v)_{i}\right\|_{\Omega}+\left\| (f\circ v)_{i}-(g\circ v)_{i}\right\|_{\Omega}\] \[\leq\left(C_{1}\operatorname{Lip}(f)+1\right)\left\|f-g\right\|_{ \mathbb{R}}.\]

The partial derivatives of \(f\circ u\) are given by \(\partial_{j}(f\circ u)_{i}=(f^{\prime}\circ u_{i})\partial_{j}u_{i}\). The triangle inequality leads to

\[\left\|\partial_{j}(f\circ u)_{i}-\partial_{j}(g\circ v)_{i}\right\| _{\Omega} \leq\left\|(f^{\prime}\circ u_{i})\partial_{j}u_{i}-(f^{\prime} \circ u_{i})\partial_{j}v_{i}\right\|_{\Omega}+\left\|(f^{\prime}\circ u_{i}) \partial_{j}v_{i}-(f^{\prime}\circ v_{i})\partial_{j}v_{i}\right\|_{\Omega}\] \[+\left\|(f^{\prime}\circ v_{i})\partial_{j}v_{i}-(g^{\prime}\circ v _{i})\partial_{j}v_{i}\right\|_{\Omega}.\]

We bound the first term by factorising by \(f^{\prime}\circ u_{i}\) and the induction hypothesis. We use the fact that \(f^{\prime}\) is Lipschitz and the induction hypothesis to bound the second term. The third term is readily bounded by \(\left\|\nabla v\right\|_{\Omega}\left\|f^{\prime}-g^{\prime}\right\|_{\mathbb{R}}\). Altogether, we obtain the following bound

\[\left\|\partial_{j}(f\circ u)_{i}-\partial_{j}(g\circ v)_{i}\right\| _{\Omega} \leq\left(C_{2}\left\|f^{\prime}\right\|_{\mathbb{R}}+C_{1}C_{v} \operatorname{Lip}(f^{\prime})\left\|g^{\prime}\right\|_{\mathbb{R}}^{\alpha} \right)\left\|f-g\right\|_{\mathbb{R}}\] \[+\left(C_{3}\left\|f^{\prime}\right\|_{\mathbb{R}}+C_{v}\left\|g^{ \prime}\right\|_{\mathbb{R}}^{\alpha}\right)\left\|f^{\prime}-g^{\prime}\right\|_{ \mathbb{R}}.\]

We observe that \(\left\|\partial_{j}u_{i}\right\|_{\Omega}\leq C_{u}\left\|f^{\prime}\right\|_{ \mathbb{R}}^{\alpha+1}\) and \(\partial_{j}v_{i}\) is bounded by a similar constant. We conclude by taking the maximum on \(1\leq i\leq m\) and \(1\leq j\leq d\). We can take \(C_{1}^{\prime}=C_{1}\operatorname{Lip}(f)+1\), \(C_{2}^{\prime}=C_{2}\left\|f^{\prime}\right\|_{\mathbb{R}}+C_{1}C_{v} \operatorname{Lip}(f^{\prime})\left\|g^{\prime}\right\|_{\mathbb{R}}^{\alpha}\), \(C_{3}^{\prime}=C_{3}\left\|f^{\prime}\right\|_{\mathbb{R}}+C_{v}\left\|g^{ \prime}\right\|_{\mathbb{R}}^{\alpha}\), \(C_{u}^{\prime}=C_{u}\), \(C_{v}^{\prime}=C_{v}\) and \(\alpha^{\prime}=\alpha+1\). 

## Appendix B Algorithms

### Mesh extraction

In this section, we provide more detail on the way to construct an adapted mesh to a neural network based on a CPWL approximation of its activation function,
from \(\mathbb{R}\) to \(\mathbb{R}\) that we apply element-wise on the coordinates of \(\pi[u_{k}]_{|R}\), and let \((\xi_{j})_{1\leq j\leq K}\) denote the breakpoints of \(\pi[\rho]\). We also introduce the local coefficients \((\alpha_{j},\beta_{j})\) of \(\pi[\rho]\) such that the restriction of \(\pi[\rho]\) to \([\xi_{j},\xi_{j+1}]\) has the expression \(x\mapsto\alpha_{j}x+\beta_{j}\).

To identify the regions where \(\pi[\rho]\circ\pi[u_{k}]_{|R}\) is linear, we need to consider the equations \(\boldsymbol{W}_{R,i}\boldsymbol{x}+\boldsymbol{b}_{R,i}=\xi_{j}\) for all \(1\leq i\leq m\) and \(1\leq j\leq K\). These hyperplanes intersect in \(R\) and define subregions adapted to \(\pi[\rho]\circ\pi[u_{k}]_{|R}\). This is the step that corresponds to cutRegion in Alg. 1. We now present specific algorithms to define these subregions in dimensions zero, one and two.

#### b.1.1. Dimensions zero and one

In dimension zero, a region is reduced to one single point so nothing needs to be done. In the one-dimensional case, \(R\) is a segment that we write \([p_{-},p_{+}]\). We first compute the range of each linear map \(h_{i}:x\mapsto w_{i}x+b_{i}\) on \(R\) by computing their values at \(p_{-}\) and \(p_{+}\). We then find which breakpoints \(\xi_{j}\) belong to that range via a binary search, and we compute the corresponding antecedents \(p_{i,j}\) of \(\xi_{j}\) by \(h_{i}\) with normalised coordinates as follows

\[p_{i,j}=\frac{1}{2}(p_{-}+p_{+})+\frac{t_{i,j}}{2}(p_{+}-p_{-}),\qquad t_{i,j} =\frac{2\xi_{j}-(h_{i}(p_{+})+h_{i}(p_{-}))}{h_{i}(p_{+})-h_{i}(p_{-})}\in[-1,+ 1]\,.\]

We found this formulation more numerically stable compared to expressing \(p_{i,j}\) in terms of \(w_{i}\), \(b_{i}\) and \(\xi_{j}\). We sort the normalised coordinates \(t_{i,j}\) and obtain the corresponding sorted abscissae \((p_{l})\). For each segment \([p_{l},p_{l+1}]\), we evaluate the linear maps \(h_{i}\) at \(\frac{1}{2}(p_{l}+p_{l+1})\) and retrieve the coefficients \((\alpha_{l,i},\beta_{l,i})\) of \(\pi[\rho]\) at this location to update the local expression \((\boldsymbol{W}_{R},\boldsymbol{b}_{R})\) into \((\mathrm{diag}(\boldsymbol{\alpha}_{l})\boldsymbol{W}_{R},\mathrm{diag}( \boldsymbol{\alpha}_{l})\boldsymbol{b}_{R}+\boldsymbol{\beta}_{l})\).

#### b.1.2. Dimension two

The two-dimensional case is much more involved as the hyperplanes are lines that can intersect with the boundary of \(R\) and among one another. We consider four steps.

_(i) Intersection of the hyperplanes with \(\partial R\)_ Since \(R\) is a convex polygon, each hyperplane \(h_{i}:\boldsymbol{x}\mapsto\boldsymbol{w}_{i}\boldsymbol{x}+b_{i}\) can intersect with \(\partial R\) at two points at most. If the number of intersections is zero or one, then \(R\) lies entirely on one side of the hyperplane thus we only need to handle the case when there are two intersection points. One notable corner case occurs when the hyperplane coincides with an edge of the polygon. In this degenerate case, the region is not cut.

The clipping of a line by a convex polygon can be done in logarithmic time, using a method akin to binary search. We refer the reader to [25] for the full algorithm. We write \((P_{i,j}^{-},P_{i,j}^{+})\) the clipping of \(h_{i}\) by \(R\), and \((P_{l})\) the collection of these vertices.

_(ii) Intersection of the hyperplanes among themselves._ Next, we check the intersection of pairwise combinations of the lines obtained above. It is important to note that for a fixed row \(i\), the lines \((P_{i,j}^{-},P_{i,j}^{+})\) are parallel. Indeed, they are solutions to \(\boldsymbol{w}_{i}\boldsymbol{x}+b_{i}=\xi\) for two different values of \(\xi\). This means that we can skip the tests corresponding to lines arising from the same row. Let \(i_{1}\neq i_{2}\) be two different row indices, \(j_{1}\) and \(j_{2}\) be two breakpoint indices, and suppose that \((P_{i_{1},j_{1}}^{-},P_{i_{1},j_{1}}^{+})\) and \((P_{i_{2},j_{2}}^{-},P_{i_{2},j_{2}}^{+})\) are two non-parallel lines. We compute their intersection \(Q\) via the following parametrisation and solve for \(s\) and \(t\):

\[Q=\overline{P}_{i_{1},j_{1}}+s\Delta P_{i_{1},j_{1}}=\overline{P}_{i_{2},j_{2} }+t\Delta P_{i_{2},j_{2}},\qquad\left[\Delta P_{i_{1},j_{1}}\quad-\Delta P_{i_{2 },j_{2}}\right]\begin{bmatrix}s\\ t\end{bmatrix}=\left[\overline{P}_{1}-\overline{P}_{2}\right],\]

where \(\overline{P}_{i,j}=\frac{1}{2}(P_{i,j}^{-}+P_{i,j}^{+})\) is the middle of the segment and \(\Delta P_{i,j}=\frac{1}{2}(P_{i,j}^{+}-P_{i,j}^{-})\) is half the direction vector. We need to check that \(s\) and \(t\) are between \(-1\) and \(+1\).

_(iii) Representation of these intersections as an adjacency graph._ We initialise a connectivity graph with the edges of the region \(R\). During step _(i)_, for each line \((P_{i,j}^{-},P_{i,j}^{+})\) that cuts the polygon, we add the points \(P_{i,j}^{-}\) and \(P_{i,j}^{+}\) to the edges they belong to, and add an arc in the graph, that connects these two points. During step _(ii)_, we only need to add new vertices to the graph.

_(iv) Extraction of the smallest cycles in this graph._ Finding the partition of \(R\) where the composition by \(\pi[\rho]\) is CPWL amounts to finding the fundamental cycle basis of the undirected graph defined above. We refer to [26] for a definition of the fundamental cycle basis of a planar graph and for an algorithm to build it.

#### b.1.3. Notes on the algorithm

In Alg. 1, \(\mathtt{weight}(\ell)\) and \(\mathtt{bias}(\ell)\) denote the parameters of the network at layer \(\ell\). The instruction coefficients is just a lookup operation to retrieve the local coefficients of \(\pi[\rho]\). Here its output are vectors because the local map on a region is vector-valued. Finally if \(\boldsymbol{\alpha}\) is a vector, \(\mathtt{diag}(\boldsymbol{\alpha})\) is the diagonal matrix whose diagonal coefficients are \(\boldsymbol{\alpha}\).
### Decomposition of convex polygons

In the general two-dimensional case, the cells of the mesh adapted to \(u\) are arbitrary convex polygons. Numerical quadrature rules are known for triangles and convex quadrangles, so we decide to split the cells into a collection of these two reference shapes. We encode hard rules to decompose convex polygons with up to ten vertices into triangles and convex quadrangles. For example, an octagon is split into one triangle and two quadrangles by considering the three cycles \((1,2,3,4)\), \((4,5,6,7)\), \((7,8,1,4)\). Larger polygons are recursively split in half until they have fewer than ten vertices.
