Autoregressive PDE Surrogates

• Autoregressive PDE surrogates are learned time-advance operators that emulate the evolution of fields governed by PDEs through recursive one-step predictions and explicit incorporation of physical structure.
• They employ strategies such as residual increment prediction and learned-derivative integration, adapting methods like CNNs, transformers, and graph neural networks to suit different spatial discretizations.
• Practical challenges include error accumulation in rollouts, which is mitigated through techniques such as curriculum learning, ensembling, and noise robustness to enhance long-horizon prediction stability.

Autoregressive PDE surrogates are learned time-advance operators that emulate the evolution of fields governed by partial differential equations by recursively applying a one-step predictor. In the standard formulation, the surrogate is trained to map the current state, or a short history of states, to the next state, and then rolled out by feeding predictions back as inputs. This paradigm appears in regular-grid CNN and transformer surrogates, graph neural simulators on irregular meshes, equivariant message-passing models, and coarse-to-fine generative solvers. Its appeal is procedural simplicity and explicit causal time marching; its central technical difficulty is that rollout quality is determined not only by one-step accuracy but also by how prediction errors accumulate, how physical structure is enforced, and how the surrogate interacts with spatial discretization and temporal resolution (Zhou et al., 2024, Bånkestad et al., 2024).

1. Formalism and problem setting

The canonical autoregressive formulation treats PDE evolution as repeated next-step prediction,

$\mathbf{u}_{n+1} = F_\theta(\mathbf{u}_n),$

with training performed on one-step pairs and inference performed by recursive rollout. On irregular geometric domains, the same idea is written as a time-marching surrogate of the form

$X^{t+1} = f_\eta(X^t; \ldots),$

with experiments explicitly distinguishing one-step error from rollout error accumulated over multiple autonomous steps (Zhou et al., 2024, Bånkestad et al., 2024).

A common refinement is residual or increment prediction rather than direct next-state regression. In small-data periodic 2D PDE forecasting, one implementation predicts

$$\hat{u}^{n+1} = \hat{u}^n + \Delta \hat{u}^{n+1}, \qquad \Delta \hat{u}^{n+1} = f_\theta(u^{n-L+1}, \dots, u^n),$$

using a fixed history length $L=7$. The training target is the increment $\Delta u^{n+1} = u^{n+1} - u^n$, not the full next field, and the paper states that this residual formulation empirically stabilizes long rollouts and encourages the model to focus on local-in-time dynamics (Nguyen et al., 13 Jan 2026).

A distinct but closely related formulation learns instantaneous time derivatives and couples them to an explicit numerical integrator. In graph neural simulators this appears as

$u^{t+1} = u^t + \Delta t \cdot \mathrm{GNS}(u^t),$

while derivative-prediction alternatives more generally train

$$\frac{\partial \mathbf{u}}{\partial t}\Big|_{t=t_n} = F_\theta(\mathbf{u}_n),$$

followed by Euler, Heun, RK4, or related ODE integration during inference. This does not remove autoregressive rollout at deployment, but it shifts the learned object from a discrete step map to a continuous-time or semi-discrete dynamical law (Nayak et al., 7 Sep 2025, Zhou et al., 2024).

The dominant misconception is that autoregressive PDE surrogates are synonymous with direct state prediction. The recent literature instead shows at least three technically distinct autoregressive regimes: direct state prediction, residual increment prediction, and learned-derivative prediction with explicit integration. This suggests that “autoregressive” is best understood as a rollout protocol rather than a unique training target.

2. Spatial representations and architectural families

Autoregressive PDE surrogates are strongly shaped by spatial discretization. On regular periodic grids, a notable design is the multi-channel U-Net with enforced periodicity, or me-UNet, which uses multi-channel temporal input, circular padding in all convolutions, average pooling, and a residual output for $\Delta u^{n+1}$. Under a common training setup across five PDE families on periodic $64 \times 64$ domains, me-UNet matched or outperformed ViT, AFNO, PDE-Transformer, and KAN-UNet in field-space error, spectral similarity, and physics-based metrics, while requiring substantially less training time (Nguyen et al., 13 Jan 2026).

On non-gridded domains, graph-based autoregressive surrogates replace stencil operations with message passing. “Flexible SE(2) graph neural networks with applications to PDE surrogates” constructs SE(2)-equivariant GNNs by aligning local node features to the $x$-axis, applying arbitrary nonlinear functions in the aligned frame, and rotating outputs back:

