PG-VPINNs: Petrov–Galerkin Neural PDE Solvers

• PG-VPINNs are variational neural methods that combine neural network trial functions with independently chosen finite-dimensional test spaces to solve PDEs via weak formulations.
• They utilize weighted residual minimization and inf-sup stability principles to lower derivative requirements and enhance numerical accuracy, particularly in challenging regimes.
• Advanced techniques such as quadrature integration, finite-element interpolation, and dual-norm residuals are employed to ensure stability, robust error control, and computational efficiency.

Petrov–Galerkin Variational Physics-Informed Neural Networks (PG-VPINNs) are variational neural methods for solving partial differential equations in which the trial space is represented by neural networks and the test space is chosen independently, typically from finite-dimensional polynomial, spectral, or finite-element spaces. In contrast to strong-form PINNs, which penalize pointwise residuals, PG-VPINNs encode the PDE through weighted residuals or minimum-residual principles in weak form, so that the loss depends on projections of the residual against prescribed test functions rather than on Dirac-like collocation alone (Kharazmi et al., 2019, Berrone et al., 2022). Across the literature, the term covers several closely related formulations: classical variational residual minimization, least-squares weak forms, residual minimization in discrete dual norms, finite-element-interpolated and discontinuous-Galerkin realizations, and problem-specific stabilized constructions for singular perturbations, advection-dominated transport, homogenization, and parabolic evolution (Mardal et al., 2023, Grekas et al., 4 Jul 2025).

1. Emergence from PINNs and variational weighted-residual methods

A canonical PG-VPINN formulation was introduced by recasting physics-informed training in a Petrov–Galerkin setting: the approximate solution $\tilde u$ is a neural network, while the test space is chosen as a linear span of analytic functions such as Legendre polynomials or sine functions (Kharazmi et al., 2019). In that framework, standard PINNs were interpreted as a limiting strong-form construction in which Dirac delta functions act as test functions, whereas VPINNs replace pointwise enforcement by integral residuals over a prescribed test space (Kharazmi et al., 2019).

This variational reinterpretation matters because it changes both the regularity requirements and the numerical structure of training. Strong-form least-squares PINNs for elliptic problems minimize losses such as

$$L_1(v; f, g) = \|-\Delta v - f\|_{L^2(\Omega)}^2 + \alpha^2 \|v - g\|_{L^2(\partial \Omega)}^2,$$

and this least-squares minimization can be recast as a Petrov–Galerkin-type variational problem,

$$(\Delta u, \Delta v)_\Omega + \alpha^2 (u, v)_{\partial\Omega} = (f, -\Delta v)_\Omega + \alpha^2 (g, v)_{\partial\Omega} \qquad \forall v \in X,$$

which makes explicit the hidden $H^2$-type regularity demands imposed by losses involving second derivatives (Mardal et al., 2023). The same paper emphasizes that the corresponding strong form is a biharmonic PDE with mixed boundary conditions, so the weak interpretation is not merely cosmetic: it clarifies what class of functions the optimization is implicitly searching over (Mardal et al., 2023).

A central historical feature of PG-VPINNs is therefore the transfer of established Petrov–Galerkin ideas—distinct trial and test spaces, inf-sup reasoning, weak consistency, and weighted residual minimization—into neural PDE solvers. This also explains why the literature frequently combines neural-network trial spaces with finite-element, spline, or spectral test spaces rather than attempting a purely neural Galerkin construction (Kharazmi et al., 2019, Berrone et al., 2022).

2. Mathematical structure of the PG-VPINN loss

For a variational problem written as

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

a standard PG-VPINN construction selects a neural-network manifold for the trial space and a finite-dimensional test space $V_h\subset V$ with basis $\{\varphi_i\}$. The discrete residuals are then

$r_{h,i}(w)=F_h(\varphi_i)-a_h(w,\varphi_i),$

and the loss takes the form

$R_h^2(w)=\sum_{i\in I_h} r_{h,i}^2(w),$

or weighted variants of this expression when norm-equivalence constants are included (Berrone et al., 2022). The resulting optimizer is the neural-network function minimizing the squared weak residual over the chosen test basis (Berrone et al., 2022).

In the original VPINN development, the variational residual for a linear operator $\mathcal L^q$ was written as

