PaperlayerView exact revision on arXiv

Source-linked AI summary

A Constitutive Markov Physics-Informed Neural Operator (MPNO) for Autoregressive Stability in Transient Dynamics

Wenpu Du, Peng Zhou, Yunlong Xia, Sinuo Xin, Congcong Zhang, Boyang Zhang, Yi Zhang, Wenzheng Xu

arXiv:2608.25744v1cs.LG

TL;DR

Transient-dynamics neural operators can become unstable during autoregressive rollout because their propagation operators lack a structural spectral-radius constraint. MPNO constructs a constitutive Markov propagation operator with spectral radius at most 1, achieving bounded stable rollouts across three PDE scenarios while using about one quarter of FNO’s parameters.

  • Problem

    Existing neural operators lack structural spectral-radius constraints, leaving autoregressive stability dependent on learned dynamics and allowing rollout errors to diverge or collapse.

  • Method

    MPNO embeds constitutive material information in a graph-based Markov propagator whose normalized Laplacian construction constrains its spectral radius to ρ(P) ≤ 1.

  • Results

    MPNO achieves 0.7304±0.0008 single-step Rel-L2 error and stable bounded-error autoregressive rollouts across Burgers, Darcy, and concrete penetration, using about 20K parameters.

  • Takeaways & Limitations

    The results support architecturally enforced spectral stability for graph-topology systems governed by local conservation laws whose evolution can be approximated as a Markov process.

  • Takeaways & Limitations

    The spectral-radius guarantee applies only to the linear propagation operator, so full-model long-horizon stability is not directly guaranteed and is confirmed empirically.

Abstract

from arXiv · show

Neural operators applied to transient-dynamics PDEs with strong discontinuities exhibit autoregressive instability: in concrete-penetration stress-field prediction, the wavelet neural operator (WNO) diverges in autoregressive rollout, while MeshGraphNets collapse to zero predictions. WNO's instability stems from the lack of a structural constraint on the spectral radius of its propagation operator; the Fourier neural operator (FNO) is stable in these measurements but only emergently, not by construction. We propose a constitutive Markov physics-informed neural operator (MPNO) modeling one-step evolution as a Markov (row-stochastic) propagation operator. Physics-coupled edge weights (acoustic-impedance harmonic mean, contact area, and traction amplitude) encode material-interface constitutive information into a nonnegative symmetric adjacency matrix W; after normalizing the graph Laplacian L = D - W by lambda_max, the propagator P = I - alpha*L~ is constructively constrained to spectral radius rho(P) <= 1, suppressing exponential amplification of autoregressive errors. Stability is thus a designable architectural property, not an optimized loss objective. On three PDEs (Burgers and two-dimensional transverse-section concrete penetration), MPNO rolls out stably with bounded error on all test seeds at 100/135/165 m/s; the single-step relative L2 error is 0.7304 +/- 0.0008, better than WNO and comparable to FNO at about one quarter of FNO's parameters. The edge-weight formula transfers across scenarios by replacing material-property variables. With about 20K parameters, MPNO delivers roughly 10^5x inference speedup over LS-DYNA.

1. Introduction

MPNO addresses autoregressive instability by constructing spectral-radius control into a physics-coupled Markov propagation architecture rather than optimizing stability through the loss. Its constitutive edge weights and online spectral normalization support bounded rollouts and transfer across PDE scenarios.

  • Motivation: WNO diverges during autoregressive concrete-penetration rollout despite convergence, while FNO remains stable without a spectral guarantee and MeshGraphNets fail to learn.The failure persisted after varying more than thirty engineering factors, including learning rate, batch size, depth, width, losses, normalization, and scheduled sampling.
  • Core principle: MPNO makes bounded autoregressive error a constructive architectural property by modeling one-step evolution as P = I − αL̃ with controlled spectral radius.The approach separates stability from loss optimization and suppresses exponential error amplification through the propagation operator’s spectral structure.
  • Spectral control: Online Rayleigh-quotient power iteration estimates λmax at O(E) cost, enabling L̃ = L/λmax and constraining the Markov propagator to ρ(P) ≤ 1.The dependency chain is physics-built W → graph-Laplacian properties → normalized L̃ → controlled spectral radius.
  • Architecture: The architecture builds W from acoustic-impedance harmonic mean, contact area, traction amplitude, and a zero-initialized MLP correction, yielding a nonnegative symmetric graph structure.This constitutive weighting embeds material-interface impedance information into propagation as a physical skeleton transferable by replacing material-property variables.
  • Validation and scope: MPNO is validated on Burgers, Darcy, and concrete-penetration PDEs for autoregressive stability and cross-scenario transferability, using about 20K parameters.Its spectral constraint depends on nonnegative symmetry of W and graph connectivity rather than a specific PDE family.