$X^{t+1} = f_\eta(X^t; \ldots),$0

The model separates scalar features $X^{t+1} = f_\eta(X^t; \ldots),$1 from rotational features $X^{t+1} = f_\eta(X^t; \ldots),$2, uses only relative positions $X^{t+1} = f_\eta(X^t; \ldots),$3, and builds edges on irregular domains by constructions such as Delaunay triangulation. Message functions can be instantiated either as SE2Conv-MLP or SE2Conv-Trans (Bånkestad et al., 2024).

A second graph-based family explicitly learns the time derivative and integrates it. In the GNS framework, edge messages and node updates take the form

$X^{t+1} = f_\eta(X^t; \ldots),$4

and the learned derivative is advanced by explicit Euler. In a strict low-data regime across 2D Burgers’ scalar, 2D coupled Burgers’ vector, and 2D Allen-Cahn systems, this framework achieved under $X^{t+1} = f_\eta(X^t; \ldots),$5 relative $X^{t+1} = f_\eta(X^t; \ldots),$6 errors with only 30 training samples out of 1000 across all three PDE systems (Nayak et al., 7 Sep 2025).

These results indicate that no single spatial architecture dominates across regimes. Periodic grids reward locality and boundary-condition alignment; irregular domains reward graph formulations and relative geometry; and equivariant constructions are particularly relevant when rotational structure is exact rather than approximate.

3. Training regimes and rollout stabilization

Autoregressive deployment creates a training–inference mismatch whenever optimization uses teacher forcing but inference uses self-fed predictions. “Learning Neural PDE Solvers with Parameter-Guided Channel Attention” makes this distinction explicit: teacher forcing uses ground-truth previous states as inputs, whereas autoregressive rollout feeds the model’s own predictions back into future steps. The paper proposes a curriculum learning strategy in which the first $X^{t+1} = f_\eta(X^t; \ldots),$7 steps of a sequence are run autoregressively and later steps switch to teacher forcing, with $X^{t+1} = f_\eta(X^t; \ldots),$8 increasing during training according to a smooth schedule. On multi-dimensional, non-linear systems like 2D Navier-Stokes, shifting to AR-only or TF-only training worsened accuracy by factors of $X^{t+1} = f_\eta(X^t; \ldots),$9 compared to curriculum learning (Takamoto et al., 2023).

Stabilization can also be achieved by altering what is aggregated at inference. “Enhanced accuracy through ensembling of randomly initialized auto-regressive models for time-dependent PDEs” trains multiple surrogates with identical architectures but different random initializations and averages their outputs at each step,

$$\hat{u}^{n+1} = \hat{u}^n + \Delta \hat{u}^{n+1}, \qquad \Delta \hat{u}^{n+1} = f_\theta(u^{n-L+1}, \dots, u^n),$$0

Across stress evolution in heterogeneous microstructures, Gray-Scott reaction–diffusion, and planetary-scale shallow water, the ensemble reported consistent reductions of $$\hat{u}^{n+1} = \hat{u}^n + \Delta \hat{u}^{n+1}, \qquad \Delta \hat{u}^{n+1} = f_\theta(u^{n-L+1}, \dots, u^n),$$1 in long-term relative error and mean absolute error compared to standalone models, with gains saturating after about 8 ensemble members (Khurjekar et al., 5 Jul 2025).

A different stabilization mechanism is noise robustness. In autoregressive deep surrogates for dendritic microstructure dynamics, training uses artificially perturbed inputs via additive Gaussian noise specifically to improve robustness to accumulated prediction errors. The same framework reports that noise injection is critical to stabilize long autoregressive rollouts (Ji et al., 5 Nov 2025).

In coarse-to-fine generative settings, recent work argues that long-horizon stability is governed by conditional-law errors rather than only by one-step reconstruction fidelity. “Memory-Conditioned Flow-Matching for Stable Autoregressive PDE Rollouts” introduces a compact online memory state $$\hat{u}^{n+1} = \hat{u}^n + \Delta \hat{u}^{n+1}, \qquad \Delta \hat{u}^{n+1} = f_\theta(u^{n-L+1}, \dots, u^n),$$2 into the conditional kernel,

$$\hat{u}^{n+1} = \hat{u}^n + \Delta \hat{u}^{n+1}, \qquad \Delta \hat{u}^{n+1} = f_\theta(u^{n-L+1}, \dots, u^n),$$3

and derives rollout bounds of the form

