Variational Physics Informed Neural Networks: the role of quadratures and test functions

Exploiting the recent advances in artificial intelligence and, in particular, in deep learning, several innovative numerical techniques have been developed in the last few years to compute numerical solutions of partial differential equations (PDEs). In such methods, the solution is approximated by a neural network that is trained by taking advantage of the knowledge of the underlying differential equation. One of the earliest models involving a neural network was described in [25]: it is based on the concept of Physics Informed Neural Networks (PINN) and it inspired further works such as e.g. [31] or [34], until the recent paper [19] which presents a very general framework for the solution of operator equations by deep neural networks.

In such papers, given an arbitrary PDE coupled with proper boundary conditions, the training of the PINN aims at finding the weights ${\mathbf{w}}$ of a neural network such that the associated function $u^{{\mathcal{NN}}}({\mathbf{x}};{\mathbf{w}})$ minimizes some functional of the equation residual while satisfying as much as possible the imposed boundary conditions. To do so, the neural network is trained to minimize the residual only at a finite set of collocation points and additional terms are added to the loss function in order to force the network to approximately satisfy the boundary conditions. Thanks to the good approximation properties of neural networks, formally proved e.g. in [8], [11], [23], [18], [24] and [10] under suitable assumptions, the PINN approach looks very promising because it is able to efficiently and accurately compute approximate solutions of arbitrary PDEs encoding their structures in the loss function.

Subsequently, the PINN paradigm has been further developed in [14] to obtain the so-called Variational Physics Informed Neural Networks (VPINN). The main differences with respect to the PINN are that the weak formulation of the PDE is exploited, the collocation points are replaced by test functions, and quadrature points are used to compute the integrals involved in the variational residuals. In such a method the solution is still approximated by a neural network, but the test functions are represented by a finite set of known functions or by a second neural network (see [35]); therefore, the technique can be seen as a Petrov-Galerkin method. The method is more flexible than the standard PINN because the integration by parts, involved in the weak formulation, decreases the required regularity of the approximate solution. Furthermore, the fact that the dataset used in the training phase consists of quadrature points significantly reduces the computational cost of the training phase. Indeed, quadrature points are, in general, much fewer than collocation points.

Combining the VPINN with the Finite Element Method (FEM), the authors of [16] developed VarNet, a VPINN that exploits the test functions of the $\mathbb{P}_{1}$-FEM. Such a work has been then extended in [15] to consider arbitrary high-order polynomials as test functions, as in the $hp$ version of the FEM. Although the authors of the cited works empirically observed that both PINNs and VPINNs are able to efficiently approximate the desired solution, no proof of a priori error estimates with convergence rates is provided for VPINNs. On the contrary, rigorous a posteriori error analyses are already available (see, for instance, [20]). Recently (see [3]) we derived a posteriori error estimates for the discretization setting considered in this paper.

The purpose of this paper is to investigate how the choice of piecewise polynomial test functions and quadrature formulas influence the accuracy of the resulting VPINN approximation of second-order elliptic problems. One might think that test functions of high polynomial degree are needed to get a high order of accuracy; we prove that this is not the case, actually we indicate that precisely the opposite is true: it is more convenient to keep the degree of the test functions as low as possible, while using quadrature formulas of precision as high as possible. Indeed, for sufficiently smooth solutions, the error decay rate is given by

$$q+2-k_{\text{test}}\,,$$

where $q$ is the precision of the quadrature formula and $k_{\text{test}}$ is the degree of the test functions.

Using a Petrov-Galerkin framework, we derive an a priori error estimate in the energy norm between the exact solution and a suitable piecewise polynomial interpolant of the computed neural network; we assume that the architecture of the neural network is fixed and sufficiently rich, and we explore the behaviour of the error versus the size of the mesh supporting the test functions. Our analysis relies upon the validity of an inf-sup condition between the spaces of test functions and the space in which the neural network is interpolated.

Numerical experiments confirm the theoretical prediction. Interestingly, in our experiments the error between the exact solution and the computed neural network decays asymptotically with the same rate as predicted by our theory for the interpolant of the network; however, this behaviour cannot be rigorously guaranteed, since in general the minimization problem which defines the computed neural network is underdetermined, and the computed neural network may be affected by spurious components. Indeed, we show that for a problem with zero data the minimization of the loss function may yield non-vanishing neural networks. With the method proposed in this paper, we combine the efficiency of the VPINN approach with the availability of a sound and certified convergence analysis.