$$L_1(v; f, g) = \|-\Delta v - f\|_{L^2(\Omega)}^2 + \alpha^2 \|v - g\|_{L^2(\partial \Omega)}^2,$$0

with a corresponding loss

$$L_1(v; f, g) = \|-\Delta v - f\|_{L^2(\Omega)}^2 + \alpha^2 \|v - g\|_{L^2(\partial \Omega)}^2,$$1

A major practical consequence is that integration by parts can move derivatives from the neural network to the test functions, thereby lowering the differentiation order required of the trial function (Kharazmi et al., 2019). For the 1D Poisson problem, one, two, or zero integrations by parts lead to equivalent variational residuals with different derivative placement, and the paper explicitly argues that this order reduction lowers training cost while improving accuracy relative to strong-form PINNs (Kharazmi et al., 2019).

Later work sharpened the same idea by formulating PG-VPINN training as residual minimization in a discrete dual norm rather than as a simple basis-dependent sum of squares. In the robust formulation,

$$L_1(v; f, g) = \|-\Delta v - f\|_{L^2(\Omega)}^2 + \alpha^2 \|v - g\|_{L^2(\partial \Omega)}^2,$$2

where $$L_1(v; f, g) = \|-\Delta v - f\|_{L^2(\Omega)}^2 + \alpha^2 \|v - g\|_{L^2(\partial \Omega)}^2,$$3 is the Gram matrix of the test basis in the test norm. This removes dependence on the chosen basis scaling and interprets the loss as the squared discrete dual norm of the residual (Rojas et al., 2023). That shift is technically important because it aligns the optimization objective with residual norms that admit reliable and efficient error estimators under Fortin-type assumptions (Rojas et al., 2023).

3. Trial spaces, test spaces, quadrature, and boundary treatment

The defining Petrov–Galerkin feature of PG-VPINNs is that the test space is not the neural trial space. The literature exhibits several recurring choices. Legendre polynomials and sine functions appear in the original spectral VPINN constructions (Kharazmi et al., 2019). Continuous, piecewise linear finite-element functions on triangulations were used in an a posteriori analysis of VPINNs for elliptic boundary-value problems (Berrone et al., 2022). More general piecewise polynomial spaces of degree $$L_1(v; f, g) = \|-\Delta v - f\|_{L^2(\Omega)}^2 + \alpha^2 \|v - g\|_{L^2(\partial \Omega)}^2,$$4 were analyzed in a quadrature-focused Petrov–Galerkin framework, where the interpolant degree was tied to quadrature precision by $$L_1(v; f, g) = \|-\Delta v - f\|_{L^2(\Omega)}^2 + \alpha^2 \|v - g\|_{L^2(\partial \Omega)}^2,$$5 (Berrone et al., 2021). Singularly perturbed one- and two-dimensional problems were treated with localized hat functions and tensor-product hat functions, respectively (Kumar et al., 13 Sep 2025, Kumar et al., 15 Jun 2026). In multiphase elasticity, local piecewise linear polynomials supported on adaptive integration cells were employed as test functions (Rodopoulos et al., 29 Jun 2025).

These choices are not merely implementational. The 2021 quadrature study showed that the convergence rate with mesh refinement depends explicitly on the precision $$L_1(v; f, g) = \|-\Delta v - f\|_{L^2(\Omega)}^2 + \alpha^2 \|v - g\|_{L^2(\partial \Omega)}^2,$$6 of the quadrature rule and on the degree of the test functions, with the estimate

$$L_1(v; f, g) = \|-\Delta v - f\|_{L^2(\Omega)}^2 + \alpha^2 \|v - g\|_{L^2(\partial \Omega)}^2,$$7

for the decay rate of the energy error between the exact solution and a suitable interpolant of the computed neural network (Berrone et al., 2021). The same work reported the counterintuitive conclusion that, for smooth solutions, the best strategy is to choose test functions of the lowest polynomial degree while using quadrature formulas of suitably high precision (Berrone et al., 2021).

A separate line of work replaces pointwise strong residual evaluation by finite-element interpolation and elementwise integration. PINN-DG introduces a DG-inspired loss in which the neural network is interpolated into a finite-element space and the residual is integrated over elements and faces rather than sampled by Monte Carlo quadrature. For Poisson’s equation, the discrete energy is

$$L_1(v; f, g) = \|-\Delta v - f\|_{L^2(\Omega)}^2 + \alpha^2 \|v - g\|_{L^2(\partial \Omega)}^2,$$8

