arXiv CSopen access

Beyond PINNs: A Unified Gauss--Newton and Petrov--Galerkin Framework for Neural and Hybrid PDE Solvers

BEYOND PINNS: A UNIFIED GAUSS–NEWTON AND PETROV–GALERKIN FRAMEWORK FOR NEURAL AND HYBRID PDE SOLVERS

arXiv:2609.20641v1 [math.NA] 17 Sep 2026

NILO SCHWENCKE∗ AND ROLAND MAIER†

Abstract. Physics-informed neural networks and finite element methods provide two different paradigms for the numerical approximation of partial differential equations: the former are commonly trained by minimizing pointwise strong residuals, whereas the latter are naturally built from weak variational formulations and the finite-dimensional systems obtained after discretization. In this work, we introduce a common framework based on the discretization of functional Gauss–Newton problems by finite families of linear measurements. We show that, through an appropriate duality pairing, the linear measurements can be represented by test functions. The resulting Gauss–Newton system is then precisely a Petrov–Galerkin discretization of the linearized functional problem. This perspective recovers pointwise collocation and natural-gradient constructions as particular cases, while making the choice of test functions an explicit algorithmic design choice. We specialize this framework to elliptic problems, where it naturally leads to weak residual formulations and to a hybrid finite element–neural construction acting on complementary approximation spaces. Numerical experiments support the proposed framework and demonstrate the effectiveness of weak Gauss–Newton formulations and hybrid finite element–neural approximations.

1. Introduction Partial differential equations (PDEs) constitute one of the principal mathematical tools for modeling physical phenomena across science and engineering. Their numerical approximation has consequently motivated a wide range of methods, among which Galerkin and Petrov–Galerkin discretizations play a central role: the solution is approximated in a prescribed finite-dimensional approximation space, while the governing equation is enforced against a suitable family of test functions through its variational formulation. Finite element methods are one of their most successful realizations, typically using piecewise-polynomial approximation spaces constructed on a mesh; see, e.g., [7, 4, 5]. Neural-network-based PDE solvers provide a rather different approximation paradigm. In particular, physics-informed neural networks (PINNs) [11, 23, 37] represent the solution by a nonlinear parametric model and determine its parameters by minimizing residuals of the governing equation and boundary conditions. Their mesh-free character and the flexibility of the neural approximation model make them attractive in settings where classical discretizations may be difficult to construct. Although early PINN formulations often exhibited limited accuracy and difficult optimization, substantial progress has been obtained through adaptive sampling [33, 44, 28, 9, 34, 24] and, in particular, through Gauss–Newton and natural-gradient-based optimization [32, 39, 19, 30, 40, 18, 20, 35, 31, 43]. Despite these developments, the classical finite element and PINN viewpoints remain rather different. Standard PINNs usually start from a strong formulation and construct a finite residual vector by evaluating the differential equation at collocation points. For a second-order elliptic operator, this requires the strong residual and its parameter derivatives to possess sufficient pointwise regularity. In contrast, the natural formulation used by finite element methods is typically 2020 Mathematics Subject Classification. 65N30, 65K10, 68T07. Key words and phrases. PINNs, Gauss–Newton methods, Petrov–Galerkin methods, variational techniques, finite element methods, hybrid FEM-NNs methods. 1

2

N. SCHWENCKE, R. MAIER

weak: the elliptic residual belongs to H −1 (Ω) and is evaluated through its action on test functions in H01 (Ω). This distinction becomes particularly relevant for nonsmooth solutions, distributional right-hand sides, or FE functions themselves, for which a pointwise strong residual may be unnatural or even undefined. Several neural approaches have investigated variational or weak formulations. VPINNs and hp-VPINNs [21, 22] construct variational losses from prescribed test spaces, while related Petrov– Galerkin methods allow more general finite element or neural test families [41]. Other approaches adapt the test spaces themselves, through adversarial optimization [45] or minimum-residual formulations [38] connected to the classical theory of optimal Petrov–Galerkin test spaces [10]. Closer to the present finite-measurement viewpoint, a related construction has recently been proposed in [3] in the context of nonlinear dynamical approximation for time-dependent PDEs. Therein, the parameters of a nonlinear decoder are evolved by projecting the PDE dynamics onto its tangent space, with this projection approximated through finite families of linear observations. The analysis focuses in particular on the stability of the resulting reconstruction and on the adaptive selection of the observations. Our setting differs in that the finite measurements arise from the discretization of a functional Gauss–Newton problem, a perspective that will lead to the Petrov– Galerkin interpretation developed below. At the same time, hybrid methods combining finite element and neural approximations have attracted increasing attention, see, e.g., [13, 1, 29]. These developments raise a common question: can the finite residual systems used by modern Gauss–Newton neural solvers, weak Petrov–Galerkin formulations, and hybrid finite element–neural approximations be cast within a single construction and handled by the same optimization machinery? The starting point of this work is the observation that Gauss–Newton, formulated at the functional level, provides such a common framework. By applying a finite family of linear measurements to the linearized functional problem, one obtains an ordinary finite residual vector and Jacobian to which standard Gauss–Newton algorithms can be applied. When these measurements are represented through pairings with test functions, the resulting system is precisely a Petrov–Galerkin discretization of the linearized problem. This interpretation recovers natural-gradient and pointwise collocation constructions as particular cases, while making the choice of measurements, and hence of test functions, an independent component of the numerical method. For elliptic PDEs, this flexibility naturally leads to weak residuals in H −1 (Ω) tested against functions in H01 (Ω), allowing the same Gauss–Newton machinery to operate directly at the variational level with a broad choice of test families, including but not limited to operator-adapted Green sections. The same approximation–test viewpoint also suggests a hybrid finite element– neural construction in which a FE component captures a prescribed approximation space, while the neural component and the test functions act on its energy-orthogonal complement. This yields a unified weak Gauss–Newton formulation without requiring alternating optimization of the finite element and neural components. The main contributions of this work are therefore threefold: (1) we establish a Petrov–Galerkin interpretation of functional Gauss–Newton problems discretized by finite linear measurements; (2) we extend this construction to weak elliptic residuals with general test families; (3) we derive a hybrid finite element–neural formulation in which the finite element and neural components act on complementary energy subspaces. Numerical experiments validate these constructions on both smooth and low-regularity problems. They show that standard Gauss–Newton-based solvers can be applied directly to weak residuals, reaching accuracies comparable to or substantially higher than specialized energy-based methods, and that the hybrid finite element–neural formulation can substantially improve the robustness of

BEYOND PINNS: A UNIFIED FRAMEWORK FOR NEURAL AND HYBRID PDE SOLVERS

3