2. Problem Setting and Data Generation

The section defines three PDE benchmarks and their dataset-generation procedures: analytic or high-accuracy numerical solutions for Burgers and Darcy, and LS-DYNA simulations for concrete penetration. It specifies the sampling, discretization, and evaluation setups used to validate MPNO across these problems.

  • Concrete penetration: The concrete-penetration dataset uses LS-DYNA explicit-dynamics simulations with random-aggregate placement in a ϕ500 × 200 mm cylindrical domain.The aggregate distribution is three-dimensional crushed granite with 42% volume fraction and 5–25 mm sizes, generated using Fuller gradation with n = 0.5.
  • Burgers equation: 100 Burgers initial conditions use random sine-wave superpositions on a 128-point periodic grid, with 30 training cases and 70 test cases.The viscosity is ν = 0.01, and t ∈[0, 1] is discretized into 50 steps; five rollout-stability seeds come from the test set.
  • Darcy equation: 100 random Darcy permeability fields form the training set, while 20 fields form the test set on a 64 × 64 grid.K(x) is an FFT-generated Matérn Gaussian random field with τ = 3 and values in [0.1, 10], and pressure is solved by a five-point finite-difference scheme.
  • Concrete penetration: 100/135/165 m/s are the concrete-penetration evaluation velocities, while 200 m/s is used only for training-data augmentation.The 200 m/s samples enrich training with higher-amplitude, steeper-gradient regimes but are excluded from evaluation because adjacent frames are highly similar.

3. The Constitutive Markov Physics-Informed Neural Operator (MPNO) · 3.1. Problem formulation: from continuum to discrete graph

MPNO reformulates transient-dynamics prediction as a Markov process on a spatial graph rather than a map between continuous function spaces. It constructs graph-based node states, physics-coupled adjacency weights, and a propagation framework whose spectral guarantees depend on nonnegative symmetry and connectivity rather than a specific mesh topology.

  • 3. The Constitutive Markov Physics-Informed Neural Operator (MPNO): MPNO models transient dynamics as a Markov process on a spatially discretized graph, defining topology, node-state vectors, adjacency weights, and propagation.This departs from the conventional neural-operator formulation as a map between continuous function spaces.
  • 3.1. Problem formulation: from continuum to discrete graph: The general framework covers both 3D hexahedral and 2D transverse-section configurations, while the numerical experiments instantiate the latter.The 2D configuration uses ±x and ±y adjacency, with z as the out-of-plane penetration direction and single-frame 10-dimensional node features.
  • 3.1.1. Mesh-to-graph topology mapping: A hexahedral mesh becomes G = (V, E, W), with one node per cell, edges for shared contact faces, and nonnegative weights on connected pairs.The graph has N = X × Y × Z nodes; interior cells have six neighbors and the undirected graph has E ≈3N edges.
  • 3.1.1. Mesh-to-graph topology mapping: Six-face adjacency is dictated by structured hexahedral meshes and the six independent components of the Cauchy stress tensor, rather than chosen as a free architecture parameter.The six face tractions jointly represent the stress tensor, while fewer faces underdetermine it and more add redundancy.
  • 3.1.1. Mesh-to-graph topology mapping: The graph construction yields a Laplacian whose diagonal entries measure incident coupling strength and whose analyzed properties include symmetry and positive semidefiniteness.The Laplacian is formed from the adjacency weight matrix through D and L = D − W.
  • 3.1.2. Cauchy traction decomposition and node states: Each cell stores six Voigt stress components, which map through t = σ·n to 18 face-traction components; delayed states and coordinates augment the node representation.The general 3D augmented state is 56-dimensional, whereas the evaluation configuration uses single-frame 10-dimensional features.
  • 3.1.3. Edge orientation and face pairing: Face-pairing indices encode the reverse contacting faces of each edge and are required because traction-amplitude coupling depends on the face carrying the force.Each edge therefore carries face-to-face contact information used by the physics-coupled edge weights.
  • 3.1.4. Generality of the framework: The framework’s spectral conclusions require only nonnegative symmetry of W and graph connectivity, extending the mathematical foundation to tetrahedral and polyhedral meshes.Other mesh types still require adapted graph topology or weight interpolation.