The paper is organized as follows. In Sect. 2 we introduce the elliptic problem we are focusing on, and we also present the way in which the Dirichlet boundary conditions are exactly imposed, which is uncommon in PINNs and VPINNs but can be generalized as in [30]. In Sect. 3 we focus on the numerical discretization; in particular, the involved neural network architecture is described in Sect. 3.1, while the problem discretization and the corresponding loss function are described in Sect. 3.2. Here we also introduce an interpolation operator $\mathcal{I}_{H}$ applied to the neural networks. Sect. 4 is the key theoretical section: through a series of preliminary results, we formally derive the a priori error estimate, the main result being Theorem 4.9. In Sect. 5 we specify the parameters of the neural network used for the numerical tests and the training phase details. We also analyse the consequences of fulfilling the inf-sup condition in connection with the VPINN efficiency. Various numerical tests are presented and discussed in Sect. 6 for two-dimensional elliptic problems. Such tests empirically confirm the validity of the a priori estimate in different scenarios. Furthermore, we compare the accuracy of the proposed method with that of a standard PINN and the non-interpolated VPINN. We also analyse the relationship between the neural network hyperparameters and the VPINN accuracy, and we highlight, with numerical experiments and analytical examples, the importance of the inf-sup condition. In Sect. 7, we show that our VPINN can be easily adapted to solve a parametric nonlinear PDE, with accurate results for the whole range of parameters. Finally, in Sect. 8, we draw some conclusions and highlight the future perspective of the current work.

Let $\Omega\subset\mathbb{R}^{n}$ be a bounded polygonal/polyhedral domain with boundary $\Gamma=\partial\Omega$, partitioned into $\Gamma=\Gamma_{D}\cup\Gamma_{N}$ with $\Gamma_{D}\cap\Gamma_{N}=\emptyset$ and $\text{meas}_{n-1}(\Gamma_{D})>0$.

Let us consider the model elliptic boundary-value problem

$$\begin{cases}Lu:=-\nabla\cdot(\mu\nabla u)+\bm{\beta}\cdot\nabla u+\sigma u=f&\text{in \ }\Omega\,,\\ u=g&\text{on \ }\Gamma_{D}\,,\\ \mu\frac{\partial u}{\partial n}=\psi&\text{on \ }\Gamma_{N}\,,\end{cases}$$

(2.1)

where $\mu,\sigma\in{\rm L}^{\infty}(\Omega)$, $\bm{\beta}\in({\rm W}^{1,\infty}(\Omega))^{n}$ satisfy $\mu\geq\mu_{0}$, $\sigma-\frac{1}{2}\nabla\cdot\bm{\beta}\geq 0$ in $\Omega$ for some constant $\mu_{0}>0$, whereas $f\in L^{2}(\Omega)$, $g=\bar{u}_{|\Gamma_{D}}$ for some $\bar{u}\in{\rm H}^{1}(\Omega)$, and $\psi\in{\rm L}^{2}(\Gamma_{N})$.

Define the spaces $U={\rm H}^{1}(\Omega)$, $V={\rm H}^{1}_{0,\Gamma_{D}}(\Omega):=\{v\in U:v_{|\Gamma_{D}}=0\}$, the bilinear form $a:U\times V\to\mathbb{R}$ and the linear form $F:V\to\mathbb{R}$ such that

$$a(w,v)=\int_{\Omega}\mu\nabla w\cdot\nabla v+\bm{\beta}\cdot\nabla w\,v+\sigma w\,v\,,\qquad F(v)=\int_{\Omega}f\,v+\int_{\Gamma_{N}}\psi\,v\,;$$

(2.2)

denote by $\alpha\geq\mu_{0}$ the coercivity constant of the form $a$, and by $\|a\|$, $\|F\|$ the continuity constants of the forms $a$ and $F$. Problem (2.1) is formulated variationally as follows: Find $u\in\bar{u}+V$ such that

$$a(u,v)=F(v)\qquad\forall v\in V\,.$$