generic weak test spaces while remaining competitive with standalone finite element discretizations with a comparable number of degrees of freedom. The remainder of the paper is organized as follows. Section 2 introduces the model PDE, parametric approximation models, Galerkin and Petrov–Galerkin methods, neural networks, and PINNs. Section 3 develops finite-dimensional and functional Gauss–Newton frameworks, introduces discretization by linear measurements, and derives its Petrov–Galerkin interpretation together with the weak-test formulation. Section 4 presents the hybrid finite element–neural model and its associated weak-test Gauss–Newton update. Section 5 contains the numerical validation. Additional connections with natural gradient, kernel approximation, pointwise collocation, sketching, and regularization are presented in Section A. We also refer to a companion blog (see https://nilo.schwencke.me/tutorials/beyond-pinns-companion/) for more detailed experimental protocols and diagnostics as well as additional numerical examples. 2. Approximation of PDE solutions 2.1. Model problem. We consider a Hilbert space H ⊂ H 1 (Ω), where Ω ⊂ Rd (d ∈ N) is a bounded Lipschitz domain.1 We define operators (2.1)

D : H → L2 (Ω)

B : H → L2 (∂Ω),

and

where D is a linear differential operator and B a linear boundary condition. We define a scalar product ⟨•, •⟩ : H × H → R by (2.2)

⟨v, w⟩ := (D[v], D[w])L2 (Ω) + (B[v], B[w])L2 (∂Ω) ,

possibly factoring out the nullspace of the pair (D, B). Given f ∈ L2 (Ω) and gbc ∈ L2 (∂Ω), an abstract PDE formulation now seeks a function u ∈ H such that (2.3)

D[u] = f

in L2 (Ω),

B[u] = gbc

in L2 (∂Ω).

We refer to (2.3) as a boundary value problem. Example 2.1 (second-order elliptic PDE). A classical example is the boundary value problem (2.4)

− div(A∇u) = f u = gbc

in Ω, on ∂Ω,

where A ∈ L∞ (Ω, Rd×d ) is symmetric and, for constants 0 < α ≤ β, satisfies (2.5)

α|η|2 ≤ (A(x)η) · η ≤ β|η|2

for almost all x ∈ Ω and every η ∈ Rd . Here, | • | denotes the Euclidean norm in Rd . In the strong formulation (2.3), one takes D = − div(A∇•) and B = γ on a sufficiently regular space Hstr ⊂ H 1 (Ω) such that D[u] ∈ L2 (Ω) and the trace γu is well defined. For homogeneous Dirichlet data, the corresponding weak formulation is naturally posed on H01 (Ω), where the elliptic operator is instead understood as a map into H −1 (Ω) = (H01 (Ω))′ . 2.2. Parametric models. Many approximation methods can be formulated in terms of a finite number of degrees of freedom. This motivates the introduction of an abstract notion that encompasses both classical discretization techniques and modern machine learning models. Let p ∈ N. A parametric model is a differentiable map (2.6)

v : Rp → H,

θ 7→ vθ .

1More generally, one could consider an abstract Hilbert space H of functions from Ω to R, but the present setting

is sufficient for the purposes of this article.

4

N. SCHWENCKE, R. MAIER

The vector θ ∈ Rp collects the degrees of freedom of the approximation. As the parameters vary, the model generates a family of functions (2.7)

Mv := {vθ | θ ∈ Rp } ⊂ H,

which may be viewed, under mild assumptions, as a finite-dimensional manifold embedded in H. Understanding how the model can evolve therefore amounts to studying the local geometry of Mv . The differential of the parametrization at θ is then given by (2.8)

dvθ : Rp → H,

ξ 7→

p X

ξi ∂i vθ .

i=1

Its image defines the tangent space of the manifold at vθ , given by (2.9)

Tθ Mv := Im(dvθ ) = span{∂i vθ | i = 1, . . . , p}.

The tangent space describes the admissible infinitesimal variations of the model around vθ . Before considering particular examples, it is useful to introduce one further construction. Let v : Rp → H be a parametric model and let F : H → Y be a differentiable map between functional spaces. We call the composition (2.10)

ΓF := F ◦ v : Rp → Y,

ΓF θ = F(vθ ),

a compound parametric model. Note that ΓF is itself a parametric model, now taking values in Y. Its admissible first-order variations are described by
F Tθ MΓF = Im dΓF θ = span{∂i Γθ | 1 ≤ i ≤ p} ⊂ Y. Two important examples of parametric models are finite element approximations and neural networks, which we introduce next. Compound models will in particular arise naturally in the PINN formulation below. 2.3. Galerkin and Petrov–Galerkin methods. To simplify the presentation, within this subsection we stick to the exemplary model problem of Example 2.1 with gbc ≡ 0. We emphasize, however, that similar constructions are possible for much more general problems. As a first step, we transform the PDE (2.4) into a weak formulation by multiplying with an appropriate test function w ∈ H01 (Ω) and integrating by parts. We then seek a function u ∈ H01 (Ω) such that (2.11)

for all w ∈ H01 (Ω).

⟨A∇u, ∇w⟩L2 (Ω) = ⟨f, w⟩L2 (Ω) ,

For convenience, we introduce the energy inner product (2.12)

a(u, w) := ⟨A∇u, ∇w⟩L2 (Ω) .

The weak formulation (2.11) can then be written as (2.13)

a(u, w) = ⟨f, w⟩L2 (Ω) ,

for all w ∈ H01 (Ω).

An advantage of the weak formulation is that it reduces the regularity assumptions on the solution: second derivatives are not required and the first derivative needs to exist in a weak sense only. Nonetheless, (2.11) cannot be solved directly, as it is still posed in the infinite-dimensional space H01 (Ω). Therefore, we solve (2.11) in a finite-dimensional subspace Vm ⊂ H01 (Ω) with dim Vm = m, that is, we seek um ∈ Vm such that (2.14)

a(um , wm ) = ⟨f, wm ⟩L2 (Ω) ,

This approach is called the Galerkin method.

for all wm ∈ Vm .

BEYOND PINNS: A UNIFIED FRAMEWORK FOR NEURAL AND HYBRID PDE SOLVERS

5

Since Vm is finite-dimensional, the Galerkin approximation can itself be viewed as a particular instance of a parametric model. Let (φi )1≤i≤m be a basis of Vm and consider the linear parametrization (2.15)

c ∈ Rm 7−→ vc :=

m X

ci φi ∈ Vm .

i=1

Its image is exactly Vm , and, since the parametrization is linear, its tangent space is independent of the coefficient vector. That is, (2.16)

Tc Mv = Vm

for all c ∈ Rm .

Thus, Galerkin methods fit naturally into the parametric-model framework: they correspond to a linear, finite-dimensional model whose tangent space coincides everywhere with the approximation space itself. Writing the Galerkin solution as (2.17)

um =

m X

ci φi ,

i=1

the variational problem (2.14) therefore reduces to a finite-dimensional linear system for the pam rameter vector c = (ci )m i=1 ∈ R . Typical examples for Vm are so-called finite element spaces that consist of piece-wise polynomial and continuous functions, see, e.g., the works [7, 4, 5] on finite element methods. More generally, one may choose a test space Wq ⊂ H01 (Ω) distinct from the approximation space Vm . The associated Petrov–Galerkin method seeks um ∈ Vm such that (2.18)

a(um , wq ) = ⟨f, wq ⟩L2 (Ω)

for all wq ∈ Wq ,

where the Galerkin method is recovered for the choice Wq = Vm . Note that this formulation allows for dim Vm = m ̸= q = dim Wq . Classical square Petrov–Galerkin discretizations usually take m = q; when m ̸= q, the discrete system is rectangular and can instead be interpreted in a least-squares or minimum-norm sense when appropriate. This will become important in Section 3. 2.4. Neural networks. Before introducing the concept of so-called physics-informed neural networks, we first provide a functional-analytic definition of a neural network that fits within the abstract framework of a parametric model, which has been introduced in Section 2.2 above. Definition 2.2 (multi-layer perceptron architecture). Based on the hyperparameters • L ∈ N, called the depth of the network, • (li )0≤i≤L ∈ NL+1 , called the widths of the layers, QL 1,∞ • (χi )1≤i≤L ∈ i=1 Wloc (Rli ; Rli ), called the activation functions, PL • p ∈ N with ≤ p ≤ i=1 li (li−1 + 1), called the number of parameters,  1Q QL L 1 p • ψ ∈ C R ; i=1 Rli ×li−1 × i=1 Rli , called the parameterization map, a multi-layer perceptron (MLP) architecture is a mapping v : Rp → C 0 (Ω; RlL ) for some subdomain Ω ⊂ Rl0 , defined by
(2.19) θ 7→ ◦L k=1 χk ψk (θ) • +ψL+k (θ) , where ψk (θ) ∈ Rlk ×lk−1 and ψL+k (θ) ∈ Rlk denote the weight matrix and bias vector of the kth layer. In view of Definition 2.2, we emphasize that activation functions are typically scalar-valued functions that are applied component-wise. That is, for a given k ∈ N and an activation function χ : Rk → Rk , there exists σ : R → R such that for x ∈ Rk , we have
χ(x) = σ(xi ) 1≤i≤k ∈ Rk .

6

N. SCHWENCKE, R. MAIER

Therefore, it is common to abuse the notation and to refer to χ through its scalar component function σ. Furthermore, we note that the parameterization map ψ encodes both structural constraints on the weights and biases, including parameter sharing as in convolutional or tied-weight architectures, as well as the distinction between learnable and fixed parameters. In this paper, we will mainly view ψ as a convenient re-parameterization of a real vector of dimension (2.20)

p=

L X

li (li−1 + 1)

i=1

QL QL li ×li−1 li to the corresponding collection of (Wi )L and (bi )L i=1 ∈ i=1 ∈ i=1 R i=1 R . The MLP architecture of Definition 2.2 therefore provides a particular instance of the parametric-model framework introduced in Section 2.2. Indeed, once the architecture (that is, the depth, layer widths, activation functions, and parameterization map ψ) has been fixed, each parameter vector θ ∈ Rp determines a unique realized function vθ . The network can thus be viewed as a map v : Rp → H,

θ 7→ vθ ,

provided that the realizations belong to the chosen Hilbert space H. In contrast with the linear parametrization associated with a Galerkin space, the dependence of vθ on θ is generally nonlinear because of the successive compositions with the activation functions. Consequently, the set of realizations Mv = {vθ | θ ∈ Rp } ⊂ H is in general a nonlinear parametric manifold. Its tangent space at vθ is generated by the parameter derivatives, Tθ Mv = span{∂i vθ | 1 ≤ i ≤ p}, and therefore generally depends on the current parameter value θ. This is the main distinction with the Galerkin parametrization above, whose tangent space is the fixed approximation space Vm at every parameter value. Remark 2.3 (Other architectures). Note that we restrict our presentation to MLPs. Many other neural network architectures have been proposed, such as convolutional neural networks [25, 14], recurrent neural networks, variants thereof [12, 17, 6], and transformers [42]. From the perspective adopted in this work, these architectures can often be viewed as particular parameterizations obtained by imposing structural constraints on a general neural network model. Since our analysis only relies on the resulting parametric map (2.19), focusing on MLPs entails no significant loss of generality. 2.5. Physics-informed neural networks. Physics-informed neural networks (PINNs) provide a natural example of the compound parametric models introduced in Section 2.2. For the boundary value problem (2.3), we introduce Y := L2 (Ω) × L2 (∂Ω) and define (2.21)

ΓD,B := (D, B) ◦ v : Rp → Y,

ΓD,B = D[vθ ], B[vθ ] . θ

With the target
y D,B := f, gbc , the associated functional residual is rθD,B := ΓD,B − y D,B . θ The PINN formulation consists of minimizing its squared norm, ˆ ˆ 2 (2.22) ℓ(θ) = rθD,B = (D[vθ ](x) − f (x))2 dx + (B[vθ ](s) − gbc (s))2 ds. Y

Ω

∂Ω

BEYOND PINNS: A UNIFIED FRAMEWORK FOR NEURAL AND HYBRID PDE SOLVERS

7

Ω In practice, this functional residual is discretized by pointwise evaluations. Let XΩ = {xi }N i=1 ⊂ Ω N∂Ω and X∂Ω = {si }i=1 ⊂ ∂Ω be families of interior and boundary collocation points. The resulting PINN loss is given by

(2.23)

ℓpinn (θ) =

N∂Ω NΩ 1 X 1 X (D[vθ ](xi ) − f (xi ))2 + (B[vθ ](si ) − gbc (si ))2 . NΩ i=1 N∂Ω i=1

Thus, standard collocation replaces the functional residual by a finite collection of pointwise measurements. These pointwise evaluations must of course be well defined. In particular, the neuralnetwork realization must possess the regularity required by the operators D and B. This typically requires sufficiently smooth activation functions. Automatic differentiation [26, 2] is then used to evaluate these quantities efficiently [37]. We emphasize that this pointwise requirement is stronger than merely assuming rθD,B ∈ Y; this distinction will become important when more general measurements of the functional residual are introduced below. For the specific problem (2.4) in Example 2.1 with gbc ≡ 0, the discrete loss becomes (2.24)

N∂Ω NΩ 1 X 1 X 2 (− div(A(xi )∇vθ (xi )) − f (xi )) + (γvθ (si ))2 , ℓpinn (θ) = NΩ i=1 N∂Ω i=1

where γ denotes the trace operator as before. 3. Connecting Gauss–Newton with Petrov–Galerkin discretizations We now introduce the second main ingredient of this work: the Gauss–Newton (GN) method and its interpretation through Petrov–Galerkin discretizations. We first recall the classical finitedimensional construction. We then show that a functional problem can be reduced to exactly the same setting by choosing a finite family of measurements of the functional residual. Different choices of approximation models and measurements recover, among others, natural-gradient methods, kernel methods, regularization and sketching strategies, both for PINNs and for direct regression problems. 3.1. Gauss–Newton in finite-dimensional optimization. Let ϱ : Rp → RN ,

θ 7→ ϱθ ,

N

be a differentiable discrete model and let y ∈ R be a target. We consider the nonlinear leastsquares problem 1 2 (3.1) ℓ(θ) = ∥ϱθ − y∥RN . 2 We denote its residual by bθ := ϱθ − y and the Jacobian of ϱ at θ by Jθ := dϱθ ∈ RN ×p . p Let θ̄ = θ − ξ, where ξ ∈ R denotes the correction to be subtracted. The first-order approximation (3.2)

bθ−ξ = bθ − Jθ ξ + o(∥ξ∥)

leads to the local least-squares problem (3.3)

ξθGN ∈ arg min ξ∈Rp

1 2 ∥bθ − Jθ ξ∥RN . 2

This formulation makes the basic idea of Gauss–Newton transparent. The residual bθ is the correction that would ideally bring ϱθ to the target, whereas parameter variations can only produce, as first-order approximations, corrections in Im(Jθ ) = span{Jθ ej | 1 ≤ j ≤ p}.

8

N. SCHWENCKE, R. MAIER

Gauss–Newton therefore selects the accessible correction that best approximates bθ . Its optimality conditions are (3.4)

Jθ ej , bθ − Jθ ξθGN RN = 0,

1 ≤ j ≤ p,

or, equivalently, Jθ⊤ Jθ ξθGN = Jθ⊤ bθ . Thus, Jθ ξθGN = ΠIm(Jθ ) bθ ,

(3.5)

where ΠIm(Jθ ) denotes the Euclidean orthogonal projection onto the range of the Jacobian Jθ . With the Moore–Penrose convention and denoting the pseudo-inverse with (•)† , the minimumnorm parameter correction is ξθGN = Jθ† bθ .

(3.6)

3.2. Discretized Gauss–Newton for functional models. We now consider a parametric model as introduced in Section 2.2, v : Rp → H, θ 7→ vθ , where H is a Hilbert space, together with a target y ∈ H. We denote the corresponding functional residual by rθ := vθ − y.

(3.7) For a parameter correction ξ ∈ Rp ,

rθ−ξ = rθ − dvθ [ξ] + o(∥ξ∥). At the linearized level, the ideal correction would therefore satisfy (3.8)

dvθ [ξ] = rθ .

The admissible corrections on the left-hand side belong to the finite-dimensional tangent space Tθ Mv = Im(dvθ ) = span{∂j vθ | 1 ≤ j ≤ p} ⊂ H. Thus, although (3.8) is an equation between functional objects, the space of admissible first-order corrections is finite-dimensional. To discretize this equation, we fix the current parameter θ and introduce the extended tangent space (3.9)

Tθ Mv := span{rθ } + Tθ Mv ⊂ H.

Its dimension is at most p + 1. Let Λθ = (λθi )N i=1 ,

′ λθi ∈ Tθ Mv ,

be a family of linear measurements on Tθ Mv . Applying these measurements to (3.8) gives the finite-dimensional system (3.10)

θ JθΛθ ξ = bΛ θ ,

θ θ bΛ θ,i := λi (rθ ),

Λθ := λθi (∂j vθ ). Jθ,ij

The linearized problem in the function space H has therefore been reduced to an ordinary finitedimensional problem of precisely the form considered in Section 3.1. The corresponding discretized Gauss–Newton problem is defined by 2 1 Λθ (3.11) ξθΛθ ∈ arg min bθ − JθΛθ ξ N . R ξ∈Rp 2 As before, the minimum-norm correction is (3.12)

† θ ξθΛθ = JθΛθ bΛ θ .

BEYOND PINNS: A UNIFIED FRAMEWORK FOR NEURAL AND HYBRID PDE SOLVERS

9

This construction has an immediate Petrov–Galerkin interpretation. More generally, suppose that the measurements admit a representation through a pairing with a test space Z, namely, for 1 ≤ i ≤ N, (3.13)

λθi (w) = ⟨w , zθi ⟩H×Z ,

for all w ∈ Tθ Mv ,

for some test function zθi ∈ Z. Consequently, (3.10) can be written equivalently as (3.14)

⟨dvθ [ξ] − rθ , zθi ⟩H×Z = 0,

1 ≤ i ≤ N.

Thus, at a fixed parameter θ ∈ Rp , choosing measurements of the linearized problem in H is equivalent, whenever such a representation is available, to choosing the Petrov–Galerkin test functions against which this problem is enforced. In particular, when Z = H and the pairing is the inner product of H, the finite dimensionality of Tθ Mv implies that every linear measurement is continuous and admits a unique Riesz representer ζiθ ∈ Tθ Mv , so that one may take zθi = ζiθ . When the resulting system is rectangular or rank deficient, Gauss–Newton provides the corresponding least-squares or minimum-norm solution through (3.11). Finally, when the measurements extend to a common linear subspace D ⊂ H containing the image of v and the target y, and which are independent of θ, they define the observation map
N EΛ : D → RN , w 7→ λi (w) i=1 , and hence a global discrete model (3.15)

ϱΛ := EΛ ◦ v : Rp → RN ,

y Λ := EΛ (y).

The construction above is then exactly the ordinary Gauss–Newton linearization of the discrete least-squares problem associated with ϱΛ and the target y Λ . The local formulation is more general, though, since the measurements only need to be defined on the current extended tangent space. 3.3. Recovering weak formulations. We now show how classical weak formulations arise naturally within the functional Gauss–Newton framework through suitable choices of linear measurements and associated test functions. In particular, operator-adapted Green sections and general weak test functions allow the functional residual equation to be discretized without requiring a pointwise strong-form residual. More broadly, the same framework recovers a variety of classical optimization and discretization methods by varying the parametric approximation model and the measurement family. Examples including natural gradient, kernel methods, empirical natural gradient, pointwise Gauss–Newton, sketching, and regularization are briefly discussed in Section A. 3.3.1. Green sections. We consider the elliptic problem of Example 2.1 with homogeneous Dirichlet boundary conditions, and denote with LA : H01 (Ω) → H −1 (Ω)

(3.16) the weak elliptic operator defined by

⟨LA w , z⟩H −1 (Ω)×H01 (Ω) = a(w, z). Given a functional parametric model v : Rp → H01 (Ω), we consider the compound model (3.17)

ΓA := LA ◦ v : Rp → H −1 (Ω)

with residual rθA = LA vθ − f = LA (vθ − u),

10

N. SCHWENCKE, R. MAIER

where u denotes the exact solution to LA u = f . At the current parameter θ, we consider the local solution space Tθ Mv := span{vθ − u} + Tθ Mv ⊂ H01 (Ω). Assume that point evaluation at x ∈ Ω is well defined on this space. Since Tθ Mv is finitedimensional, point evaluation is then continuous and admits a corresponding local Green section gθ,x ∈ Tθ Mv satisfying (3.18)

a(w, gθ,x ) = w(x),

for all w ∈ Tθ Mv .

Choosing gθ,xi as test functions in the Petrov–Galerkin formulation gives A ⟨dΓA θ [ξ] − rθ , gθ,xi ⟩H −1 (Ω)×H01 (Ω) = 0.

Since dΓA θ [ξ] = LA dvθ [ξ], the reproducing property yields dvθ [ξ](xi ) = vθ (xi ) − u(xi ),

(3.19)

1 ≤ i ≤ N.

Thus, Green sections transform weak measurements of the PDE residual into pointwise measurements of the desired correction in the solution space. Importantly, the differential operator is transferred from the approximation model to the test functions, so that it does not need to be evaluated directly on the model. 3.3.2. General weak test functions. Green sections are only one particular choice of weak test functions. We consider again the compound model ΓA := LA ◦ v : Rp → H −1 (Ω), with functional residual rθA := ΓA θ − f = LA vθ − f,

f ∈ H −1 (Ω).

1 Let (zi )N i=1 be arbitrary functions in H0 (Ω). At the current parameter θ, they define the local measurements λθi (r) := ⟨r , zi ⟩H −1 (Ω)×H01 (Ω)

on the extended tangent space of ΓA . The discretized functional Gauss–Newton equation therefore becomes A ⟨dΓA θ [ξ] − rθ , zi ⟩H −1 (Ω)×H01 (Ω) = 0,

(3.20)

1 ≤ i ≤ N.

Using the definition of LA , this is equivalent to a(dvθ [ξ], zi ) = a(vθ , zi ) − ⟨f , zi ⟩H −1 (Ω)×H01 (Ω) ,

(3.21) 2

1 ≤ i ≤ N. 2

When f ∈ L (Ω), the last duality pairing reduces to the usual L (Ω) inner product. The corresponding finite residual and Jacobian components in (3.10) are therefore (3.22)

bZ θ,i := a(vθ , zi ) − ⟨f , zi ⟩H −1 (Ω)×H01 (Ω) ,

and (3.23)

Z := a(∂j vθ , zi ). Jθ,ij

The parameter correction is then obtained by applying ordinary finite-dimensional Gauss–Newton to these weak residual components. This formulation makes explicit the regularity advantage of weak measurements. Although the compound residual takes values in H −1 (Ω), it is tested against functions in H01 (Ω). For a secondorder elliptic operator, the bilinear form a(·, ·) involves only weak spatial derivatives of first order, so that the strong differential operator never needs to be explicitly evaluated in the approximation model. In this sense, derivatives are transferred to the weak form and the test functions, allowing residuals of substantially weaker regularity than in a pointwise strong-form discretization.

BEYOND PINNS: A UNIFIED FRAMEWORK FOR NEURAL AND HYBRID PDE SOLVERS

11

The Green-section construction of the previous paragraph is recovered as a particular operatoradapted choice of these weak test functions. By the Lax–Milgram theorem, the weak elliptic operator LA : H01 (Ω) → H −1 (Ω) is an isomorphism and therefore admits a bounded inverse. Hence, for any λ ∈ H −1 (Ω), we may define its Green lift 1 gλ := L−1 A [λ] ∈ H0 (Ω).

By construction, for all w ∈ H01 (Ω).

a(gλ , w) = ⟨λ , w⟩H −1 (Ω)×H01 (Ω)

Consequently, choosing zi = gλi in (3.20) transforms the weak residual measurements into the corresponding measurements of the desired correction in the solution space, namely

(3.24) λi dvθ [ξ] = λi vθ − u , 1 ≤ i ≤ N. Arbitrary weak test functions therefore extend the same principle beyond point evaluations while retaining the low-regularity H −1 formulation. The test functions can consequently be selected according to the variational structure, regularity, or computational properties of the problem. In the next section, we combine this freedom with a hybrid finite element–neural approximation model to obtain a hybrid weak-test method. 4. Hybrid finite element–neural method The preceding framework shows that optimization and discretization methods can be constructed by varying the approximation model and the test functions in the common Petrov– Galerkin formulation (3.14). We now exploit this freedom to construct a hybrid finite element– neural method in which a finite element component accounts exactly for a prescribed finitedimensional space, while the neural component is restricted to its a-orthogonal complement. 4.1. Projected hybrid approximation model. Let Vn = span{φi }ni=1 ⊂ H01 (Ω) be a finite-dimensional approximation space. In the following, we primarily consider finite element spaces of moderate dimension, although the construction does not depend on this particular choice. We denote by Pan : H01 (Ω) → Vn the a-orthogonal projection onto Vn , and set Qan := I − Pan . Hence,  Im(Qan ) = Vn⊥,a := v ∈ H01 (Ω) | a(v, wn ) = 0 for all wn ∈ Vn . Let un ∈ Vn denote the Galerkin approximation associated with Vn , i.e., (4.1)

a(un , wn ) = ⟨f , wn ⟩H −1 (Ω)×H01 (Ω)

for all wn ∈ Vn .

Starting from the functional neural model v : Rp → H01 (Ω),

θ 7→ vθ ,

we define the projected hybrid model (4.2)

vnhyb : Rp → H01 (Ω),

θ 7→ vθ,n := un + Qan vθ .

Thus, the finite element component accounts for the directions in Vn , while the trainable contribution is restricted to Vn⊥,a . Equivalently, (4.3)

vθ,n = vθ + zn,θ ,

zn,θ := un − Pan vθ ∈ Vn .

12

N. SCHWENCKE, R. MAIER

The finite element term zn,θ therefore acts as a compensator, rather than as an additional set of trainable parameters. By construction, a(vθ,n , wn ) = ⟨f , wn ⟩H −1 (Ω)×H01 (Ω)

(4.4)

for all wn ∈ Vn .

For each θ, computing the compensator amounts to solving a finite-dimensional Galerkin system whose matrix depends only on Vn and a, and can therefore be factorized once and reused throughout training. The differential of the hybrid model with respect to θ is
(4.5) d vnhyb θ = Qan dvθ , so that its tangent directions are themselves restricted to Vn⊥,a . For the elliptic problem, we introduce the associated compound model := LA ◦ vnhyb : Rp → H −1 (Ω), ΓA,hyb n

(4.6) with residual

A := ΓA,hyb rθ,n − f = LA vθ,n − f. θ,n

(4.7)

Equation (4.4) immediately implies A ⟨rθ,n , wn ⟩H −1 (Ω)×H01 (Ω) = 0

(4.8)

for all wn ∈ Vn .

The weak equation is therefore satisfied exactly on Vn , independently of the neural parameters. 1 4.2. Projected Petrov–Galerkin update. Let (zi )N i=1 be arbitrary test functions in H0 (Ω). A,hyb Inserting the compound model Γn into (3.14) gives
A (4.9) ⟨d ΓA,hyb [ξ] − rθ,n , zi ⟩H −1 (Ω)×H01 (Ω) = 0, 1 ≤ i ≤ N, n θ

where ξ ∈ Rp denotes, as before, the correction subtracted from the current parameter. Both the residual and the tangent directions vanish when tested against Vn . Consequently, only the complementary components of the test functions matter, and we may replace each zi by ezi := Qan zi ∈ Vn⊥,a .

(4.10) Since ezi ∈ Vn⊥,a ,

a(Qan dvθ [ξ],ezi ) = a(dvθ [ξ],ezi ),

and, since un ∈ Vn , A ⟨rθ,n , ezi ⟩H −1 (Ω)×H01 (Ω) = a(vθ ,ezi ) − ⟨f , ezi ⟩H −1 (Ω)×H01 (Ω) .

Therefore, (4.9) reduces to (4.11)

a(dvθ [ξ],ezi ) = a(vθ ,ezi ) − ⟨f , ezi ⟩H −1 (Ω)×H01 (Ω) ,

1 ≤ i ≤ N.

Thus, the proposed method is again precisely of the form of (3.14), now with the approximation model ΓA,hyb and projected test functions ezi = Qan zi . n In the notation of (3.10), the corresponding finite residual and Jacobian are (4.12)

bnθ,i := a(vθ ,ezi ) − ⟨f , ezi ⟩H −1 (Ω)×H01 (Ω) ,

and (4.13)

n := a(∂j vθ ,ezi ). Jθ,ij

The parameter correction is therefore obtained from a similar finite-dimensional Gauss–Newton problem as before and reads (4.14)

Jθn ξ = bnθ ,

understood in the least-squares sense when the system is rectangular or rank deficient. The whole construction can be summarized as follows. The finite element component is treated exactly in Vn , while Qan = I − Pan removes the corresponding finite element directions from both

BEYOND PINNS: A UNIFIED FRAMEWORK FOR NEURAL AND HYBRID PDE SOLVERS

13

the neural correction and the test functions. Consequently, the remaining optimization problem is posed entirely on the a-orthogonal complement Vn⊥,a . 4.3. Numerical realization. The hybrid formulation does not require an alternating optimization between finite element and neural variables. The Galerkin solution un and a factorization of the finite element stiffness matrix are computed once. For a given parameter θ, evaluating the hybrid model only requires computing the projection Pan vθ , while projecting a test function through Qan uses the same prefactorized finite-dimensional system. At each neural iteration, a family of test functions is selected and projected onto Vn⊥,a . The residual vector bnθ and Jacobian Jθn are then assembled from (4.12) and (4.13), respectively, defining an ordinary finite-dimensional Gauss–Newton problem. After updating the neural parameters, the new approximation is simply evaluated through vθ,n = un + Qan vθ . Recomputing the projection at the updated parameter is therefore part of evaluating the hybrid approximation model itself, and should not be interpreted as a block-alternating optimization step. 5. Numerical experiments The experiments are organized around two questions. First, we verify that the functional discretization developed above allows Gauss–Newton algorithms designed for finite residual vectors to be applied directly to weak PDE formulations. Second, we examine whether this construction remains effective beyond Green sections and smooth solutions, and whether the hybrid finite element–neural model improves the robustness of weak neural solvers. Throughout the section, we monitor both relative L2 and full relative H 1 errors. Complementary implementation details, optimizer and activation ablations, additional figures, and complete method-level results are provided in the companion blog ↗ . 5.1. Gauss–Newton optimization of weak residuals. 5.1.1. Benchmark problems and weak-residual discretization. We first consider two smooth onedimensional problems on Ω = (−1, 1) with homogeneous Neumann boundary conditions. Both may be written as (5.1)

−u′′ + ρ(u) = f,

u′ (−1) = u′ (1) = 0,

with ρ(s) = s

or

ρ(s) = s3 ,

corresponding respectively to the linear and nonlinear benchmarks. In both cases, the solution is u(x) = cos(πx), so that
f (x) = π 2 cos(πx) + ρ cos(πx) . The natural energy space is H 1 (−1, 1), equipped with ˆ 1
a(v, w) = v ′ (x)w′ (x) + v(x)w(x) dx. −1

For t ∈ [−1, 1], let gt denote the Riesz representer of point evaluation at t with respect to a, given by (5.2)

gt (x) =

cosh(min{x, t} + 1) cosh(1 − max{x, t}) . sinh(2)

14

N. SCHWENCKE, R. MAIER

For t ∈ (−1, 1), this is the Neumann Green section associated with −∂xx + 1. Testing the residual against these functions and integrating by parts, while retaining the Neumann conditions as separate residual coordinates, gives (5.3)

rtρ (v) = v(t) + v ′ (−1)gt (−1) − v ′ (1)gt (1) ˆ 1 ˆ 1
f (x)gt (x) dx. ρ(v(x)) − v(x) gt (x) dx − + −1

−1

For the linear problem, the first integral in the second line vanishes, whereas for the nonlinear problem it contains v 3 − v. We discretize the weak PDE residual with 200 fixed sections gt centered at equispaced points of [−1, 1], including both endpoints, and append the two Neumann residuals v ′ (−1) and v ′ (1), resulting in a residual vector with 202 scalar components. For interior centers t, gt′ has a jump at x = t. The integrals in (5.3) are therefore split at t and evaluated with a 64-point Gauss–Legendre rule on each subinterval. These quadrature points are used only to evaluate the corresponding linear measurements entering the residual vector and its Jacobian; they do not define additional measurements and therefore do not increase the dimension of either. 5.1.2. Energy references. We compare ordinary Gauss–Newton (GN) solvers applied to this weak residual with specialized energy-based methods. For the linear problem, the reference is the Gauss– Newton Deep Ritz method of [16], associated with ˆ ˆ 1
1 1 ′ |v (x)|2 + v(x)2 dx − f (x)v(x) dx. (5.4) Elin (v) = 2 −1 −1 Following the reference implementation, its integrals are evaluated by two-point Gauss–Legendre quadrature on uniform cells with htrain = 1/3000 and htest = 1/4000, corresponding to 12000 and 16000 quadrature nodes, respectively, and we use the same geometric step-length search. For the nonlinear problem, the reference is the Energy Natural Gradient (NG) method of [32], based on ˆ ˆ ˆ 1 1 1 ′ 1 1 (5.5) Enl (v) = |v (x)|2 dx + v(x)4 dx − f (x)v(x) dx, 2 −1 4 −1 −1 whose second variation, (5.6)

ˆ 1

ˆ 1 w′ (x)z ′ (x) dx + 3

D2 Enl (v)[w, z] = −1

v(x)2 w(x)z(x) dx, −1

defines the corresponding state-dependent natural-gradient metric. We report both the publishedprotocol reference, using a width-32 tanh network and trapezoidal integration with 20000 training and 200000 evaluation points, and a matched version using the same metric with the controlled width-64 protocol. 5.1.3. Residual Gauss–Newton solvers. Let rθ ∈ RN denote a strong or weak residual vector and let Jθ = U ΣV ⊤ , Σ = diag(σ1 , . . . , σq ), q = min{N, p}, be a thin singular-value decomposition of its Jacobian, with singular values in non-increasing order. With the convention that the parameter update is θ̄ = θ − ξ, we consider the cutoff and ridge Gauss–Newton corrections   1{σi ≥τ } (5.7) U ⊤ rθ , ξcut = V diag σi   σi (5.8) ξridge = V diag U ⊤ rθ . σi2 + λ

BEYOND PINNS: A UNIFIED FRAMEWORK FOR NEURAL AND HYBRID PDE SOLVERS

15

Table 1. Median final relative errors and optimization wall-clock times over ten seeds for the smooth tanh benchmarks. Energy NG (baseline) follows the published protocol of [32], whereas Energy NG (matched) uses the controlled width64 protocol. Wall-clock times reflect the stated protocol of each method and are therefore not all intended as matched computational-cost comparisons. Problem

Method

rel. L2

rel. H 1

time [s]

1.12 × 10 5.61 × 10−8 6.16 × 10−8 3.22 × 10−15 6.70 × 10−15

−7

Linear

GN Deep Ritz (baseline) ridge GN strong ridge GN weak DSGNAR strong DSGNAR weak

−8

1.04 × 10 3.89 × 10−8 5.15 × 10−7 1.36 × 10−15 4.33 × 10−14

127.3 140.5 76.9 24.6 21.0

Nonlinear

Energy NG (baseline) Energy NG (matched) ridge GN strong ridge GN weak DSGNAR strong DSGNAR weak

1.27 × 10−8 1.74 × 10−8 1.39 × 10−8 2.07 × 10−8 6.30 × 10−15 3.31 × 10−15

1.23 × 10−7 1.77 × 10−7 1.37 × 10−8 1.31 × 10−7 5.46 × 10−15 2.81 × 10−14

3.9 133.6 143.8 91.8 28.5 24.3

For ridge Gauss–Newton, the regularization parameter is held fixed throughout each run. We use λ = 10−12 for the strong residuals, and λ = 10−13 and 10−14 for the weak residuals in the linear and nonlinear problems, respectively. We also apply AMStramGRAM [40] and DSGNAR [43] directly to the same finite residual vectors. No modification of these algorithms is required for the weak formulation: once the test functionals have been evaluated, the optimizer receives an ordinary residual vector together with its parameter Jacobian. For the smooth benchmarks, we additionally consider the corresponding strong-residual formulations as controls. Since the MLP trial functions use smooth tanh activations, their pointwise second-order derivatives are well defined, allowing us to distinguish the effect of the nonlinear least-squares geometry from that of weak testing. 5.1.4. Experimental protocol and results. All matched experiments use the same width-64 tanh architecture and are repeated over ten seeds, with identical initial neural parameters across methods for each seed. The strong formulation uses 200 pointwise PDE residuals and the weak formulation 200 Green-section residuals, evaluated at the same equispaced points of [−1, 1]; the two Neumann boundary residuals are appended in both cases. Relative L2 and full relative H 1 errors with respect to the solution u are evaluated using an independent high-resolution quadrature. We report medians over the ten seeds, with interquartile bands in the convergence plots; reported wall-clock times exclude compilation and warm-up. The principal weak-residual comparisons use ridge Gauss–Newton and DSGNAR. Final errors and wall-clock times for the corresponding strong-residual controls are also reported below. Additional results for cutoff Gauss–Newton and AMStramGRAM, together with complete convergence curves and further optimizer diagnostics, are provided in the companion blog ↗ . The linear benchmark exhibits essentially the same qualitative behavior as the nonlinear one; its complete convergence curves are provided in the companion blog ↗ . We therefore show only the nonlinear convergence curves in Figure 5.1, while Table 1 reports the final results for both problems. Table 1 shows first that no specialized energy metric is required to optimize the weak formulation. Ridge Gauss–Newton applied directly to the weak residual reaches an accuracy comparable to the dedicated energy methods [16, 32], whereas DSGNAR drives both weak problems close to machine precision.

16

N. SCHWENCKE, R. MAIER

100 10 2

10 3

relative H 1 error

relative L 2 error

10 1

10 5 10 7 10 9 10 11

10 4 10 6 10 8 10 10 10 12

10 13 0

200

400

600

800

iteration

1000

0

200

400

600

800

iteration

1000

100 10 2

10 3

relative H 1 error

relative L 2 error

10 1

10 5 10 7 10 9 10 11

10 4 10 6 10 8 10 10 10 12

10 13 0

20

40

60

80

100

optimization wall-clock [s]

120

Energy NG (baseline) Energy NG (matched)

140

Ridge GN (strong) Ridge GN (weak)

0

20

40

60

80

100

optimization wall-clock [s]

120

140

DSGNAR (strong) DSGNAR (weak)

Figure 5.1. Nonlinear smooth benchmark with tanh networks. Median relative L2 and full relative H 1 errors over ten seeds, with interquartile bands, against iterations and optimization wall-clock. Energy NG (baseline) follows the protocol of [32], whereas Energy NG (matched) uses the controlled width-64 protocol. The strong-residual controls reach comparable or better accuracies, with DSGNAR reaching near-machine precision in both formulations. Hence, on these smooth problems, the results indicate that the dominant improvement is associated with the residual least-squares formulation and its Gauss–Newton treatment, rather than with weak testing itself. The significance of the weak formulation is that this same optimization structure survives after integration by parts, without requiring the strong differential operator to be evaluated pointwise on the approximation model. The nonlinear Energy-NG baseline additionally exhibits substantial seed dependence, as visible in Figure 5.1. Six of the ten runs reach the high-accuracy regime, with final relative H 1 errors between approximately 8 × 10−8 and 1.3 × 10−7 , whereas four remain close to unit relative error. Its median is therefore accurate while its interquartile band remains broad. Under the matched protocol, all ten runs reach the accurate branch. 5.2. General weak formulations and hybrid finite element–neural models. We next investigate whether the framework remains effective when operator-adapted Green sections are unavailable or suboptimal, and when the right-hand side or solution has limited regularity. We consider four benchmarks.

BEYOND PINNS: A UNIFIED FRAMEWORK FOR NEURAL AND HYBRID PDE SOLVERS

17

Multiscale diffusion (MS). On Ω = (0, 1), we consider the one-dimensional instance of the elliptic problem of Example 2.1,   d d (5.9) − Aε (x) u(x) = 1, u(0) = u(1) = 0, dx dx with an oscillating diffusion coefficient (5.10)

Aε (x) =

1 , 2 + cos(2πx/ε)

ε=

1 . 16

The corresponding solution involves ε-dependent oscillations; see also [27, Ch. 2]. Jump forcing (JF). On Ω = (0, 1)2 , we consider the homogeneous Dirichlet Poisson problem with manufactured solution  2   x − 3x , x < 1 , 2 2 8 u(x, y) = p(x)y(1 − y), p(x) = x − 1   , x ≥ 12 . 8 The corresponding right-hand side f = −∆u is discontinuous across x = 1/2, while the solution has no distributional singularity. Line source (LS). Again on Ω = (0, 1)2 , we prescribe u(x, y) = min(x, 1 − x)y(1 − y) +

(5.11)

1 sin(2πx) sin(2πy), 2

which solves (5.12)

−∆u = freg + 2y(1 − y) δ{x=1/2} ,

freg = 2 min(x, 1 − x) + 4π 2 sin(2πx) sin(2πy).

Here

 ˆ 1  1 , y dy, δ{x=1/2} , z = z 2 0 so the singular source is treated directly as a one-dimensional functional and is never replaced by a regularized volumetric source. Reentrant corner (RC). Finally, on the L-shaped domain
(5.13) ΩL = (−1, 1)2 \ [0, 1) × (−1, 0] , we consider the homogeneous Dirichlet Poisson problem with manufactured solution   2ϑ (5.14) u(x, y) = (1 − x2 )(1 − y 2 )r2/3 sin , 3 where (r, ϑ) are polar coordinates centered at the reentrant corner. The factor r2/3 sin(2ϑ/3) gives the classical corner singularity, with |∇u| ∼ r−1/3 near the origin. 5.2.1. Trial and test spaces. For MS, the pure neural model contains 1255 trainable parameters and the hybrid model combines a 375-parameter neural component with 15 interior P1 finite element degrees of freedom, giving a total approximation dimension of 390. For JF and LS, the corresponding dimensions are 6561 for the pure neural model and 1605 + 165 = 1770 for the hybrid model. For RC, the hybrid model combines the same 1605-parameter neural component with 403 graded finite element degrees of freedom, giving a total approximation dimension of 2008. The neural parametrizations use smooth tanh activation functions; MS, JF, and LS additionally use fixed Fourier features. The finite element compensation spaces are deliberately coarse. The MS compensation space consists of P1 functions on 16 one-dimensional elements. For JF and LS, the compensation space is built on a 16×12 interface-fitted triangular grid, so that x = 1/2 is represented exactly. For RC, the mesh is graded toward the reentrant corner. The corresponding compensation meshes are shown in Figure 5.2.

18

N. SCHWENCKE, R. MAIER

1.0

0.8

y

0.6

0.4

ε = 1/16 0.2

0.0 0.0

0.2

0.4

0.6

0.8

0.0

1.0

0.2

0.4

0.6

0.8

1.0

x

x

(a) MS: P1 mesh with 15 interior degrees of freedom on 16 elements.

(b) JF: interface-fitted triangular mesh, with x = 1/2 represented exactly.

1.0

1.00 0.75

0.8

0.50

0.6

y

y

0.25 0.00

0.4 −0.25 −0.50

0.2

−0.75

0.0

−1.00

0.0

0.2

0.4

0.6

0.8

1.0

−1.0

x

−0.5

0.0

0.5

1.0

x

(c) LS: interface-fitted triangular mesh, with the line source supported on x = 1/2.

(d) RC: graded P1 mesh toward the reentrant corner.

Figure 5.2. Finite element compensation meshes. The deliberately coarse finite element meshes used inside the hybrid models for MS, JF, LS, and RC. The JF and LS meshes resolve the line x = 1/2 exactly, while the RC mesh is graded toward the reentrant corner. These meshes are distinct from the standalone finite element reference discretizations used in Table 2.

We compare several weak test families. In one dimension, besides the operator-adapted Aε Green sections, we use the H01 (0, 1) point-evaluation representers (5.15)

gt (x) = min(x, t) − xt.

BEYOND PINNS: A UNIFIED FRAMEWORK FOR NEURAL AND HYBRID PDE SOLVERS

19

Table 2. Median final relative errors over five seeds. “Best pure weak” and “best hybrid” select the test family with the smallest relative H 1 error within the corresponding model class; the reported L2 value is that of the same selected method. No strong-residual result is reported for LS because its line source is distributional. The JF finite element value is exceptional because the interfacefitted P4 space contains the manufactured solution. Method

Metric

MS

JF

LS −3

RC

2

GN Deep Ritz

L H1

1

−2

4.48 × 10 5.25 × 102

1.71 × 10 2.53 × 10−2

1.31 × 10 1.14 × 10−1

2.28 × 10−1 1.15

Strong PINN

L2 H1

9.98 × 10−1 9.98 × 10−1

5.97 × 10−4 1.60 × 10−3

— —

1.60 × 10−3 1.26 × 10−2

Best pure weak

L2 H1

2.13 × 10−15 1.06 × 10−13

1.45 × 10−1 2.04

1.16 × 10−4 5.61 × 10−3

1.05 1.40

Best hybrid

L2 H1

1.49 × 10−13 1.08 × 10−12

4.26 × 10−5 4.26 × 10−4

6.72 × 10−4 6.66 × 10−3

5.04 × 10−4 1.20 × 10−2

FE near hybrid budget

L2 H1

2.19 × 10−7 7.31 × 10−5

2.95 × 10−14 5.44 × 10−14

1.36 × 10−3 1.48 × 10−2

1.23 × 10−3 1.60 × 10−2

For JF and LS, if (λk , φk ) denote discrete Dirichlet-Laplacian eigenpairs, we use truncated eigenGreen sections (5.16)

gz(K) (x) =

K X φk (z)φk (x) k=1

λk

,

as well as tensor products of the one-dimensional functions in (5.15). For RC, we additionally −1 consider the shifted Sobolev family obtained by replacing λ−1 . In all four k in (5.16) with (λk + 1) problems, we also consider compact piecewise-affine hat functions with randomly chosen locations. The global test families use 4096 test locations in one dimension and 8192 in two dimensions, while the local-hat experiments use 1000 sampled test functions. Fixed test functions are normalized in the corresponding energy norm. The Fourier features belong only to the trial parametrization and are not used as test functions. Weak residuals are evaluated directly as variational functionals. For JF, quadrature is split at the interface x = 1/2 and, when necessary, at the breakpoints of the test functions. For LS, the singular contribution is evaluated separately as   ˆ 1 1 2 y(1 − y)z , y dy. 2 0 For RC, both residual and error quadratures are graded toward the reentrant corner. As in the smooth benchmarks, the quadrature points used to evaluate a weak functional do not constitute additional residual measurements. 5.2.2. Optimization and evaluation. All trainable neural components are optimized with DSGNAR [43] in double precision, with at most 300 iterations over five seeds; we report median relative L2 and full relative H 1 errors evaluated by quadratures independent of the training residual, with interquartile bands in convergence plots. Ablations comparing the exact projected hybrid formulation of Section 4 with lagged-Jacobian and genuinely alternating finite element–neural variants are provided in the companion blog ↗ .

20

N. SCHWENCKE, R. MAIER

100

100

relative error

relative error

10−2 10−4 10−6 10−8 10−10 10

10−1 10−2 10−3 10−4

−12

10−5 0

50

100

150

0

200

25

50

(a) MS: multiscale diffusion.

100

125

150

(b) JF: jump forcing.

100

100

relative error

relative error

75

iteration

iteration

10−1

10−1

10−2

10−2

10−3

0

50

100

150

200

0

50

(c) LS: line source. Pure neural – L 2

100

150

200

iteration

iteration

(d) RC: reentrant corner. Pure neural – H 1

Hybrid FE–neural – L 2

Hybrid FE–neural – H 1

Figure 5.3. Effect of the hybrid finite element–neural construction for generic weak test functions. Median relative L2 and full relative H 1 errors over five seeds, with interquartile bands, for pure neural and hybrid finite element– neural models using the same compact random-hat test family. The hybrid construction is essentially neutral for MS, where the pure neural model already reaches very high accuracy, but substantially improves JF, LS, and RC. The pure neural runs for JF, LS, and RC terminate after 15 iterations under the optimizer stopping criterion.

5.2.3. Results. The results first show that no single test construction is uniformly preferable. On JF, the best hybrid method uses tensor products of the one-dimensional H01 (0, 1) sections and reaches relative errors 4.26 × 10−5 in L2 and 4.26 × 10−4 in H 1 . On RC, compact local hats are best, reaching 5.04×10−4 and 1.20×10−2 , respectively. On LS, by contrast, the eigen-Green family gives the best pure neural result, with errors 1.16 × 10−4 in L2 and 5.61 × 10−3 in H 1 . Finally,

BEYOND PINNS: A UNIFIED FRAMEWORK FOR NEURAL AND HYBRID PDE SOLVERS

21

Table 3. Median final relative errors over five seeds for pure neural and hybrid finite element–neural solvers using the same random-hat test family. Method

Metric

MS

JF −13

LS

RC −1

2

Pure

L H1

2.19 × 10 8.26 × 10−12

2.05 4.43

9.93 × 10 9.89 × 10−1

1.052 1.40

Hybrid

L2 H1

9.14 × 10−13 2.16 × 10−11

1.71 × 10−4 1.34 × 10−3

4.07 × 10−3 3.32 × 10−2

5.04 × 10−4 1.20 × 10−2

on MS the pure H01 family reaches 2.13 × 10−15 and 1.06 × 10−13 . The test space is therefore a numerical design choice rather than a fixed consequence of the operator. The LS experiment illustrates a more structural advantage of the weak formulation. Its line source defines a bounded functional on the variational test space but cannot be represented by an ordinary pointwise residual. Nevertheless, the weak solver reaches errors of order 10−4 in L2 and 10−3 in H 1 . Together with the multiscale and reentrant-corner benchmarks, this shows that the finite-dimensional Gauss–Newton construction extends naturally to problems for which the strong residual is inconvenient or not defined pointwise. To isolate the effect of the finite element component in the hybrid construction, we compare pure neural and hybrid finite element–neural solvers using the same random-hat test family. The corresponding final errors are reported in Table 3. Hybridization is not a monotone improvement in final accuracy. On MS, the pure neural solver already resolves the solution to very high accuracy, and the hybrid construction provides essentially no further improvement at convergence. On JF, LS, and RC, however, it converts the same generic weak test construction into an accurate solver. The relative H 1 errors decrease from 4.43 to 1.34 × 10−3 on JF, from 9.89 × 10−1 to 3.32 × 10−2 on LS, and from 1.40 to 1.20 × 10−2 on RC, with analogous improvements in L2 . The standalone finite element (FE) references in Table 2 are selected from degree and mesh sweeps with approximation dimensions between approximately 0.8 and 1.2 times the corresponding hybrid budget. At these comparable dimensions, the hybrid models substantially outperform the selected finite element reference on MS, improve upon it on LS, and are moderately more accurate on RC in both reported norms. JF is exceptional: once the interface x = 1/2 is fitted, the manufactured solution belongs to the tested P4 space and is therefore recovered to machine precision. For RC, the selected finite element reference has 2409 degrees of freedom, compared with a total approximation dimension of 2008 for the hybrid model. These standalone references are independent of the coarse finite element compensation spaces based on the meshes shown in Figure 5.2. To complement the comparison at fixed representation dimension, we also report indicative wallclock timings for the selected standalone finite element references and for the corresponding pure neural and hybrid methods. For the neural methods, we report both the median optimization time Toptim and the median total runtime Ttotal over the available seeds; for the finite element references, we report the median warm solve time TFE . These timings should be interpreted cautiously: they come from the present implementations and hardware setup, and they measure computational cost rather than approximation power. In particular, the finite element timings concern standalone solves, whereas the neural timings correspond to iterative optimization procedures. The timings in Table 4 show that hybridization can substantially reduce the computational cost of the neural solve. Relative to the corresponding pure neural method, the total runtime is reduced by factors of approximately 10.5, 2.1, and 5.5 on MS, JF, and LS, respectively. The apparently shorter pure neural runtime on RC should not be interpreted as an efficiency advantage:

22

N. SCHWENCKE, R. MAIER

Table 4. Indicative wall-clock timings for the pure neural, hybrid, and standalone finite element methods. For the neural methods, Toptim denotes the median optimization time and Ttotal the median total runtime. For the finite element references, TFE denotes the median warm solve time. Benchmark

Method

Toptim (s)

Ttotal (s)

TFE (s)

MS

Pure neural Hybrid

299.5 26.4

304.2 28.9

0.040

JF

Pure neural Hybrid

75.7 36.2

80.0 38.8

0.309

LS

Pure neural Hybrid

134.5 22.5

137.3 25.1

0.345

RC

Pure neural† Hybrid

9.7 36.2

12.3 38.4

2.172

† The pure neural RC run stops after 15 iterations without reaching an accurate solution; the corresponding

runtime is therefore not representative of successful convergence.

the optimizer stalls after only 15 iterations and does not reach an accurate solution, whereas the hybrid method continues to convergence. Standalone finite element solves remain considerably faster in absolute wall-clock time for the discretizations considered here. This comparison, however, should be distinguished from representation efficiency. At comparable approximation dimensions, the hybrid models are substantially more accurate than the selected finite element reference on MS, improve upon it on LS, and are slightly more accurate on RC, while JF is a favorable special case for interface-fitted finite elements. Thus, the hybrid construction can both reduce the optimization cost relative to a pure neural solver and provide a more compact high-accuracy approximation than a finite element space of comparable dimension. The comparison is also relevant in light of [15], where standard PINN formulations were found to compare unfavorably with finite element methods in both accuracy and computational cost. Subsequent natural-gradient approaches [32, 39, 19, 40, 18, 20, 30, 43] have already shown that the accuracy attainable by PINNs can depend strongly on the optimization method. The present experiments extend this observation to weak formulations: once the functional residual is discretized and optimized within the Gauss–Newton framework developed above, weak neural solvers can attain high accuracy even for problems for which pointwise residual formulations are poorly suited or not defined. The hybrid finite element–neural construction further improves the robustness of such formulations for generic test functions. A direct quantitative comparison with the benchmark study of [15], under matched problem, implementation, and hardware conditions, is left for future work. 5.2.4. Representative solution decompositions. To complement the aggregate error measures, we examine representative pointwise diagnostics for MS and RC. For each problem, we select the seed whose final relative H 1 error is closest to the five-seed median of the corresponding hybrid method; no best-seed selection is performed. For MS, we also compare with an independently computed finite element approximation near the hybrid budget. For MS, the representative diagnostics are shown in Figure 5.4. Panels a–c correspond to the representative hybrid run, while d shows the independent finite element reference computed with a number of free parameters close to the hybrid model. The diagnostics in Figure 5.4 show that the hybrid and neural errors are essentially indistinguishable, while the finite element compensator itself remains numerically negligible. Thus the hybrid approximation is carried almost entirely by the neural component.

BEYOND PINNS: A UNIFIED FRAMEWORK FOR NEURAL AND HYBRID PDE SOLVERS

|vθ, n − u|

23

|vθ − u| 10−13

10−14 10−14

10−15

10−15

10−16

10−16 0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

x

0.6

0.8

1.0

x

(a) Hybrid error |vθ,n − u|.

(b) Neural-component error |vθ − u|.

|wn|

|uh − u| 10

10−13

−7

10−8 10−9

10−14

10−10

10−15

10−11 10−12

10

−16

0.0

0.2

0.4

0.6

0.8

1.0

0.0

x

0.2

0.4

0.6

0.8

1.0

x

(c) Magnitude of finite element compensator |wn |.

(d) Finite element near hybrid budget error |uh −u|.

Figure 5.4. MS: diagnostics for the hybrid decomposition. The first three panels correspond to the representative hybrid run: the hybrid approximation error |vθ,n − u|, the neural component error |vθ − u|, and the magnitude |wn | of the finite element compensator. The fourth panel shows the pointwise absolute error of the independent finite element approximation near the hybrid budget. The selected finite element reference uses degree-4 elements on 96 uniform elements, with 383 degrees of freedom, compared with a total approximation dimension of 390 for the hybrid model. This behavior, while initially counter-intuitive, has a simple structural explanation. Writing vθ,n = vθ + wn ,

wn ∈ Vn ,

exact recovery of the smooth solution would imply wn = u − vθ . For the present experiment, both u and vθ are smooth, whereas wn belongs to the continuous piecewise-affine finite element space Vn . Hence exact recovery would require wn ∈ Vn ∩ C 1 ([0, 1]).

24

N. SCHWENCKE, R. MAIER exact solution u

1.00

1.00

0.75

hybrid approximation vθ, n 1.00

0.75 0.4

0.4 0.50

0.006

0.50

0.004

0.25

0.002

0.00

0.000

−0.25

−0.002

−0.50

−0.004

0.25

−0.25 −0.50

0.1

0.00 0.2

−0.25 −0.50

0.1

−0.75

−0.75 −1.00 −1.0

0.0 −0.5

0.0

0.5

1.0

−1.00 −1.0

y

0.2

0.3

vθ, n

u

0.00

y

0.3

−0.75

0.0 −0.5

0.0

0.5

1.0

x

x

(a) Exact solution u.

(b) Hybrid approximation vθ,n .

−1.00 −1.0

wn

0.50

0.25

y

FE compensator wn

0.75

−0.006

−0.5

0.0

0.5

1.0

x

(c) Finite element compensator wn .

Figure 5.5. RC: hybrid solution decomposition. Exact solution u, hybrid approximation vθ,n , and finite element compensator wn for the representative hybrid run. The exact and hybrid solutions are displayed on the same color scale, whereas the finite element compensator uses its own scale to reveal its much smaller, localized contribution near the reentrant corner. A globally C 1 piecewise-affine function is necessarily globally affine; combined with the homogeneous Dirichlet conditions, this gives Vn ∩ C 1 ([0, 1]) = {0}. Thus, in the exact-recovery limit, necessarily wn = 0 and vθ = u. The numerically negligible finite element contribution observed in Figure 5.4 is therefore consistent with the regularity mismatch: introducing a nonzero piecewise-affine compensator would force the neural component to approximate the less regular difference u − wn , rather than the smooth solution u itself. This argument, however, concerns the converged decomposition and not the optimization path leading to it. During training, the finite element component can make a non-negligible contribution before becoming negligible near convergence, consistently with the behavior observed in Section 5.2.3. Thus, the nearly vanishing final compensator should not be interpreted as indicating that hybridization does not play a role in the optimization. More detailed trajectory diagnostics are provided in the companion blog ↗ . The finite element approximation with a budget of free parameters close to the hybrid model (shown in Figure 5.4) serves a different purpose. It is an independent finite element reference at approximately the same total approximation budget and should not be confused with the much coarser finite element compensation space with 15 degrees of freedom used inside the hybrid model. For RC, the representative decomposition is shown in Figure 5.5. In contrast to MS, the finite element compensator is non-negligible at convergence and is strongly localized near the reentrant corner. To examine this localization more closely, Figure 5.6 shows the signed hybrid difference vθ,n − u together with the signed finite element compensator wn over successively smaller neighborhoods of the corner. Unlike the MS case, the finite element contribution does not vanish in the converged RC approximation. Instead, Figure 5.6b shows a localized, sign-changing correction concentrated around the reentrant corner, where the exact solution has reduced regularity. The successive zooms show that this structure persists down to the smallest displayed scale; the [−0.05, 0.05]2 window still contains several layers of the graded finite element mesh. The signed hybrid difference in Figure 5.6a provides the corresponding local approximation diagnostic on the same spatial and color scales. Detailed test-family comparisons, finite element mesh and degree sweeps, per-seed results, hybrid-update ablations, and further numerical diagnostics are provided in the companion blog ↗ .

BEYOND PINNS: A UNIFIED FRAMEWORK FOR NEURAL AND HYBRID PDE SOLVERS

25 0.008

1.0

vθ, n − u on [−1, 1]2

0.35

vθ, n − u on [−0.35, 0.35]2

0.12

vθ, n − u on [−0.12, 0.12]2

0.05

vθ, n − u on [−0.05, 0.05]2 0.004

0.000

0.00

vθ, n − u

0.00

y

0.00

y

0.0

y

y

0.5

−0.5

−0.004 −1.0 −1.0

−0.5

0.0

0.5

−0.35 −0.35

1.0

x

0.00

−0.12 −0.12

0.35

0.00

x

−0.05 −0.05

0.12

x

0.00

0.05

x −0.008

(a) Signed hybrid difference vθ,n − u. 0.008

1.0

wn on [−1, 1]2

0.35

wn on [−0.35, 0.35]2

0.12

wn on [−0.12, 0.12]2

0.05

wn on [−0.05, 0.05]2 0.004

0.000

0.00

wn

0.00

y

0.00

y

0.0

y

y

0.5

−0.5

−0.004 −1.0 −1.0

−0.5

0.0

x

0.5

1.0

−0.35 −0.35

0.00

0.35

−0.12 −0.12

x

0.00

0.12

x

−0.05 −0.05

0.00

0.05

x −0.008

(b) Signed finite element compensator wn .

Figure 5.6. RC: localization near the reentrant corner. Signed hybrid difference vθ,n − u and finite element compensator wn for the representative hybrid run. From left to right, the panels show the full domain and the nested windows [−0.35, 0.35]2 , [−0.12, 0.12]2 , and [−0.05, 0.05]2 , restricted to the L-shaped domain. Both quantities are displayed with the same symmetric color scale, allowing their magnitudes and spatial localization to be compared directly. The multilevel zoom layout and associated plotting implementation are adapted from [8].

References [1] H. Barucq, M. Duprez, F. Faucher, E. Franck, F. Lecourtier, V. Lleras, V. Michel-Dansac, and N. Victorion. Enriching continuous Lagrange finite element approximation spaces using neural networks. ArXiv Preprint, 2502.04947, 2025. [2] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, and J. M. Siskind. Automatic differentiation in machine learning: A survey. Journal of Marchine Learning Research, 18:1–43, 2018. [3] D. Bon, B. Caris, and O. Mula. Stable Nonlinear Dynamical Approximation with Dynamical Sampling, May 2025. [4] D. Braess. Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics. Cambridge University Press, Cambridge, 2007. [5] S. C. Brenner and L. R. Scott. The Mathematical Theory of Finite Element Methods, volume 15 of Texts in Applied Mathematics. Springer, New York, 3 edition, 2008. doi: 10.1007/978-0-387-75934-0. [6] K. Cho, B. van Merrienboer, D. Bahdanau, and Y. Bengio. On the Properties of Neural Machine Translation: Encoder-Decoder Approaches. ArXiv Preprint, 1409.1259, 2014. doi: 10.48550/arXiv.1409.1259.

26

N. SCHWENCKE, R. MAIER

[7] P. G. Ciarlet. The Finite Element Method for Elliptic Problems, volume 4 of Studies in Mathematics and its Applications. North-Holland, Amsterdam, 1978. [8] A. Combette, A. Venaille, and N. Pustelnik. A new initialisation to Control Gradients in Sinusoidal Neural network. ArXiv Preprint, 2512.06427, 2025. [9] A. Daw, J. Bu, S. Wang, P. Perdikaris, and A. Karpatne. Mitigating Propagation Failures in Physics-informed Neural Networks using Retain-Resample-Release (R3) Sampling, 2022. [10] L. Demkowicz and J. Gopalakrishnan. A class of discontinuous Petrov–Galerkin methods. II. Optimal test functions. Numerical Methods for Partial Differential Equations, 27(1):70–105, 2011. [11] M. W. M. G. Dissanayake and N. Phan-Thien. Neural-network-based approximations for solving partial differential equations. Communications in Numerical Methods in Engineering, 10(3):195–201, 1994. [12] J. L. Elman. Finding Structure in Time. Cognitive Science, 14(2):179–211, Mar. 1990. [13] E. Franck, V. Michel-Dansac, L. Navoret, and V. Vigon. Neural semi-Lagrangian method for high-dimensional advection-diffusion problems. Computer Methods in Applied Mechanics and Engineering, 448(B):118481, 2026. [14] I. Goodfellow, Y. Bengio, and A. Courville. Deep Learning. MIT press, 2016. [15] T. G. Grossmann, U. J. Komorowska, J. Latz, and C.-B. Schönlieb. Can physics-informed neural networks beat the finite element method? IMA Journal of Applied Mathematics, 89 (1):143–174, 2024. [16] W. Hao, Q. Hong, and X. Jin. Gauss Newton method for solving variational problems of PDEs with neural network discretizaitons. ArXiv Preprint, 2306.08727, 2023. [17] S. Hochreiter and J. Schmidhuber. Long short-term memory. Neural computation, 9(8):1735– 1780, 1997. [18] A. Jnini and F. Vella. Dual Natural Gradient Descent for Scalable Training of PhysicsInformed Neural Networks. ArXiv Preprint, 2505.21404, 2025. [19] A. Jnini, F. Vella, and M. Zeinhofer. Gauss-Newton natural gradient descent for physicsinformed computational fluid dynamics. Computers & Fluids, 307:106955, 2025. [20] A. Jnini, E. Kiyani, K. Shukla, J. F. Urban, N. A. Daryakenari, J. Muller, M. Zeinhofer, and G. E. Karniadakis. Curvature-Aware Optimization for High-Accuracy Physics-Informed Neural Networks. ArXiv Preprint, 2604.05230, 2026. [21] E. Kharazmi, Z. Zhang, and G. E. Karniadakis. Variational Physics-Informed Neural Networks For Solving Partial Differential Equations. ArXiv Preprint, 1912.00873, 2019. [22] E. Kharazmi, Z. Zhang, and G. E. Karniadakis. Hp-VPINNs: Variational physics-informed neural networks with domain decomposition. Computer Methods in Applied Mechanics and Engineering, 374:113547, 2021. [23] I. E. Lagaris, A. Likas, and D. I. Fotiadis. Artificial neural networks for solving ordinary and partial differential equations. IEEE transactions on neural networks, 9(5):987–1000, 1998. [24] G. K. R. Lau, A. Hemachandra, S.-K. Ng, and B. K. H. Low. PINNACLE: PINN adaptive ColLocation and experimental points selection. In The Twelfth International Conference on Learning Representations, 2024. [25] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner. Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(11):2278–2324, 1998. [26] S. Linnainmaa. Taylor expansion of the accumulated rounding error. BIT Numerical Mathematics, 16(2):146–160, 1976. [27] A. Målqvist and D. Peterseim. Numerical homogenization by localized orthogonal decomposition, volume 5 of SIAM Spotlights. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2021. [28] Z. Mao and X. Meng. Physics-informed neural networks with residual/gradient-based adaptive sampling methods for solving partial differential equations with sharp solutions. Applied Mathematics and Mechanics, 44(7):1069–1084, 2023.

BEYOND PINNS: A UNIFIED FRAMEWORK FOR NEURAL AND HYBRID PDE SOLVERS

27

[29] N. Margenberg, R. Jendersie, C. Lessig, and T. Richter. DNN-MG: A hybrid neural network/finite element method with applications to 3D simulations of the Navier–Stokes equations. Computer Methods in Applied Mechanics and Engineering, 420:116692, 2024. [30] M. B. McKay, A. Kaur, C. Greif, and B. Wetton. Near-optimal Sketchy Natural Gradients for Physics-Informed Neural Networks, 2025. URL https://openreview.net/forum?id= bKsZomnmqn. [31] M. B. McKay, N. P. Lawrence, B. Wetton, and R. B. Gopaluni. Error whitening: Why Gauss-Newton outperforms Newton. ArXiv Preprint, 2605.11316, 2026. [32] J. Müller and M. Zeinhofer. Achieving high accuracy with PINNs via energy natural gradient descent. In International Conference on Machine Learning, pages 25471–25485. PMLR, 2023. [33] M. A. Nabian, R. J. Gladstone, and H. Meidani. Efficient training of physics-informed neural networks via importance sampling. Computer-Aided Civil and Infrastructure Engineering, 36 (8):962–977, 2021. [34] T. N. K. Nguyen, T. Dairay, R. Meunier, C. Millet, and M. Mougeot. Fixed-Budget Online Adaptive Learning for Physics-Informed Neural Networks. Towards Parameterized Problem Inference. In J. Mikyška, C. de Mulatier, M. Paszynski, V. V. Krzhizhanovskaya, J. J. Dongarra, and P. M. Sloot, editors, Computational Science – ICCS 2023, pages 453–468, Cham, 2023. Springer Nature Switzerland. [35] A. Nouy and A. Somacal. Natural gradient descent with momentum. ArXiv Preprint, 2604.15554, 2026. [36] V. I. Paulsen and M. Raghupathi. An Introduction to the Theory of Reproducing Kernel Hilbert Spaces, volume 152. Cambridge University Press, 2016. [37] M. Raissi, P. Perdikaris, and G. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686–707, 2019. [38] S. Rojas, P. Maczuga, J. Muñoz-Matute, D. Pardo, and M. Paszyński. Robust Variational Physics-Informed Neural Networks. Computer Methods in Applied Mechanics and Engineering, 425:116904, 2024. [39] N. Schwencke and C. Furtlehner. ANaGRAM: A natural gradient relative to adapted model for efficient PINNs learning. In The Thirteenth International Conference on Learning Representations, 2025. [40] N. Schwencke, C. Rousselot, A. Shilova, and C. Furtlehner. AMStramGRAM: Adaptive multicutoff strategy modification for ANaGRAM. ArXiv Preprint, 2510.15998, 2025. [41] Y. Shang, F. Wang, and J. Sun. Deep Petrov-Galerkin Method for Solving Partial Differential Equations. ArXiv Preprint, 2201.12995, 2022. [42] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. Polosukhin. Attention Is All You Need. ArXiv Preprint, 1706.03762, 2017. [43] J. Webb, S. Jerad, and C. Cartis. An Optimisation Framework for the Well-Conditioned Training of Physics-Informed Neural Networks. ArXiv Preprint, 2607.02194, 2026. [44] C. Wu, M. Zhu, Q. Tan, Y. Kartha, and L. Lu. A comprehensive study of non-adaptive and residual-based adaptive sampling for physics-informed neural networks. Computer Methods in Applied Mechanics and Engineering, 403:115671, 2023. [45] Y. Zang, G. Bao, X. Ye, and H. Zhou. Weak adversarial networks for high-dimensional partial differential equations. Journal of Computational Physics, 411:109409, 2020.

Appendix A. Classical methods as instances of the Petrov–Galerkin framework The framework of Section 3.2 separates two ingredients: the local approximation space given by the tangent space of the parametric model and the measurements, or equivalently test functions, used to discretize the linearized residual equation. We briefly record how several familiar

28

N. SCHWENCKE, R. MAIER

constructions arise from particular choices of these ingredients. Throughout, we use the notation introduced in Section 3.2. For more detailed developments and complementary explanations, particularly from a machine-learning perspective, we refer to the companion blog ↗ . A.1. Natural gradient. Let v : Rp → H be a parametric model with residual rθ = vθ − y. Choosing the tangent directions ∂i vθ themselves as test functions in (3.14) gives (A.1)

p X

⟨∂j vθ , ∂i vθ ⟩H ξj = ⟨rθ , ∂i vθ ⟩H ,

1 ≤ i ≤ p.

j=1

Writing Gθ,ij = ⟨∂i vθ , ∂j vθ ⟩H , this is Gθ ξ = ∇θ ℓ(θ) for ℓ(θ) = 12 ∥rθ ∥2H . Hence the corresponding correction is the natural-gradient direction [32, 19, 31], and its functional action satisfies (A.2)

dvθ [ξ] = ΠTθ Mv rθ .

Thus, natural gradient is the Galerkin choice in which the approximation and test spaces both coincide with the tangent space. A.2. Kernel and empirical natural-gradient methods. Assume that point evaluation is well defined and continuous on a Hilbert space H. Its Riesz representer k(x, •) ∈ H satisfies (A.3)

w(x) = ⟨w , k(x, •)⟩H .

N b For points X = {xi }N i=1 , choosing HX = span{k(xi , •)}i=1 as both approximation and test space b therefore yields the orthogonal projection of y onto HX , equivalently the standard kernel interpolant associated with the points X [36]. The empirical natural gradient of [39] is the local version of this construction. Point evaluation on the tangent space Tθ Mv admits representers kθ (xi , •), and the corresponding empirical tangent space is

(A.4)

Tbθ,X Mv = span{kθ (xi , •) | 1 ≤ i ≤ N } ⊂ Tθ Mv .

The empirical natural-gradient correction is consequently characterized by (A.5)

vbθ,X (c) = ΠTbθ,X Mv rθ .

It thus replaces the full tangent-space projection of natural gradient by an orthogonal projection onto the kernel space generated by the chosen measurements. A.3. Pointwise Gauss–Newton and PINN collocation. For a compound parametric model Γ : Rp → Y, let point evaluation at X = {xi }N i=1 be well defined on the corresponding extended tangent space. Choosing (A.6)

λθi (w) = w(xi )

in (3.10) gives (A.7)

bX θ,i = rθ (xi ),

X Jθ,ij = ∂j Γθ (xi ).

The resulting least-squares problem is precisely the standard pointwise Gauss–Newton discretization. Taking Γ to be the PDE compound model therefore recovers the usual collocation-based Gauss–Newton formulation for PINNs. Since point evaluation admits a Riesz representer on the finite-dimensional extended tangent space, this construction is itself a particular Petrov–Galerkin discretization of the form (3.14). Observational data are incorporated in the same way by augmenting the PDE compound model with the identity, (D, B, IH ) ◦ v, and applying point-evaluation measurements to this additional component, thereby recovering the usual supervised data term in PINNs.

BEYOND PINNS: A UNIFIED FRAMEWORK FOR NEURAL AND HYBRID PDE SOLVERS

29

A.4. Sketching and regularization. The same viewpoint also gives concise interpretations of common algebraic modifications of Gauss–Newton systems. If P ∈ Rp×sp restricts parameter corrections and S ∈ Rst ×N combines the measurements, the discretized system becomes (A.8)

θ SJθΛθ P η = SbΛ θ .

Hence parameter-side sketching restricts the local approximation space, whereas residual-side sketching replaces the test family by linear combinations of its elements [30, 43]. Likewise, ridge-regularized Gauss–Newton follows by augmenting the functional model with the scaled identity, vθα = (vθ , αθ). After recentering the augmented target at the current iterate, the local problem becomes 1 α2 minp ∥rθ − dvθ [ξ]∥2H + ∥ξ∥2Rp , ξ∈R 2 2 which yields the usual ridge-regularized Gauss–Newton correction [40]. (A.9)

∗ ENS de Lyon, CNRS, Université Claude Bernard Lyon 1, Inria, LIP, UMR 5668, 69342, Lyon cedex 07, France Email address: [email protected] † Institute for Applied and Numerical Mathematics, Karlsruhe Institute of Technology, Engler-

str. 2, 76131 Karlsruhe, Germany Email address: [email protected]