Numerical analysis of physics-informed neural networks and related models in physics-informed machine learning

1. Introduction

• Physics-informed machine learning uses governing equations to solve PDEs without requiring data about their underlying solutions.This motivates approaches that need only small amounts of data, given the potential expense of generating or accessing solution data.Source: §1 · p. 2
• PINNs minimize a neural-network approximation's strong-form PDE residual using (stochastic) gradient descent.The modern PINN framework was introduced by Raissi and Karniadakis (2018) and Raissi, Perdikaris and Karniadakis (2019).Source: §1 · pp. 2–3
• Related models use variational or weak PDE residuals, energy minimization, or alternative ansatz spaces such as Gaussian processes and neural operators.Examples include VPINNs, wPINNs, and the DeepRitz method.Source: §1 · p. 3
• The review critically analyzes numerical guarantees and limits through approximation, stability, generalization, and training error analyses.It also introduces PDE formulations, residuals, and approximating ansatz spaces.Source: §1 · p. 3

2. Physics-informed machine learning

• The article introduces an abstract PDE setting and examples, model classes, PDE formulations and residuals, then discretizes and optimizes the resulting physics-informed loss.Section 2.6 summarizes the method.Source: §2 · p. 4

2.1. PDE setting

• The abstract setting places PDE solutions in Banach spaces and specifies a differential operator, source, boundary operator, boundary data, and domain.The framework assumes a unique solution for each source and identifies the domain as spatial or space-time.Source: §2.1 · p. 4

2.1.1. Forward problems

• Forward problems map known PDE inputs, such as source, boundary, or initial data, to solutions; learning this whole mapping is operator learning.For time-dependent equations, an example maps initial conditions to the solution or its value at a later time.Source: §2.1.1 · pp. 4–5
• Forward-problem examples include semilinear heat equations, Navier–Stokes equations, scalar conservation laws, and Poisson's equation.The semilinear heat example assumes a globally Lipschitz nonlinearity; inviscid conservation laws may develop shocks even from smooth initial data.Source: §2.1.1 · pp. 5–6

2.1.2. Inverse problems

• Inverse problems recover unknown PDE inputs from solution observations when incomplete inputs prevent a unique forward solution.The review focuses on unique continuation, also called data assimilation for time-dependent problems, using potentially noisy observations.Source: §2.1.2 · pp. 6–7
• When measurements cover only an interior observation domain, the inverse task can infer the solution beyond that domain; the examples include heat and Poisson equations.For the heat equation, the task includes inferring initial conditions and the later-time solution from measurements in time.Source: §2.1.2 · pp. 7–8

2.2. Ansatz spaces

• Ansatz spaces are classes of parameterized functions used to approximate the true PDE solution.The parameter space is a subset of R^n, where n counts the real parameters.Source: §2.2 · p. 8

2.2.1. Linear models

• Linear models combine fixed functions, with spectral methods choosing global bases and FEM using local polynomials on subdomains.FEM produces a sparse linear system; accuracy can be improved by reducing mesh size, increasing polynomial degree, or both.Source: §2.2.1 · pp. 8–9
• RBFs provide mesh-free basis functions, while random feature models draw basis functions from a chosen distribution.RBF nodes can be randomly generated, which is useful in high-dimensional settings; random neural features connect to networks that train only the outer layer.Source: §2.2.1 · pp. 9–10

2.2.2. Nonlinear models Neural networks.

• Neural networks are nonlinear parametric models; an MLP composes affine maps and activation functions, whose properties affect gradient-based training.Sigmoidal activations can have vanishing gradients away from zero, while ReLU reduces that problem but can have dying neurons.Source: §2.2.2 · pp. 10–11
• DeepONets and neural operators represent maps between infinite-dimensional spaces, while Gaussian process regression models functions using a Gaussian measure.Neural operators use learned operator layers; Fourier Neural Operators are a special case, implemented through finite-dimensional mappings after discretization.Source: §2.2.2 · pp. 11–13

2.3.1. Classical formulation

• The classical formulation minimizes PDE and boundary residuals, adding an initial-condition residual for time-dependent problems; boundary conditions may be enforced softly or exactly.The objective can also include data residuals for measurements and parameter regularization.Source: §2.3.1 · pp. 15–16

2.3.2. Weak formulation