(2.3)

We assume that we can represent $u$ in the form

$$u=\bar{u}+\Phi\tilde{u}\,,$$

(2.4)

for some (known) smooth function $\Phi\in V$ and some $\tilde{u}\in U$ having the same smoothness of $u$. Let us introduce the affine mapping

$$B:U\to\bar{u}+V\qquad\text{ such that}\quad Bw=\bar{u}+\Phi w$$

(2.5)

which enforces the given Dirichlet boundary condition. Then, Problem (2.3) can be equivalently formulated as follows: Find $\tilde{u}\in U$ such that

$$a(B\tilde{u},v)=F(v)\qquad\forall v\in V\,.$$

(2.6)

In this section, we first introduce the class of neural networks used in this paper, then we describe the numerical discretization of the boundary-value problem (2.6), which uses neural networks to represent the discrete solution and piecewise polynomial functions to enforce the variational equations. An inf-sup stable Petrov-Galerkin formulation is introduced which guarantees stability and convergence, as indicated in Sect. 4; this is the main difference between the proposed method and other formulations, such as [14, 15].

Let $u^{\cal N\!N}\in U^{\cal N\!N}$ be any solution of the minimization problem (3.12); let us set

$$u_{H}^{{\cal N\!N}}={\cal I}_{H}Bu^{\cal N\!N}\in U_{H}\,.$$

(4.1)

Recalling the definition (2.4) of the affine mapping $B$, it holds

$$u_{H}^{{\cal N\!N}}=\bar{u}_{H}+u_{H}^{{\cal N\!N},0}\,,\qquad\text{with}\quad\bar{u}_{H}={\cal I}_{H}\bar{u}\quad\text{and}\quad u_{H,0}^{\cal N\!N}={\cal I}_{H}(\Phi u^{\cal N\!N})\in U_{H,0}\,;$$

(4.2)

note that $\bar{u}_{H}$ is a discrete lifting in $U_{H}$ of the Dirichlet data $g$.

We aim at estimating the error between $u$ and $u_{H}^{\cal N\!N}$. To accomplish this task, we need several definitions, assumptions, and technical results.

Next, we introduce the consistency errors due to numerical quadratures

$$E_{h}^{a}(w_{H},v_{h})=a(w_{H},v_{h})-a_{h}(w_{H},v_{h})\qquad\forall w_{H}\in U_{H},\ \forall v_{h}\in V_{h}\,,$$

(4.4)

$$E_{h}^{F}(v_{h})=F(v_{h})-F_{h}(v_{h})\qquad\forall v_{h}\in V_{h}\,,$$

(4.5)

and we provide a bound on these errors. To this end, let us assume that the quadrature rules used in the elements in ${\cal T}_{h}$ are obtained by affine transformations from a quadrature rule $\{(\hat{\xi}_{\iota},\hat{\omega}_{\iota}):\iota\in\hat{I}\}$ on a reference element $\hat{E}\subset\mathbb{R}^{n}$; similarly, let us assume that the quadrature rules used in the edges on $\partial{\cal T}_{h}(\Gamma_{N})$ are obtained by affine transformations from a quadrature rule $\{(\check{\xi}_{\iota},\check{\omega}_{\iota}):\iota\in\check{I}\}$ on a reference element $\check{e}\subset\mathbb{R}^{n-1}$.

Consequently, let us introduce the following notation

	$\displaystyle{\cal N}_{k}(\mu,\bm{\beta},\sigma)$	$\displaystyle=$	$\displaystyle\|\mu\|_{{\rm W}^{k,\infty}(\Omega)}+\|\bm{\beta}\|_{({\rm W}^{k,\infty}(\Omega))^{n}}+\|\sigma\|_{{\rm W}^{k,\infty}(\Omega)}\,,$		(4.8)
	$\displaystyle{\cal N}_{k}(f,\psi)$	$\displaystyle=$	$\displaystyle\|f\|_{{\rm W}^{k,\infty}(\Omega)}+\|\psi\|_{{\rm W}^{k,\infty}(\Gamma_{N})}\,,$		(4.9)
	$\displaystyle\|w_{H}\|_{k,{\cal T}_{H}}$	$\displaystyle=$	$\displaystyle\left(\sum_{G\in{\cal T}_{H}}\|w_{H|G}\|_{{\rm H}^{k}(G)}\right)^{1/2}\quad\forall w_{H}\in U_{H}\,.$		(4.10)

