Preconditioned Physics-Informed Neural Operator Training

1 INTRODUCTION

• Solver-generated datasets are a substantial offline cost, while finite datasets cap attainable accuracy and must be revisited during training.Physics-informed training can instead draw fresh samples each step and adapt to a new input distribution or geometry without new data.Source: §1 · pp. 1–2
• Differential operators amplify high frequencies, and the physics-informed loss Hessian has condition number growing like h^-4, making training harder as the mesh is refined.The supervised loss Hessian is the identity, while the physics-informed least-squares Hessian is A⊤A.Source: §1 · p. 2

1.1 MAIN CONTRIBUTIONS

• The proposed preconditioned loss reduces conditioning from O(h^-4) to O(1) for elliptic problems and to O(h^-2) for saddle point problems.Multigrid implementation adds negligible training cost; the loss is architecture agnostic and adds no inference cost.Source: §1.1 · p. 3
• Across Poisson, Allen–Cahn, and stationary Stokes, label-free training closely matches or surpasses supervised training; PINO and physics-informed DeepONet are four to twenty-five times less accurate.For Poisson, fresh samples at every step match supervised training on 4096 labelled solutions.Source: §1.1 · p. 3

1.2 RELATED LITERATURE

• Earlier neural-operator physics-informed methods used strong-form residuals, reweighting, sampling, architecture changes, labelled data, or finite-element residuals, but did not address PDE-operator ill-conditioning.The cited approaches include Li et al. (2024), Wang et al. (2021b), Chen et al. (2026), Medvedev et al. (2026), Yamazaki et al. (2025), and Wen et al. (2026).Source: §1.2 · p. 3

2 METHODS

• A neural operator approximates a map between function spaces, and its interpolated output is evaluated through a discretized PDE residual.The discrete problem seeks a finite-element solution whose residual is zero for each PDE input.Source: §2 · pp. 3–4
• The interpolation space and operator are modeling choices; interpolation belongs to the loss, not the model, and does not tie the operator to a fixed resolution.The applications use Lagrange interpolation and Lagrange finite-element spaces.Source: §2 · p. 4

2.1 PHYSICS-INFORMED LOSS FUNCTIONS

• Preconditioning changes the loss landscape but preserves its minimizers, and the loss-only construction is architecture agnostic with no additional inference cost.Invertibility of P and positive definiteness of B ensure the preconditioned and unpreconditioned integrands vanish at the same solutions.Source: §2.1 · p. 4
• For Poisson, the unpreconditioned choice (P, B) = (I, I) has O(h^-4) conditioning, while (P, B) ≈ (A^-1, M) yields O(1).The alternative (I, A^-1) gives O(h^-2) conditioning and requires A to be symmetric positive definite.Source: §2.1 · p. 4

3 RESULTS

• Preconditioned training matches supervised training on Poisson, trails it slightly on Allen–Cahn, and surpasses it on Stokes; physics-informed baselines are less accurate on Poisson and Stokes.PINO comes within a factor of two on Allen–Cahn, while physics-informed DeepONet reaches 90% error; per-epoch overhead ranges from under 10% to about a third.Source: §3 · p. 5

3.1 POISSON EQUATION

• 2.1% relative L2 error with infinite-data Poisson training beats the 2.8% error from the largest finite dataset, while finite datasets through 2048 are indistinguishable at this budget.In a preconditioner sweep, geometric multigrid performs virtually identically to data-driven training; diminishing returns appear below conditioning of about 10^3–10^4.Source: §3.1 · pp. 6–7, Figure 3 · p. 6, Figure 2 · p. 6

3.2 ALLEN-CAHN EQUATION

• On Allen–Cahn, preconditioning accelerates convergence and improves accuracy over the baselines, nearly matching data-driven training; physics-informed roll-outs remain stable beyond the training horizon.The training horizon is t = 0.1, and the roll-out is evaluated through T = 1.Source: §3.2 · p. 7, Figure 4 · p. 8

3.3 STOKES EQUATION

• 1.6% velocity and 0.85% pressure error with preconditioned training beat supervised errors of 3.3% and 1.8%, respectively, on Stokes.The unpreconditioned least-squares loss stalls at 30.6% velocity and 11.4% pressure error; the preconditioner leaves O(h^-2) conditioning rather than O(h^-4).Source: §3.3 · pp. 8–9

4 CONCLUSION AND LIMITATIONS

• The method matches or closely approaches supervised training on Poisson, Allen–Cahn, and stationary Stokes equations across structured and unstructured grids.The authors identify efficient preconditioners for the problem class as the main limitation, particularly for convection-dominated and high-frequency wave problems.Source: §4 · p. 9

AI USE STATEMENT

• The work uses generative AI for coding, writing, and mathematical reasoning assistance, with authors reviewing and verifying its output.The authors take responsibility for the final content, including AI-assisted text, claims, and artifacts.Source: AI USE STATEMENT · p. 9