3.2. Graph-Laplacian propagation and spectral stability

MPNO constructs graph propagation from a normalized Laplacian so its Markov operator has bounded spectral radius and nonexpansive rollouts. Online power iteration estimates the normalization scale, while residual propagation preserves local features across multiple rounds.

  • Online spectral estimation: K = 20 yields a relative error of about 5 × 10−6, corresponding to an underestimate of 0.0005% in the online λmax estimate.Power iteration uses only matrix–vector products and estimates λmax at O(E) cost per training batch.
  • Spectral stability: P = I − αL~ has eigenvalues in [0, 1], giving ρ(P) = 1 and the nonexpansive rollout bound ∥v_t∥_2 ≤ ∥v_0∥_2.The bound follows from α ∈ [0, 1] and the normalized Laplacian spectrum L~λ_i ∈ [0, 1].
  • Spectral stability: ρ(P) ≤1 is carried by constructing P = I − αL~, rather than relying on a spectral-radius loss penalty.The architecture normalizes against an online estimate of λmax, making stability a structural property.
  • Multi-round propagation: R = 8 propagation rounds cost O(8 · E · d_hidden), about 45% of a single-step inference through sparse matrix–vector products.The propagation uses node-degree weighting and sparse edge aggregation.
  • Multi-round propagation: β = 0.3 initial residual connections preserve local node features after 8 rounds while supporting long-range wave propagation and local material scattering.Each round applies an APPNP-style residual to the initial encoder embedding.

3.3. Spectral-radius criterion for autoregressive collapse

Autoregressive error is governed by the spectral radius of the effective propagation operator: values above one cause divergence, while MPNO enforces the critical stable regime architecturally rather than relying on training. This explains why MPNO remains bounded in rollout while WNO diverges and FNO is stable without a structural guarantee.

  • Teacher forcing and autoregressive rollout: Teacher-forced training avoids accumulated errors, but autoregressive inference feeds predictions back into the model, causing errors to compound through nonlinear dynamics.This training–inference distribution shift is known as exposure bias.
  • Spectral-radius dichotomy: ρ(P) > 1 causes catastrophic divergence for almost every initial error, whereas ρ(P) < 1 drives errors to zero; ρ(P) = 1 is the critical boundary.The controlled regime requires architectural control because uncontrolled dissipation or amplification eventually causes rollout failure.
  • Constructive stability of MPNO: MPNO’s P = I − α ˜L satisfies ρ(P) = 1 exactly, placing its propagation operator on the critical boundary independently of training.The unit eigenvalue corresponds to the translation-invariant mode of global conservation.
  • Baseline behavior: 0.7422 Rel-L2: WNO diverges on all seeds at the 10^12 scale, while 0.7210 Rel-L2: FNO remains stable on clean test seeds without a structural spectral-radius guarantee.Neither Fourier nor wavelet convolution layers constrain the spectral radius of the implicit Jacobian.
  • Experimental evidence: 0.7304 ± 0.0008 single-step Rel-L2: MPNO rolls out across 29 steps with bounded error of about 1.0, while λmax normalization keeps its propagator within the stability bound.Removing normalization yields αλmax values of 4–7, exceeding the stability bound of 2.

3.4. Physics-coupled edge weights