Finally, we pose a fundamental assumption.

This assumption together with Property 4.3 yields the following result.

We are ready to estimate the error $\|u-u_{H}^{\cal N\!N}\|_{1,\Omega}$. Recalling the decomposition (4.2), we use the triangle inequality

$$\|u-u_{H}^{\cal N\!N}\|_{1,\Omega}\leq\|u-u_{H}\|_{1,\Omega}+\|u_{H}-u_{H}^{\cal N\!N}\|_{1,\Omega}\,,$$

(4.15)

where $u_{H}$ is a suitable element in the affine subspace $\bar{u}_{H}+U_{H,0}\subset U_{H}$. Writing $u_{H}=\bar{u}_{H}+u_{H,0}$ with $u_{H,0}\in U_{H,0}$, one has $u_{H}-u_{H}^{\cal N\!N}=u_{H,0}-u_{H,0}^{\cal N\!N}\in U_{H,0}$; hence, we can apply (4.14) to get

$$\|u_{H}-u_{H}^{\cal N\!N}\|_{1,\Omega}\leq\frac{1}{\tilde{\alpha}_{\star}}\sup_{v_{h}\in V_{h}}\frac{a_{h}(u_{H}-u_{H}^{\cal N\!N},v_{h})}{\|v_{h}\|_{1,\Omega}}\,.$$

(4.16)

Recalling the definitions (4.4) and (4.5), it holds

	$\displaystyle a_{h}(u_{H},v_{h})$	$\displaystyle=$	$\displaystyle a(u_{H},v_{h})-E_{h}^{a}(u_{H},v_{h})$
		$\displaystyle=$	$\displaystyle a(u,v_{h})-a(u-u_{H},v_{h})-E_{h}^{a}(u_{H},v_{h})$
		$\displaystyle=$	$\displaystyle F(v_{h})-a(u-u_{H},v_{h})-E_{h}^{a}(u_{H},v_{h})$
		$\displaystyle=$	$\displaystyle F_{h}(v_{h})+E_{h}^{F}(v_{h})-a(u-u_{H},v_{h})-E_{h}^{a}(u_{H},v_{h})\,.$

Thus, the numerator in (4.16) is given by

$$a_{h}(u_{H}-u_{H}^{\cal N\!N},v_{h})=F_{h}(v_{h})-a_{h}(u_{H}^{\cal N\!N},v_{h})-a(u-u_{H},v_{h})-E_{h}^{a}(u_{H},v_{h})+E_{h}^{F}(v_{h})\,.$$

On the other hand, recalling (3.10) we have

$$F_{h}(v_{h})-a_{h}(u_{H}^{\cal N\!N},v_{h})=F_{h}(v_{h})-a_{h}({\cal I}_{H}Bu^{\cal N\!N},v_{h})=\sum_{i\in I_{h}}r_{h,i}(u^{\cal N\!N})\,v_{i}\,,$$

(4.17)

hence, by (3.11) and (4.3),

$$|F_{h}(v_{h})-a_{h}(u_{H}^{\cal N\!N},v_{h})|\leq R_{h}(u^{\cal N\!N})\|\bm{v}\|_{\gamma}\leq C_{h}R_{h}(u^{\cal N\!N})\|v_{h}\|_{1,\Omega}\,.$$

Using the bounds (4.11) and (4.12), we obtain the following inequality

$$\begin{split}\|u-u_{H}^{\cal N\!N}\|_{1,\Omega}&\lesssim\left(1+\frac{1}{\tilde{\alpha}_{\star}}\right)\Bigl{(}\ \inf_{u_{H}\in\bar{u}_{H}+U_{H,0}}\left(\|u-u_{H}\|_{1,\Omega}+h^{k}{\cal N}_{k}(\mu,\bm{\beta},\sigma)\|u_{H}\|_{k,{\cal T}_{H}}\right)\\ &\qquad\qquad\qquad\qquad+C_{h}R_{h}(u^{\cal N\!N})+h^{k}{\cal N}_{k}(f,\psi)\ \Bigr{)}\,.\end{split}$$