with a penalty enforcing weak continuity and boundary data (Grekas et al., 4 Jul 2025). This construction avoids repeated pointwise second-derivative evaluations of the neural network and instead differentiates the finite-element basis functions analytically (Grekas et al., 4 Jul 2025).

Boundary treatment is another major design axis. Four strategies were compared systematically for PINNs and VPINNs: penalization in the loss, exact imposition through a modified output using an approximate distance function, a product-type distance-function variant, and Nitsche’s method in the variational residual (Berrone et al., 2022). The reported conclusion was that the best results are achieved by exactly enforcing the Dirichlet boundary conditions by means of an approximate distance function, while variational imposition through Nitsche’s method leads to suboptimal solvers (Berrone et al., 2022). Hard enforcement also appears in singularly perturbed PG-VPINNs through ansätze such as $$L_1(v; f, g) = \|-\Delta v - f\|_{L^2(\Omega)}^2 + \alpha^2 \|v - g\|_{L^2(\partial \Omega)}^2,$$9 in one dimension and $$(\Delta u, \Delta v)_\Omega + \alpha^2 (u, v)_{\partial\Omega} = (f, -\Delta v)_\Omega + \alpha^2 (g, v)_{\partial\Omega} \qquad \forall v \in X,$$0 in two dimensions (Kumar et al., 13 Sep 2025, Kumar et al., 15 Jun 2026).

4. Stability, convergence, and error control

The theoretical literature on PG-VPINNs is dominated by inf-sup stability, residual-based a posteriori estimation, coercivity, and variational convergence. In the quadrature-and-test-function analysis, the key stability hypothesis is a discrete inf-sup condition of the form

$$(\Delta u, \Delta v)_\Omega + \alpha^2 (u, v)_{\partial\Omega} = (f, -\Delta v)_\Omega + \alpha^2 (g, v)_{\partial\Omega} \qquad \forall v \in X,$$1

which underpins the a priori error estimate between the exact solution and the interpolated neural approximation (Berrone et al., 2021). The same study documented that non-interpolated formulations or insufficiently rich test sets can admit spurious neural solutions for zero-data problems, illustrating that Petrov–Galerkin structure alone does not eliminate underdetermination (Berrone et al., 2021).

A posteriori theory for elliptic VPINNs provides a local error estimator composed of three parts: a residual-type term, a loss-function term, and data oscillation terms (Berrone et al., 2022). With local indicators $$(\Delta u, \Delta v)_\Omega + \alpha^2 (u, v)_{\partial\Omega} = (f, -\Delta v)_\Omega + \alpha^2 (g, v)_{\partial\Omega} \qquad \forall v \in X,$$2, the global estimator

$$(\Delta u, \Delta v)_\Omega + \alpha^2 (u, v)_{\partial\Omega} = (f, -\Delta v)_\Omega + \alpha^2 (g, v)_{\partial\Omega} \qquad \forall v \in X,$$3

was proved reliable and efficient for the energy norm error, and numerical tests showed close agreement between estimator and true $$(\Delta u, \Delta v)_\Omega + \alpha^2 (u, v)_{\partial\Omega} = (f, -\Delta v)_\Omega + \alpha^2 (g, v)_{\partial\Omega} \qquad \forall v \in X,$$4-error under mesh refinement (Berrone et al., 2022). This is a major distinction from strong-form PINNs, for which comparable finite-element-style localized estimators are not part of the standard formulation.

Robust VPINNs reformulate the loss as the discrete dual norm of the residual and show that, under the existence of a local Fortin operator, this norm is a reliable and efficient estimator of the true error in the energy norm (Rojas et al., 2023). The paper also stresses that the classical VPINN loss depends heavily on the chosen basis functions of the test space, whereas the Gram-matrix-corrected dual-norm loss is basis independent (Rojas et al., 2023).