MPNO constructs edge weights from a physics-motivated product skeleton with a lightweight neural correction, rather than learning them purely from data. The design incorporates interface coupling, stress transfer, symmetry, and zero-initialized corrections to preserve a physics-driven starting point.

  • Physics-coupled construction: MPNO defines each edge weight as a physical formula skeleton multiplied by a lightweight neural correction, providing an inductive bias absent from purely data-driven edge functions.The skeleton is inspired by acoustic-impedance coupling, traction decomposition, and contact geometry.
  • Symmetrization: Symmetrizing opposing directional weights makes W nonnegative and symmetric, ensuring L = D - W is positive semidefinite and preserves the spectral conclusions.For boundary edges, the original one-sided weight is retained without halving.
  • Physics-coupled construction: The construction combines acoustic-impedance coupling, contact area, traction magnitude, and an MLP correction to represent force-transfer strength across interfaces.The MLP captures residual effects including multiple scattering, local plastic-deformation impedance changes, and mesostructural irregularities.
  • Initialization: Zero initialization makes the correction exp(MLP(·)) equal 1 initially, so training begins with purely physics-driven edge weights before learning residual corrections.This prevents a randomly initialized neural component from overriding the physical prior early in training.
  • Ablation: Ablation E1 replaces the physics-coupled weights with a randomly initialized pure-MLP alternative to test whether material-property and traction coupling is necessary under sparse penetration data.The stated expectation is that pure-MLP weights overfit and fail to generalize to unseen material distributions.

3.5. Row stochasticity and total-mass conservation

The propagation operators provide algebraic unweighted total-mass conservation in addition to stability, with row-stochastic and Laplacian-filter forms offering complementary viewpoints. This conservation is operator-level rather than physical field conservation and remains valid through shock discontinuities because it requires no spatial differentiation.

  • Stability and conservation: The conservation property complements the non-divergence guarantee by making 1 an eigenvector with λ = 1 for the propagation operator.This algebraic property applies to the diffusion-form propagation dynamics and preserves the unweighted total mass.
  • Algebraic conservation: Both propagation forms conserve unweighted total mass: D^-1W satisfies P1 = 1, while I − αL~ satisfies 1^T P = 1^T.The two forms are complementary: the Markov matrix is generally asymmetric and row-stochastic, whereas the Laplacian-filter form is real symmetric.
  • Algebraic conservation: 1.2×10−8 operator-level verification confirms the row-stochastic conservation identity numerically.The invariant concerns the latent state and is not transmitted pointwise through the nonlinear encoder–decoder map.
  • Shock-discontinuity immunity: Every-step invariance is immune to shock discontinuities because constructing P as D^-1W or I − αL~ requires no spatial differentiation.The conserved quantity is the unweighted sum, not weighted physical momentum or mass.
  • Power-iteration convergence: 0.0005% mean underestimate is achieved by the Rayleigh-quotient power iteration for λmax, whose convergence yields the realized constraint ρ(P) ≤ 1.The iteration converges at rate O((λN−1/λmax)2k) and always underestimates λmax.

3.6. Temporal evolution

MPNO advances temporal states by encoding the current frame, applying a Markov propagation operator in latent space, and decoding the next frame’s absolute stress directly. Its 8-iteration propagation and skip-1 autoregression avoid drift accumulation associated with incremental prediction schemes.

  • Temporal evolution: MPNO evolves each time step by encoding the current frame, propagating its latent state with P, and decoding the next absolute stress field.The decoder directly outputs absolute stress at the next time step.
  • Temporal evolution: 8 propagation iterations are used in practice before decoding the next frame’s six-component Voigt stress at each node.The decoder uses two linear layers, 32 →32 →6, with a ReLU between them.
  • Temporal evolution: Skip-1 autoregression feeds the predicted frame back as input while predicting absolute stress σt+1 rather than increment ∆σ, avoiding incremental-scheme drift accumulation.Velocity information is incorporated into the edge weights through the physics-coupled propagation mechanism.

3.7. Loss function design

MPNO uses a hard-masked MSE to address sparse, discontinuous stress fields, focusing training on active stress-wave regions rather than near-zero nodes. The joint loss consists solely of this data loss.

  • Loss motivation: Spatial sparsity makes standard all-node MSE vulnerable to learning the trivial output-zero solution.Most grid points have normalized stress |σ| < 10^-3, while stress-concentration zones occupy only a few points.
  • Hard-masked MSE: MPNO computes MSE only over active nodes whose ground-truth stress norm exceeds τ = 0.01.This hard mask prevents zero-stress nodes from dominating the gradient.
  • Hard-masked MSE: About 30% of the trimmed-impact-zone subgraph is active, concentrating the training signal on stress-wave-dominated regions.The same masking convention is used during evaluation, with frames included only when they contain at least 5 active nodes.
  • Joint loss: The joint loss is simply the masked data loss, L = Ldata.The method abandons pointwise PDE residuals that fail at strong discontinuities.