$$\hat{u}^{n+1} = \hat{u}^n + \Delta \hat{u}^{n+1}, \qquad \Delta \hat{u}^{n+1} = f_\theta(u^{n-L+1}, \dots, u^n),$$4

The paper reports improved accuracy, markedly more stable long-horizon rollouts, and better fine-scale spectral and statistical fidelity on compressible flows with shocks and multiscale mixing (Armegioiu, 6 Feb 2026).

4. Geometric and physical structure

A central theme in autoregressive PDE surrogates is that hard-wired geometric and physical structure can reduce rollout drift more effectively than purely data-driven next-state regression. In the SE(2)-equivariant GNN framework, translation equivariance is enforced globally by shifting the centroid of node positions to the origin and locally by using only relative positions; rotation equivariance is enforced by aligning vector features to a common axis, applying nonlinear transformations, and rotating back. The paper reports that this sidesteps many constraints while preserving SE(2) equivariance and enables arbitrary nonlinear functions inside the equivariant block (Bånkestad et al., 2024).

For CFD-style staggered grids, “Geometric and Physical Constraints Synergistically Enhance Neural PDE Surrogates” introduces new input and output layers that map between staggered-grid field representations and equivariant regular representations. The networks are constructed to satisfy

$$\hat{u}^{n+1} = \hat{u}^n + \Delta \hat{u}^{n+1}, \qquad \Delta \hat{u}^{n+1} = f_\theta(u^{n-L+1}, \dots, u^n),$$5

for discrete symmetry groups $$\hat{u}^{n+1} = \hat{u}^n + \Delta \hat{u}^{n+1}, \qquad \Delta \hat{u}^{n+1} = f_\theta(u^{n-L+1}, \dots, u^n),$$6, $$\hat{u}^{n+1} = \hat{u}^n + \Delta \hat{u}^{n+1}, \qquad \Delta \hat{u}^{n+1} = f_\theta(u^{n-L+1}, \dots, u^n),$$7, and $$\hat{u}^{n+1} = \hat{u}^n + \Delta \hat{u}^{n+1}, \qquad \Delta \hat{u}^{n+1} = f_\theta(u^{n-L+1}, \dots, u^n),$$8. Physical constraints are enforced as hard output-layer constructions rather than as loss penalties. For example, mass conservation is imposed by subtracting the mean of the predicted increment, and incompressibility can be imposed by learning a vector potential and constructing velocities as the discrete curl. Across shallow water equations with closed boundaries and decaying incompressible turbulence, symmetries were more effective than physical constraints, but surrogates with both performed best, even compared to baselines with data augmentation or pushforward training (Huang et al., 5 Jun 2025).

A still more structural approach is to predict transport operators rather than next states. “FluxNet: Learning Capacity-Constrained Local Transport Operators for Conservative and Bounded PDE Surrogates” advances each cell through neighborhood exchanges,

$$\hat{u}^{n+1} = \hat{u}^n + \Delta \hat{u}^{n+1}, \qquad \Delta \hat{u}^{n+1} = f_\theta(u^{n-L+1}, \dots, u^n),$$9

which guarantees exact global conservation on periodic grids. Lower- and upper-bound constraints are enforced by capacity-limited L-head and U-head parameterizations, and dual-bounded quantities use a D-head with a dual consistency loss. On shallow-water equations, FluxNet-LAP reported lowest rollout error by $L=7$0 over projection baselines, zero bound violations, and conservation error at $L=7$1; on traffic flow, FluxNet-D achieved machine-precision conservation with negligible bound-violation magnitudes below $L=7$2 (Lan et al., 2 Feb 2026).

The broader implication is not merely that constraints help, but that their implementation matters. The literature distinguishes soft penalties, post-hoc projection, and architectural enforcement, and the strongest reported rollout gains generally come from the latter two categories, especially when combined with symmetry.

5. Empirical regimes, benchmarks, and applications

Autoregressive PDE surrogates now span regular-grid, irregular-graph, multiscale, and materials applications. In irregular-graph fluid surrogates, SE(2) GNNs were trained on 2D Navier-Stokes plus advected scalar fields sampled as Delaunay graphs with 1024 random points from a $L=7$3 grid. Against non-equivariant models with otherwise identical architectures, the SE(2) models reached the one-step error of the baseline with one-eighth the data and showed much reduced error accumulation in autonomous rollouts, including obstacle and inlet configurations in which non-equivariant predictions degraded substantially (Bånkestad et al., 2024).