(4.18)

From now on, we assume that $u\in{\rm H}^{k+1}(\Omega)$. Then, assumption (3.8) yields the inequalities

$$\|u-{\cal I}_{H}u\|_{1,\Omega}\lesssim H^{k}|u|_{k+1,\Omega}\lesssim h^{k}|u|_{k+1,\Omega}$$

(4.19)

and

$$\|{\cal I}_{H}u\|_{k,{\cal T}_{H}}\leq\|u\|_{k,\Omega}+\|u-{\cal I}_{H}u\|_{k,{\cal T}_{H}}\lesssim\|u\|_{k,\Omega}+H|u|_{k+1,\Omega}\lesssim\|u\|_{k+1,\Omega}\,.$$

(4.20)

Choosing $u_{H}={\cal I}_{H}u\in\bar{u}_{H}+U_{H,0}$ in (4.18) and using these estimates, we arrive at the following intermediate result, which can be viewed as a mixed a priori/a posteriori error estimate.

Our next task will be bounding the term $R_{h}(u^{\cal N\!N})$. To this end, we use the minimality condition (3.12) to get

$$R_{h}(u^{\cal N\!N})\leq R_{h}(w^{\cal N\!N})\qquad\forall w^{\cal N\!N}\in U^{\cal N\!N}\,.$$

(4.21)

On the other hand, since $R_{h}(w^{\cal N\!N})$ is a weighted $\ell_{2}$-norm in $\mathbb{R}^{|I_{h}|}$, we can write

$$R_{h}(w^{\cal N\!N})=\sup_{{\mathbf{z}}\in\mathbb{R}^{|I_{h}|}}\frac{1}{\|{\mathbf{z}}\|_{\gamma}}\sum_{i\in I_{h}}r_{h,i}(w^{\cal N\!N})z_{i}\,,$$

where, similarly to (4.17),

$$\sum_{i\in I_{h}}r_{h,i}(w^{\cal N\!N})z_{i}=F_{h}(z_{h})-a_{h}({\cal I}_{H}Bw^{\cal N\!N},z_{h})\qquad\text{with \ }z_{h}=\sum_{i\in I_{h}}z_{i}\varphi_{i}\in V_{h}\,.$$

For convenience, in analogy with (4.1), let us set

$$w_{H}^{{\cal N\!N}}={\cal I}_{H}Bw^{\cal N\!N}\in U_{H}\,.$$

(4.22)

Thus, recalling (4.3), we obtain

$$R_{h}(w^{\cal N\!N})\leq\frac{1}{c_{h}}\sup_{z_{h}\in V_{h}}\frac{F_{h}(z_{h})-a_{h}(w_{H}^{\cal N\!N},z_{h})}{\|z_{h}\|_{1,\Omega}}\qquad\forall w^{\cal N\!N}\in U^{\cal N\!N}\,.$$

(4.23)

The numerator can be manipulated as above, using

$$F_{h}(z_{h})=F(z_{h})-E_{h}^{F}(z_{h})=a(u,z_{h})-E_{h}^{F}(z_{h})$$

and

$$a_{h}(w_{H}^{\cal N\!N},z_{h})=a(w_{H}^{\cal N\!N},z_{h})-E_{h}^{a}(w_{H}^{\cal N\!N},z_{h})\,,$$

whence, using once more Property 4.3, we get

$$R_{h}(w^{\cal N\!N})\lesssim\frac{1}{c_{h}}\left(\|u-w_{H}^{\cal N\!N}\|_{1,\Omega}+h^{k}{\cal N}_{k}(\mu,\bm{\beta},\sigma)\|w_{H}^{\cal N\!N}\|_{k,{\cal T}_{H}}+h^{k}{\cal N}_{k}(f,\psi)\right)\,.$$

(4.24)

In order to bound the terms containing $w_{H}^{\cal N\!N}$, we introduce the quantity

$$e^{{\cal N\!N}}=u-Bw^{\cal N\!N}\,,$$

(4.25)

which, recalling the definitions (2.4) and (2.5), can be written as

$$e^{{\cal N\!N}}=\Phi(\tilde{u}-w^{\cal N\!N})\,,$$