3.8. Gradient propagation and training stability … 3.11. Architecture overview

MPNO is designed for stable gradient propagation and efficient inference through a spectrally constrained, sparse graph propagator with about 20K parameters. Its modular ablations and architecture explicitly test and enforce the spectral constraint underlying the method.

  • 3.8. Gradient propagation and training stability: ∥P^t∥2 = 1 prevents exponential gradient decay or amplification during backpropagation through time by algebraic construction rather than training.Six neighbors per node localize gradient flow along edges, preserving physical locality and smoothing the loss landscape.
  • 3.8. Gradient propagation and training stability: 500 training epochs require no gradient clipping for MPNO, while PINO’s gradient variance is two orders of magnitude higher.The Rayleigh-quotient power iteration runs inside torch.no_grad(), excluding spectral-estimate sampling error from the training signal.
  • 3.9. Complexity and efficiency: 20K trainable parameters result from zero-parameter physics-coupled priors, a parameter-free graph-Laplacian propagator, and compact encoder, post-processing, and decoder components.The MLP correction has only a few hundred parameters.
  • 3.9. Complexity and efficiency: O(E · dhidden) = O(N · dhidden) single-step inference scales linearly with node count, while eight sparse propagation rounds form the main bottleneck.Power iteration estimates λmax in O(5E), below the O(N3) cost of full eigendecomposition.
  • 3.10. Component separability and ablation verification: Each component can be disabled independently, including physics-coupled edge weighting, spectral normalization, propagation post-processing, velocity coupling, and temporal-output format.This predefined scheme is intended to ensure falsifiable component comparisons rather than favorable post hoc selection.
  • 3.10. Component separability and ablation verification: E3 removes λmax scaling from P = I −αL under controlled settings to test whether the predicted spectral-radius dichotomy produces collapse or divergence.All ablations use SEED=42, Adam lr=3 × 10−3, hidden=32, R=8, hard-masked MSE, and absolute-stress targets.
  • 3.11. Architecture overview: ρ(P) ≤1 bounds the online λmax estimate and completes normalization with ˜L = L/λmax in the full encoding-to-propagation-to-decoding flow.The architecture takes current-frame stress σt, aggregate volume fraction Vf, and spatial coordinates x on the impact-zone trimmed subgraph.

4. Numerical experiments

Numerical experiments show that MPNO provides constructively constrained propagation and stable autoregressive prediction for transient impact dynamics, while unconstrained baselines can diverge or collapse. Across benchmark settings, MPNO achieves competitive single-step accuracy, bounded long-horizon error, and verified conservation and spectral-radius properties.

  • Scope and limitations: On the steady-state Darcy problem, MPNO, FNO, and WNO achieve test Rel-L2 of 0.558, 0.053, and 0.434, respectively, with MPNO weakest.The experiments frame MPNO’s core value as cross-scenario generalization of its physics-coupling formula rather than steady-state accuracy.
  • Spectral verification: ρ(P) = 1 ≤1 is confirmed numerically, supporting the intended spectral-radius constraint of the propagation operator.Power iteration estimates λmax at O(E) cost on the trimmed impact-zone graph; offline analysis uses K = 20 for a stricter estimate.
  • Conservation verification: 1.2 × 10−8 is the mean of ∥P1 −1∥ over 10 random aggregate distributions, indicating row sums remain 1 to floating-point precision.The construction therefore approximately conserves global total mass at every autoregressive step.
  • Penetration benchmarks: WNO diverges to infinity on all penetration-data seeds, whereas MeshGraphNets collapse to zero prediction; FNO remains stable in measurement but lacks a constructive spectral guarantee.The observed baseline behaviors match the predicted spectral-radius dichotomy, while MeshGraphNets’ decay is attributed to underfitting without learning.
  • Penetration benchmarks: MPNO remains stable over 29 autoregressive steps across all six test seeds and 100/135/165 m/s cases, with error about 0.98 at step 29.Its single-step accuracy is 0.7278 ± 0.0030 with about 37K parameters.
  • Single-step accuracy: 0.7304 ± 0.0008 is MPNO’s three-case single-step Rel-L2 over 50 test seeds, improving over Persistence at 0.8917 and Zero at 1.0000.The reported improvements are 18.1% over Persistence and 27.0% over Zero.