In strict small-data periodic 2D forecasting, me-UNet generalized qualitatively to unseen initial conditions with as few as approximately 20 training simulations. Across advection, diffusion, dislocation dynamics, Kolmogorov flow, and Gray-Scott reaction–diffusion, the reported metrics included RMSE, azimuthally averaged power-spectral-density cosine similarity, and physics-based invariant measures such as mass, energy, or dislocation density. The paper emphasizes that convolutional architectures with locality and periodic padding remain strong contenders in such regimes (Nguyen et al., 13 Jan 2026).

In data-efficient graph-based time stepping, the GNS framework reported average autoregressive error reductions of $L=7$4 relative to FNO AR and $L=7$5 relative to DON AR across three canonical PDE systems. The same study introduced a PCA+KMeans trajectory-selection strategy for low-data training, arguing that representative trajectory choice materially affects generalization under sample scarcity (Nayak et al., 7 Sep 2025).

Autoregressive surrogates have also been pushed into phase-field microstructure evolution. “Scalable Autoregressive Deep Surrogates for Dendritic Microstructure Dynamics” trains an SI-ConvNeXt surrogate on short trajectories from quantitative phase-field simulations in limited domains and then rolls it out over much larger spatial and temporal scales. The paper reports more than two orders of magnitude speed-up, normalized MSE of $L=7$6 for long rollouts, tip-dynamics deviations within about $L=7$7, and $L=7$8 for primary-spacing correlation. The authors state that no bottleneck or U-Net-style spatial compression is crucial for spatial extrapolation to large domains (Ji et al., 5 Nov 2025).

Together, these benchmarks show that empirical success depends on matching inductive bias to regime. Equivariance is particularly valuable when rotational symmetry is exact; periodic convolutions are effective on translation-invariant toroidal domains; graph simulators are highly data-efficient when local interactions dominate; and coarse-grained autoregressive surrogates can scale to long-horizon materials prediction when trained on short, high-fidelity trajectories.

6. Limitations, critiques, diagnostics, and adjacent formulations

The most persistent limitation of autoregressive PDE surrogates is error accumulation. “Predicting Change, Not States” argues that standard state prediction entangles spatial and temporal dynamics, fixes the temporal discretization to the training step size, encourages time-axis downsampling, and can suffer from flat, redundant loss landscapes when $L=7$9 is close to $\Delta u^{n+1} = u^{n+1} - u^n$0. The same study reports that derivative prediction, combined with numerical integration at inference, can improve stability and long-term accuracy in finely discretized regimes, while also permitting flexible time stepping and higher-order integrators (Zhou et al., 2024).

A more radical critique is that discrete autoregressive step maps are not the only way to learn time-dependent PDE dynamics. “CFO: Learning Continuous-Time PDE Dynamics via Flow-Matched Neural Operators” trains a neural operator to match analytic velocity fields obtained from spline-fitted trajectories, avoiding backpropagation through ODE solvers during training and then integrating the learned vector field at inference. CFO is described as time-resolution invariant, able to train on arbitrary non-uniform time grids, and capable of reverse-time inference and arbitrary temporal querying. Across Lorenz, 1D Burgers, 2D diffusion-reaction, and 2D shallow water, CFO trained on only $\Delta u^{n+1} = u^{n+1} - u^n$1 of irregularly subsampled time points outperformed autoregressive baselines trained on complete data, with relative error reductions up to $\Delta u^{n+1} = u^{n+1} - u^n$2 (Hou et al., 4 Dec 2025).

Another limitation concerns genuine generalization. “Generalization vs. Memorization in Autoregressive Deep Learning” applies influence-function analysis and reports rapid temporal decay of off-diagonal influence, near-complete absence of gradient overlap between different initial-condition classes, and similar qualitative failures under both standard and physics-informed losses. Its central claim is that validation loss can mask fragmentation and memorization in autoregressive PDE surrogates, and it proposes time-aware influence kernels and irreducible-representation summaries as diagnostics of trust horizon and cross-scenario coherence (Amarel et al., 18 Aug 2025).

These critiques do not imply that autoregression is obsolete. Rather, they clarify its operating envelope. A plausible synthesis is that autoregressive PDE surrogates remain highly effective when their rollout protocol is paired with suitable structure: geometry-aware architectures, conservation-aware outputs, explicit or implicit control of training–inference mismatch, memory or ensemble mechanisms for long horizons, and diagnostics that test more than pointwise one-step accuracy.