• Weak formulations allow solutions that need only be integrable, rather than differentiable, but the available PDE formulations are not universal.A finite test-function set yields VPINNs, while maximizing over a parametric test function avoids choosing a basis but creates a more challenging min-max optimization.Source: §2.3 · p. 14, §2.3.2 · pp. 16–17
• For scalar conservation laws, weak solutions may be nonunique, so entropy conditions are added to recover uniqueness.The cited account states that entropy solutions are unique and continuous in time.Source: §2.3.2 · pp. 17–18

2.3.3. Variational formulation

• Variational formulations replace solving a PDE with minimizing its associated energy functional when the differential operator derives from that functional.The Deep Ritz method minimizes the energy over deep neural networks; an extra term addresses boundary conditions not automatically satisfied by the networks.Source: §2.3.3 · pp. 18–19

• Quadrature approximates each loss integral by a weighted sum of values at selected points, and the resulting discretized objective is minimized to produce the trained model.For the time-dependent classical loss, separate quadratures cover the PDE interior, spatial boundary over time, and initial domain; their points form the training set.Source: §2.4 · pp. 20–21
• For d ≤ 4, grid-based midpoint or composite Gauss rules are suitable; for 4 ≤ d ≤ 20, Sobol or Halton points are suggested, with error rate (log(N))^d N^-1 for bounded Hardy-Krause variation.Source: §2.4 · p. 20
• For dimensions d >> 20, Monte Carlo quadrature is the stated method of choice, using randomly chosen independent, identically distributed points.Source: §2.4 · p. 21

2.5. Optimization

• The high-dimensional parameter space and generally non-convex objective make minimization challenging, although its almost-everywhere differentiability permits gradient-based iterative methods.Source: §2.5 · p. 21
• Full-gradient evaluation requires N = |S| gradients, making it slow and memory-intensive for large training sets; SGD evaluates one data point per iteration at lower per-iteration cost but with noisier convergence.Source: §2.5 · p. 22
• Mini-batch gradient descent evaluates a subset of training data per iteration, interpolating between full-batch gradient descent and SGD as batch size changes.Source: §2.5 · p. 22
• Adam uses moving averages of gradients and adjusts its step size when gradients vary substantially, reducing noise relative to mini-batch gradient descent.Source: §2.5 · pp. 22–23
• L-BFGS is given as a popular quasi-Newton method for scientific computing, while Hessian computation and inversion can be costly.First-order optimizers remain dominant in many other application areas.Source: §2.5 · p. 23

2.5.1. Parameter initialization

• Xavier initialization was proposed with supervised learning in mind, and different schemes may be better suited to the loss functions used in physics-informed learning.A poor initialization can cause exploding or vanishing gradients.Source: §2.5.1 · p. 24

• Choose a model class and a loss based on the classical, weak, or variational PDE formulation.Source: §2.6 · p. 24
• Generate a training set with a suitable quadrature rule to make the loss feasibly computable, then initialize parameters and optimize to an approximate local minimum.Source: §2.6 · p. 24
• The resulting trained function is the model used to approximate the PDE solution.Source: §2.6 · p. 24

3. Analysis

• Q1 asks whether the hypothesis class contains a model with small generalization error, while Q2 asks whether small generalization error implies small total error.The analysis relates Q1 to approximation error and seeks stability bounds for Q2.Source: §3 · p. 25
• Q3 asks whether small training error and a sufficiently large training set yield small generalization error, a question tied to the accuracy of the chosen quadrature rule.The review uses standard arguments to address this question.Source: §3 · p. 26
• Q4 asks whether training can bring the training error close to the minimum of the empirical loss, recognizing that this optimization can be difficult or infeasible.The analysis reviews settings associated with optimization difficulty and potential ways to address training issues.Source: §3 · p. 26
• The total error has approximation, generalization-gap, and optimization components; the optimization component vanishes when the trained parameter equals the empirical minimizer.Affirmative answers to Q1–Q4 lead to a small total error.Source: §3 · p. 27
• The general results are applied to Navier–Stokes, heat, inviscid scalar conservation, and Poisson equations to illustrate different PDE settings.These examples span low- and potentially high-dimensional problems, potentially discontinuous solutions, and a variational formulation.Source: §3 · p. 27

4. Approximation error

• Because physics-informed generalization error measures PDE loss, uniform function approximation alone may leave inaccurate derivatives and large residuals, motivating Sobolev-norm approximation.Source: §4 · pp. 27–28

4.1. First general result