5. Discussion and limitations

MPNO’s stability guarantee applies to its linear propagation operator and depends on assumptions about state representation, data, resolution, and physics coupling. The methodology therefore defines an applicability boundary while offering an architectural alternative when stability must be designed rather than optimized as a loss objective.

  • Applicability assumptions: The Markov approximation is most natural for short-range-correlated systems, while high-strain-rate penetration compensates for rate-dependent constitutive memory through state-implicit strain-rate and acceleration information.The stated penetration regime is 10^3–10^5 s−1.
  • Limitations: MPNO’s spectral guarantee applies to the linear propagation operator P, not the full nonlinear model, and does not imply arbitrary accuracy.The reported accuracy is a measured post-convergence value rather than a direct consequence of the spectral constraint.
  • Limitations: Extending the spectral constraint to encoder and decoder nonlinearities, such as through spectral normalization, remains future work.This could tighten the model’s Lipschitz constant.
  • Limitations: 400 samples, four initial velocities (100/135/165/200 m/s), and a 10mm mesh limit the penetration experiments to macroscopic low-frequency impact response.Higher-resolution phenomena including material fracture and local damage patterns lie outside the considered data range.
  • Contribution boundary: MPNO’s contribution is a design path using graph-propagation spectral structure when stability must be an architectural property rather than a loss-function objective.The paper does not present MPNO as a replacement for all neural-operator paradigms.

6. Conclusions

The conclusions present MPNO as a graph-structured neural operator that reformulates transient-dynamics PDEs as Markov propagation on a spatially discretized graph. Experiments support its accuracy, efficiency, stable autoregressive behavior, and applicability under defined conservation-law conditions.

  • Core methodology: MPNO reformulates transient-dynamics PDEs as a Markov propagation process on a spatially discretized graph.Its methodology chains physical adjacency construction with graph-Laplacian spectral analysis and online estimation of λmax.
  • Experimental results: 0.7304±0.0008 single-step Rel-L2 is achieved on concrete penetration with about 20K parameters, better than WNO and comparable to FNO at about one quarter its parameters.The conclusions also report stable 29-step autoregressive rollout across the entire active horizon.
  • Applicability: MPNO applies to graph-topology systems governed by local conservation laws whose evolution can be approximated as a Markov process.The applicability conditions are delimited in §5.1.
  • Design principle: Constructively embedding physics priors in the architecture allows the data-driven components to remain lightweight.This is presented as a broader design principle derived from MPNO.

Data availability statement

All datasets are generated by independently reproducible procedures described in Section 2, with generation scripts accompanying the code.

  • Data generation: The datasets are generated reproducibly from analytic, stochastic, finite-difference, and LS-DYNA/Explicit simulation procedures described in Section 2.Burgers uses the Cole–Hopf analytic solution; Darcy uses an FFT Gaussian random field and five-point finite differences; concrete penetration uses LS-DYNA/Explicit simulations.
  • Reproducibility: Generation scripts accompany the code, supporting independent reproduction of the datasets.

Code availability and reproducibility

The PyTorch implementation specifies the data, seeding, evaluation, training, and deterministic long-horizon settings needed for reproducibility.

  • Reproducibility settings: The reproducibility protocol fixes a 70:15:15 data split with SEED=42, N = 3 baseline initializations, single-seed ablations, and a deterministic long-horizon seed rule.Evaluation uses single-step Rel-L2 with τ = 0.01 and ≥5 active nodes per frame, alongside per-scenario training hyperparameters.
  • Implementation: MPNO is implemented in PyTorch with code covering data generation, resampling, and training/evaluation entry points.The passage also identifies the torch.no_gr… component of the codebase, but the supplied text is truncated.
Loading 2608.25744v1…