Robust Variational Physics-Informed Neural Networks
We introduce a Robust version of the Variational Physics-Informed Neural
Networks method (RVPINNs). As in VPINNs, we define the quadratic loss
functional in terms of a Petrov-Galerkin-type variational formulation of the
PDE problem: the trial space is a (Deep) Neural Network (DNN) manifold, while
the test space is a finite-dimensional vector space. Whereas the VPINN's loss
depends upon the selected basis functions of a given test space, herein, we
minimize a loss based on the discrete dual norm of the residual. The main
advantage of such a loss definition is that it provides a reliable and
efficient estimator of the true error in the energy norm under the assumption
of the existence of a local Fortin operator. We test the performance and
robustness of our algorithm in several advection-diffusion problems. These
numerical results perfectly align with our theoretical findings, showing that
our estimates are sharp.
Sergio Rojas, Paweł Maczuga, Judit Muñoz-Matute, David Pardo, Maciej Paszynski
# Robust Variational Physics-Informed Neural Networks

###### Abstract

We introduce a Robust version of the Variational Physics-Informed Neural Networks (RVPINNs) to approximate the Partial Differential Equations (PDEs) solution. We start from a weak Petrov-Galerkin formulation of the problem, select a discrete test space, and define a quadratic loss functional as in VPINNs. Whereas in VPINNs the loss depends upon the selected basis functions of a given test space, herein we minimize a loss based on the residual in the discrete dual norm, which is independent of the test space's choice of test basis functions. We demonstrate that this loss is a reliable and efficient estimator of the true error in the energy norm. The proposed loss function requires computation of the Gram matrix inverse, similar to what occurs in traditional residual minimization methods. To validate our theoretical findings, we test the performance and robustness of our algorithm in several advection-dominated-diffusion problems in one spatial dimension. We conclude that RVPINNs is a robust method.

_Keywords_: Robustness, Variational Physics-Informed Neural Networks, Petrov-Galerkin formulation, Riesz representation, Minimum Residual principle, a posteriori error estimation.

###### Contents

* 1 Introduction
* 2 Preliminaries
	* 2.1 Abstract framework
	* 2.2 Neural Network framework
	* 2.3 Variational Physics-Informed Neural Networks
	* 2.4 An alternative definition of VPINNs
* 3 Robust Variational Physics-Informed Neural Networks
	* 3.1 Orthonormal discrete basis and relation with other VPINNs
* 4 Error estimates for RVPINNs
	* 4.1 A posteriori error estimates for RVPINNs in the sense of equivalence classes
	* 4.2 Energy norm error estimates based on local semi-discrete inf-sup condition
* 5 Numerical examples
	* 5.1 Diffusion-advection model problem
	* 5.2 Discrete setting
	* 5.3 A smooth diffusion problem
* 6
5.4 Delta source problem * 5.5 Advection-dominated diffusion problem
* 6 Conclusions
* 7 Acknowledgments

## 1 Introduction

The remarkable success of Deep Learning (DL) algorithms across different scientific areas [22, 30, 21] in the last decade has recently led to explore the potential of this discipline to address classical problems in physics and mathematics. These problems include approximating the solutions of Partial Differential Equations (PDEs) employing (deep) Neural Networks (NN), which can accurately approximate continuous functions. The exponential growth of interest in these techniques started with the Physics Informed Neural Networks (PINNs) ([40]). This method incorporates the governing physical laws the PDE describes in the learning process. Then, the network is trained on a dataset that consists of a random selection of points in the physical domain and its boundary. PINNs have been successfully applied to solve a wide range of problems in scientific computing, including fluid mechanics [5, 35], wave propagation [41], or inverse problems [10, 36], among many others. However, the loss function in PINNs is given by the strong form of the PDE residual, and it is known that in some problems with low regularity of the data, the solution only makes sense in a variational form. Therefore, PINNs fails to provide accurate solutions in those cases [29].

A natural continuation to PINNs addressing the aforementioned limitation is the so-called Variational PINNs (VPNNs) [25]. Here, the authors introduce a variational formulation of the underlying PDE within a Petrov-Galerkin framework. The solution is then approximated by a (deep) NN, whereas the test functions belong to linear vector spaces. Finally, the authors define a variational loss function to minimize during the training process. However, this methodology was not as popular as its predecessor. The main reason might be that the approach is sensible to the choice of the basis test functions. That is, given a discrete test space, the method's loss, stability, and robustness heavily depends upon the choice of the basis functions of the given test space. In [2], authors present an a posteriori error analysis for discretizing elliptic boundary-value problems with VPINNs employing piecewise polynomials for the test space. As the loss function in VPINNs is, in general, not robust with respect to the true error (for example, the loss function can tend to zero even if the true error does not), they provide an error estimator employing classical techniques in finite element analysis to obtain practical information on the quality of the approximation. They also study in [3] the importance of selecting appropriate quadrature rules and the degree of piecewise polynomial test functions to guarantee optimal convergence rates for VPINNs when performing mesh refinements. Finally, in [26], the same authors from VPINNs introduced _hp_-VPINNs. In the latter, they employ piecewise polynomials for testing and allow for _hp_-refinements via domain decomposition by selecting non-overlapping test functions. However, the aforementioned issue of selecting an appropriate loss function for each problem remains also in this approach.

The strategy we propose in this article overcomes the aforementioned limitation by following the core ideas introduced in Minimum Residual (MinRes) methods. The latter is a class of numerical methods for solving PDEs that guarantee stability. The spirit of MinRes methods is to minimize the dual norm of the residual of the PDE given in variational form. Many strategies have been developed based on this idea over the last five decades, including the families of Galerkin Least-Squares methods [4, 24, 23], First Order Least-squares [6, 7], residual minimization methods on dual discontinuous Galerkin norms [8, 12, 11, 43, 31], isogeometric residual minimization methods [9, 33, 32, 34], and Discontinuous Petrov-Galerkin (DPG) methods [14, 15, 16, 37, 42]. As the residual operator lives in the dual of the test space and the dual norm is difficult to compute, the use of the Riesz representation theorem [38] is natural in this context. The latter maps elements of the dual space into elements on the test space. Therefore, in many of the aforementioned methods, instead of minimizing the residual in the dual norm, they minimize its Riesz representative (which is an element of the test space) in a given test norm.