• Fourier-basis and tanh-network results provide derivative-aware approximation bounds for sufficiently smooth functions; the tanh construction uses two hidden layers of width N.The Fourier result assumes s>d/2 and s≥k; the tanh result covers derivative orders k<s.Source: §4.1 · pp. 28–29
• For Navier–Stokes, sufficiently regular divergence-free initial data yield solution regularity that supports two-hidden-layer tanh approximants for velocity and pressure (De Ryck et al. 2023b).Source: §4.1 · pp. 29–30

4.2. Second general result

• Under the stated smoothness and domain assumptions, finite differences can turn Lq approximation into bounds on PDE residuals by estimating derivatives of the model error.The domain assumption enables piecewise finite-difference operators that are well-defined across the domain.Source: §4.2 · p. 31
• Residual bounds become small when βℓ < α(r−ℓ); for the Fourier basis, the resulting convergence rate is N^−(r−ℓ)/d, matching Theorem 4.3.The condition favors fast Lq convergence, slowly growing model derivatives, or greater solution smoothness; Theorem 4.11 need not yield optimal rates with loose derivative bounds.Source: §4.2 · pp. 31–32
• For discontinuous scalar conservation-law solutions, BV approximation and piecewise-smooth tanh networks support small entropy-residual loss without assuming Sobolev regularity.Lemma 4.16 gives two-hidden-layer networks with O(ε^−2) neurons for BV functions; the piecewise-smooth result assumes a smooth discontinuity surface.Source: §4.2 · pp. 32–33

4.3. Neural networks overcome curse of dimensionality in approximation error

• A time-dependent Sobolev-regularity framework yields space-time tanh networks with small physics-informed loss from efficient fixed-time approximators, and can extend to operator learning.The framework assumes controlled approximation error and Sobolev norms for the fixed-time networks.Source: §4.3 · pp. 33–34
• A second route avoids the efficient-initial-condition assumption by restricting PDE solutions to Barron space, where functions admit approximation without the curse of dimensionality.The remaining challenge is proving that a given PDE solution belongs to Barron space.Source: §4.3 · p. 33
• Generic DeepONet approximation bounds use Z^d sensors and specify branch and trunk network width and depth under operator stability assumptions (De Ryck and Mishra 2022b).The stated widths are O(M(Z^d + N^d width(U_ε))) for the branch network and O(MN^d(N + ln(N))) for the trunk; trunk depth is 3.Source: §4.3 · pp. 35–36
• For linear Kolmogorov PDEs, PINNs can overcome the curse of dimensionality with size O(poly(dρ_d)ε^-β), and ρ_d is polynomial in d for a subclass (De Ryck and Mishra 2022a, 2022b).The result is demonstrated under assumptions of polynomial growth of solution regularity and neural-network approximability of the initial condition.Source: §4.3 · p. 36
• PINNs also overcome the curse of dimensionality for nonlinear parabolic PDEs, while sufficiently smooth solutions yield dimension-independent convergence rates for physics-informed DeepONets (De Ryck and Mishra 2022b).The DeepONet result assumes C^(2,r) functions; assuming greater regularity did not necessarily improve the convergence rate, and it improves on a prior result that incurred the curse of dimensionality (Lanthaler, Mishra and Karniadakis 2022).Source: §4.3 · p. 37

4.3.2. Results based on Barron regularity

• Barron-space results characterize PDE solution regularity, including Poisson solutions with ∥u∥_B^s(Ω) ≤ d∥f∥_B^(s+2)(Ω) and static Schrödinger solutions with ∥u∥_B^s(Ω) ≤ C∥f∥_B^(s+2)(Ω).The Poisson result assumes zero-mean forcing and zero Neumann boundary conditions; the Schrödinger result assumes V ∈ B^(s+2)(Ω) and V(x) ≥ V_min > 0.Source: §4.3.2 · pp. 38–39
• Barron norms are independent of p, Barron spaces are algebras under conditions on the activation, and neural networks approximate Barron functions in Sobolev norm without the curse of dimensionality (Chen et al. 2021; E et al. 2022).The Sobolev approximation result is presented as an adaptation of Chen et al. (2021, Theorem 2.5).Source: §4.3.2 · p. 39
• If ∥f∥_B^2(Ω) is known to scale with d, the Barron approximation framework can establish neural networks overcoming the curse of dimensionality for physics-informed loss functions under Assumption 4.1.The argument follows the approach of Section 4.1.Source: §4.3.2 · p. 41

5.1. Stability for forward problems