4.1 REPRODUCIBILITY

• The authors provide code to reproduce the experiments in the TensorPILS GitHub repository.The passage identifies the repository as the location of the reproduction code.Source: §4.1 · p. 10

• The interpolated operator can be evaluated at any resolution because interpolation is part of the loss, not the model.The training resolution still limits accuracy because features smaller than h cannot be resolved; the passage notes a similar limitation for data-driven labels.Source: A DETAILS ON INTERPOLATED NEURAL… · p. 12
• A triangle-inequality analysis separates the interpolated operator’s error into network error and interpolation error.The discussion defines the true operator, neural surrogate, and interpolated neural operator before deriving this restriction on achievable accuracy.Source: A DETAILS ON INTERPOLATED NEURAL… · p. 13
• For the 65 × 65 grids, the interpolation error floor is of order 10^-4 in L2, below all measured errors.For nodal interpolation into the stated finite-element spaces, the interpolation error scales with h^2; the passage gives h^2 = 64^-2 ≈ 2.4 · 10^-4.Source: A DETAILS ON INTERPOLATED NEURAL… · p. 13
• Preconditioned and supervised operators trained at 65 × 65 degrade alike when evaluated from 33 × 33 to 257 × 257.The passage attributes the resolution-related accuracy boundary to the training resolution, rather than either method being uniquely tied to that grid.Source: A DETAILS ON INTERPOLATED NEURAL… · p. 13

• Least-squares preconditioning adds no inference cost because the preconditioner is part of the loss, while allowing a choice of preconditioner and applying to general PDEs.The choice remains PDE-dependent; the section introduces these advantages while comparing available approaches for nonlinear neural-operator outputs.Source: B THE DIFFERENT OPTIONS OF… · p. 13
• The discussion specializes to quadratic losses and learning a solution operator, with a positive-definite symmetric A additionally assumed for the variational loss.The losses are considered on a single sample, and u(θ) is interpreted as an interpolated neural operator.Source: B THE DIFFERENT OPTIONS OF… · p. 13

• Up to the Jacobian, the first-order optimality conditions are the normal equations of the left-preconditioned system PAu = Pf in the B-inner product.The passage defines J(θ) as the Jacobian of θ ↦ u(θ) and gives the associated Gauss–Newton matrix.Source: §B.1 · p. 14
• The lemma quantifies the Gauss–Newton matrix’s conditioning up to the influence of the Jacobians, assuming B is positive definite and symmetric and A and P are invertible.A further estimate is stated when A is also symmetric.Source: §B.1 · p. 14
• For the symmetric-A case, the proof uses a matrix similar to BA; its symmetry lets the authors express its condition number through its eigenvalues.The proof also invokes submultiplicativity and invariance under transposition for the earlier estimate.Source: §B.1 · p. 14

B.2 PRECONDITIONING THROUGH ARCHITECTURE ADAPTION

• A non-trainable final layer maps the network output through a bijective linear preconditioner P, which must also be applied at inference.This architecture modification is used to precondition the residual loss.Source: §B.2 · p. 14
• For least squares, P ≈ A^-1 mitigates ill-conditioning, while the variational formulation calls for P ≈ A^-1/2.The Gauss–Newton matrix takes the form J(θ)⊤[P⊤AP]J(θ), linking the choice of P to the conditioning of the loss.Source: §B.2 · p. 15
• Related approaches precondition PINN-type parametric problems with BPX folded into the ansatz or Cholesky-derived block factors for Stokes and Navier–Stokes.The cited examples concern low-dimensional parameter spaces.Source: §B.2 · p. 15

B.3 PRECONDITIONING THROUGH CHANGE OF INNER PRODUCT

• Changing the inner product gives a left-preconditioned variational condition, and a good P reduces the condition number of PA.For example, P ≈ A^-1; the construction requires A and P to be symmetric positive definite and to commute.Source: §B.3 · p. 15
• The commutation requirement AP = PA is restrictive: polynomial preconditioners in A satisfy it, but multigrid methods do not.Source: §B.3 · p. 15

B.4 PRECONDITIONING IN PARAMETER SPACE

• Parameter-space preconditioning updates parameters using Pk∇L(θk), with Newton, natural gradient, and Adam as examples.Scaling such techniques to large networks requires specialized methods, and the paper considers this approach orthogonal to its own.Source: §B.4 · p. 15

Data Set

• Source–solution pairs are constructed from Laplacian eigenfunctions, so they satisfy the Poisson equation exactly without requiring a PDE solve.The frequency parameter K controls the number of frequencies; larger K makes the learning problem more challenging.Source: C DETAILS OF THE POISSON EQUATION… · p. 16

Discretization

• The problem is discretized with continuous piecewise bilinear finite elements on a uniform grid, using nx = 65 and n = 3969 interior degrees of freedom.The domain is the unit square with homogeneous Dirichlet boundary conditions.Source: C DETAILS OF THE POISSON EQUATION… · pp. 15–16