In this work, we revisit the initial work on VPINNs from Kharazmi et al. in [25], and we introduce a Robust version of Variational Physics-Informed Neural Networks (RVPINNs). As in VPINNs, we consider a Petrov-Galerkin formulation of the PDE, approximate the solution employing NNs, and select a finite-dimensional test space. Our goal is to minimize the residual in the discrete dual norm. Therefore, we select a single test function that is the Riesz representation of the weak residual over the selected discrete test space. We select an inner product and span the Riesz representative on a discrete basis so the resulting loss function will include the inverse of the Gram matrix corresponding to the selected inner product. We prove that the norm of such test function is efficient and reliable, i.e., the norm of the true error is bounded from below and above (up to an oscillation
term) by the norm of the residual representative. Moreover, as we show in the numerical results, our strategy is insensitive (unlike VPINNs) to the choice of the discrete basis spanning the test space (assuming that the numerical integration and the inversion of the Gram matrix are sufficiently accurate). It is easy to see that minimizing the proposed loss functional is equivalent to minimizing the test norm of the Riesz representation of the residual, which is the idea behind classical MinRes methods. Summarizing, we provide a general mathematical framework to define robust loss functionals in VPINNs. Our strategy relies on two ideas: (a) the appropriate selection of the inner product in the test space, ensuring the stability of the variational formulation, and (b) the selection of a single test function, that is the Riesz representation of the weak residual.

In particular, if we select an orthonormal discrete basis with respect to the inner product in the test space, the Gram matrix becomes the identity, and we recover the original definition of the loss functional in VPINNs. The difference with classical VPINNs is that our strategy is robust and independent of the choice of the basis functions. Other works employ similar ideas closely related to minimum residual methods. In [45], the authors minimize the dual norm of the weak residual for symmetric variational formulations in \(H^{1}\) and rectangular domains via a spectral decomposition. In [44], the same authors extend this method to time-harmonic Maxwell equations in the context of symmetric \(H(curl)\)-formulations. Also, in [46], in the context of Petrov-Galerkin formulations, the authors minimize the dual norm of the residual based on the concept of optimal testing form [16]. They approximate both the solution of the variational problem and the corresponding optimal test functions employing NNs by solving two nested deep Ritz problems.

Finally, apart from the usual limitations of DL-based technologies (optimizer, integration, etc.), the largest bottleneck of RVPINNs at the moment is that in certain configurations, we need to invert the Gram matrix corresponding to the selected basis and inner product on the test space. We will explore how to optimize this step in the future. However, there are particular cases where the inversion of the Gram matrix is trivial (as in the spectral case) or easy to compute. For example, if we consider a strong variational formulation, the test inner product is \(L^{2}\), and selecting piecewise discontinuous polynomials for testing the Gram matrix becomes block diagonal. On the other hand, for parametric problems, the Gram matrix may remain the same for different values of the parameters. Therefore, the inversion can be done offline, being valid for a large class of problems. Parametric problems are essential to solve inverse problems [20], so RVPINNs could be of great interest in this area.

The article is organized as follows: In Section 2, we introduce the variational formulation of the model problem we consider in this article, the NN framework for approximation of the solution, and a brief overview of VPINNs. We also provide an alternative definition of VPINNs employing a single test function in the loss. Section 3 presents the methodology of RVPINNs and the connection with other methods based on VPINNs. Section 4 is devoted to the derivation of robust error estimates. We first prove that the residual representative is a reliable and efficient a posteriori error estimator of the true error in the sense of equivalent classes. Then, we demonstrate under the assumption of a local Fortin operator's existence that the true error is equivalent to the residual error estimator up to an oscillation term. In Section 5 we test our method in several 1D advection-diffusion problems, showing the robustness of the approach. Finally, Section 6 summarizes the conclusions and future research lines.

## 2 Preliminaries

### Abstract framework

Let \(U\) and \(V\) denote Hilbert spaces with norms \(\|\cdot\|_{U}\) and \(\|\cdot\|_{V}\). Assume that we are interested in obtaining an approximation of a PDE problem admitting a variational formulation of the form:

\[\text{Find}\;u\in U,\;\text{such that:}\;\;r(u,v):=l(v)-a(u,v)=0,\,\forall\,v \in V, \tag{1}\]

where \(l(\cdot)\in V^{\prime}\) is a bounded linear form, with \(V^{\prime}\) denoting the dual space of \(V\); \(a(\cdot,\cdot)\) is a bounded inf-sup stable bilinear form in \(U\times V\). That is, there exist constants \(\mu,\alpha>0\), respectively, such that:

\[a(w,v)\leq\mu\|w\|_{U}\|v\|_{V},\quad\forall\,w\in U,v\in V, \tag{2}\]

and

\[\sup_{0\neq v\in V}\frac{a(w,v)}{\|v\|_{V}}\geq\alpha\|w\|_{U},\quad\forall\,w \in U. \tag{3}\]

We also assume that, for all \(v\in V\), the operator \(a(\cdot,v)\in U^{\prime}\) satisfies:

\[a(w,v)=0,\,\forall\,w\in U\Longrightarrow v=0. \tag{4}\]
Problem (1) admits a unique solution by means of the Banach-Necas-Baboska Theorem, and the following a priori estimate holds (see [17, Theorem 1.1]):

\[\|u\|_{U}\leq\frac{1}{\alpha}\|I(\cdot)\|_{V^{\prime}}, \tag{5}\]

with

\[\|I(\cdot)\|_{V^{\prime}}:=\sup_{0\neq v\in V}\frac{I(v)}{\|v\|_{V}}. \tag{6}\]