(4.26)

and we formulate a final assumption.

Note that (4.27) implies in particular $u\in{\rm H}^{k+1}(\Omega)$ with $\|u\|_{k+1,\Omega}\lesssim\|\tilde{u}\|_{k+1,\Omega}\,\|\Phi\|_{k+1,\infty,\Omega}$; on the other hand, (4.28) implies $e^{\cal N\!N}\in{\rm H}^{2}(\Omega)$. (We refer to Remark 4.11 for another set of assumptions on the neural network.)

Recalling (4.22) and using the identity

$$u-w_{H}^{\cal N\!N}=(u-{\cal I}_{H}u)+{\cal I}_{H}e^{\cal N\!N}=(u-{\cal I}_{H}u)-e^{\cal N\!N}+(I-{\cal I}_{H})e^{\cal N\!N}\,,$$

(4.29)

we can write

$$\|u-w_{H}^{\cal N\!N}\|_{1,\Omega}\lesssim\|u-{\cal I}_{H}u\|_{1,\Omega}+\|e^{\cal N\!N}\|_{1,\Omega}+H|e^{\cal N\!N}|_{2,\Omega}$$

(4.30)

and, using a standard inverse inequality in $\mathbb{P}_{k}(G)$ for any $G\in{\cal T}_{H}$,

$$\begin{split}\|w_{H}^{\cal N\!N}\|_{k,{\cal T}_{H}}&\lesssim\|{\cal I}_{H}u\|_{k,{\cal T}_{H}}+H^{1-k}\|{\cal I}_{H}e^{\cal N\!N}\|_{1,\Omega}\\ &\lesssim\|{\cal I}_{H}u\|_{k,{\cal T}_{H}}+H^{1-k}\left(\|e^{\cal N\!N}\|_{1,\Omega}+H|e^{\cal N\!N}|_{2,\Omega}\right)\,.\end{split}$$

(4.31)

Keeping into account (4.19) and (4.20), in order to conclude we need to identify a function $\tilde{w}^{\cal N\!N}\in U^{\cal N\!N}$ for which a bound of the type

$$|e^{\cal N\!N}|_{m,\Omega}\lesssim|\tilde{u}-\tilde{w}^{\cal N\!N}|_{m,\Omega}\lesssim H^{k+1-m}|\tilde{u}|_{k+1,\Omega}$$

(4.32)

holds true for $m=1,2$. The existence of such a function is guaranteed by one of the available results on the approximation of functions in Sobolev spaces by neural networks (see [7, Theorem 5.1, Remark 5.2]; see also [24]), provided the number of layers $L$ and the widths of the layers in the chosen ${\cal N\!N}$ satisfy suitable conditions depending on the target accuracy (hence, in our case depending on $H^{k}$). Indeed, suppose one is interested in using meshes with meshsize as small as $H_{\min}$ in the domain $\Omega$ (here assumed to satisfy $\Omega\subset[0,1]^{n}$ for the sake of simplicity), and let $N\in\mathbb{N}$ be such that

$$3^{n}(1+\delta)(2(m+1))^{3m}\max\{R^{m},\text{ln}^{m}(\beta N^{k+n+3})\}\dfrac{C(n,m,k,\tilde{u})}{N^{k+1-m}}\leq H_{\min}^{k+1-m}|\tilde{u}|_{k+1,\Omega}\,,$$

(4.33)

where $\delta$, $R$, $\beta$ and $C$ are constants not depending on $N$ defined in [7]. Then, a function $\tilde{w}^{\cal N\!N}$ exists which fulfils (4.32) and is represented as a neural network with the hyperbolic tangent as activation function and two hidden layers with $N_{1}$ and $N_{2}$ neurons respectively, satisfying

$$N_{1}\leq 3\left\lceil\frac{k+1}{2}\right\rceil\left|P_{k,n+1}\right|+n(N-1),\hskip 8.5359ptN_{2}\leq 3\left\lceil\frac{n+2}{2}\right\rceil\left|P_{n+1,n+1}\right|N^{n},$$

(4.34)

where

$$\left|P_{a,b}\right|=\left(a+b-1\atop a\right),\hskip 28.45274pt\forall a,b\in\mathbb{N},\ b\geq 2.$$