Neural Operator

• The FNO output grid coincides with the finite-element nodes, and projecting its output to zero on the boundary enforces the boundary condition exactly.Thus the model output lies in the finite-element space without a boundary penalty.Source: C DETAILS OF THE POISSON EQUATION… · p. 16

Baselines

• The PINO baseline minimizes a strong-form residual with second-order central differences, while PI-DeepONet uses automatic differentiation and a soft boundary penalty.PI-DeepONet runs remain descending at the end of the fixed budget, so its reported values do not represent a converged state.Source: C DETAILS OF THE POISSON EQUATION… · pp. 16–17

Multigrid Preconditioner

• The Poisson preconditioner is one geometric multigrid V-cycle with weighted Jacobi smoothing, bilinear prolongation, and a direct solve on the coarsest grid.It uses dyadic coarsening on grids of fineness 2k + 1 to ensure nested grids.Source: C DETAILS OF THE POISSON EQUATION… · p. 17

Optimization

• Training uses Adam and cosine annealing, with learning rates selected per loss by validation error over a 500-epoch budget.Final configurations are repeated over three seeds, and test error is reported at the checkpoint with the lowest validation error.Source: C DETAILS OF THE POISSON EQUATION… · p. 17, Table 2 · p. 17, Table 3 · p. 18

Spatial Discretization

• The experiment uses homogeneous Dirichlet conditions on the unit square and bilinear finite elements on a 129 × 129 grid.The initial conditions use the Poisson dataset with K = 4.Source: D DETAILS OF THE ALLEN–CAHN EXPERIMENT · p. 18

Time-stepping operator

• The learned operator advances the Allen–Cahn solution by a fixed step τ = 0.01 using a first-order, unconditionally stable convex–concave splitting scheme.The solution is represented through iterates approximating the solution semigroup.Source: D DETAILS OF THE ALLEN–CAHN EXPERIMENT · p. 18

Preconditioner

• The preconditioned least-squares loss applies P to the residual, and its Gauss–Newton matrix in w is J(w)⊤P⊤PJ(w).The residual Jacobian depends on the state and is not symmetric, so the method freezes it at w² = 1 to obtain a symmetric positive-definite, state-independent matrix.Source: D DETAILS OF THE ALLEN–CAHN EXPERIMENT · p. 19
• The frozen-Jacobian multigrid hierarchy is built once before training, and the unpreconditioned least-squares loss sees κ(J)².The analysis bounds the frozen matrix using the mass and stiffness matrix scalings and uses hε = 1/4 to resolve interfaces.Source: D DETAILS OF THE ALLEN–CAHN EXPERIMENT · p. 19

Conditioning

• Conditioning of O(h^-2) remains in theory, but the time-step size τ strongly counteracts it in practice.Source: D DETAILS OF THE ALLEN–CAHN EXPERIMENT · p. 19

Operator and discretization

• The learned map sends body force to velocity-pressure pairs on a unit square with an off-centre circular hole, using an unstructured Taylor–Hood P2/P1 mesh.The mesh has 9031 degrees of freedom, with 568 velocity components constrained.Source: E DETAILS OF THE STOKES EQUATION… · pp. 19, 21

Data set

• Body-force components are truncated sine series; reference fields are discrete Taylor–Hood solutions, with pressure pinned at one node and then regauged to zero mean.The reference factorization is shared across samples, and its relative residual is 4.4 · 10^-6.Source: E DETAILS OF THE STOKES EQUATION… · pp. 19, 21

Neural operator and baseline

• The experiment uses GAOT for Stokes; PI-DeepONet is the only tested literature baseline that transfers to the unstructured mesh, but its second-derivative evaluation takes 6.0 seconds per epoch.The stochastic-collocation variant takes 1.5 seconds per epoch for accuracy within the seed spread.Source: E DETAILS OF THE STOKES EQUATION… · p. 21

Preconditioner

• The loss uses block-diagonal P = diag(K^-1, S^-1): one algebraic multigrid V(2, 2) cycle per velocity component approximates A^-1, while the lumped pressure mass matrix approximates S^-1.Source: E DETAILS OF THE STOKES EQUATION… · p. 21

Preconditioner ablation

• Removing P worsens velocity accuracy by a factor of 20; at ωS = 256, pressure error is 3.9% while velocity error is 22.9%.The unpreconditioned residual has κ(A2) = 1.2·1011, five orders of magnitude worse than the best preconditioned value.Source: E DETAILS OF THE STOKES EQUATION… · p. 22, Table 6 · p. 22

Qualitative results

• LPLS and supervised outputs visually match reference velocity and pressure; the unpreconditioned control loses recirculation amplitude, while PI-DeepONet captures only the coarsest structures.The unpreconditioned pressure error forms a ring around the obstacle.Source: Figure 9 · p. 23