• For forward problems, a bounded inverse of a linear operator or of a nonlinear operator’s linearization ensures that small physics-informed loss implies small total error.The section gives explicit stability examples for semilinear heat, Navier–Stokes, scalar conservation, and Poisson equations, and reports an affirmative result for radiative transfer.Source: §5.1 · pp. 42–43
• For scalar conservation laws, strong-residual bounds apply when viscosity is positive but can fail at low viscosity; weak-residual wPINNs retain stability for shock-containing solutions.The section therefore gives an affirmative answer using either the classical residual for ν > 0 or the weak residual for ν = 0.Source: §5.1 · pp. 45–46
• For high-dimensional HJB equations, an L2-based classical-residual loss can yield solutions arbitrarily far from the truth, while stability requires an Lp residual with p scaling with dimension.The cited result states that the linear dependence of p on dimension d cannot be relaxed.Source: §5.1 · p. 47
• cPINNs have affirmative stability results for Navier–Stokes, advection-diffusion, and scalar conservation equations, whereas XPINNs were not stable for Poisson or heat equations; regularized PINNs also have consistency results for elliptic and parabolic PDEs.The consistency result requires a specific Ck,α-regularization term and gives no quantitative error estimates; a related advection-equation result is more quantitative.Source: §5.1 · p. 48

5.2. Stability for inverse problems

• For inverse problems, the general error estimate bounds error using PDE and data residuals, provided a conditional stability estimate holds for sufficiently regular solutions.The estimate assumes sufficiently smooth models and square-integrable residuals.Source: §5.2 · pp. 48–49
• For inverse Poisson problems, the three-balls inequality supplies conditional stability, and additional results bound error across the whole domain, including in Sobolev norm.For heat-equation data assimilation, the section likewise presents an error bound based on classical stability results and Carleman estimates.Source: §5.2 · pp. 49–50
• The inverse data-assimilation problem seeks the Stokes velocity and pressure from the forcing, observed velocity data, and boundary data.The analysis considers estimating model error on balls in the domain whose radii exceed that of the largest ball contained in the observation domain.Source: §5.2 · p. 50
• For sufficiently smooth models and the stated regularity of the Stokes data and solution, Theorem 5.23 bounds generalization error on balls larger than the observation ball.The theorem assumes k ≥ 2; the supplied passage refers to the bound as (5.42) but does not include its formula, and states that C depends on u and uθ.Source: §5.2 · p. 51

General a posteriori bound

• Theorem 6.2 gives a computable a posteriori bound for expected generalization error using a randomly generated training set, under boundedness and Lipschitz assumptions.For neural networks, the sizes of those constants may be unclear, and bounds on c can depend exponentially on depth.Source: §6 · pp. 52–53

Navier-Stokes

• With a sufficiently large training set, the Navier–Stokes total error scales as E(θ) ∼ √E_T(θ, S), and the bound requires N_int ≫ N_s ≫ N_t.The figure plots total and training error against N_int, alongside the Theorem 6.5 bound, and against each other.Source: §6 · p. 54, Figure 6.1 · p. 55

High-dimensional heat equation

• For high-dimensional heat equations, the training-set size need grow at most polynomially with input dimension to achieve a specified accuracy.The comparison with Navier–Stokes also shows that fewer interior points are needed: N_s ≫ N_t ≈ N_int is sufficient.Source: §6 · pp. 54–55

Heat equation

• For the heat equation, the total error also scales as E(θ) ∼ √E_T(θ, S) with sufficiently many training points, while N_s ≫ N_t ≈ N_int is sufficient.The one-dimensional experiment reports decreasing generalization and training errors as quadrature points increase.Source: §6 · p. 55

Parametric heat equation

• The heat equation’s relative generalization error grows sub-quadratically, rather than exponentially, with parameter dimension d and remains below 2% even for very large d.Figure 6.3 plots relative generalization error against parameter dimension.Source: §6 · p. 56, Figure 6.3 · p. 56

7. Training

• PINNs can be slow or infeasible to train, and training may fail to converge or may end at unacceptably large loss values.The section identifies principled analysis of gradient-descent training as limited.Source: §7 · pp. 57–58

7.1. Global minimum of the loss approximates PDE solution well

• A priori error estimates follow when a global minimum of the physics-informed loss can be found, by combining approximation, stability, and generalization results.The Navier–Stokes estimate requires a sufficiently smooth solution and sufficiently large network and training set.Source: §7.1 · pp. 58–60
• For the Poisson example, the convergence rate is independent of dimension, mitigating the curse of dimensionality, although the dimension dependence of the constant is unclear.The constant may depend on the solution’s Barron norm.Source: §7.1 · p. 60

7.2. Characterization of gradient descent