PINN-DG extends the variational theory to a DG-inspired energy on finite-element interpolants. The discrete functional is shown coercive: $$(\Delta u, \Delta v)_\Omega + \alpha^2 (u, v)_{\partial\Omega} = (f, -\Delta v)_\Omega + \alpha^2 (g, v)_{\partial\Omega} \qquad \forall v \in X,$$5 and a $$(\Delta u, \Delta v)_\Omega + \alpha^2 (u, v)_{\partial\Omega} = (f, -\Delta v)_\Omega + \alpha^2 (g, v)_{\partial\Omega} \qquad \forall v \in X,$$6-convergence-type result is established: $$(\Delta u, \Delta v)_\Omega + \alpha^2 (u, v)_{\partial\Omega} = (f, -\Delta v)_\Omega + \alpha^2 (g, v)_{\partial\Omega} \qquad \forall v \in X,$$7 for discrete minimizers under refinement (Grekas et al., 4 Jul 2025). The same analysis introduces a lifted discrete Laplacian

$$(\Delta u, \Delta v)_\Omega + \alpha^2 (u, v)_{\partial\Omega} = (f, -\Delta v)_\Omega + \alpha^2 (g, v)_{\partial\Omega} \qquad \forall v \in X,$$8

which combines finite-element projection with a lifting of jump terms and is central to both stable evaluation and convergence proofs (Grekas et al., 4 Jul 2025).

5. Stabilized formulations and problem-specific variants

PG-VPINNs have been developed most aggressively for PDE classes where strong-form PINNs are known to be fragile. Advection-dominated diffusion is a recurrent example. A comparative study of strong and weak residual neural methods showed that standard Galerkin FEM cannot solve a model advection-dominated problem without oscillations, discussed stabilization through Petrov–Galerkin FEM and SUPG, and found that PINN and VPINN formulations on uniform meshes fail for sufficiently small diffusion parameters, while adaptive point or test-function placement can recover accurate solutions (Sikora et al., 2023). The same study explicitly noted that no fully explicit Petrov–Galerkin trick was implemented in the VPINN itself; rather, richer and adapted test functions act in the spirit of PG stabilization (Sikora et al., 2023).

For singularly perturbed boundary-value problems and parabolic PDEs, a 2025 PG-VPINN framework employed neural-network trial functions and localized hat functions as test functions, introduced interface penalty terms to enhance numerical stability and accurately capture boundary layers, imposed Dirichlet data through hard constraints, and reported significantly improved accuracy in both $$(\Delta u, \Delta v)_\Omega + \alpha^2 (u, v)_{\partial\Omega} = (f, -\Delta v)_\Omega + \alpha^2 (g, v)_{\partial\Omega} \qquad \forall v \in X,$$9 and maximum norms compared to the standard VPINN approach for one-dimensional singularly perturbed differential equations (Kumar et al., 13 Sep 2025). A two-dimensional extension used tensor-product hat functions as test functions and direct imposition of Dirichlet conditions; the reported computational experiments showed high accuracy in both the maximum and $H^2$0 norms and emphasized robustness in resolving sharp boundary layers for one- and two-parameter singularly perturbed problems (Kumar et al., 15 Jun 2026).

Weak-form neural methods have also been adapted to discontinuous coefficients and multiscale media. In periodic homogenization of thermo-conductive composites, variational PINNs with spectral and neural-network-based test functions were compared against strong-form PINNs. The reported conclusion was that strong-form PINNs may outperform VPINNs in controlled settings but are sensitive to material discontinuities and may fail without clear diagnostics, whereas VPINNs accommodate piecewise constant material parameters directly but require careful selection of test functions to avoid instability (Gaynutdinova et al., 9 Sep 2025). The same work introduced a dual formulation yielding guaranteed upper and lower bounds and used the primal-dual gap as an internal convergence diagnostic (Gaynutdinova et al., 9 Sep 2025).

In multiphase architected materials, a Petrov–Galerkin VPINN with local polynomial test functions was combined with a Galerkin Boundary Element Method data generator and a quadtree-based adaptive integration algorithm (Rodopoulos et al., 29 Jun 2025). There the weak form was used to handle displacement discontinuities and stress jumps at interfaces, and the tree-based integration scheme was reported to capture inner material discontinuities robustly while reducing computational cost relative to uniform integration (Rodopoulos et al., 29 Jun 2025).

Parabolic evolution problems have also been treated variationally. A time-discrete VPINN for heat conduction used backward Euler in time and, at each time step, minimized the dual norm of the weak residual. In that formulation the trial and test spaces coincide, so the method reduces to a Galerkin or Bubnov–Galerkin style rather than a strict Petrov–Galerkin one, but it remains part of the same weak-form neural lineage (Olivares et al., 5 Mar 2026).