Additionally, it is well-known that (3) and (4) imply that the following adjoint inf-sup condition holds (see, e.g., [13, Theorem 1]:

\[\sup_{0\neq w\in U}\frac{a(w,v)}{\|w\|_{U}}\geq\alpha\|v\|_{V},\quad\forall\,v\in V. \tag{7}\]

**Remark 1** (The weak residual).: _For a given \(w\in U\), we refer to \(r(w,\cdot)\in V^{\prime}\) (see (1)) as the weak residual._

### Neural Network framework

To numerically approximate (1), we consider a (Deep) Neural Network (DNN) function with input \(\mathbf{x}=(x_{1},\ldots,x_{d})\) and output \(u_{\theta}(\mathbf{x})\), where \(\theta\in\mathbb{R}^{S}\) represents the trainable parameters. For simplicity, in this work, we employ a simple fully-connected feedforward Neural Network structure composed of \(L\) layers. Each layer \(l\) in the Neural Network consists of a set of neurons, which compute a weighted sum of their inputs plus a bias followed by a nonlinear activation function in the first \(L-1\) layers and the identity function as the activation in the last layer. More precisely, the output of layer \(l\), with \(l=1,\ldots L-1\), is given by:

\[\mathbf{z}^{(l)}=\sigma(\mathbf{w}^{(l)}\mathbf{z}^{(l-1)}+\mathbf{b}^{(l)}), \tag{8}\]

where \(\sigma\) is a nonlinear activation function (e.g., the tanh activation function), \(\mathbf{w}^{(l)}\), \(\mathbf{b}^{(l)}\) are the weights and biases respectively associated with the layer \(l\), and \(\mathbf{z}^{(0)}=\mathbf{x}\) is the input to the first layer. The final layer \(L\) is given by:

\[u_{\theta}=\mathbf{w}^{(L)}\mathbf{z}^{(L-1)}+\mathbf{b}^{(L)}. \tag{9}\]

Using an optimization algorithm, the Neural Network weights and biases are learned from a training set by minimizing a loss functional of input-output pairs.

We denote by \(U_{NN}\) the manifold consisting of all possible realizations for a given DNN architecture belonging to \(U\), where \(U\) is defined in Section 2.1. That is,

\[U_{NN}:=\{u_{\theta},\,\forall\,\theta\in\mathbb{R}^{S}\}\cap U. \tag{10}\]

Finally, for a given discrete space \(V_{M}\subseteq V\), we consider the following Petrov-Galerkin discretization of (1):

\[\text{Find }u_{NN}\in U_{NN},\text{ such that: }r(u_{NN},v_{M})=0,\,\forall\,v_{M} \in V_{M}. \tag{11}\]

### Variational Physics-Informed Neural Networks

Variational Physics Informed Neural Networks (VPINNs) approximate the solution \(u\) of problem (1) by minimizing a loss functional defined in terms of the weak residual (see Remark 1), and a discrete space \(V_{M}\subseteq V\) with basis \(\{\varphi_{m}\}_{m=1}^{M}\). That is, it approximates \(u\) by solving the following minimization problem:

\[\text{Find }u_{\theta^{*}},\text{ such that }\theta^{*}=\underset{\theta\in \mathbb{R}^{S}}{\text{argmin}}\,\mathcal{L}_{r}\left(u_{\theta}\right), \tag{12}\]

with \(\mathcal{L}_{r}(u_{\theta})\) being a loss functional defined in terms of the weak residual and the discrete basis. This approach allows to consider scenarios where the solution belongs to less regular spaces than in collocation PINNs. For example, elliptic problems with solutions in \(H^{2-\varepsilon}(\Omega)\), with \(0<\varepsilon\leq 1\).

A typical definition for the loss functional in (12) found in the literature is the following (see, e.g., [25, 27, 2, 3]):

\[\mathcal{L}_{r}(u_{\theta}):=\sum_{m=1}^{M}r(u_{\theta},\varphi_{m})^{2}+C(u_ {\theta}), \tag{13}\]

where \(C(\cdot)\) is a \(V_{M}\)-independent quadratic functional employed to impose boundary conditions for \(u_{NN}\), that can be taken as \(C(u_{\theta})=0\) if boundary conditions are strongly imposed in the DNN structure.
**Remark 2** (Classical approach).: _We will refer to problem (12), with the loss functional defined by (13), as the classical VPINNs approach._

If we assume sufficient approximability of the DNN structure and an ideal optimizer (solver), a minimizer for (13) for which the loss functional vanishes is also a solution to the Petrov-Galerkin problem (11). However, adopting such a loss functional may be computationally inefficient. Indeed, given the test space \(V_{M}\), Eq. (13) depends heavily upon the choice of its basis. In particular, a simple re-scaling of one basis function (e.g., by making it sufficiently large) may easily lead to catastrophic results (cf. [45]). Consequently, there is a natural demand for robust loss definitions in practical VPINNs.

### An alternative definition of VPINNs

To motivate how robust losses can be constructed, we first notice that, for a given set of trainable parameters \(\theta\), we can define the following function in \(V_{M}\):

\[\widetilde{\varphi}=\sum_{m=1}^{M}r(u_{\theta}\,,\varphi_{m})\,\varphi_{m}. \tag{14}\]

Then, as a consequence of the linearity of the weak residual, it holds (cf. [25]):

\[\mathcal{L}_{r}\,(u_{\theta})=\sum_{m=1}^{M}r(u_{\theta}\,,\varphi_{m})^{2}+C (u_{\theta})=r(u_{\theta}\,,\widetilde{\varphi})+C(u_{\theta}). \tag{15}\]

The last equation suggests that VPINNs can be generalized in the following way. For a given \(\widetilde{\varphi}\in V_{M}\), a VPINNs loss functional can be defined as:

\[\mathcal{L}_{r}\,(u_{\theta})=r(u_{\theta}\,,\widetilde{\varphi})+C(u_{ \theta}). \tag{16}\]

The goal now is to define the loss functional as (16) in terms of a particular discrete test function \(\widetilde{\varphi}\in V_{M}\) such that VPINNs becomes robust. In the following Section, we present a strategy based on defining such a discrete function as the Riesz residual representative in \(V_{M}\) with respect to the norm inducing the inf-sup stability (3).

## 3 Robust Variational Physics-Informed Neural Networks

In Robust Variational Physics-Informed Neural Networks (RVPINNs), we first introduce an inner product \((\cdot,\cdot)_{V}\)1, whose associated norm satisfies the inf-sup stability condition (3). Next, for a given trainable parameter \(\theta\), we compute \(\phi:=\phi(\theta)\in V_{M}\) being the solution of the following Galerkin problem:

Footnote 1: i.e., \(\|v\|_{V}^{2}:=(v,v)v\), for all \(v\in V\).

\[(\phi,\varphi_{n})_{V}=r(u_{\theta},\varphi_{n}),\text{ for }n=1,\ldots,M, \tag{17}\]

and define the loss functional as:

\[\mathcal{L}_{r}^{\phi}\,(u_{\theta}):=r(u_{\theta},\phi)+C(u_{\theta}), \tag{18}\]

where \(C(\cdot)\) is a \(V_{M}\)-independent quadratic functional to impose boundary conditions as in (13). Finally, we obtain the approximation of the solution of problem (1) by solving the following minimization problem (cf., [25]):

\[\text{Find }u_{\theta^{*}},\text{ such that }\theta^{*}=\underset{\theta\in \mathbb{R}^{S}}{\text{argmin}}\,\mathcal{L}_{r}^{\phi}\,(u_{\theta})\,, \tag{19}\]

**Remark 3** (Riesz representative of the weak residual).: _The solution \(\phi\) of problem (17) is the Riesz representative in \(V_{M}\) of the residual._

As a consequence of the linearity of the weak residual and the linearity of the inner product with respect to the second variable, defining:

\[\phi:=\sum_{m=1}^{M}\eta_{m}(\theta)\varphi_{m}, \tag{20}\]

problem (17) leads to the resolution of the following problem written in matrix form:

\[G\eta(\theta)=\mathcal{R}(\theta), \tag{21}\]
with \(\eta(\theta)\) the vector of coefficients \(\eta_{m}(\theta)\), \(G\) the (\(\theta\)-independent) symmetric and positive definite Gram matrix of coefficients \(G_{nm}=(\varphi_{m},\varphi_{n})_{V}\), and \(\mathcal{R}(\theta)\) the vector of coefficients \(\mathcal{R}_{n}(\theta)=r(u_{\theta},\varphi_{n})\), with \(n=1,\ldots,M\). Thus,

\[r(u_{\theta}\,,\phi)=\sum_{m=1}^{M}\eta_{m}(\theta)r(u_{\theta}\,,\phi_{m})= \sum_{m=1}^{M}\mathcal{R}_{m}(\theta)\eta_{m}(\theta)=\mathcal{R}(\theta)^{T }G^{-1}\mathcal{R}(\theta), \tag{22}\]

implying that the loss functional (18) is equivalently written as:

\[\mathcal{L}_{r}^{\phi}(u_{\theta})=\mathcal{R}(\theta)^{T}G^{-1}\mathcal{R}( \theta)+C(u_{\theta}). \tag{23}\]

We also notice that, by definition, to minimize the loss functional (23) is equivalent, up to the constraint \(C(u_{\theta})\), to minimize the quantity \(\|\phi\|_{V}^{2}\). Indeed, from (17), we deduce:

\[\|\phi\|_{V}^{2}=(\phi,\phi)_{V}=r(u_{\theta},\phi). \tag{24}\]

In Section 4.2, we prove that \(\|\phi\|_{V}\) is, up to an oscillation term, a local robust error estimation for the error \(\|u-u_{\theta}\|_{U}\) under the assumption of the existence of a local Fortin's operator.

### Orthonormal discrete basis and relation with other VPINNs

When a test space \(V_{M}\) is the span of an orthonormal set \(\{\varphi_{m}\}_{m=1}^{M}\) with respect to the \(\|\cdot\|_{V}\)-norm, the Gram matrix \(G\) coincides with the identity matrix; therefore, the corresponding residual representative has the form:

\[\phi=\sum_{m=1}^{M}r(u_{\theta},\varphi_{m})\varphi_{m}, \tag{25}\]

and the loss functional is explicitly written as

\[\mathcal{L}_{r}^{\phi}(u_{\theta})=\sum_{m=1}^{M}r(u_{\theta},\varphi_{m})^{2 }+C(u_{\theta}). \tag{26}\]

Thus, the classical VPINNs loss functional definition (13) is recovered for this case. This is also the case of the recently proposed Deep Fourier Residual method (see [45, 44]), where authors first define the loss functional as the continuous dual norm of the weak residual, that later is approximated by considering a truncation of a series in terms of an orthonormal basis for the test space. In particular, in [45], authors consider diffusion problems and approximate the norm of the weak residual in \(H^{-1}\). This turns out to be equivalent to (26) considering the standard \(H^{1}\)-norm for the space \(V\), and the functions \(\varphi_{m}\) as orthonormalized sinusoidal (spectral) functions with respect to the \(H^{1}\)-norm.

## 4 Error estimates for RVPINNs

One of the main complexities for proving the robustness of the residual estimator \(\|\phi\|_{V}\) is that the solution of the Petrov-Galerkin problem (11) may not have a solution or, if there exists, it may be non-unique since the space of all possible realizations for the NN structure defines a manifold instead of a finite-dimensional space (see, e.g., [3, Section 6.3]). Thus, standard FEM arguments based on a discrete inf-sup condition cannot be applied in this context.

Nevertheless, we can derive a posteriori error estimates using a different strategy. First, we introduce in Section 4.1 an equivalence class that allows us to neglect the part of the error that is \(a\)-orthogonal to \(V_{M}\). For that equivalence class, we prove that the residual representative always defines a reliable and efficient a posteriori estimator for the error. Then, in Section 4.2, we consider the case of the full error, for which we demonstrate its equivalence to the residual error estimator up to an oscillation term and under the assumption of the existence of a local Fortin operator.

### A posteriori error estimates for RVPINNs in the sense of equivalence classes

We start by introducing the following Null space of the operator \(A:U\mapsto V_{M}^{\prime}\), defined in terms of the bilinear form \(a(\cdot\,,\cdot)\) associated with problem (1) (cf. [13]):

\[U_{M}^{0}:=\left\{w\in U\,:\,\langle A(w)\,,v_{M}\rangle:=a(w\,,\,v_{M})=0,\, \forall\,v_{M}\in V_{M}\right\}, \tag{27}\]
and the following norm for the quotient space \(U/U_{M}^{0}\):

\[\|[w]\|_{U/U_{M}^{0}}:=\inf_{w_{0}\in U_{M}^{0}}\|w+w_{0}\|_{U}. \tag{28}\]

We extend the definition of the bilinear form \(a(\cdot\,,\cdot)\) to the product space \(U/U_{M}^{0}\times V\) as:

\[a([w],v)=a(w,v),\quad\text{with }w\in[w]\text{ being any arbitrary representative of the equivalence class.} \tag{29}\]

The following result follows as a consequence of the boundedness and inf-sup stability properties of the bilinear form \(a(\cdot\,,\cdot)\) (see equations (2) and (3)), respectively:

**Proposition 1**.: _The following boundedness and semi-discrete inf-sup conditions, respectively, are satisfied by the bilinear form (29):_

\[a([w],v_{M})\leq\mu\|[w]\|_{U/U_{M}^{0}}\|v_{M}\|_{V},\quad\forall[w]\in U/U_{ M}^{0},\,v_{M}\in V_{M}, \tag{30}\]

\[\sup_{0\neq v_{M}\in V_{M}}\frac{a([w],v_{M})}{\|v_{M}\|_{V}}\geq\alpha\|[w] \|_{U/U_{M}^{0}},\quad\forall[w]\in U/U_{M}^{0}. \tag{31}\]

Proof.: First, notice that (30) is direct from the definition (28) and the boundedness property of the bilinear form \(a(\cdot\,,\cdot)\). Indeed, for all \(w\in U\) and \(w_{0}\in U_{M}^{0}\), it holds:

\[a(w,v_{M})=a(w+w_{0}\,,v_{M})\leq\mu\|w+w_{0}\|_{U}\|v_{M}\|_{V},\qquad\forall v _{M}\in V_{M}.\]

Thus, (30) follows from taking the infimum. To prove (31), first notice that, for all \(v_{M}\in V_{M}\), (7) implies:

\[\alpha\|v_{M}\|_{V}\leq\sup_{0\neq w\in U}\frac{a(w\,,v_{M})}{\|w\|_{U}}=\sup_ {[0]\neq[w]\in U/U_{M}^{0}}\frac{a([w]\,,v_{M})}{\|w\|_{U}}\leq\sup_{[0]\neq[ w]\in U/U_{M}^{0}}\frac{a([w]\,,v_{M})}{\|[w]\|_{U/U_{M}^{0}}},\]

where the last inequality follows from the fact that \(\|[w]\|_{U/U_{M}^{0}}\leq\|w\|_{U}\), since \(w_{0}=0\in U_{M}^{0}\). Additionally, by construction, \(\forall[w]\in U/U_{M}^{0}\),

\[a([w]\,,v_{M})=0,\,\forall\,v_{M}\in V_{M}\Longrightarrow[w]=[0]. \tag{32}\]

Therefore, (31) follows as a consequence of Theorem 1 in [13]. 

Using the previous proposition, we can establish the main result of this section:

**Theorem 2** (Error class bounds in terms of the residual representative).: _Let \(u\in U\) be the solution of the continuous problem (1); \(u_{\theta}\in U_{NN}\) denote a DNN structure with trainable parameters \(\theta\in\mathbb{R}^{S}\); \(V_{M}\subseteq V\) denote a finite dimensional space, dotted with norm \(\|\cdot\|_{V}\); and \(\phi\in V_{M}\) be the solution of problem (17). It holds:_

\[\frac{1}{\mu}\|\phi\|_{V}\leq\|[u-u_{\theta}]\|_{U/U_{M}^{0}}\leq\frac{1}{ \alpha}\|\phi\|_{V}. \tag{33}\]

Proof.: First notice that, by definition of \(\phi\) and consistency of the analytical solution \(u\), it holds:

\[\|\phi\|_{V}=\sup_{0\neq v_{M}\in V_{M}}\frac{(\phi,v_{M})_{V}}{\|v_{M}\|_{V} }=\sup_{0\neq v_{M}\in V_{M}}\frac{r(u_{\theta},v_{M})}{\|v_{M}\|_{V}}=\sup_{ 0\neq v_{M}\in V_{M}}\frac{a([u-u_{\theta}],v_{M})}{\|v_{M}\|_{V}}.\]

Therefore, the left inequality in (33) is obtained by using the boundedness property (30). On another side, as a consequence of the inf-sup stability (31), it holds:

\[\alpha\left\|[u-u_{\theta}]\right\|_{U/U_{M}^{0}}\leq\sup_{0\neq v _{M}\in V_{M}}\frac{a([u-u_{\theta}],v_{M})}{\|v_{M}\|_{V}} =\sup_{0\neq v_{M}\in V_{M}}\frac{a(u-u_{\theta},v_{M})}{\|v_{M}\| _{V}}\] \[=\sup_{0\neq v_{M}\in V_{M}}\frac{r(u_{\theta},v_{M})}{\|v_{M}\|_{V}}\] \[=\sup_{0\neq v_{M}\in V_{M}}\frac{(\phi,v_{M})_{V}}{\|v_{M}\|_{V}}\] \[=\|\phi\|_{V},\]

proving the right inequality in (33).
**Remark 4** (Robustness of the residual representative in the sense of classes).: _Inequalities (33) show that, for any \(\theta\), \(\|\phi\|_{V}\) always, up to some multiplicative constants, defines an efficient (lower) and a reliable (upper) bound (thus, robust) estimation for the error \(\|[u-u_{\theta}]\|_{U/U_{M}^{0}}\). This, in practice, implies that to minimize \(\|\phi\|_{V}\) is equivalent to minimize \(\|[u-u_{\theta}]\|_{U/U_{M}^{0}}\). Moreover, if we assume that problem (11) admits a solution, when \(\|\phi\|_{V}\to 0^{+}\), \(u-u_{\theta}\) converges to a function that belongs to \(U_{M}^{0}\) as a consequence of the robustness. Thus, since \(u\not\in U_{M}^{0}\), we conclude that \(u_{\theta}\) converges to a function \(u_{\theta}\), satisfying:_

\[0=a(u-u_{\theta^{*}},v_{M})=r(u_{\theta^{*}},v_{M}),\quad\forall v_{M}\in V_{ M}.\]

_Therefore, \(u_{\theta}\) converges to a solution of the Petrov-Galerkin problem (11), as expected._

### Energy norm error estimates based on local semi-discrete inf-sup condition

Even if the a posteriori error estimates of the previous section only ensure that the RVPINNs definition of the loss functional is robust in the sense of equivalence classes, numerical experiments confirm that it can also be robust concerning the true error in the energy norm.

We first notice that, as a consequence of Theorem 2 and the quotient norm definition (28), the following result can be deduced:

**Corollary 3** (Lower bound for the error in the energy norm).: _Under the same hypothesis of Theorem 2, it holds:_

\[\frac{1}{\mu}\|\phi\|_{V}\leq\|u-u_{\theta}\|_{U}. \tag{34}\]

Therefore, for a given \(\theta\), \(\|\phi\|_{V}\) always defines an efficient bound for the error in the energy norm.

If we assume that the NN structure allows for solutions of the Petrov-Galerkin problem (11), we can also obtain a local reliable bound for \(\|u-u_{\theta}\|_{U}\) through the following Assumption:

**Assumption 1** (Local Fortin's condition).: _There exists a parameter \(\theta^{*}\in\mathbb{R}^{S}\), such that \(u_{\theta^{*}}\in U_{NN}\) is a solution of the Petrov-Galerkin problem (11). Additionally, there exists \(R>0\) such that, for all \(\theta\in B(\theta^{*},R)\), there is an operator \(\Pi_{\theta}:V\mapsto V_{M}\), and a \(\theta\)-independent constant \(C_{\Pi}>0\), verifying:_

1. \(a(u_{\theta},v-\Pi_{\theta V})=0,\quad\forall v\in V\)_,_
2. \(\|\Pi_{\theta}v\|_{V}\leq C_{\Pi}\|v\|_{V},\quad\forall v\in V\)_,_

_where \(B(\theta^{*},R)\) denotes an open ball of center \(\theta^{*}\) and radius \(R\), with respect to a given norm of \(\mathbb{R}^{S}\)._

**Proposition 4** (Upper bound of the error in the energy norm).: _Under the same hypothesis of Theorem 2. If Assumption 1 is satisfied, it holds:_

\[\|u-u_{\theta}\|_{U}\leq\frac{1}{\alpha}\operatorname{osc}(u)+\frac{1}{C_{ \Pi}\alpha}\|\phi\|_{V},\quad\forall\,\theta\in B(\theta^{*},R), \tag{35}\]

_with_

\[\operatorname{osc}(u):=\sup_{0\neq v\in V}\frac{a(u,v-\Pi_{\theta}v)}{\|v\|_ {V}}.\]

Proof.: If Assumption 1 is satisfied, it holds:

\[\alpha\|u-u_{\theta}\|_{U}\leq\sup_{0\neq v\in V}\frac{a(u-u_{ \theta},v)}{\|v\|_{V}} =\sup_{0\neq v\in V}\frac{a(u-u_{\theta},v-\Pi_{\theta}v)}{\|v\|_ {V}}+\sup_{0\neq v\in V}\frac{a(u-u_{\theta},\Pi_{\theta}v)}{\|v\|_{V}}\] \[\leq\operatorname{osc}(u)+\frac{1}{C_{\Pi}}\sup_{0\neq v\in V} \frac{a(u-u_{\theta},\Pi_{\theta}v)}{\|v\|_{V}}\] \[\leq\operatorname{osc}(u)+\frac{1}{C_{\Pi}}\sup_{0\neq v_{M}\in V _{M}}\frac{a(u-u_{\theta},v_{M})}{\|v_{M}\|_{V}}\] \[=\operatorname{osc}(u)+\frac{1}{C_{\Pi}}\sup_{0\neq v_{M}\in V_{ M}}\frac{r(u_{\theta},v_{M})}{\|v_{M}\|_{V}}\] \[=\operatorname{osc}(u)+\frac{1}{C_{\Pi}}\sup_{0\neq v_{M}\in V_{ M}}\frac{(\theta,v_{M})_{V}}{\|v_{M}\|_{V}}\] \[=\operatorname{osc}(u)+\frac{1}{C_{\Pi}}\|\phi\|_{V},\]

proving (35).
**Remark 5** (Existence of a local Fortin's operator).: _First, notice that, for a given \(\theta\), a Fortin operator's existence is only possible if \(u_{\theta}\neq 0\) does not belong to the space \(U^{0}_{M}\). Indeed, if \(u_{\theta}\in U^{0}_{M}\), we have_

\[a(u_{\theta},\Pi_{\theta}v)=0,\text{ for all Fortin operator }\Pi_{\theta},\]

_since \(\Pi_{\theta^{\prime}}\in V_{M}\). Thus, in such a case, the residual representative \(\|\phi\|_{V}\) is only a reliable estimate in the spirit of Theorem 2. Nevertheless, it is expected that, in the minimization procedure, the NN solution approaches a solution of the Petrov-Galerkin problem (11) (see Remark 4), implying that the NN solution will not belong to the space \(U^{0}_{M}\) in a further iteration. Therefore, it is meaningful to assume that, after a sufficient number of iterations, the nonlinear solver will reach a parameter \(\theta\) belonging to a neighborhood of a parameter \(\theta^{*}\), in which \(\theta^{*}\) is the local minimum, implying that a local semi-discrete inf-sup condition is satisfied, as stated in the following Lemma (cf. [19])_

**Lemma 5** (Local Fortin's lemma).: _If Assumption 1 is satisfied, then the following local semi-discrete inf-sup condition is satisfied:_

\[\frac{\alpha}{C_{\Pi}}\|w_{NN}\|_{U}\leq\sup_{0\neq v_{M}\in V_{M}}\frac{a(w_ {NN},v_{M})}{\|v_{M}\|_{V}},\quad\forall\,w_{NN}\in U^{\theta^{*}}_{NN}:=\{u_ {\theta}:\theta\in B(\theta^{*},R)\}\,. \tag{36}\]

Proof.: Using the properties of Fortin's operator and the linearity of the bilinear form \(a(\cdot,\cdot)\) with respect to the first variable, it holds:

\[\sup_{0\neq v_{M}\in V_{M}}\frac{a(w_{NN},v_{M})}{\|v_{M}\|_{V}}\geq\sup_{0 \neq v\in V}\frac{a(w_{NN},\Pi_{\theta}v)}{\|\Pi_{\theta}v\|_{V}}=\sup_{0\neq v \in V}\frac{a(w_{NN},v)}{\|\Pi_{\theta}v\|_{V}}\geq\frac{1}{C_{\Pi}}\sup_{0 \neq v\in V}\frac{a(w_{NN},v)}{\|v\|_{V}}\geq\frac{\alpha}{C_{\Pi}}\|w_{NN}\|_ {U}.\]

Finally, with the help of the previous Lemma, we can prove the following local a priori error estimate:

**Proposition 6** (Local a priori error estimate).: _If Assumption 1 is satisfied, it holds:_

\[\|u-u_{\theta^{*}}\|_{U}\leq\left(1+\frac{\mu C_{\Pi}}{\alpha}\right)\inf_{w_ {NN}\in U^{\theta^{*}}_{NN}}\|u-w_{NN}\|_{U}. \tag{37}\]

Proof.: First notice that, for all \(w_{NN}\in U^{\theta^{*}}_{NN}\), it holds

\[\frac{\alpha}{C_{\Pi}}\|u_{\theta^{*}}-w_{NN}\|_{U}\leq\sup_{0\neq v_{M}\in V _{M}}\frac{a(u_{\theta^{*}}-w_{NN},v_{M})}{\|v_{M}\|_{V}}=\sup_{0\neq v_{M} \in V_{M}}\frac{a(u-w_{NN},v_{M})}{\|v_{M}\|_{V}}\leq\mu\|u-w_{NN}\|_{U}.\]

Therefore, the result (37) follows from the triangular inequality

\[\|u-u_{\theta^{*}}\|_{U}\leq\|u-w_{NN}\|_{U}+\|u_{\theta^{*}}-w_{NN}\|_{U}.\]

**Remark 6** (Fortin's lemma converse).: _Following [18], it can be proved that the converse of Lemma 5 is also valid. This could be useful in scenarios where proving the local semi-discrete inf-sup condition is simpler than proving the existence of a Fortin operator._

**Remark 7** (Norm for the discrete space).: _For the sake of simplicity, we assume that the norm for the discrete space \(V_{M}\) coincides with the norm \(\|\cdot\|_{V}\), of the space \(V\). However, under certain conditions, it could be convenient to consider a different norm for the space \(V_{M}\) being equivalent to \(\|\cdot\|_{V}\). For instance, when properties of the discrete test allow obtaining a better inf-sup constant for the local semi-discrete inf-sup condition. In such a case, previous results are straightforwardly extended to this scenario._

## 5 Numerical examples

To illustrate the proposed strategy's performance, we consider several diffusion-advection problems subject to Dirichlet-type boundary conditions. We focus on problems that either are challenging for the standard finite element method (advection-dominated problems) or do not admit a solution belonging to \(H^{2}(\Omega)\). To construct the loss functionals, we consider two distinct discrete test spaces: one employing standard FE piece-wise polynomials and another employing spectral orthonormal test functions. Here, we do not intend to compare the performance between different test spaces, but rather our aim is to numerically validate the robustness of RVPINNs.
### Diffusion-advection model problem

For \(d\geq 1\), given a bounded and open Lipschitz polyhedra \(\Omega\in\mathbb{R}^{d}\) with boundary \(\partial\Omega\), this work considers the following linear diffusion-advection problem:

Find \(u\), such that:

\[-\nabla\cdot(\varepsilon\nabla u)+\beta\cdot\nabla u =f,\quad\text{ in }\Omega,\] \[u =0,\quad\text{ on }\partial\Omega, \tag{38}\]

where \(\varepsilon>0\) denotes a diffusive coefficient term, \(\beta\in[L^{m}(\Omega)]^{d}\) is a divergence-free (i.e., \(\text{div}\beta\equiv 0\)) advective coefficient function. Finally, \(f\in V^{\prime}\) is a given source term, with \(V^{\prime}:=H^{-1}(\Omega)\) being the dual space of \(V:=H^{1}_{0}(\Omega)\). Problem (38) admits the following continuous variational formulation:

\[\text{Find }u\in U:=V\,:\,r(u,v):=l(v)-a(u,v),\,\forall\,v\in V, \tag{39}\]

with

\[a(u,v):=\left(\varepsilon\nabla u-\beta\,u,\nabla v\right)_{0},\quad\text{and }\quad l(v)=\langle f,v\rangle, \tag{40}\]

where \((\cdot,\cdot)_{0}\) denotes the \(L^{2}\)-inner product2, and \(\langle\cdot,\cdot\rangle\) denotes the duality map between \(V^{\prime}\) and \(V\), coinciding with \((f,v)_{0}\) when \(f\) belongs to \(L^{2}(\Omega)\).

Footnote 2: We adopt the same notation for the scalar and vectorial \(L^{2}\)-inner products for the sake of simplicity.

Let us equip the Hilbert spaces \(U,V\) with the norms

\[\|\cdot\|_{U}^{2}:=\|\cdot\|_{V}^{2}:=\varepsilon\,(\nabla\cdot,\nabla\cdot)_ {0}. \tag{41}\]

Using the Cauchy-Schwarz inequality, the following boundedness estimate is proved (cf. (2)):

\[a(w,v)\leq\left(1+\frac{C_{\Omega}\,\|\beta\|_{w}}{\varepsilon}\right)\|w\|_{ U}\|v\|_{V},\quad\forall\,w\in U,v\in V, \tag{42}\]

where \(C_{\Omega}\) denotes the Poincare constant (see [39, 1]). Moreover, as a consequence of the fact that (recall that \(\text{div}\beta\equiv 0\)):

\[(\beta\cdot\nabla v,v)_{0}=0,\quad\forall\,v\in V, \tag{43}\]

a straightforward verification shows that the following coercive estimate holds:

\[a(v,v)\geq\|v\|_{V}^{2},\quad\forall\,v\in V. \tag{44}\]

Notice also that the coercive property (44) implies that the following inf-sup condition is satisfied (cf. (3)):

\[\sup_{0\neq v\in V}\frac{a(w,v)}{\|v\|_{V}}\geq\|w\|_{U},\quad\forall\,w\in U. \tag{45}\]

### Discrete setting

For simplicity, we consider only one-dimensional problems. We set \(\Omega:=(-1,1)\) and define two kinds of test spaces. First, we consider \(V_{M}\) as FE space of the standard globally continuous and piece-wise linear functions that belong to \(H^{1}_{0}(\Omega)\), defined over a uniform mesh partition of \(\Omega\) into \(M+1\) elements \([x_{i},x_{i+1}]\), with \(x_{0}=-1<x_{1}<\dots<x_{M}<x_{M+1}=1\). The corresponding Gram matrix is constructed using the inner product, inducing the norm \(\|\cdot\|_{V}\), defined in (41). The second test space is intended to show examples involving an orthonormal basis for the discrete space \(V_{M}\). Here, we normalize with respect to the inner product (41) the first \(M\) eigenfunctions of the Laplacian satisfying the Dirichlet homogeneous boundary conditions. The \(m\)-th eigenfunction and the corresponding orthonormal set, respectively, are given by

\[s_{m}:=\sin\left(m\pi\frac{x+1}{2}\right),\qquad E_{M}=\left\{\varphi_{m}:= \frac{2s_{m}}{\sqrt{\varepsilon}\pi m}\right\}_{m=1}^{M}. \tag{46}\]

Therefore, the loss functional (26) is explicitly given, up to the constraint \(C(u_{\theta})\), as:

\[\mathcal{L}_{r}^{\phi}(u_{\theta})=\frac{4}{\varepsilon\pi^{2}}\sum_{m=1}^{M} \frac{1}{m^{2}}\,r\,(u_{\theta}\,,s_{m})^{2}\,. \tag{47}\]

For the architecture of the DNN approximator \(u_{\theta}\), we adopt the following structure: a fixed feed-forward fully connected network comprising five layers, each with 25 neurons. Throughout all layers, the activation function
employed is tanh. When strongly enforcing the homogeneous Dirichlet boundary conditions, we multiply the last layer's output by \((x+1)(x-1)\). Another version includes the constrained imposition of homogeneous Dirichlet boundary conditions by selecting:

\[C(u_{\theta})=|u_{\theta}(-1)|^{2}+|u_{\theta}(1)|^{2}. \tag{48}\]

As our choice for the nonlinear solver, we opt for the ADAM optimizer (refer to [28]) and initialize it with a learning rate of \(0.0005\). Our nonlinear solver undergoes up to \(6000\) iterations (epochs). To approximate the integral terms appearing in the loss functional definitions, we employ (a) a fixed Gaussian quadrature rule with five nodes per element when considering FE functions, and (b) a trapezoidal rule with \(4000\) equally distributed nodes when considering spectral functions. Lastly, the error \(|u-u_{\theta}|_{U}\) is numerically estimated using a trapezoidal rule, with \(10000\) equally distributed nodes.

### A smooth diffusion problem

We first consider a simple diffusion problem with a smooth solution. More precisely, setting \(\Omega=(-1,1)\), we consider the variational problem (39) with \(\epsilon=1\), \(\beta=0\) and forcing term \(f\) defined in such a way that the analytical solution is given by

\[u(x)=x\sin(\pi(x+1)).\]

Figure 1 shows the numerical results obtained with RVPINNs with \(50\) sinusoidal test functions and strong imposition of the BCs. We observe a good performance of the RVPINNs strategy, obtaining an accurate approximation of the analytical solution, as shown in Figure 0(a). Moreover, we observe a perfect match between the quantities \(\|u-u_{\theta}\|_{U}\), \(\|\phi\|_{V}\), and \(\sqrt{\mathcal{L}_{\nu}^{\phi}(u_{\theta})}\) as shown in figures 0(b) and 0(c) where, in particular, the last one shows a perfect correlation between the true error and the square root of the loss functional. Figure 2 shows an almost identical behavior when considering \(100\) FE test functions and strong imposition of BCs. We obtain similar results when considering the constrained imposition of BCs, but we omit them here for brevity.

Note that in the case of pure diffusive problems, the continuity constant \(\mu\) and the inf-sup stability constant \(\alpha\) are equal to one when considering the \(H^{1}\)-seminorm (cf. (42) and (45)). Therefore, \(\|\phi\|_{V}\) always defines a lower bound for the error (cf. (34)) and an upper bound up to a constant close to one if the oscillation term in (35) is sufficiently small. These theoretical considerations are confirmed in figures 1 and 2

Finally, to show the effect of the oscillation term (35), Figure 2(b) shows the results obtained with only five test FE basis functions. Figure 2(b) reveals that the estimation is not robust after (approximately) \(100\) iterations and only becomes a lower bound, as the DNN has reached (up to implementation precision) a solution of the Petrov-Galerkin problem (11).

### Delta source problem

We consider again the pure \(1\)D diffusion problem (i.e., eq. (39) with \(\epsilon=1\) and \(\beta=0\)), defined over \(\Omega=(-1,1)\) and subject to homogeneous Dirichlet BCs, but now with a Dirac delta source \(\delta_{1/2}\in H^{-1}(\Omega)\), that is explicitly defined through the action:

\[l(v):=\langle\delta_{1/2}\,,\,v\rangle:=v(1/2),\quad\forall v\in H^{1}_{0}( \Omega). \tag{49}\]

Figure 1: RVPINNs approximation of the smooth diffusion problem; strong BCs imposition, \(50\) spectral test functions.
The analytical solution to this problem is explicitly given as:

\[u(x)=\left\{\begin{array}{ll}\frac{1}{4}(x+1),&-1\leq x\leq\frac{1}{2},\\ \frac{3}{4}(1-x),&\frac{1}{2}<x\leq 1,\end{array}\right.\]

which is only a \(C^{0}(\Omega)\) solution, while the activation function tanh is smooth. Therefore, a major dominance of the oscillation term is expected, even when considering a sufficiently large number of test functions, as appreciated in Figure 4 where we have considered 100 FE test functions and strong imposition of the BCs. Figure 5 shows similar results for the same configuration but with the constrained imposition of the BCs. In both cases, a good representation of the analytical solution is obtained despite being approximated with a smooth function. Additionally, figures 4 and 5 show that the loss functional defines a robust estimation of the error in the pre-asymptotic regime (before the dominance of the oscillation term). A similar behavior is obtained when considering the loss functional with spectral test functions- omitted here for brevity.

Figure 4: RVPINNs approximation of the delta source diffusion problem; strong BCs imposition, 100 FE test functions.

Figure 3: RVPINNs approximation of the smooth diffusion problem; strong BCs imposition, 5 FE test functions.

Figure 2: RVPINNs approximation of the smooth diffusion problem; strong BCs imposition, 100 FE test functions.
### Advection-dominated diffusion problem

As our last example, we consider problem (39) in an advection-dominant regime, with the spirit of testing the performance and robustness of the estimates for RVPINNs in near unstable scenarios. For this, we set \(\beta=1\), the source term as \(f=1\), and consider small values of \(\epsilon\). For a given \(\epsilon\), the analytical solution of (39) in this case is given by:

\[u(x)=\frac{2(1-e^{\frac{\epsilon-1}{\epsilon}})}{1-e^{-\frac{ \epsilon}{2}}}+x-1,\]

exhibiting a strong gradient near the right boundary when \(\epsilon\to 0\) (cf. figures (a)a and (a)a). We restrict the discussion to spectral test functions, as we observe similar behaviors when considering FE test functions.

Figure 6 shows the results obtained for \(\epsilon=0.1\), considering 50 test functions and strong imposition of the BCs. Contrary to the pure diffusive case, Figure (b)b reflects that the lines representing the error and the square root of the loss functional (respectively, the norm of the discrete Riesz representative) do not overlap. This behavior is expected as, in this case, the continuity constant \(\mu\) is not equal to one. It deteriorates when \(\epsilon\to 0\) (cf. (42)). We also observe a decay in the expected optimal correlation between the error and the square root of the loss starting near iteration 3000, apparently due to stability issues. This suboptimality is improved when considering the constrained imposition of the boundary conditions, as shown in Figure 7. This constrained imposition of BCs could thus serve as a stabilization technique. This is further confirmed in Figure 8, where we show the results obtained for \(\epsilon=0.005\), with 200 test basis functions and constrained BCs. When considering instead the strong imposition of the BCs and the same number of test functions, we do not observe convergence within the first 6000 iterations--omitted here for brevity. We observe that when considering the constrained BCs, the loss functional exhibits some differences with the norm of the Riesz representative. However, this also agrees with our findings, as the constraint does not depend on the test space. Thus, it is not represented by the Riesz representative, as the boundary is not taken into consideration in the norm definition of the test space. Nevertheless, in this case, we still observe a good correlation between the loss and the true error.

Figure 5: RVPINNs approximation of the delta diffusion problem; constrained BCs imposition, 100 FE test functions.

Figure 6: RVPINNs approximation of the diffusion-advection problem with \(\epsilon=0.1\); strong BCs imposition, 50 spectral test functions.
## 6 Conclusions

In this article, we provide a general mathematical framework to construct robust loss functionals based on VPINNs. For that, we first generalize the definition of the loss functional in VPINNs to the choice of a single test function. Then, following a minimum residual principle, we select such test function as the Riesz representative of the weak residual over a given discrete test space. We expand the Riesz representation over the discrete test space so the loss functional includes the inversion of the Gram matrix corresponding to the inner product. We prove that the true error in the energy norm is equivalent to the test norm of the residual error estimator. To prove the robustness of the method, we need to select a test norm that induces the inf-sup stability at the continuous level. We numerically show the robustness of our a posteriori error estimator and also that our methodology is insensitive to the choice of the basis functions in the selected discrete test space.

Possible future research lines include: (a) the application of the framework to other families of PDEs like wave propagation, nonlinear or time-dependent problems, (b) the extension to nonconforming formulations, (c) the application of RVPINNs to parametric problems with an off-line inversion of the Gram matrix, and (d) an efficiency study of the method for a given variational formulation, NN configuration, and test space.

## 7 Acknowledgments

The work of Sergio Rojas was done in the framework of the Chilean grant ANID FONDECYT No. 3210009. Judit Munoz-Matute has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie individual fellowship grant agreement No. 101017984 (GEODPG). David Pardo has received funding from: the Spanish Ministry of Science and Innovation projects with references TED2021-132783B-I00, PID2019-1081118B-I00 (FEDER/AEI) and PDC2021-121093-I00 (MCIN / AEI / 10.13039/501100011033/Next Generation EU), the "BCAM Severo Ochoa" accreditation of excellence CEX2021-001142-S / MICIN / AEI / 1.13039/501100011033; the Spanish Ministry of Economic and Digital Transformation with Misiones Project IA4TES (MIA.2021.M04.008 / NextGenerationEU PRTR); and the Basque Government through the BERC 2022-2025 program, the Elkartek project SIGZE (KK-2021/00095), and the Consolidated Research Group MATHMODE (IT1456-22) given by the Department of Education, and the Euskampus Foundation through the ORLEG-IA project. The work of Maciej Paszynski and Pawel Maczuga was supported by

Figure 8: RVPINNs approximation of the diffusion-advection problem with \(\varepsilon=0.005\); constrained BCs imposition, 200 spectral test functions.

Figure 7: RVPINNs approximation of the diffusion-advection problem with \(\varepsilon=0.1\); constrained BCs imposition, 50 spectral test functions.
the program "Excellence initiative - research university" for the AGH University of Science and Technology.