• When the linearization error stays small, simplified gradient descent approximates full gradient descent; wide networks make that error small as their neural tangent kernels remain approximately constant.The small-error condition is ensured for linear models and can follow from approximate linearity near initialization.Source: §7.2 · pp. 62–63
• For the simplified dynamics, convergence speed depends on the operator’s condition number, and a large condition number requires many more gradient-descent steps.The analysis identifies preconditioning as a way to alleviate training issues.Source: §7 · p. 58, §7.2 · p. 63

7.3. Training and operator preconditioning

• Ill-conditioning of the Hermitian-square/kernel operator implies ill-conditioning of the training matrix and can make gradient descent very slow.With an invertible Gram matrix, the theorem gives equality between the matrix and operator condition numbers; the relation also extends to losses with boundary conditions.Source: §7.3 · p. 64

7.4.1. General framework for preconditioning.

• Preconditioning seeks a model kernel operator approximating the inverse Hermitian square, implemented by a linear parameter transform that changes A to P⊤AP.The transformed parameters are equivalent to preconditioning the loss gradient by a positive definite matrix.Source: §7.4.1 · p. 65

7.4.2. Preconditioning linear physics-informed machine learning models.

• For Poisson Fourier features, the unpreconditioned condition number grows as K^4, while preconditioning keeps it near 1 and yields frequency-independent exponential loss decay.The text reports analogous results for the Helmholtz equation; for linear advection, the condition number grows as β^2 without preconditioning and remains constant with it, alongside faster training.Source: §7.4.2 · pp. 66–68, Figure 7.1 · p. 67, Figure 7.2 · p. 67

7.4.3. Preconditioning nonlinear physics-informed machine learning models.

• For nonlinear MLPs, a spread-out spectrum accompanies slow training; Fourier-feature rescaling makes the spectrum more even and the loss decays more rapidly.The reported rescaling uses α_k = 1/k^2 for k ≠ 0, though near-zero eigenvalues remain.Source: §7.4.3 · pp. 68–69, Figure 7.3 · p. 68

7.5. Strategies for improving training in physics-informed machine learning

• Different strategies for choosing the boundary-loss weight λ produce different scalings with Fourier frequency K, yet the resulting condition numbers are reported as similar.The excerpt introduces an optimal λ and an annealing approach, but its comparison statement is truncated.Source: §7.5.1 · pp. 69–70

7.5.2. Hard boundary conditions.

• In this example, hard-boundary variant 2 has condition number 1, and hard-boundary conditions are strictly better conditioned than soft-boundary conditions; different choices of the boundary-enforcing function can change conditioning.Source: §7.5.2 · p. 71

7.5.3. Second-order optimizers.

• For linear physics-informed models, (quasi-)Newton methods approximately invert the loss Hessian, thereby preconditioning the training matrix A.The stated Hessian identity applies to linear models or neural networks in the NTK regime.Source: §7.5.3 · p. 71

7.5.4. Domain decomposition.

• For linear advection, splitting the time domain into N subintervals reduces each submodel’s condition number by a factor of N^2 because conditioning scales quadratically with advection speed β.The passage connects this domain-decomposition interpretation to causal-learning training of PINNs.Source: §7.5.4 · p. 72

8. Conclusion

• The total error can be estimated through approximation error, generalization gap, and optimization error when the PDE solution is suitably stable.The review analyzes these components generally and for prototypical PDEs.Source: §8 · p. 72
• PDE stability under perturbations and sufficient solution regularity are important to error estimates, and the required regularity can depend on the PDE formulation.The review notes that stability requirements can be PDE-specific.Source: §8 · p. 72
• Training is a key bottleneck because it depends on the Hermitian-square spectrum of the PDE operator, although preconditioning might ameliorate training difficulties.The review says theory guides use: regularity informs strong versus weak formulations, and PINNs may be attractive for high-dimensional, complex-geometry, or inverse problems.Source: §8 · pp. 72–73

A. Sobolev spaces

• The appendix defines classical or distributional (weak) derivatives for functions on open subsets of R^d.It introduces the derivative notation before defining function spaces.Source: A. Sobolev spaces · p. 74
• Sobolev spaces W^{k,p}(Ω) are defined using seminorms based on derivatives, with a norm specified for p < ∞.The excerpt does not state the displayed formula for that norm.Source: A. Sobolev spaces · p. 74
• C^k(Ω) denotes functions that are k times continuously differentiable and is equipped with the W^{k,∞}(Ω) norm.Source: A. Sobolev spaces · p. 74