6. Optimization, automatic differentiation, and computational realizations

Optimization in PG-VPINNs is not restricted to end-to-end stochastic gradient descent. A 2024 least-squares formulation observed that the residual is linear in the weights of the last neural-network layer, so those weights can be updated by solving a least-squares problem while the hidden parameters are optimized by gradient descent (Uriarte et al., 2024). With

$H^2$1

the last-layer update solves

$H^2$2

followed by a gradient step in the hidden weights (Uriarte et al., 2024). The same work showed that a traditional backward-mode automatic-differentiation implementation is prohibitively expensive, and proposed either forward-mode automatic differentiation or an ultraweak-type scheme that avoids differentiation of trial functions in the discrete weak formulation; the proposed alternatives were reported to be up to one hundred times faster than the traditional one (Uriarte et al., 2024).

Finite-element interpolation can further reduce differentiation costs. PINN-DG reported 3–11x speedups and substantially reduced GPU memory usage compared to collocation-based PINNs, precisely because second derivatives are computed analytically on finite-element basis functions rather than repeatedly by automatic differentiation on the raw network (Grekas et al., 4 Jul 2025).

A neighboring Petrov–Galerkin family replaces fully trainable networks with randomized feature models. Randomized Neural Networks with Petrov-Galerkin Methods approximate the solution by randomized neural networks, keep the hidden layers fixed after random initialization, and determine the output layer by linear least squares against finite-element-type test functions (Shang et al., 2023). Although these methods are not usually labeled VPINNs, they occupy the same trial-space/test-space design niche and show how Petrov–Galerkin neural solvers can be realized without nonlinear optimization of all weights (Shang et al., 2023).

7. Misconceptions, limitations, and extensions

A recurrent misconception is that increasing the polynomial degree of test functions necessarily improves PG-VPINN accuracy. The quadrature analysis showed the opposite for smooth elliptic problems: the best strategy is to use the lowest polynomial degree for the test space and increase quadrature precision instead (Berrone et al., 2021). Another misconception is that weak-form training automatically guarantees stability. Multiple papers state more qualified conditions: correct test/trial containment matters in inverse formulations, discrete inf-sup conditions matter in elliptic formulations, Fortin operators matter for reliable dual-norm residual estimates, and test-function selection remains critical in discontinuous or advection-dominated settings (Mardal et al., 2023, Rojas et al., 2023, Gaynutdinova et al., 9 Sep 2025).

The literature also distinguishes between explicit and implicit uses of Petrov–Galerkin stabilization. In advection-dominated diffusion, one study emphasized that its VPINN did not implement an explicit Petrov–Galerkin stabilization analogous to SUPG or optimal test functions; instead, adapted test-function placement mimicked some PG effects (Sikora et al., 2023). This objective caveat is important because it separates genuinely Petrov–Galerkin neural formulations from weak residual methods that merely use different trial and test spaces without tailored stability mechanisms.

More recent work pushes the PG viewpoint beyond solution networks. PG-VarMiON, an operator-learning framework, embeds the optimal Petrov–Galerkin weak form into the network architecture and trains a neural operator to learn the optimal weighting functions implicitly (Charles et al., 6 Mar 2025). This suggests a broader continuation of PG-VPINN ideas from instance-wise PDE solvers to supervised operator approximation.

A unified generalization analysis covering PINNs and VPINNs showed that nonlinear differential operators can be represented by Taylor expansion on a higher-dimensional space and that the nonlinearity of the differential operator exponentially enlarges the bound, highlighting its significant impact on generalization (Hashimoto et al., 13 May 2026). The same analysis accommodates arbitrary test functions in weak-form neural methods, so it is directly relevant to Petrov–Galerkin VPINNs and reinforces a theme visible throughout the literature: weak formulations improve regularity and stability properties, but they do not remove the dependence of accuracy and generalization on operator structure, test-space design, and computational realization (Hashimoto et al., 13 May 2026).

In aggregate, PG-VPINNs form a technically diverse but conceptually coherent class of neural PDE solvers. Their unifying principle is not a single architecture or optimizer, but the use of neural trial spaces within a Petrov–Galerkin or minimum-residual variational framework that imports quadrature, inf-sup theory, stabilization, dual norms, and finite-element machinery into physics-informed learning (Kharazmi et al., 2019, Grekas et al., 4 Jul 2025).