Structure-preserving uncertainty quantification for GENERIC dynamics

# Structure-preserving uncertainty quantification for GENERIC dynamics
Source: [https://arxiv.org/html/2608.12624](https://arxiv.org/html/2608.12624)
Zequn HeEmail:[hezequn@engineering\.upenn\.edu](mailto:[email protected])Affiliation:Department of Mechanical Engineering and Applied Mechanics, University of Pennsylvania, Philadelphia, 19104, PA, United States of AmericaCorresponding author:Corresponding authors\.Celia ReinaEmail:[creina@engineering\.upenn\.edu](mailto:[email protected])Affiliation:Department of Mechanical Engineering and Applied Mechanics, University of Pennsylvania, Philadelphia, 19104, PA, United States of AmericaCorresponding author:Corresponding authors\.

###### Abstract

Structure\-preserving machine learning embeds physical structure directly into model architectures, yet uncertainty quantification \(UQ\) for such hard\-constrained models remains limited because standard UQ methods may violate the encoded admissibility conditions, require architectural modifications, or impose substantial computational costs\. In this work, we propose Structure\-Preserving Epistemic Neural Networks \(S\-PENNs\), a general framework for UQ in scientific machine learning models with hard architectural constraints, and instantiate it for GENERIC \(General Equation for Non\-Equilibrium Reversible\-Irreversible Coupling\) dynamics\. S\-PENNs preserve the structural constraints of a pretrained model by attaching lightweight epinets to its constrained components, ensuring that every sampled realization remains physically admissible by construction\. When applied to GENERIC dynamics, such a proposed framework yields thermodynamically consistent rollouts that preserve the first and second laws\. Furthermore, we combine S\-PENNs with split conformal prediction as a post\-hoc calibration method to produce prediction intervals with finite\-sample marginal coverage guarantees\. We validate S\-PENNs on three numerical examples: a harmonic oscillator coupled to a heat bath and an idealized chemical motor, both governed by ODEs, and a one\-dimensional viscoplastic model governed by PDEs\. Across all three examples, S\-PENNs produce thermodynamically consistent stochastic realizations and well\-calibrated prediction intervals while reducing the computational cost by about1−31\-3orders of magnitude compared to deep ensembles\. Although the present study focuses on GENERIC dynamics, S\-PENNs can be extended more broadly to scientific machine learning models in computational mechanics with either hard or soft constraints\.

###### Keywords:

Structure\-preserving machine learning , Uncertainty quantification , Thermodynamic consistency , Coverage guarantee , GENERIC , Non\-equilibrium thermodynamics , Generalized gradient flows , Constitutive modeling

## 1Introduction

Scientific machine learning has become an important computational tool in applied mechanics, where learned models are used to approximate governing equations, constitutive laws, and solution operators from data and physical knowledge\. Representative subclasses include neural ODEs for continuous\-time dynamics[7](https://arxiv.org/html/2608.12624#bib.bib73), physics\-informed neural networks \(PINNs\) for differential\-equation\-constrained learning[59](https://arxiv.org/html/2608.12624#bib.bib56),[35](https://arxiv.org/html/2608.12624#bib.bib50), and neural operators for maps between function spaces[45](https://arxiv.org/html/2608.12624#bib.bib51),[43](https://arxiv.org/html/2608.12624#bib.bib52),[36](https://arxiv.org/html/2608.12624#bib.bib57)\. A central distinction among these methods is how physical constraints are imposed\. In soft\-constrained formulations, the governing equations or admissibility conditions enter the training objective as residual or penalty terms\. In hard\-constrained formulations, often referred to as structure\-preserving machine learning, the model class is restricted through the architecture or parameterization itself, so selected invariants, stability properties, or thermodynamic laws are satisfied by construction rather than only encouraged during optimization[6](https://arxiv.org/html/2608.12624#bib.bib49),[32](https://arxiv.org/html/2608.12624#bib.bib48),[74](https://arxiv.org/html/2608.12624#bib.bib46),[78](https://arxiv.org/html/2608.12624#bib.bib30)\.

Within this structure\-preserving paradigm, reversible dynamics have served as the initial proving ground\. Hamiltonian Neural Networks[20](https://arxiv.org/html/2608.12624#bib.bib74)and Lagrangian Neural Networks[9](https://arxiv.org/html/2608.12624#bib.bib75)infer the dynamics from learned energy or action principles, whereas SympNets[32](https://arxiv.org/html/2608.12624#bib.bib48), Symplectic ODE\-Net[80](https://arxiv.org/html/2608.12624#bib.bib76), and Poisson neural networks[31](https://arxiv.org/html/2608.12624#bib.bib47)encode symplectic or Poisson geometry directly in the model class\. However, many mechanical systems are governed by coupled reversible and irreversible dynamics\. The GENERIC \(General Equation for Non\-Equilibrium Reversible\-Irreversible Coupling\) formalism[21](https://arxiv.org/html/2608.12624#bib.bib54),[54](https://arxiv.org/html/2608.12624#bib.bib55),[55](https://arxiv.org/html/2608.12624#bib.bib4),[57](https://arxiv.org/html/2608.12624#bib.bib3), which is tightly connected to metriplectic evolution[49](https://arxiv.org/html/2608.12624#bib.bib77), provides a natural foundation for such systems by decomposing the evolution into reversible and irreversible contributions, subject to structural conditions that imply energy conservation and nonnegative entropy production\. These thermodynamic guarantees make GENERIC a natural inductive bias for machine learning where neural networks provide flexible parameterizations of the thermodynamic building blocks, while the GENERIC structure controls their physical admissibility\. Early work in this direction adopted soft constraints, as in Structure\-Preserving Neural Networks, which learn conservative\-dissipative dynamics by adding the so\-called degeneracy conditions to the training objective[27](https://arxiv.org/html/2608.12624#bib.bib28)\. In contrast, hard\-constrained formulations build the structural conditions into the parameterization itself and have focused largely on the linear GENERIC/metriplectic setting, where the irreversible contribution is represented by a friction or dissipative operator applied to the entropy gradient\. In this setting, bracket\-based models parameterize admissible dissipative brackets[40](https://arxiv.org/html/2608.12624#bib.bib29)\. GFINNs reparameterize the reversible and irreversible operators so the GENERIC degeneracy conditions are satisfied for deterministic and stochastic systems[78](https://arxiv.org/html/2608.12624#bib.bib30)\. Neural Metriplectic Systems provide superior scalability of metriplectic dynamics that retains the energy\-entropy structure[23](https://arxiv.org/html/2608.12624#bib.bib32)\. Nonlinear GENERIC\-Embedded Neural Networks \(N\-GENNs\) move beyond this linear setting by learning nonlinear GENERIC dynamics with general dissipation potentials, including non\-quadratic ones[69](https://arxiv.org/html/2608.12624#bib.bib31), combining GFINN\-style projections for the reversible degeneracy condition with a convex dissipation\-potential parameterization as introduced in VONNs[29](https://arxiv.org/html/2608.12624#bib.bib40)\. These models secure thermodynamic admissibility by architectural design, but admissibility is not the same as predictive reliability\. The predictions may satisfy the thermodynamic laws while still being poorly supported by sparse data, noisy measurements, model misspecification, or extrapolation beyond the training regime\. To use these learned models in practice, one crucial requirement is to obtain uncertainty estimates that preserve the imposed structure while quantifying predictive reliability\.

Uncertainty quantification \(UQ\) provides the standard framework for characterizing such predictive uncertainty and has received substantial attention in scientific machine learning\. In broad terms, uncertainty in learned surrogates is commonly separated into aleatoric uncertainty, which captures irreducible variability from noise or unresolved stochasticity, and epistemic uncertainty, which captures lack of knowledge due to limited data, model misspecification, or unobserved regimes[58](https://arxiv.org/html/2608.12624#bib.bib67)\. Standard UQ methods include Bayesian neural networks \(BNNs\)[50](https://arxiv.org/html/2608.12624#bib.bib43),[18](https://arxiv.org/html/2608.12624#bib.bib2), which place priors over network parameters and infer a posterior predictive distribution\. Deep ensembles[39](https://arxiv.org/html/2608.12624#bib.bib17)estimate uncertainty from the spread of predictions across independently initialized and trained models\. Monte Carlo dropout[11](https://arxiv.org/html/2608.12624#bib.bib18)keeps dropout active at test time so that repeated stochastic forward passes are interpreted as samples from an approximate variational posterior\. In physics\-informed machine learning, existing UQ methods can be organized by where uncertainty is introduced\. Bayesian PINNs \(B\-PINNs\) place priors on network parameters \(or unknown physical parameters in inverse problems\) and infer posterior predictive distributions, for example, with Hamiltonian Monte Carlo or variational inference[72](https://arxiv.org/html/2608.12624#bib.bib14)\. NN\-aPC uses deep NNs to learn modal functions of arbitrary Polynomial Chaos \(aPC\) expansion of the solution or parameters in stochastic PDEs, and employs dropout to quantify their uncertainty[77](https://arxiv.org/html/2608.12624#bib.bib19)\. Adversarial UQ PINNs use latent\-variable generative models to represent PDE solution distributions under random physical inputs \(e\.g\., boundary or initial data\) or noisy measurements, while WGAN\-PINNs learn uncertainty in initial or boundary data and propagate it to interior solution fields through physics\-informed constraints[73](https://arxiv.org/html/2608.12624#bib.bib22),[13](https://arxiv.org/html/2608.12624#bib.bib23)\. Physics\-informed variational autoencoders \(PI\-VAE\) and physics\-informed variational inference introduce latent variables and physics\-based constraints to solve forward and inverse problems involving stochastic differential equations[79](https://arxiv.org/html/2608.12624#bib.bib20),[63](https://arxiv.org/html/2608.12624#bib.bib21)\. Physics\-informed polynomial chaos expansions augment conventional PCE surrogate construction with differential\-equation and boundary\-condition constraints[51](https://arxiv.org/html/2608.12624#bib.bib16)\. More recently, conformal prediction \(CP\) has been used to calibrate existing PINN uncertainty estimators into intervals with finite\-sample coverage guarantees and spatial adaptivity[75](https://arxiv.org/html/2608.12624#bib.bib15)\. These approaches target soft\-penalized or physics\-regularized surrogates, rather than models with hard\-encoded physical constraints in which every sampled realization must satisfy structural admissibility conditions\.

By contrast, UQ for models with hard\-encoded physical constraints must preserve admissibility rather than merely estimate variation around a physics\-regularized predictor\. Existing methods have demonstrated such a requirement in several specialized settings\. For Hamiltonian systems with additive dissipation, Symplectic Spectrum Gaussian Processes \(SSGP\) use a symplectic Gaussian\-process prior with random Fourier features to infer Hamiltonians from noisy, sparse data while retaining conservative or dissipative structure[68](https://arxiv.org/html/2608.12624#bib.bib33)\. Similarly, Hamiltonian Gaussian Processes \(HGP\) use a decoupled Gaussian\-process Hamiltonian and energy\-conserving shooting for inference from long noisy trajectories[60](https://arxiv.org/html/2608.12624#bib.bib34), while Bayesian identification of nonseparable Hamiltonians combines a structure\-preserving Hamiltonian parameterization, dependent additive and multiplicative noise models, and reduced\-order Bayesian inference[12](https://arxiv.org/html/2608.12624#bib.bib35)\. Thermodynamic and constitutive examples use analogous admissible representations\. Thermodynamically constrained Gaussian\-process equations of state capture model and data uncertainty subject to consistency and stability constraints[62](https://arxiv.org/html/2608.12624#bib.bib38)\. Bayesian\-EUCLID performs Bayesian sparse discovery over a hyperelastic feature library from displacement and reaction\-force data, quantifying uncertainty in active constitutive terms and data noise through a momentum\-balance likelihood and spike\-slab prior[33](https://arxiv.org/html/2608.12624#bib.bib36)\. Bayesian constitutive artificial neural networks learn probability distributions for the weights of free\-energy\-based constitutive networks, yielding credible intervals for the used model terms and associated constitutive responses[44](https://arxiv.org/html/2608.12624#bib.bib37)\. More recently, conformal quantile regression adopts a frequentist approach to probabilistic modeling and has calibrated tensor\-valued quantile predictions from a strain\-invariant, polyconvex, and thermodynamically consistent constitutive network, with the resulting intervals targeting aleatoric uncertainty from data variability rather than epistemic uncertainty[3](https://arxiv.org/html/2608.12624#bib.bib39)\. This body of work shows that UQ for hard\-encoded physics constraints is both desirable and feasible, but it also indicates that generic UQ methods do not automatically yield hard\-constrained UQ easily\.

The core limitation is that hard\-encoded constraints restrict the admissible model class itself, so sampled UQ realizations must remain in a space that is closed under those constraints\. Randomizing unconstrained weights or perturbing intermediate activations can move a realization outside that space, while preserving admissibility by independent replication requires training complete constrained models\. BNNs can preserve structure only when the prior, posterior approximation, and inference procedure are defined over the constrained surrogate itself; otherwise, generic weight\-space uncertainty does not enforce the required relations, and posterior inference over a constrained architecture can be expensive at scale[34](https://arxiv.org/html/2608.12624#bib.bib80)\. Deep ensembles can preserve constraints by training several full constrained models with cost growing essentially linearly in the ensemble size[39](https://arxiv.org/html/2608.12624#bib.bib17),[14](https://arxiv.org/html/2608.12624#bib.bib81)\. MC dropout uses stochastic internal masks, so hard\-encoded properties may be lost unless every masked subnetwork still satisfies the same constraints[11](https://arxiv.org/html/2608.12624#bib.bib18),[14](https://arxiv.org/html/2608.12624#bib.bib81)\. This difficulty increases in structure\-preserving models whose admissibility usually depends on several interacting constrained building blocks rather than on a single output map\. Therefore, a practical UQ framework is expected to introduce uncertainty inside the constrained parameterization, preserve the structural relations under each realization, propagate uncertainty coherently through the coupled constrained blocks, and keep the added computational cost tractable\.

To address the gap and meet the requirements, we propose Structure\-Preserving Epistemic Neural Networks \(S\-PENNs\), an epistemic UQ framework that injects uncertainty into constrained scientific machine learning parameterizations without violating their encoded physical and structural relations\. S\-PENNs build on the Epistemic Neural Networks \(ENNs\) introduced in[53](https://arxiv.org/html/2608.12624#bib.bib44), in which a conventional base network is augmented by a lightweight auxiliary network called the epinet\. Recent work has applied this idea to operator learning through NEON[24](https://arxiv.org/html/2608.12624#bib.bib24), to PINNs through E\-PINNs[30](https://arxiv.org/html/2608.12624#bib.bib27), and to thermodynamics\-informed diffusion models through EVODMs[26](https://arxiv.org/html/2608.12624#bib.bib25)and SPIEDiff[25](https://arxiv.org/html/2608.12624#bib.bib26)\. Nevertheless, none of these works address structure\-preserving architectures formally in settings in which each UQ realization must satisfy coupled admissibility conditions by construction\. To demonstrate the feasibility of S\-PENNs, we develop and instantiate the framework for GENERIC dynamics, specifically, on top of the N\-GENNs[69](https://arxiv.org/html/2608.12624#bib.bib31)\. This is a particularly demanding test bed for structure\-preserving UQ because the GENERIC formalism couples four thermodynamic building blocks through skew\-symmetry, convexity, and degeneracy conditions that jointly enforce the first and second laws of thermodynamics, which will be detailed in Section[2\.1](https://arxiv.org/html/2608.12624#S2.SS1)\. Any UQ method equipped on N\-GENNs must preserve these constraints jointly rather than one block at a time, and the perturbations must propagate coherently across blocks to yield consistent uncertainty estimates for the resulting dynamics\. The key idea here for S\-PENNs is to treat the deterministic N\-GENNs as base networks and attach a separate epinet to each building block of N\-GENNs so that each epinet inherits the architectural properties of its corresponding block\. Since the perturbed blocks are then assembled through the same structure\-preserving reparameterizations as in N\-GENNs, every sampled realization remains thermodynamically admissible by construction\. The component\-wise perturbations are driven by a common source of randomness\. Each epinet is conditioned on stop\-gradient features from its corresponding deterministic block and on the same inputs as that block, while the shared epistemic index couples the block\-wise perturbations across the assembled GENERIC vector field\. One uncertainty realization therefore perturbs the entire GENERIC parameterization coherently and propagates uncertainty consistently across blocks without breaking the thermodynamic structure\. The nonlinear GENERIC instantiation of S\-PENNs also requires uncertainty representations for global, input\-independent quantities, such as the state\-independent terms in the GENERIC reparameterization\. This setting is analogous to inverse problems with unknown physical parameters, where the parameters are inferred jointly with the solution field but are not themselves functions of the input variables\. Prior epinet\-based work in scientific machine learning has either not considered such global parameters[24](https://arxiv.org/html/2608.12624#bib.bib24),[26](https://arxiv.org/html/2608.12624#bib.bib25),[25](https://arxiv.org/html/2608.12624#bib.bib26)or handled them indirectly by softly penalizing parameter variance during optimization[30](https://arxiv.org/html/2608.12624#bib.bib27)\. To have a hard\-constrained solution for this problem, we introduce an epinet construction tailored to global parameters, which enforces input\-independent uncertainty by construction with no soft\-penalized variance term\. Lastly, beyond preserving thermodynamic consistency, we use split conformal prediction[2](https://arxiv.org/html/2608.12624#bib.bib58)on S\-PENNs as a post\-hoc calibration step to obtain distribution\-free finite\-sample marginal coverage guarantees on the resulting prediction intervals\. The proposed framework is evaluated on three numerical examples: a harmonic oscillator coupled to a heat bath, an idealized chemical motor, and a one\-dimensional viscoplastic model\. Across all three examples, S\-PENNs yield thermodynamically admissible uncertainty realizations and well\-calibrated prediction intervals while reducing training cost by about1−31\-3orders of magnitude compared to deep ensembles\. Although the presented work focuses on GENERIC dynamics, the proposed idea extends to broader classes of scientific machine learning models that encode physical structure through either hard architectural constraints or soft physics\-informed penalties\.

The remainder of the paper is organized as follows\. Section[2](https://arxiv.org/html/2608.12624#S2)reviews the GENERIC formalism and the N\-GENNs framework of[69](https://arxiv.org/html/2608.12624#bib.bib31), introduces the S\-PENNs construction, and presents the conformal calibration procedure\. Section[3](https://arxiv.org/html/2608.12624#S3)reports the main numerical results for a harmonic oscillator coupled to a heat bath \(Section[3\.1](https://arxiv.org/html/2608.12624#S3.SS1)\), an idealized chemical motor \(Section[3\.2](https://arxiv.org/html/2608.12624#S3.SS2)\), and a one\-dimensional viscoplastic model \(Section[3\.3](https://arxiv.org/html/2608.12624#S3.SS3)\)\. Section[4](https://arxiv.org/html/2608.12624#S4)discusses limitations and future directions\. Additional inverse\-problem benchmarks demonstrating the proposed epinet construction are provided in[D](https://arxiv.org/html/2608.12624#A4)\.

## 2Methodology

This section begins by reviewing the structure\-preserving N\-GENNs framework for nonlinear GENERIC dynamics in Section[2\.1](https://arxiv.org/html/2608.12624#S2.SS1)\. Section[2\.2](https://arxiv.org/html/2608.12624#S2.SS2)reviews the epinet construction, and Section[2\.3](https://arxiv.org/html/2608.12624#S2.SS3)specializes it to N\-GENNs to obtain the proposed S\-PENNs architecture\. Finally, Section[2\.4](https://arxiv.org/html/2608.12624#S2.SS4)presents the split conformal calibration procedure used to obtain finite\-sample marginal coverage guarantees for the resulting prediction intervals\.

### 2\.1Nonlinear GENERIC\-embedded neural networks \(N\-GENNs\)

We first recall the nonlinear GENERIC formalism, which decomposes the evolution into a reversible Hamiltonian part and an irreversible generalized gradient flow\. For a finite\-dimensional state vector𝐱∈ℝd\\mathbf\{x\}\\in\\mathbb\{R\}^\{d\}and its conjugate variable𝐱∗\\mathbf\{x\}^\{\*\}, the nonlinear GENERIC dynamics are written as

𝐱˙=L⁡\(𝐱\)​D​E​\(𝐱\)\+D𝐱∗​Ξ​\(𝐱,𝐱∗\)\|𝐱∗=D​S​\(𝐱\)\.\\dot\{\\mathbf\{x\}\}=L\(\\mathbf\{x\}\)\\mathrm\{D\}E\(\\mathbf\{x\}\)\+\\left\.\\mathrm\{D\}\_\{\\mathbf\{x\}^\{\*\}\}\\Xi\(\\mathbf\{x\},\\mathbf\{x\}^\{\*\}\)\\right\|\_\{\\mathbf\{x\}^\{\*\}=\\mathrm\{D\}S\(\\mathbf\{x\}\)\}\.\(1\)HereLLis the Poisson operator describing reversible dynamics,EEis the total energy,SSis the entropy, andΞ\\Xiis the dissipation potential\. In finite\-dimensional settings,D\\mathrm\{D\}andD𝐱∗\\mathrm\{D\}\_\{\\mathbf\{x\}^\{\*\}\}denote ordinary partial derivatives; in infinite\-dimensional settings, where the state variables are fields, they denote the corresponding functional derivatives\. Among the thermodynamic constraints imposed on the GENERIC building blocks, the two degeneracy conditions are

reversible degeneracy condition\\displaystyle\\text\{reversible degeneracy condition\}L⁡\(𝐱\)​D​S​\(𝐱\)=𝟎,∀𝐱,\\displaystyle L\(\\mathbf\{x\}\)\\mathrm\{D\}S\(\\mathbf\{x\}\)=\\mathbf\{0\},\\quad\\forall\\mathbf\{x\},\(2\)irreversible degeneracy condition\\displaystyle\\text\{irreversible degeneracy condition\}Ξ⁡\(𝐱,𝐱∗\+λ​D​E​\(𝐱\)\)=Ξ⁡\(𝐱,𝐱∗\),∀𝐱,𝐱∗,∀λ∈ℝ\.\\displaystyle\\Xi\(\\mathbf\{x\},\\mathbf\{x\}^\{\*\}\+\\lambda\\mathrm\{D\}E\(\\mathbf\{x\}\)\)=\\Xi\(\\mathbf\{x\},\\mathbf\{x\}^\{\*\}\),\\quad\\forall\\mathbf\{x\},\\mathbf\{x\}^\{\*\},\\quad\\forall\\lambda\\in\\mathbb\{R\}\.When combined with the skew\-symmetry of the Poisson operator,L⁡\(𝐱\)=−L​\(𝐱\)⊤L\(\\mathbf\{x\}\)=\-L\(\\mathbf\{x\}\)^\{\\top\}, and the admissibility conditionsΞ⁡\(𝐱,𝟎\)=0\\Xi\(\\mathbf\{x\},\\mathbf\{0\}\)=0,Ξ⁡\(𝐱,𝐱∗\)≥0\\Xi\(\\mathbf\{x\},\\mathbf\{x\}^\{\*\}\)\\geq 0, and convexity ofΞ\\Xiwith respect to𝐱∗\\mathbf\{x\}^\{\*\}, these degeneracy conditions ensure conservation of total energy and nonnegative entropy production[22](https://arxiv.org/html/2608.12624#bib.bib41),[37](https://arxiv.org/html/2608.12624#bib.bib42)\.

Building on the above GENERIC structure, N\-GENNs[69](https://arxiv.org/html/2608.12624#bib.bib31)represent the thermodynamic building blocks with neural networks and enforce the constraints through structure\-preserving reparameterizations\. With the deterministic trainable parameter set denoted by𝝍\\bm\{\\psi\}111In this subsection,𝝍\\bm\{\\psi\}denotes the aggregate trainable parameter set of the deterministic N\-GENNs backbone\. It includes neural\-network weights and biases as well as finite\-dimensional trainable variables used in the structure\-preserving reparameterizations\. A shared subscript𝝍\\bm\{\\psi\}is used here only for clarity and does not imply weight sharing among distinct networks\., the deterministic backbone uses scalar networksE𝝍E\_\{\\bm\{\\psi\}\}andS𝝍S\_\{\\bm\{\\psi\}\}, a reversible operator networkL~𝝍\\tilde\{L\}\_\{\\bm\{\\psi\}\}, and a raw dissipation\-potential networkΞ~𝝍\\tilde\{\\Xi\}\_\{\\bm\{\\psi\}\}\. Let𝑩𝝍=\(𝑩1​𝝍,…,𝑩d​𝝍\)∈ℝd×d×d\\bm\{B\}\_\{\\bm\{\\psi\}\}=\(\\bm\{B\}\_\{1\\bm\{\\psi\}\},\\dots,\\bm\{B\}\_\{d\\bm\{\\psi\}\}\)\\in\\mathbb\{R\}^\{d\\times d\\times d\}denote the trainable matrices included in𝝍\\bm\{\\psi\}, and define

𝑨i​𝝍=𝑩i​𝝍−𝑩i​𝝍⊤,i=1,…,d\.\\bm\{A\}\_\{i\\bm\{\\psi\}\}=\\bm\{B\}\_\{i\\bm\{\\psi\}\}\-\\bm\{B\}\_\{i\\bm\{\\psi\}\}^\{\\top\},\\qquad i=1,\\dots,d\.\(3\)The matrixQS𝝍​\(𝐱\)∈ℝd×dQ\_\{S\_\{\\bm\{\\psi\}\}\}\(\\mathbf\{x\}\)\\in\\mathbb\{R\}^\{d\\times d\}is constructed so that itsiith row is\(𝑨i​𝝍​D​S𝝍​\(𝐱\)\)⊤\(\\bm\{A\}\_\{i\\bm\{\\psi\}\}\\mathrm\{D\}S\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)\)^\{\\top\}, and the reparameterized Poisson operator is defined as

L𝝍​\(𝐱\)=QS𝝍​\(𝐱\)⊤​L~𝝍​\(𝐱\)​QS𝝍​\(𝐱\)\.L\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)=Q\_\{S\_\{\\bm\{\\psi\}\}\}\(\\mathbf\{x\}\)^\{\\top\}\\tilde\{L\}\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)Q\_\{S\_\{\\bm\{\\psi\}\}\}\(\\mathbf\{x\}\)\.\(4\)Because each𝑨i​𝝍\\bm\{A\}\_\{i\\bm\{\\psi\}\}is skew\-symmetric,QS𝝍​\(𝐱\)​D​S𝝍​\(𝐱\)=0Q\_\{S\_\{\\bm\{\\psi\}\}\}\(\\mathbf\{x\}\)\\mathrm\{D\}S\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)=0, which givesL𝝍​\(𝐱\)​D​S𝝍​\(𝐱\)=0L\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)\\mathrm\{D\}S\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)=0\. This enforces the reversible degeneracy condition by construction\.

For the irreversible part, the raw dissipation potentialΞ~𝝍​\(𝐱,𝐱∗\)\\tilde\{\\Xi\}\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\},\\mathbf\{x\}^\{\*\}\)is represented by a partially input\-convex neural network \(PICNN\)[1](https://arxiv.org/html/2608.12624#bib.bib82),[29](https://arxiv.org/html/2608.12624#bib.bib40), treating𝐱∗\\mathbf\{x\}^\{\*\}as the convex variable and𝐱\\mathbf\{x\}as a conditioning input\. LetPE𝝍P\_\{E\_\{\\bm\{\\psi\}\}\}be the projection matrix onto the orthogonal complement ofD​E𝝍\\mathrm\{D\}E\_\{\\bm\{\\psi\}\},

PE𝝍​\(𝐱\)=I−D​E𝝍​\(𝐱\)​D​E𝝍​\(𝐱\)⊤‖D​E𝝍​\(𝐱\)‖22\.P\_\{E\_\{\\bm\{\\psi\}\}\}\(\\mathbf\{x\}\)=I\-\\frac\{\\mathrm\{D\}E\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)\\mathrm\{D\}E\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)^\{\\top\}\}\{\\\|\\mathrm\{D\}E\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)\\\|\_\{2\}^\{2\}\}\.\(5\)The dissipation potential used in the dynamics is obtained by the reparameterization as follows

Ξ𝝍​\(𝐱,𝐱∗\)\\displaystyle\\Xi\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\},\\mathbf\{x\}^\{\*\}\)=Ξ~𝝍​\(𝐱,PE𝝍​\(𝐱\)​𝐱∗\)−Ξ~𝝍​\(𝐱,𝟎\)−\(D𝐱∗​Ξ~𝝍​\(𝐱,𝐱∗\)\|𝐱∗=𝟎\)⊤​PE𝝍​\(𝐱\)​𝐱∗\.\\displaystyle=\\tilde\{\\Xi\}\_\{\\bm\{\\psi\}\}\\\!\\left\(\\mathbf\{x\},P\_\{E\_\{\\bm\{\\psi\}\}\}\(\\mathbf\{x\}\)\\mathbf\{x\}^\{\*\}\\right\)\-\\tilde\{\\Xi\}\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\},\\mathbf\{0\}\)\-\\left\(\\left\.\\mathrm\{D\}\_\{\\mathbf\{x\}^\{\*\}\}\\tilde\{\\Xi\}\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\},\\mathbf\{x\}^\{\*\}\)\\right\|\_\{\\mathbf\{x\}^\{\*\}=\\mathbf\{0\}\}\\right\)^\{\\top\}P\_\{E\_\{\\bm\{\\psi\}\}\}\(\\mathbf\{x\}\)\\mathbf\{x\}^\{\*\}\.\(6\)Here,𝐱∗\\mathbf\{x\}^\{\*\}enters only throughPE𝝍​\(𝐱\)​𝐱∗P\_\{E\_\{\\bm\{\\psi\}\}\}\(\\mathbf\{x\}\)\\mathbf\{x\}^\{\*\}\. SincePE𝝍​\(𝐱\)​D​E𝝍​\(𝐱\)=0P\_\{E\_\{\\bm\{\\psi\}\}\}\(\\mathbf\{x\}\)\\mathrm\{D\}E\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)=0, replacing𝐱∗\\mathbf\{x\}^\{\*\}by𝐱∗\+λ​D​E𝝍​\(𝐱\)\\mathbf\{x\}^\{\*\}\+\\lambda\\mathrm\{D\}E\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)leaves this projected argument unchanged, and thereforeΞ𝝍​\(𝐱,𝐱∗\+λ​D​E𝝍​\(𝐱\)\)=Ξ𝝍​\(𝐱,𝐱∗\)\\Xi\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\},\\mathbf\{x\}^\{\*\}\+\\lambda\\mathrm\{D\}E\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)\)=\\Xi\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\},\\mathbf\{x\}^\{\*\}\)\. As a result, the irreversible degeneracy condition is satisfied\. The reparameterization preserves convexity in𝐱∗\\mathbf\{x\}^\{\*\}because it combines a linear composition of the raw convex potential with an affine subtraction\. The affine correction further givesΞ𝝍​\(𝐱,𝟎\)=0\\Xi\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\},\\mathbf\{0\}\)=0andD𝐱∗​Ξ𝝍​\(𝐱,𝟎\)=0\\mathrm\{D\}\_\{\\mathbf\{x\}^\{\*\}\}\\Xi\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\},\\mathbf\{0\}\)=0, so by convexity𝐱∗=𝟎\\mathbf\{x\}^\{\*\}=\\mathbf\{0\}is a global minimizer andΞ𝝍​\(𝐱,𝐱∗\)≥0\\Xi\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\},\\mathbf\{x\}^\{\*\}\)\\geq 0\. These properties makeΞ𝝍\\Xi\_\{\\bm\{\\psi\}\}admissible by construction\.

The constrained building blocks define the deterministic N\-GENNs vector field

g𝝍​\(𝐱\)=L𝝍​\(𝐱\)​D​E𝝍​\(𝐱\)\+D𝐱∗​Ξ𝝍​\(𝐱,𝐱∗\)\|𝐱∗=D​S𝝍​\(𝐱\)\.g\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)=L\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)\\mathrm\{D\}E\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)\+\\left\.\\mathrm\{D\}\_\{\\mathbf\{x\}^\{\*\}\}\\Xi\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\},\\mathbf\{x\}^\{\*\}\)\\right\|\_\{\\mathbf\{x\}^\{\*\}=\\mathrm\{D\}S\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)\}\.\(7\)Since the reversible and irreversible constraints are built into the parameterization,g𝝍g\_\{\\bm\{\\psi\}\}is thermodynamically admissible by construction\.

### 2\.2Epistemic neural networks and the epinet

According to[53](https://arxiv.org/html/2608.12624#bib.bib44), ENNs refer to conventional neural networks augmented with a lightweight auxiliary network, the epinet\. Such models can be written as

τϑ​\(𝐱,𝐳\)=μ𝝍​\(𝐱\)\+σϕ​\(𝒉¯𝝍​\(𝐱\),𝐳\),𝒉¯𝝍​\(𝐱\)=\[sg⁡\(𝒉𝝍​\(𝐱\)\),𝐱\],\\tau\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{z\}\)=\\mu\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)\+\\sigma\_\{\\bm\{\\phi\}\}\(\\bar\{\\bm\{h\}\}\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\),\\mathbf\{z\}\),\\qquad\\bar\{\\bm\{h\}\}\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)=\[\\mathrm\{sg\}\\\!\\left\(\\bm\{h\}\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)\\right\),\\mathbf\{x\}\],\(8\)where𝝍\\bm\{\\psi\}andϕ\\bm\{\\phi\}denote the base\-network and epinet trainable parameter sets, respectively, andϑ=\(𝝍,ϕ\)\\bm\{\\vartheta\}=\(\\bm\{\\psi\},\\bm\{\\phi\}\)denotes the trainable parameter set of the augmented model\. Hereμ𝝍​\(𝐱\)\\mu\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)denotes the base network andσϕ​\(𝒉¯𝝍​\(𝐱\),𝐳\)\\sigma\_\{\\bm\{\\phi\}\}\(\\bar\{\\bm\{h\}\}\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\),\\mathbf\{z\}\)is the epinet\.𝐳∈ℝdz\\mathbf\{z\}\\in\\mathbb\{R\}^\{d\_\{z\}\}is an epistemic index drawn from a chosen reference distributionπ⁡\(𝐳\)\\pi\(\\mathbf\{z\}\)\.𝒉𝝍​\(𝐱\)\\bm\{h\}\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)represents the base\-network features, typically taken from the last hidden layer, andsg⁡\(⋅\)\\mathrm\{sg\}\(\\cdot\)is the stop\-gradient operator that prevents gradients from flowing into the base network\. Throughout the paper, the notation𝒉¯𝝍​\(𝐱\)=\[sg⁡\(𝒉𝝍​\(𝐱\)\),𝐱\]\\bar\{\\bm\{h\}\}\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)=\[\\mathrm\{sg\}\(\\bm\{h\}\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)\),\\mathbf\{x\}\]means the concatenation of stop\-gradient hidden features and the raw input variables\. This is the same typical setting used in[53](https://arxiv.org/html/2608.12624#bib.bib44), but we state it explicitly here to avoid ambiguity\.

The epinet is decomposed into a learnable networkσϕlearn\\sigma\_\{\\bm\{\\phi\}\}^\{\\mathrm\{learn\}\}and a prior networkσprior\\sigma^\{\\mathrm\{prior\}\},

σϕ​\(𝒉¯,𝐳\)=σϕlearn​\(𝒉¯,𝐳\)\+w​σprior​\(𝒉¯,𝐳\),\\sigma\_\{\\bm\{\\phi\}\}\(\\bar\{\\bm\{h\}\},\\mathbf\{z\}\)=\\sigma\_\{\\bm\{\\phi\}\}^\{\\mathrm\{learn\}\}\(\\bar\{\\bm\{h\}\},\\mathbf\{z\}\)\+w\\,\\sigma^\{\\mathrm\{prior\}\}\(\\bar\{\\bm\{h\}\},\\mathbf\{z\}\),\(9\)wherew\>0w\>0is the prior\-weight hyperparameter that scales the contribution of the fixed random prior network relative to the learnable network, thereby controlling the strength of the prior perturbation and hence the initial epistemic spread\. Since the prior network is randomly initialized and kept fixed throughout training, its parameters are omitted from the notation for conciseness\. The learnable and prior networks are written as

σϕlearn​\(𝒉¯,𝐳\)\\displaystyle\\sigma\_\{\\bm\{\\phi\}\}^\{\\mathrm\{learn\}\}\(\\bar\{\\bm\{h\}\},\\mathbf\{z\}\)=NNϕ​\(\[𝒉¯,𝐳\]\)⊤​𝐳,\\displaystyle=\\mathrm\{NN\}\_\{\\bm\{\\phi\}\}\(\[\\bar\{\\bm\{h\}\},\\mathbf\{z\}\]\)^\{\\top\}\\mathbf\{z\},\(10\)σprior​\(𝒉¯,𝐳\)\\displaystyle\\sigma^\{\\mathrm\{prior\}\}\(\\bar\{\\bm\{h\}\},\\mathbf\{z\}\)=∑n=1dzzn​σnprior​\(𝒉¯\),\\displaystyle=\\sum\_\{n=1\}^\{d\_\{z\}\}z\_\{n\}\\,\\sigma\_\{n\}^\{\\mathrm\{prior\}\}\(\\bar\{\\bm\{h\}\}\),\(11\)respectively\. HereNNϕ\\mathrm\{NN\}\_\{\\bm\{\\phi\}\}is a trainable network whose output dimension matches that of𝐳\\mathbf\{z\}, and\{σnprior\}n=1dz\\\{\\sigma\_\{n\}^\{\\mathrm\{prior\}\}\\\}\_\{n=1\}^\{d\_\{z\}\}are fixed randomly initialized networks that typically share the same architecture as the learnable network, with one prior network for each component of the epistemic index\.

### 2\.3Structure\-preserving epistemic neural networks \(S\-PENNs\)

We now introduce S\-PENNs and specialize the framework to N\-GENNs for learning GENERIC dynamics with quantified uncertainty\. S\-PENNs attach epinets to the thermodynamic building blocks of the N\-GENNs backbone, and reassemble the perturbed objects through the same structure\-preserving parameterization\. Consequently, every draw of the epistemic index yields a thermodynamically consistent dynamics sample\. Fig\.[1](https://arxiv.org/html/2608.12624#S2.F1)depicts the overall framework\.

Figure 1:Schematic of the S\-PENNs construction for GENERIC dynamics\. The state variables𝐱\\mathbf\{x\}are passed through the deterministic N\-GENNs base networks, which produce the raw thermodynamic building blocksL~𝝍\\tilde\{L\}\_\{\\bm\{\\psi\}\},E𝝍E\_\{\\bm\{\\psi\}\},S𝝍S\_\{\\bm\{\\psi\}\}, andΞ~𝝍\\tilde\{\\Xi\}\_\{\\bm\{\\psi\}\}, together with their hidden representations\. The stop\-gradient hidden representations and the relevant inputs are used to construct epinet feature vectors𝒉¯𝝍L~\\bar\{\\bm\{h\}\}\_\{\\bm\{\\psi\}\}^\{\\tilde\{L\}\},𝒉¯𝝍E\\bar\{\\bm\{h\}\}\_\{\\bm\{\\psi\}\}^\{E\},𝒉¯𝝍S\\bar\{\\bm\{h\}\}\_\{\\bm\{\\psi\}\}^\{S\}, and𝒉¯𝝍Ξ~\\bar\{\\bm\{h\}\}\_\{\\bm\{\\psi\}\}^\{\\tilde\{\\Xi\}\}\. These feature vectors, together with a shared epistemic index𝐳∼π⁡\(𝐳\)\\mathbf\{z\}\\sim\\pi\(\\mathbf\{z\}\), are passed to block\-specific epinets, which generate perturbations of the reversible operator, energy, entropy, and dissipation\-potential branches\. The epinet outputs are then assembled with the corresponding base\-network outputs to obtainL~ϑ\\tilde\{L\}\_\{\\bm\{\\vartheta\}\},EϑE\_\{\\bm\{\\vartheta\}\},SϑS\_\{\\bm\{\\vartheta\}\}, andΞ~ϑ\\tilde\{\\Xi\}\_\{\\bm\{\\vartheta\}\}; the same structure\-preserving reparameterizations used in N\-GENNs map the raw reversible operator and dissipation\-potential blocks toLϑL\_\{\\bm\{\\vartheta\}\}andΞϑ\\Xi\_\{\\bm\{\\vartheta\}\}\. Hence each draw of𝐳\\mathbf\{z\}defines one admissible GENERIC vector fieldLϑ​\(𝐱,𝐳\)​D​Eϑ​\(𝐱,𝐳\)\+D𝐱∗​Ξϑ​\(𝐱,𝐱∗,𝐳\)\|𝐱∗=D​Sϑ​\(𝐱,𝐳\)L\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{z\}\)\\mathrm\{D\}E\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{z\}\)\+\\left\.\\mathrm\{D\}\_\{\\mathbf\{x\}^\{\*\}\}\\Xi\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{x\}^\{\*\},\\mathbf\{z\}\)\\right\|\_\{\\mathbf\{x\}^\{\*\}=\\mathrm\{D\}S\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{z\}\)\}, and repeated draws form a structure\-preserving dynamics ensemble\.To preserve the thermodynamic structure under epistemic perturbations, we attach a separate epinet to each deterministic building block in Eq\. \([7](https://arxiv.org/html/2608.12624#S2.E7)\), i\.e\., the reversible operator networkL~𝝍\\tilde\{L\}\_\{\\bm\{\\psi\}\}, the energy networkE𝝍E\_\{\\bm\{\\psi\}\}, the entropy networkS𝝍S\_\{\\bm\{\\psi\}\}, and the raw dissipation\-potential networkΞ~𝝍\\tilde\{\\Xi\}\_\{\\bm\{\\psi\}\}\. The augmented blocks can be written as

L~ϑ​\(𝐱,𝐳\)\\displaystyle\\tilde\{L\}\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{z\}\)=L~𝝍​\(𝐱\)\+σϕL~​\(𝒉¯𝝍L~​\(𝐱\),𝐳\),\\displaystyle=\\tilde\{L\}\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)\+\\sigma\_\{\\bm\{\\phi\}\}^\{\\tilde\{L\}\}\(\\bar\{\\bm\{h\}\}\_\{\\bm\{\\psi\}\}^\{\\tilde\{L\}\}\(\\mathbf\{x\}\),\\mathbf\{z\}\),\(12\)Eϑ​\(𝐱,𝐳\)\\displaystyle E\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{z\}\)=E𝝍​\(𝐱\)\+σϕE​\(𝒉¯𝝍E​\(𝐱\),𝐳\),\\displaystyle=E\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)\+\\sigma\_\{\\bm\{\\phi\}\}^\{E\}\(\\bar\{\\bm\{h\}\}\_\{\\bm\{\\psi\}\}^\{E\}\(\\mathbf\{x\}\),\\mathbf\{z\}\),Sϑ​\(𝐱,𝐳\)\\displaystyle S\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{z\}\)=S𝝍​\(𝐱\)\+σϕS​\(𝒉¯𝝍S​\(𝐱\),𝐳\),\\displaystyle=S\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)\+\\sigma\_\{\\bm\{\\phi\}\}^\{S\}\(\\bar\{\\bm\{h\}\}\_\{\\bm\{\\psi\}\}^\{S\}\(\\mathbf\{x\}\),\\mathbf\{z\}\),Ξ~ϑ​\(𝐱,𝐱∗,𝐳\)\\displaystyle\\tilde\{\\Xi\}\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{x\}^\{\*\},\\mathbf\{z\}\)=Ξ~𝝍​\(𝐱,𝐱∗\)\+σϕΞ~​\(𝒉¯𝝍Ξ~​\(𝐱,𝐱∗\),𝐳\)\.\\displaystyle=\\tilde\{\\Xi\}\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\},\\mathbf\{x\}^\{\*\}\)\+\\sigma\_\{\\bm\{\\phi\}\}^\{\\tilde\{\\Xi\}\}\(\\bar\{\\bm\{h\}\}\_\{\\bm\{\\psi\}\}^\{\\tilde\{\\Xi\}\}\(\\mathbf\{x\},\\mathbf\{x\}^\{\*\}\),\\mathbf\{z\}\)\.Here and below, the shared subscripts𝝍\\bm\{\\psi\}andϕ\\bm\{\\phi\}denote aggregate parameter collections and do not imply weight sharing across thermodynamic blocks\. The inputs to the correspoding epinet branches are denoted as𝒉¯𝝍L~\\bar\{\\bm\{h\}\}\_\{\\bm\{\\psi\}\}^\{\\tilde\{L\}\},𝒉¯𝝍E\\bar\{\\bm\{h\}\}\_\{\\bm\{\\psi\}\}^\{E\},𝒉¯𝝍S\\bar\{\\bm\{h\}\}\_\{\\bm\{\\psi\}\}^\{S\}, and𝒉¯𝝍Ξ~\\bar\{\\bm\{h\}\}\_\{\\bm\{\\psi\}\}^\{\\tilde\{\\Xi\}\}, more specifically, they are constructed as

𝒉¯𝝍L~​\(𝐱\)\\displaystyle\\bar\{\\bm\{h\}\}\_\{\\bm\{\\psi\}\}^\{\\tilde\{L\}\}\(\\mathbf\{x\}\)=\[sg⁡\(𝒉𝝍L~​\(𝐱\)\),𝐱\],\\displaystyle=\[\\mathrm\{sg\}\\\!\\left\(\\bm\{h\}\_\{\\bm\{\\psi\}\}^\{\\tilde\{L\}\}\(\\mathbf\{x\}\)\\right\),\\mathbf\{x\}\],\(13\)𝒉¯𝝍E​\(𝐱\)\\displaystyle\\bar\{\\bm\{h\}\}\_\{\\bm\{\\psi\}\}^\{E\}\(\\mathbf\{x\}\)=\[sg⁡\(𝒉𝝍E​\(𝐱\)\),𝐱\],\\displaystyle=\[\\mathrm\{sg\}\\\!\\left\(\\bm\{h\}\_\{\\bm\{\\psi\}\}^\{E\}\(\\mathbf\{x\}\)\\right\),\\mathbf\{x\}\],𝒉¯𝝍S​\(𝐱\)\\displaystyle\\bar\{\\bm\{h\}\}\_\{\\bm\{\\psi\}\}^\{S\}\(\\mathbf\{x\}\)=\[sg⁡\(𝒉𝝍S​\(𝐱\)\),𝐱\],\\displaystyle=\[\\mathrm\{sg\}\\\!\\left\(\\bm\{h\}\_\{\\bm\{\\psi\}\}^\{S\}\(\\mathbf\{x\}\)\\right\),\\mathbf\{x\}\],𝒉¯𝝍Ξ~​\(𝐱,𝐱∗\)\\displaystyle\\bar\{\\bm\{h\}\}\_\{\\bm\{\\psi\}\}^\{\\tilde\{\\Xi\}\}\(\\mathbf\{x\},\\mathbf\{x\}^\{\*\}\)=\[sg⁡\(𝒉𝝍Ξ~​\(𝐱,𝐱∗\)\),𝐱,𝐱∗\]\.\\displaystyle=\[\\mathrm\{sg\}\\\!\\left\(\\bm\{h\}\_\{\\bm\{\\psi\}\}^\{\\tilde\{\\Xi\}\}\(\\mathbf\{x\},\\mathbf\{x\}^\{\*\}\)\\right\),\\mathbf\{x\},\\mathbf\{x\}^\{\*\}\]\.The same epistemic index𝐳\\mathbf\{z\}is used in all augmented blocks, so each draw selects a joint perturbation of the full GENERIC parameterization\.

#### 2\.3\.1Input\-independent epinet branch for global tensors

As introduced in Section[2\.1](https://arxiv.org/html/2608.12624#S2.SS1), the N\-GENNs backbone contains finite\-dimensional trainable quantities that are not state\-dependent neural\-network outputs\. The tensor𝑩𝝍=\(𝑩1​𝝍,…,𝑩d​𝝍\)∈ℝd×d×d\\bm\{B\}\_\{\\bm\{\\psi\}\}=\(\\bm\{B\}\_\{1\\bm\{\\psi\}\},\\dots,\\bm\{B\}\_\{d\\bm\{\\psi\}\}\)\\in\\mathbb\{R\}^\{d\\times d\\times d\}is one such quantity\. It is shared over the state space and enters the reversible operator reparameterization only through the skew matrices𝑨i​𝝍\\bm\{A\}\_\{i\\bm\{\\psi\}\}\. For each fixed epistemic index, the perturbation should be a single state\-independent tensor whereas a standard epinet does not enforce this property, and soft penalizing[30](https://arxiv.org/html/2608.12624#bib.bib27)would only encourage, rather than guarantee, input independence\. Alternatively, we propose to use a linear epinet branch driven only by the epistemic index,

𝑩ϑ​\(𝐳\)\\displaystyle\\bm\{B\}\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{z\}\)=𝑩𝝍\+𝝈ϕB​\(𝐳\),\\displaystyle=\\bm\{B\}\_\{\\bm\{\\psi\}\}\+\\bm\{\\sigma\}\_\{\\bm\{\\phi\}\}^\{B\}\(\\mathbf\{z\}\),\(14\)𝝈ϕB​\(𝐳\)\\displaystyle\\bm\{\\sigma\}\_\{\\bm\{\\phi\}\}^\{B\}\(\\mathbf\{z\}\)=𝝈ϕB,learn​\(𝐳\)\+w​𝝈B,prior​\(𝐳\),\\displaystyle=\\bm\{\\sigma\}\_\{\\bm\{\\phi\}\}^\{B,\\mathrm\{learn\}\}\(\\mathbf\{z\}\)\+w\\,\\bm\{\\sigma\}^\{B,\\mathrm\{prior\}\}\(\\mathbf\{z\}\),𝝈ϕB,learn​\(𝐳\)\\displaystyle\\bm\{\\sigma\}\_\{\\bm\{\\phi\}\}^\{B,\\mathrm\{learn\}\}\(\\mathbf\{z\}\)=∑n=1dzznϕnB,𝝈B,prior\(𝐳\)=∑n=1dzzn𝜻nB\.\\displaystyle=\\sum\_\{n=1\}^\{d\_\{z\}\}z\_\{n\}\\,\\bm\{\\phi\}\_\{n\}^\{B\},\\qquad\\bm\{\\sigma\}^\{B,\\mathrm\{prior\}\}\(\\mathbf\{z\}\)=\\sum\_\{n=1\}^\{d\_\{z\}\}z\_\{n\}\\,\\bm\{\\zeta\}\_\{n\}^\{B\}\.HereϕnB,𝜻nB∈ℝd×d×d\\bm\{\\phi\}\_\{n\}^\{B\},\\bm\{\\zeta\}\_\{n\}^\{B\}\\in\\mathbb\{R\}^\{d\\times d\\times d\}have the same shape as𝑩𝝍\\bm\{B\}\_\{\\bm\{\\psi\}\}\. The tensorsϕB=\{ϕnB\}n=1dz\\bm\{\\phi\}^\{B\}=\\\{\\bm\{\\phi\}\_\{n\}^\{B\}\\\}\_\{n=1\}^\{d\_\{z\}\}are included in the epinet parameter setϕ\\bm\{\\phi\}, whereas𝜻B=\{𝜻nB\}n=1dz\\bm\{\\zeta\}^\{B\}=\\\{\\bm\{\\zeta\}\_\{n\}^\{B\}\\\}\_\{n=1\}^\{d\_\{z\}\}are fixed prior coefficients and are not trained\. For each realization, the skew matrices used in the reversible operator reparameterization are

𝑨i​ϑ\(𝐳\)=𝑩i​ϑ\(𝐳\)−𝑩i​ϑ\(𝐳\)⊤,i=1,…,d\.\\bm\{A\}\_\{i\\bm\{\\vartheta\}\}\(\\mathbf\{z\}\)=\\bm\{B\}\_\{i\\bm\{\\vartheta\}\}\(\\mathbf\{z\}\)\-\\bm\{B\}\_\{i\\bm\{\\vartheta\}\}\(\\mathbf\{z\}\)^\{\\top\},\\qquad i=1,\\dots,d\.\(15\)
In this construction, the state\-dependent branch perturbs the reversible operatorL~𝝍​\(𝐱\)\\tilde\{L\}\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\)as a function of\(𝐱,𝐳\)\(\\mathbf\{x\},\\mathbf\{z\}\), whereas the input\-independent branch perturbs the global matrices𝑩𝝍\\bm\{B\}\_\{\\bm\{\\psi\}\}as a function of𝐳\\mathbf\{z\}alone\. Because both branches are driven by the same epistemic index, a single realization jointly perturbs the reversible operatorL~ϑ​\(𝐱,𝐳\)\\tilde\{L\}\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{z\}\)and the skew\-symmetric matrices\{𝑨i​ϑ​\(𝐳\)\}i=1d\\\{\\bm\{A\}\_\{i\\bm\{\\vartheta\}\}\(\\mathbf\{z\}\)\\\}\_\{i=1\}^\{d\}used in the reparameterization\. This preserves the state\-independent character of𝑩𝝍\\bm\{B\}\_\{\\bm\{\\psi\}\}while allowing the state\-dependent uncertainty to be expressed throughL~ϑ\\tilde\{L\}\_\{\\bm\{\\vartheta\}\}\. Although introduced here for the global matrices in the GENERIC reparameterization, the input\-independent branch can also be extended to global unknown quantities in inverse problems\. In[D](https://arxiv.org/html/2608.12624#A4), we further discuss the details, validate this extension on the Kraichnan–Orszag and Korteweg–de Vries benchmarks studied in[82](https://arxiv.org/html/2608.12624#bib.bib66), and compare the proposed method with B\-PINNs\-HMC[72](https://arxiv.org/html/2608.12624#bib.bib14),[82](https://arxiv.org/html/2608.12624#bib.bib66)\.

#### 2\.3\.2Structure\-preserving reparameterization

According to[53](https://arxiv.org/html/2608.12624#bib.bib44), the epinet framework does not prescribe a fixed architecture for the learnable and prior networks\. This flexibility is central in the present structure\-preserving setting, since the epinet attached to each thermodynamic building block can be selected to preserve the admissible function class required by the corresponding N\-GENNs component\. For the reversible contribution, the reversible operator epinet associated withL~𝝍\\tilde\{L\}\_\{\\bm\{\\psi\}\}uses unconstrained neural networks for both the learnable network and the prior network, and incorporates the two\-branch design of Section[2\.3\.1](https://arxiv.org/html/2608.12624#S2.SS3.SSS1)\. The energy and entropy epinets associated withE𝝍E\_\{\\bm\{\\psi\}\}andS𝝍S\_\{\\bm\{\\psi\}\}are also represented by unconstrained multilayer perceptrons \(MLPs\)\. Finally, the dissipation\-potential epinet must preserve the convexity conditions imposed in N\-GENNs\. We therefore use the PICNNs for the raw dissipation\-potential epinet to ensure the outputs are convex in𝐱∗\\mathbf\{x\}^\{\*\}\.

Specifically, to enforce the reversible degeneracy condition, we apply the reversible operator reparameterization at each epistemic index\. For sampled𝐳\\mathbf\{z\}, the base networks in Eq\. \([4](https://arxiv.org/html/2608.12624#S2.E4)\) are replaced bySϑS\_\{\\bm\{\\vartheta\}\},L~ϑ\\tilde\{L\}\_\{\\bm\{\\vartheta\}\}, and the matrices defined in Eq\. \([15](https://arxiv.org/html/2608.12624#S2.E15)\)\. This yields the reparameterized reversible operator

Lϑ​\(𝐱,𝐳\)=QSϑ​\(𝐱,𝐳\)⊤​L~ϑ​\(𝐱,𝐳\)​QSϑ​\(𝐱,𝐳\),L\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{z\}\)=Q\_\{S\_\{\\bm\{\\vartheta\}\}\}\(\\mathbf\{x\},\\mathbf\{z\}\)^\{\\top\}\\tilde\{L\}\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{z\}\)Q\_\{S\_\{\\bm\{\\vartheta\}\}\}\(\\mathbf\{x\},\\mathbf\{z\}\),\(16\)whereQSϑ​\(𝐱,𝐳\)Q\_\{S\_\{\\bm\{\\vartheta\}\}\}\(\\mathbf\{x\},\\mathbf\{z\}\)is constructed from\{𝑨i​ϑ​\(𝐳\)\}i=1d\\\{\\bm\{A\}\_\{i\\bm\{\\vartheta\}\}\(\\mathbf\{z\}\)\\\}\_\{i=1\}^\{d\}andD​Sϑ​\(𝐱,𝐳\)\\mathrm\{D\}S\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{z\}\)in the same manner as in the deterministic construction\.

For the irreversible contribution, we form the energy\-orthogonal projection using the augmented energy,

PEϑ​\(𝐱,𝐳\)=I−D​Eϑ​\(𝐱,𝐳\)​D​Eϑ​\(𝐱,𝐳\)⊤‖D​Eϑ​\(𝐱,𝐳\)‖22\.P\_\{E\_\{\\bm\{\\vartheta\}\}\}\(\\mathbf\{x\},\\mathbf\{z\}\)=I\-\\frac\{\\mathrm\{D\}E\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{z\}\)\\mathrm\{D\}E\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{z\}\)^\{\\top\}\}\{\\\|\\mathrm\{D\}E\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{z\}\)\\\|\_\{2\}^\{2\}\}\.\(17\)The augmented dissipation potential used in the S\-PENNs dynamics is obtained by applying the same projection\-based reparameterization to the augmented raw potential,

Ξϑ​\(𝐱,𝐱∗,𝐳\)\\displaystyle\\Xi\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{x\}^\{\*\},\\mathbf\{z\}\)=Ξ~ϑ​\(𝐱,PEϑ​\(𝐱,𝐳\)​𝐱∗,𝐳\)−Ξ~ϑ​\(𝐱,𝟎,𝐳\)−\(D𝐱∗​Ξ~ϑ​\(𝐱,𝐱∗,𝐳\)\|𝐱∗=𝟎\)⊤​PEϑ​\(𝐱,𝐳\)​𝐱∗\.\\displaystyle=\\tilde\{\\Xi\}\_\{\\bm\{\\vartheta\}\}\\\!\\left\(\\mathbf\{x\},P\_\{E\_\{\\bm\{\\vartheta\}\}\}\(\\mathbf\{x\},\\mathbf\{z\}\)\\mathbf\{x\}^\{\*\},\\mathbf\{z\}\\right\)\-\\tilde\{\\Xi\}\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{0\},\\mathbf\{z\}\)\-\\left\(\\left\.\\mathrm\{D\}\_\{\\mathbf\{x\}^\{\*\}\}\\tilde\{\\Xi\}\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{x\}^\{\*\},\\mathbf\{z\}\)\\right\|\_\{\\mathbf\{x\}^\{\*\}=\\mathbf\{0\}\}\\right\)^\{\\top\}P\_\{E\_\{\\bm\{\\vartheta\}\}\}\(\\mathbf\{x\},\\mathbf\{z\}\)\\mathbf\{x\}^\{\*\}\.\(18\)Together with the constrained epinet parameterizations described above, each epistemic realization remains in the same admissible function classes as the deterministic N\-GENNs building blocks\.

It is worth noting that the dissipation\-potential epinet preserves the required convexity only when their constrained components are combined with nonnegative coefficients\. Since the epinet outputs in Eqs\. \([10](https://arxiv.org/html/2608.12624#S2.E10)\)–\([11](https://arxiv.org/html/2608.12624#S2.E11)\) are weighted by the epistemic index, the coefficients multiplying the dissipation\-potential components must be nonnegative\. Since the reference distributionπ⁡\(𝐳\)\\pi\(\\mathbf\{z\}\)can be freely chosen as stated in[53](https://arxiv.org/html/2608.12624#bib.bib44),π⁡\(𝐳\)\\pi\(\\mathbf\{z\}\)having nonnegative support is picked here\.

The constrained augmented building blocks define the S\-PENNs vector field

gϑ​\(𝐱,𝐳\)=Lϑ​\(𝐱,𝐳\)​D​Eϑ​\(𝐱,𝐳\)\+D𝐱∗​Ξϑ​\(𝐱,𝐱∗,𝐳\)\|𝐱∗=D​Sϑ​\(𝐱,𝐳\)\.g\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{z\}\)=L\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{z\}\)\\mathrm\{D\}E\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{z\}\)\+\\left\.\\mathrm\{D\}\_\{\\mathbf\{x\}^\{\*\}\}\\Xi\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{x\}^\{\*\},\\mathbf\{z\}\)\\right\|\_\{\\mathbf\{x\}^\{\*\}=\\mathrm\{D\}S\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\},\\mathbf\{z\}\)\}\.\(19\)Since the reversible and irreversible constraints are built into the augmented parameterization,gϑ​\(⋅,𝐳\)g\_\{\\bm\{\\vartheta\}\}\(\\cdot,\\mathbf\{z\}\)is thermodynamically consistent for every fixed draw of𝐳\\mathbf\{z\}\.

#### 2\.3\.3Training and predictive sampling

Let𝒟train=\{𝐗i\}i=1Ntrain\\mathcal\{D\}\_\{\\mathrm\{train\}\}=\\\{\\mathbf\{X\}\_\{i\}\\\}\_\{i=1\}^\{N\_\{\\mathrm\{train\}\}\}denote the training dataset, where𝐗i=\(𝐱i\(t\)\)t=0Nt\\mathbf\{X\}\_\{i\}=\(\\mathbf\{x\}\_\{i\}^\{\(t\)\}\)\_\{t=0\}^\{N\_\{t\}\}is theii\-th reference rollout and𝐱i\(0\)\\mathbf\{x\}\_\{i\}^\{\(0\)\}is its initial state\. S\-PENNs are trained in two stages\. First, the deterministic base networks are fitted by minimizing

ℒS−PENNsbase​\(𝝍\)=1Ntrain​\(Nt\+1\)​∑i=1Ntrain∑t=0Nt‖g𝝍​\(𝐱i\(t\)\)−𝐱˙i\(t\)‖22,\\mathcal\{L\}^\{\\mathrm\{base\}\}\_\{\\mathrm\{S\-PENNs\}\}\(\\bm\{\\psi\}\)=\\frac\{1\}\{N\_\{\\mathrm\{train\}\}\(N\_\{t\}\+1\)\}\\sum\_\{i=1\}^\{N\_\{\\mathrm\{train\}\}\}\\sum\_\{t=0\}^\{N\_\{t\}\}\\left\\\|g\_\{\\bm\{\\psi\}\}\(\\mathbf\{x\}\_\{i\}^\{\(t\)\}\)\-\\dot\{\\mathbf\{x\}\}\_\{i\}^\{\(t\)\}\\right\\\|\_\{2\}^\{2\},\(20\)where𝐱˙i\(t\)\\dot\{\\mathbf\{x\}\}\_\{i\}^\{\(t\)\}denotes the ground\-truth time derivative\. Second, the base\-network parameters𝝍\\bm\{\\psi\}are frozen and only the epinet parametersϕ\\bm\{\\phi\}are optimized\. For each trajectoryii, an epistemic index𝐳i∼π\\mathbf\{z\}\_\{i\}\\sim\\piis sampled and kept fixed over the rollout\. The epinet objective is

ℒS−PENNsepinet​\(ϕ\)=𝔼𝐳i​∼i\.i\.d\.​π​\[1Ntrain​\(Nt\+1\)​∑i=1Ntrain∑t=0Nt‖gϑ​\(𝐱i\(t\),𝐳i\)−𝐱˙i\(t\)‖22\]\.\\mathcal\{L\}^\{\\mathrm\{epinet\}\}\_\{\\mathrm\{S\-PENNs\}\}\(\\bm\{\\phi\}\)=\\mathbb\{E\}\_\{\\mathbf\{z\}\_\{i\}\\overset\{\\mathrm\{i\.i\.d\.\}\}\{\\sim\}\\pi\}\\left\[\\frac\{1\}\{N\_\{\\mathrm\{train\}\}\(N\_\{t\}\+1\)\}\\sum\_\{i=1\}^\{N\_\{\\mathrm\{train\}\}\}\\sum\_\{t=0\}^\{N\_\{t\}\}\\left\\\|g\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\}\_\{i\}^\{\(t\)\},\\mathbf\{z\}\_\{i\}\)\-\\dot\{\\mathbf\{x\}\}\_\{i\}^\{\(t\)\}\\right\\\|\_\{2\}^\{2\}\\right\]\.\(21\)At test time, letNsN\_\{s\}denote the number of predictive samples\. For each testing trajectoryii, independent draws\{𝐳i\(r\)\}r=1Ns\\\{\\mathbf\{z\}\_\{i\}^\{\(r\)\}\\\}\_\{r=1\}^\{N\_\{s\}\}generate predictive rollouts\{𝐗^i\(r\)\}r=1Ns\\\{\\hat\{\\mathbf\{X\}\}\_\{i\}^\{\(r\)\}\\\}\_\{r=1\}^\{N\_\{s\}\}initialized at𝐱i\(0\)\\mathbf\{x\}\_\{i\}^\{\(0\)\}\. This ensemble is summarized by the empirical predictive meanμ^ϑ\\hat\{\\mu\}\_\{\\bm\{\\vartheta\}\}and componentwise predictive standard deviationσ^ϑ\\hat\{\\sigma\}\_\{\\bm\{\\vartheta\}\}\.

### 2\.4Post\-hoc calibration via split conformal prediction

The sampled S\-PENNs realizations are thermodynamically consistent by construction, but their empirical uncertainty summaries, the meanμ^ϑ\\hat\{\\mu\}\_\{\\bm\{\\vartheta\}\}and componentwise standard deviationσ^ϑ\\hat\{\\sigma\}\_\{\\bm\{\\vartheta\}\}, do not by themselves guarantee to be well\-calibrated\. Standard regression calibration techniques, including variance scaling[42](https://arxiv.org/html/2608.12624#bib.bib68)and CDF\-based corrections[38](https://arxiv.org/html/2608.12624#bib.bib69),[76](https://arxiv.org/html/2608.12624#bib.bib70), are less suitable here because they need additional model fitting and can modify the implied predictive distribution\. We instead use split conformal prediction[56](https://arxiv.org/html/2608.12624#bib.bib61),[41](https://arxiv.org/html/2608.12624#bib.bib59)to calibrate the prediction intervals\. This post\-processing step leaves the sampled trajectories unchanged and gives distribution\-free finite\-sample marginal coverage under the exchangeability assumption stated next\.

The calibration and testing samples used in the split conformal prediction are required to be jointly exchangeable[56](https://arxiv.org/html/2608.12624#bib.bib61),[41](https://arxiv.org/html/2608.12624#bib.bib59)\. For dynamical systems, the exchangeable unit is the full trajectory rather than an individual timestep, since time levels within one rollout are coupled by the evolution equations\. Following trajectory\-level conformal treatments for dynamical systems[66](https://arxiv.org/html/2608.12624#bib.bib63),[67](https://arxiv.org/html/2608.12624#bib.bib65),[19](https://arxiv.org/html/2608.12624#bib.bib60), we regard each reference rollout𝐗i\\mathbf\{X\}\_\{i\}as one sample\. Using the notation of Section[2\.3\.3](https://arxiv.org/html/2608.12624#S2.SS3.SSS3), the available trajectories are partitioned into three disjoint datasets: a training dataset𝒟train\\mathcal\{D\}\_\{\\mathrm\{train\}\}for training the S\-PENNs \(Section[2\.3\.3](https://arxiv.org/html/2608.12624#S2.SS3.SSS3)\), a calibration dataset𝒟cal\\mathcal\{D\}\_\{\\mathrm\{cal\}\}ofNcalN\_\{\\mathrm\{cal\}\}trajectories reserved for calibration, and a testing dataset𝒟test\\mathcal\{D\}\_\{\\mathrm\{test\}\}used only for evaluation\. In the numerical examples below, this assumption is satisfied by drawing the initial conditions independently from the same distribution for the calibration and testing datasets\.

For each calibration trajectory𝐗i=\(𝐱i\(t\)\)t=0Nt∈𝒟cal\\mathbf\{X\}\_\{i\}=\(\\mathbf\{x\}\_\{i\}^\{\(t\)\}\)\_\{t=0\}^\{N\_\{t\}\}\\in\\mathcal\{D\}\_\{\\mathrm\{cal\}\}, initialized at𝐱i\(0\)\\mathbf\{x\}\_\{i\}^\{\(0\)\}, letj∈ℐj\\in\\mathcal\{I\}index one scalar entry of the corresponding forecast rollout\. For ODE trajectories,jjspecifies a state variable and a noninitial time point; for field\-valued problems, it specifies a state variable and a noninitial space–time grid point\. The initial state is prescribed as an input to the predictor and is therefore excluded fromℐ\\mathcal\{I\}\. Similarly, prescribed boundary conditions are also excluded fromℐ\\mathcal\{I\}\. Calibration is applied only to the forecast portion of the rollout, namely the time pointst=1,…,Ntt=1,\\dots,N\_\{t\}for trajectory\-valued problems and the corresponding non\-prescribed space–time grid points for field\-valued problems\. We define the normalized nonconformity score[2](https://arxiv.org/html/2608.12624#bib.bib58)

si\(j\)=\|𝐗i​\(j\)−μ^ϑ​\(𝐱i\(0\),j\)\|σ^ϑ​\(𝐱i\(0\),j\),i=1,…,Ncal\.s\_\{i\}\(j\)=\\frac\{\\left\|\\mathbf\{X\}\_\{i\}\(j\)\-\\hat\{\\mu\}\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\}\_\{i\}^\{\(0\)\};j\)\\right\|\}\{\\hat\{\\sigma\}\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\}\_\{i\}^\{\(0\)\};j\)\},\\qquad i=1,\\dots,N\_\{\\mathrm\{cal\}\}\.\(22\)Then, for each fixed indexj∈ℐj\\in\\mathcal\{I\}, lets\(1\)​\(j\)≤⋯≤s\(Ncal\)​\(j\)s\_\{\(1\)\}\(j\)\\leq\\cdots\\leq s\_\{\(N\_\{\\mathrm\{cal\}\}\)\}\(j\)denote the ordered calibration scores atjj, obtained by sorting the normalized nonconformity scores\{si​\(j\)\}i=1Ncal\\\{s\_\{i\}\(j\)\\\}\_\{i=1\}^\{N\_\{\\mathrm\{cal\}\}\}\. For a target miscoverage levelα∈\(0,1\)\\alpha\\in\(0,1\), set

rα=⌈\(Ncal\+1\)​\(1−α\)⌉,q^1−α​\(j\)=\{s\(rα\)​\(j\),rα≤Ncal,\+∞,rα=Ncal\+1\.r\_\{\\alpha\}=\\left\\lceil\(N\_\{\\mathrm\{cal\}\}\+1\)\(1\-\\alpha\)\\right\\rceil,\\qquad\\hat\{q\}\_\{1\-\\alpha\}\(j\)=\\begin\{cases\}s\_\{\(r\_\{\\alpha\}\)\}\(j\),&r\_\{\\alpha\}\\leq N\_\{\\mathrm\{cal\}\},\\\\ \+\\infty,&r\_\{\\alpha\}=N\_\{\\mathrm\{cal\}\}\+1\.\\end\{cases\}\(23\)Here,rαr\_\{\\alpha\}is the finite\-sample\-adjusted integer rank andq^1−α​\(j\)\\hat\{q\}\_\{1\-\\alpha\}\(j\)is the corresponding split\-conformal score threshold atjj\. For a testing trajectory𝐗=\(𝐱\(t\)\)t=0Nt\\mathbf\{X\}=\(\\mathbf\{x\}^\{\(t\)\}\)\_\{t=0\}^\{N\_\{t\}\}, this threshold rescales the predictive standard deviation to give

C^1−α​\(𝐱\(0\),j\)=\[μ^ϑ​\(𝐱\(0\),j\)−q^1−α​\(j\)​σ^ϑ​\(𝐱\(0\),j\),μ^ϑ​\(𝐱\(0\),j\)\+q^1−α​\(j\)​σ^ϑ​\(𝐱\(0\),j\)\]\.\\widehat\{C\}\_\{1\-\\alpha\}\(\\mathbf\{x\}^\{\(0\)\};j\)=\\left\[\\hat\{\\mu\}\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\}^\{\(0\)\};j\)\-\\hat\{q\}\_\{1\-\\alpha\}\(j\)\\hat\{\\sigma\}\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\}^\{\(0\)\};j\),\\;\\hat\{\\mu\}\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\}^\{\(0\)\};j\)\+\\hat\{q\}\_\{1\-\\alpha\}\(j\)\\hat\{\\sigma\}\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{x\}^\{\(0\)\};j\)\\right\]\.\(24\)
For any testing trajectory that is jointly exchangeable with the calibration trajectories, the interval in Eq\. \([24](https://arxiv.org/html/2608.12624#S2.E24)\) satisfies, for every fixedj∈ℐj\\in\\mathcal\{I\},

ℙ⁡\(𝐗⁡\(j\)∈C^1−α​\(𝐱\(0\),j\)\)≥1−α\.\\mathbb\{P\}\\\!\\left\(\\mathbf\{X\}\(j\)\\in\\widehat\{C\}\_\{1\-\\alpha\}\(\\mathbf\{x\}^\{\(0\)\};j\)\\right\)\\geq 1\-\\alpha\.\(25\)This is the standard split\-conformal rank argument[61](https://arxiv.org/html/2608.12624#bib.bib62),[2](https://arxiv.org/html/2608.12624#bib.bib58): exchangeability of the trajectories implies exchangeability of theNcalN\_\{\\mathrm\{cal\}\}calibration scores and the corresponding test score at the same indexjj\. Therefore, the test score exceeds the empirical quantileq^1−α​\(j\)\\hat\{q\}\_\{1\-\\alpha\}\(j\)with probability at mostα\\alpha\. A proof of this result is given in Appendix D of[2](https://arxiv.org/html/2608.12624#bib.bib58), following[56](https://arxiv.org/html/2608.12624#bib.bib61)\. This finite\-sample statement is distribution\-free in the sense that it does not assume a parametric data\-generating law, a correctly specified predictive distribution, or a particular residual model\. Its only statistical requirement is the exchangeability condition above\.

## 3Numerical examples

In this section, we evaluate the proposed S\-PENNs framework against two benchmark UQ methods, deep ensembles and MC dropout, on three numerical examples of increasing complexity\. As discussed in Section[2\.3](https://arxiv.org/html/2608.12624#S2.SS3), the reference distribution for the epistemic index𝐳\\mathbf\{z\}is a modeling choice\. We therefore consider three cases: a half\-normal distribution, denoted “S\-PENNs \(half\-normal\)”; a uniform distribution on\[0,5\]\[0,5\], denoted “S\-PENNs \(uniform\)”; and a unit\-rate exponential distribution truncated to\[0,5\]\[0,5\], denoted “S\-PENNs \(exponential\)”\. For quantitative evaluation, we consider the following metrics\. Predictive accuracy for each state variable is measured by the relativeℓ2\\ell^\{2\}error \(RL2E\) of the predictive mean over the simulated time interval or over the space–time grid, and by pointwise absolute errors for field\-valued quantities\. Probabilistic performance is evaluated with the continuous ranked probability score \(CRPS\) and the energy score \(ES\), both proper scoring rules for predictive distributions[16](https://arxiv.org/html/2608.12624#bib.bib8)\. The CRPS is evaluated on scalar marginals and provides a componentwise measure expressed in the physical units of each predicted variable[47](https://arxiv.org/html/2608.12624#bib.bib10),[28](https://arxiv.org/html/2608.12624#bib.bib11)\. The ES extends this assessment to vector\-valued probabilistic forecasts and is used here to evaluate the predictive distribution of the full rollout[17](https://arxiv.org/html/2608.12624#bib.bib12)\. This is appropriate in the present setting because the learned dynamics couple the state variables, whereas scalar marginal scores do not assess their joint behavior\. CRPS and ES are evaluated on the uncalibrated predictive samples from all methods\. Meanwhile, the calibration quality of split conformal prediction is separately assessed through the empirical coverage of the calibrated S\-PENNs prediction intervals\. Unless otherwise stated, the calibrated prediction intervals shown in the trajectory and field plots useα=0\.05\\alpha=0\.05, corresponding to a target marginal coverage level of1−α=95%1\-\\alpha=95\\%\. Definitions of these metrics are given in[A](https://arxiv.org/html/2608.12624#A1)\.

### 3\.1Harmonic oscillator in a heat bath

We first consider the harmonic oscillator example used in the N\-GENNs study[69](https://arxiv.org/html/2608.12624#bib.bib31)\. This low\-dimensional problem has a closed\-form GENERIC representation, and its irreversible response is generated by a quadratic dissipation potential\. It therefore provides a controlled setting for testing whether the proposed UQ construction can perturb a learned thermodynamic model without violating its conservation and dissipation properties\.

The state of the system is given by𝐱=\(q,p,ϵ\)\\mathbf\{x\}=\(q,p,\\epsilon\), whereqqandppare the oscillator position and conjugate momentum, respectively, andϵ\\epsilonrepresents the internal energy of the heat bath\. The total energy and the canonical Poisson operator are given by

E⁡\(𝐱\)=p22​m\+12​k​q2\+ϵ,L=\(010−100000\)\.E\(\\mathbf\{x\}\)=\\frac\{p^\{2\}\}\{2m\}\+\\frac\{1\}\{2\}kq^\{2\}\+\\epsilon,\\qquad L=\\begin\{pmatrix\}0&1&0\\\\ \-1&0&0\\\\ 0&0&0\\end\{pmatrix\}\.\(26\)For a heat bath at constant temperatureTbathT\_\{\\mathrm\{bath\}\}, the entropy and the mobility matrix are

S⁡\(𝐱\)=ϵTbath,M⁡\(𝐱\)=γ​Tbath​\(00001−pm0−pmp2m2\),S\(\\mathbf\{x\}\)=\\frac\{\\epsilon\}\{T\_\{\\mathrm\{bath\}\}\},\\qquad M\(\\mathbf\{x\}\)=\\gamma T\_\{\\mathrm\{bath\}\}\\begin\{pmatrix\}0&0&0\\\\ 0&1&\-\\frac\{p\}\{m\}\\\\ 0&\-\\frac\{p\}\{m\}&\\frac\{p^\{2\}\}\{m^\{2\}\}\\end\{pmatrix\},\(27\)wheremmis the mass,kkis the stiffness of the linear spring connecting the mass to a fixed support, andγ\\gammais the damping coefficient\. The dissipation potential isΞ⁡\(𝐱,𝐱∗\)=12​\(𝐱∗\)T​M​\(𝐱\)​𝐱∗\\Xi\(\\mathbf\{x\},\\mathbf\{x\}^\{\*\}\)=\\frac\{1\}\{2\}\(\\mathbf\{x\}^\{\*\}\)^\{T\}M\(\\mathbf\{x\}\)\\mathbf\{x\}^\{\*\}\. Altogether, these building blocks give the evolution equations

q˙=pm,p˙=−k​q−γ​pm,ϵ˙=γ​p2m2\.\\dot\{q\}=\\frac\{p\}\{m\},\\qquad\\dot\{p\}=\-kq\-\\gamma\\frac\{p\}\{m\},\\qquad\\dot\{\\epsilon\}=\\gamma\\frac\{p^\{2\}\}\{m^\{2\}\}\.\(28\)
#### 3\.1\.1Data generation and model validation

For this example, the dataset is generated by using an implicit midpoint integrator for the ODE system above\. The simulations use a uniform timestepΔ​t=0\.015\\Delta t=0\.015and physical parametersm=1\.0m=1\.0,k=1\.2k=1\.2, andγ=0\.4\\gamma=0\.4\. Initial states are drawn uniformly fromq⁡\(0\)∈\[−1\.0,1\.0\]q\(0\)\\in\[\-1\.0,\\,1\.0\]andp⁡\(0\)∈\[−1\.0,1\.0\]p\(0\)\\in\[\-1\.0,\\,1\.0\], while the heat\-bath energy is initialized asϵ⁡\(0\)=0\\epsilon\(0\)=0\. The generated dataset consists of100100trajectories with10001000state records per trajectory\. The trajectory\-level partition follows the protocol described in Section[2\.4](https://arxiv.org/html/2608.12624#S2.SS4):6060trajectories are used for training,5050for calibration, and4040for testing, where the splitting is performed randomly\.

Model validation is performed by rolling out each trained model from the initial conditions of the testing trajectories and comparing the predictions with the corresponding simulated reference solutions\. All rollout predictions use a standard fourth\-order Runge–Kutta \(RK4\) integrator with the same timestep as the generated dataset,Δ​t=0\.015\\Delta t=0\.015\. S\-PENNs and MC dropout useNs=2000N\_\{s\}=2000UQ realizations for each testing trajectory, while the deep ensembles use5050independently trained N\-GENNs with different architecture choices and random initializations\. The network architectures, optimization settings, and inference parameters are summarized in Table[4](https://arxiv.org/html/2608.12624#A3.T4)of[C](https://arxiv.org/html/2608.12624#A3)\.

#### 3\.1\.2Results and discussion

![Refer to caption](https://arxiv.org/html/2608.12624v1/figs/results/harmonicOscillatorMain.png)Figure 2:Predictive trajectories and calibration diagnostics for the harmonic oscillator example\. The first three columns compare deep ensembles, MC dropout, and S\-PENNs \(half\-normal\) on a representative testing trajectory; rows correspond toqq,pp, andϵ\\epsilon\. In each panel, the predictive mean is shown as a solid line, the reference trajectory as a dashed line, and the uncalibratedμ^±2​σ^\\hat\{\\mu\}\\pm 2\\hat\{\\sigma\}band as a blue shaded region\. For S\-PENNs \(half\-normal\), the calibrated 95% conformal prediction interval is overlaid in red\. The rightmost column plots empirical coverage against target coverage for S\-PENNs \(half\-normal\), with the diagonal indicating ideal calibration\. Titles report the relativeℓ2\\ell^\{2\}error \(RL2E\) for each method and state variable\.Fig\.[2](https://arxiv.org/html/2608.12624#S3.F2)compares the predictions of deep ensembles, MC dropout, and S\-PENNs \(half\-normal\) on a representative testing trajectory\. The first three columns correspond to the different methods, while the rows show the results for the state variablesqq,pp, andϵ\\epsilon\. Each panel displays the predictive trajectory together with its associated uncertainty band\. The rightmost column reports empirical coverage for the uncalibrated and calibrated S\-PENNs \(half\-normal\) intervals\. All three methods reproduce the damped oscillator response, including the oscillatory exchange betweenqqandppand the monotonic increase in the bath energyϵ\\epsilon\. Deep ensembles give the smallest relativeℓ2\\ell^\{2\}errors on this trajectory, with2\.25×10−22\.25\\times 10^\{\-2\}forqq,1\.67×10−21\.67\\times 10^\{\-2\}forpp, and1\.54×10−31\.54\\times 10^\{\-3\}forϵ\\epsilon\. S\-PENNs \(half\-normal\) remain accurate on this displayed trajectory, with errors of3\.019×10−23\.019\\times 10^\{\-2\}forqq,2\.20×10−22\.20\\times 10^\{\-2\}forpp, and1\.24×10−21\.24\\times 10^\{\-2\}forϵ\\epsilon\. MC dropout is the worst of the three methods for all three variables, giving errors of6\.43×10−26\.43\\times 10^\{\-2\}forqq,5\.86×10−25\.86\\times 10^\{\-2\}forpp, and1\.38×10−21\.38\\times 10^\{\-2\}forϵ\\epsilon\. The uncertainty bands in the trajectory panels provide a qualitative comparison of predictive spread\. MC dropout produces the broadest bands, whereas the other two methods produce narrower bands\. The rightmost panels show that the uncalibrated S\-PENNs deviate from the target coverage, as empirical coverage can lie above or below the target depending on the state variable and coverage level\. Split conformal calibration moves the curves toward the ideal diagonal, improving coverage without changing the sampled trajectories\.

![Refer to caption](https://arxiv.org/html/2608.12624v1/figs/results/harmonicCombinedDiagnostics.png)Figure 3:Predictive accuracy and structure\-preservation diagnostics for the harmonic oscillator example\. The top panel shows the componentwise trajectory\-wise relativeℓ2\\ell^\{2\}error \(RL2E\) over the testing dataset forqq,pp, andϵ\\epsilon, comparing deep ensembles, MC dropout, and S\-PENNs using half\-normal, uniform, and exponential distributions for the epistemic index\. The box plots show the median with the solid lines and the boxes span the interquartile ranges\. The bottom\-left panel shows the energy rateE˙\\dot\{E\}over time, with the median as a solid line and the 95% prediction interval as a shaded band across sampled rollouts\. The bottom\-right panel shows the pointwise minimum entropy productionmin⁡\(S˙\)\\min\(\\dot\{S\}\)over time on a logarithmic scale, where the minimum is taken across the sampled S\-PENNs rollouts at each timestep\.As a broader assessment, Fig\.[3](https://arxiv.org/html/2608.12624#S3.F3)extends the comparison to the full testing dataset and examines the thermodynamic structure of the S\-PENNs rollouts\. The top panel reports the trajectory\-wise relativeℓ2\\ell^\{2\}errors\. All three S\-PENNs variants have lower median errors than MC dropout forqqandpp\. Among them, the half\-normal distribution gives the lowest median errors for all three state variables\. Its medians forqqandppare nearly identical to those of deep ensembles, while deep ensembles retain the lowest median error forϵ\\epsilon\. The uniform and exponential variants of S\-PENNs have higher median errors onϵ\\epsilonthan MC dropout\. The bottom panels assess the first and second laws for the sampled rollouts of the representative testing trajectory\. For all three S\-PENNs variants, the median energy rate remains centered at zero, and the 95% intervals stay at the10−710^\{\-7\}scale\. The pointwise minimum entropy\-production rate over the sampled rollouts remains positive throughout the trajectory and decreases overall as the damped oscillator approaches equilibrium\. These results indicate that, in this example, S\-PENNs preserve thermodynamic consistency across all three reference distributions, while predictive accuracy remains stable, with slightly greater variation forϵ\\epsilonthan forqqandpp\.

Beyond predictive accuracy and thermodynamic consistency, computational cost and the statistical quality of the resulting uncertainty estimates are central to the practical value of the UQ method\. Table[1](https://arxiv.org/html/2608.12624#S3.T1)compares serial wall time and proper scoring rules computed from the uncalibrated predictive samples\. Deep ensembles attain the lowest ES and CRPS values, while all S\-PENNs variants outperform MC dropout in ES and in CRPS forqqandpp; the half\-normal variant also improves CRPS forϵ\\epsilon\. Measured wall times are9\.1×1039\.1\\times 10^\{3\}s for the 50\-member deep ensemble,2\.3×1022\.3\\times 10^\{2\}s for MC dropout, and1\.2×1021\.2\\times 10^\{2\}–1\.3×1021\.3\\times 10^\{2\}s for S\-PENNs\. Thus, S\-PENNs are less costly despite using 2000 rather than 50 samples\. Scaling the deep ensemble to 2000 members gives3\.6×1053\.6\\times 10^\{5\}s, which is three orders of magnitude higher than the cost for S\-PENNs\.

Table 1:Serial wall time and proper scoring rules for the harmonic oscillator example\. The wall time includes model training and generation of predictive samples for all testing trajectories\. ES denotes the energy score over the full multivariate trajectory, and CRPS is averaged for each state variable\. Scores are computed from uncalibrated predictive samples and averaged over the testing dataset\. Lower ES and CRPS values indicate better probabilistic predictions\.∗For deep ensembles, the reported wall time is the serial cost of the 50\-member ensemble, and the scoring rules are computed from the same ensemble\.

Overall, the harmonic oscillator example shows that S\-PENNs can introduce epistemic uncertainty into a hard\-constrained learned dynamics model at a low computational cost while guaranteeing thermodynamic consistency\. Across the testing dataset, S\-PENNs produce accurate mean predictions, improve clearly over MC dropout, and remain largely insensitive to the choice of reference distribution\. The sampled rollouts preserve energy conservation and nonnegative entropy production, and split conformal prediction improves empirical coverage without modifying these admissible samples\. Although deep ensembles remain the strongest baseline in accuracy and proper scoring rules for this low\-dimensional problem, S\-PENNs achieve competitive predictive accuracy and UQ performance at substantially lower computational cost\.

### 3\.2Idealized chemical motor

The second example is the idealized chemical motor considered in the N\-GENNs study[69](https://arxiv.org/html/2608.12624#bib.bib31)\. A piston of massmmand cross\-sectional areaAAis attached to a linear spring with natural lengthq0q\_\{0\}and encloses a van der Waals mixture undergoing a chemical reactionX1\+X2⇌X3\\mathrm\{X\}\_\{1\}\+\\mathrm\{X\}\_\{2\}\\rightleftharpoons\\mathrm\{X\}\_\{3\}\. The absolute piston positionqqdetermines the chamber volume throughV⁡\(q\)=A​qV\(q\)=Aq, whileq−q0q\-q\_\{0\}is the spring extension\. The state variables and vector of species mole numbers are

𝐱=\(q,p,ϵ,n1,n2,n3\),𝒏=\(n1,n2,n3\)T,\\mathbf\{x\}=\(q,p,\\epsilon,n\_\{1\},n\_\{2\},n\_\{3\}\),\\qquad\\bm\{n\}=\(n\_\{1\},n\_\{2\},n\_\{3\}\)^\{T\},\(29\)whereppis the piston momentum,ϵ\\epsilonis the internal energy of the gas mixture, andnμn\_\{\\mu\}is the number of moles of specieXμ\\mathrm\{X\}\_\{\\mu\}\. The total energy combines the kinetic energy of the piston, the elastic energy of the spring, and the internal energy of the mixture,

E⁡\(𝐱\)=p22​m\+12​k​\(q−q0\)2\+ϵ,E\(\\mathbf\{x\}\)=\\frac\{p^\{2\}\}\{2m\}\+\\frac\{1\}\{2\}k\(q\-q\_\{0\}\)^\{2\}\+\\epsilon,\(30\)wherekkis the spring stiffness\. The Poisson operator for this system is noncanonical and couples the piston momentum to the gas internal energy as,

L⁡\(𝐱\)=\(010000−10F⁡\(𝐱\)0000−F⁡\(𝐱\)0000000000000000000000\),F⁡\(𝐱\)=\(∂S/∂q\)p,ϵ\(∂S/∂ϵ\)p,q\.L\(\\mathbf\{x\}\)=\\begin\{pmatrix\}0&1&0&0&0&0\\\\ \-1&0&F\(\\mathbf\{x\}\)&0&0&0\\\\ 0&\-F\(\\mathbf\{x\}\)&0&0&0&0\\\\ 0&0&0&0&0&0\\\\ 0&0&0&0&0&0\\\\ 0&0&0&0&0&0\\end\{pmatrix\},\\qquad F\(\\mathbf\{x\}\)=\\frac\{\(\\partial S/\\partial q\)\_\{p,\\epsilon\}\}\{\(\\partial S/\\partial\\epsilon\)\_\{p,q\}\}\.\(31\)HereF⁡\(𝐱\)F\(\\mathbf\{x\}\)is the force exerted by the gas on the piston\.

We adopt the N\-GENNs dimensionless conventionR=h=NA=1R=h=N\_\{A\}=1, whereRR,hh, andNAN\_\{A\}denote the gas constant, Planck constant, and Avogadro constant, respectively\. WithS=Sphys/RS=S\_\{\\mathrm\{phys\}\}/R, the gas\-constant\-scaled van der Waals entropy is

S⁡\(𝐱\)=∑μ=13nμ​\[52\+ln⁡\(V−∑ν=13nν​bνnμ​\[4​π​mμ3​∑β=13nβ​\(ϵ\+∑ρ,σ=13nρ​nσ​aρ​σV\)\]3/2\)\],S\(\\mathbf\{x\}\)=\\sum\_\{\\mu=1\}^\{3\}n\_\{\\mu\}\\left\[\\frac\{5\}\{2\}\+\\ln\\left\(\\frac\{V\-\\sum\_\{\\nu=1\}^\{3\}n\_\{\\nu\}b\_\{\\nu\}\}\{n\_\{\\mu\}\}\\left\[\\frac\{4\\pi m\_\{\\mu\}\}\{3\\sum\_\{\\beta=1\}^\{3\}n\_\{\\beta\}\}\\left\(\\epsilon\+\\frac\{\\sum\_\{\\rho,\\sigma=1\}^\{3\}n\_\{\\rho\}n\_\{\\sigma\}a\_\{\\rho\\sigma\}\}\{V\}\\right\)\\right\]^\{3/2\}\\right\)\\right\],\(32\)heremμm\_\{\\mu\}is the molecular mass of speciesXμ\\mathrm\{X\}\_\{\\mu\},bμb\_\{\\mu\}is its excluded\-volume parameter, andaρ​σa\_\{\\rho\\sigma\}is the interaction parameter between speciesXρ\\mathrm\{X\}\_\{\\rho\}andXσ\\mathrm\{X\}\_\{\\sigma\}\. The mass\-action kinetics are generated by the non\-quadratic dissipation potential

Ξ⁡\(𝐱,𝐱∗\)=α​n1​n2​n3​\[cosh⁡\(n1∗\+n2∗−n3∗2\)−1\],\\Xi\(\\mathbf\{x\},\\mathbf\{x\}^\{\*\}\)=\\alpha\\sqrt\{n\_\{1\}n\_\{2\}n\_\{3\}\}\\left\[\\cosh\\left\(\\frac\{n\_\{1\}^\{\*\}\+n\_\{2\}^\{\*\}\-n\_\{3\}^\{\*\}\}\{2\}\\right\)\-1\\right\],\(33\)where𝐱∗=\(q∗,p∗,ϵ∗,n1∗,n2∗,n3∗\)\\mathbf\{x\}^\{\*\}=\(q^\{\*\},p^\{\*\},\\epsilon^\{\*\},n\_\{1\}^\{\*\},n\_\{2\}^\{\*\},n\_\{3\}^\{\*\}\)is the conjugate state, andα\>0\\alpha\>0is the rate coefficient\. Evaluating the GENERIC evolution with Eqs\. \([30](https://arxiv.org/html/2608.12624#S3.E30)\)–\([33](https://arxiv.org/html/2608.12624#S3.E33)\) gives

q˙=\\displaystyle\\dot\{q\}=\{\}pm,\\displaystyle\\frac\{p\}\{m\},\(34\)p˙=\\displaystyle\\dot\{p\}=\{\}−k⁡\(q−q0\)\+23​ϵ​A​q\+∑ρ,σ=13nρ​nσ​aρ​σq⁡\(A​q−∑ν=13nν​bν\)−∑ρ,σ=13nρ​nσ​aρ​σA​q2,\\displaystyle\-k\(q\-q\_\{0\}\)\+\\frac\{2\}\{3\}\\frac\{\\epsilon Aq\+\\sum\_\{\\rho,\\sigma=1\}^\{3\}n\_\{\\rho\}n\_\{\\sigma\}a\_\{\\rho\\sigma\}\}\{q\\left\(Aq\-\\sum\_\{\\nu=1\}^\{3\}n\_\{\\nu\}b\_\{\\nu\}\\right\)\}\-\\frac\{\\sum\_\{\\rho,\\sigma=1\}^\{3\}n\_\{\\rho\}n\_\{\\sigma\}a\_\{\\rho\\sigma\}\}\{Aq^\{2\}\},ϵ˙=\\displaystyle\\dot\{\\epsilon\}=\{\}−pm​\[23​ϵ​A​q\+∑ρ,σ=13nρ​nσ​aρ​σq⁡\(A​q−∑ν=13nν​bν\)−∑ρ,σ=13nρ​nσ​aρ​σA​q2\],\\displaystyle\-\\frac\{p\}\{m\}\\left\[\\frac\{2\}\{3\}\\frac\{\\epsilon Aq\+\\sum\_\{\\rho,\\sigma=1\}^\{3\}n\_\{\\rho\}n\_\{\\sigma\}a\_\{\\rho\\sigma\}\}\{q\\left\(Aq\-\\sum\_\{\\nu=1\}^\{3\}n\_\{\\nu\}b\_\{\\nu\}\\right\)\}\-\\frac\{\\sum\_\{\\rho,\\sigma=1\}^\{3\}n\_\{\\rho\}n\_\{\\sigma\}a\_\{\\rho\\sigma\}\}\{Aq^\{2\}\}\\right\],n˙1=\\displaystyle\\dot\{n\}\_\{1\}=\{\}α4​\(K​n3−n1​n2K\),\\displaystyle\\frac\{\\alpha\}\{4\}\\left\(Kn\_\{3\}\-\\frac\{n\_\{1\}n\_\{2\}\}\{K\}\\right\),n˙2=\\displaystyle\\dot\{n\}\_\{2\}=\{\}α4​\(K​n3−n1​n2K\),\\displaystyle\\frac\{\\alpha\}\{4\}\\left\(Kn\_\{3\}\-\\frac\{n\_\{1\}n\_\{2\}\}\{K\}\\right\),n˙3=\\displaystyle\\dot\{n\}\_\{3\}=\{\}α4​\(n1​n2K−K​n3\)\.\\displaystyle\\frac\{\\alpha\}\{4\}\\left\(\\frac\{n\_\{1\}n\_\{2\}\}\{K\}\-Kn\_\{3\}\\right\)\.HereKKis the reaction constant\. Withn=∑β=13nβn=\\sum\_\{\\beta=1\}^\{3\}n\_\{\\beta\},𝒂=\(aρ​σ\)\\bm\{a\}=\(a\_\{\\rho\\sigma\}\), and𝒃=\(b1,b2,b3\)T\\bm\{b\}=\(b\_\{1\},b\_\{2\},b\_\{3\}\)^\{T\}, it is given by

K=exp⁡\[3​n​∑σ=13\(a1​σ\+a2​σ−a3​σ\)​nσ2​\(ϵ​A​q\+𝒏T​𝒂​𝒏\)−n⁡\(b1\+b2−b3\)2​\(A​q−𝒏T​𝒃\)\]​\(m1​m2m3\)3/4​\(A​q−𝒏T​𝒃\)1/2​\[4​π3​n​\(ϵ\+𝒏T​𝒂​𝒏A​q\)\]3/4\.K=\\exp\\left\[\\frac\{3n\\sum\_\{\\sigma=1\}^\{3\}\(a\_\{1\\sigma\}\+a\_\{2\\sigma\}\-a\_\{3\\sigma\}\)n\_\{\\sigma\}\}\{2\(\\epsilon Aq\+\\bm\{n\}^\{T\}\\bm\{a\}\\bm\{n\}\)\}\-\\frac\{n\(b\_\{1\}\+b\_\{2\}\-b\_\{3\}\)\}\{2\(Aq\-\\bm\{n\}^\{T\}\\bm\{b\}\)\}\\right\]\\left\(\\frac\{m\_\{1\}m\_\{2\}\}\{m\_\{3\}\}\\right\)^\{3/4\}\(Aq\-\\bm\{n\}^\{T\}\\bm\{b\}\)^\{1/2\}\\left\[\\frac\{4\\pi\}\{3n\}\\left\(\\epsilon\+\\frac\{\\bm\{n\}^\{T\}\\bm\{a\}\\bm\{n\}\}\{Aq\}\\right\)\\right\]^\{3/4\}\.\(35\)More details on this example can be found in Section 4\.2\.1 and Appendix B of[69](https://arxiv.org/html/2608.12624#bib.bib31)\.

#### 3\.2\.1Data generation and model validation

For this example, the dataset is generated using the implicit midpoint integrator\. The simulations use a uniform timestepΔ​t=0\.01\\Delta t=0\.01, and the physical parameters

m=0\.8,A=1\.0,q0=2\.0,k=31\.6,α=2\.0,𝒂=\(0\.10\.10\.10\.10\.10\.10\.10\.120\.0\),𝒃=\(0\.05,0\.05,0\.05\)T,\(m1,m2,m3\)=\(1\.0,1\.2,2\.2\)\.\\begin\{gathered\}m=0\.8,\\quad A=1\.0,\\quad q\_\{0\}=2\.0,\\quad k=31\.6,\\quad\\alpha=2\.0,\\\\ \\bm\{a\}=\\begin\{pmatrix\}0\.1&0\.1&0\.1\\\\ 0\.1&0\.1&0\.1\\\\ 0\.1&0\.1&20\.0\\end\{pmatrix\},\\quad\\bm\{b\}=\(0\.05,0\.05,0\.05\)^\{T\},\\quad\(m\_\{1\},m\_\{2\},m\_\{3\}\)=\(1\.0,1\.2,2\.2\)\.\\end\{gathered\}\(36\)The piston starts from rest,p⁡\(0\)=0p\(0\)=0, while the initial mole numbers, piston position, and dimensionless temperatureT0T\_\{0\}are sampled uniformly and independently from

q\(0\)∈\[1\.2,2\.8\],T0∈\[0\.8,1\.2\],n1\(0\)∈\[0\.8,1\.3\],n2\(0\)∈\[1\.2,2\.4\],n3\(0\)∈\[0\.002,0\.003\]\.\\begin\{gathered\}q\(0\)\\in\[1\.2,2\.8\],\\qquad T\_\{0\}\\in\[0\.8,1\.2\],\\\\ n\_\{1\}\(0\)\\in\[0\.8,1\.3\],\\qquad n\_\{2\}\(0\)\\in\[1\.2,2\.4\],\\qquad n\_\{3\}\(0\)\\in\[0\.002,0\.003\]\.\\end\{gathered\}\(37\)The initial internal energy is set toϵ⁡\(0\)=32​n​\(0\)​T0\\epsilon\(0\)=\\tfrac\{3\}\{2\}n\(0\)T\_\{0\}, which givesϵ⁡\(0\)∈\[2\.4,6\.7\]\\epsilon\(0\)\\in\[2\.4,6\.7\]to the reported precision\. The generated dataset consists of 220 trajectories with 1000 state records per trajectory ont∈\[0,9\.99\]t\\in\[0,9\.99\]\. Among these generated trajectories, 140 trajectories are used for training, 50 for calibration, and 30 for testing\.

Model validation is performed by rolling out each trained model from the initial condition of every testing trajectory and comparing the predictions with the corresponding reference solution\. RK4 integrator is used with the same timestep employed to generate the dataset,Δ​t=0\.01\\Delta t=0\.01\. For each test trajectory, S\-PENNs and MC dropout useNs=2000N\_\{s\}=2000uncertainty realizations, whereas the deep\-ensemble method comprises 50 independently trained N\-GENNs with different architecture choices and random initializations\. Further details on the network architectures, optimization settings, and inference parameters are provided in Table[4](https://arxiv.org/html/2608.12624#A3.T4)\.

#### 3\.2\.2Results and discussion

In Fig\.[4](https://arxiv.org/html/2608.12624#S3.F4), we compare the three UQ methods on a representative testing trajectory\. The six columns correspond to the state variablesqq,pp,ϵ\\epsilon,n1n\_\{1\},n2n\_\{2\}, andn3n\_\{3\}, while the first three rows show deep ensembles, MC dropout, and S\-PENNs \(half\-normal\), respectively\. Deep ensembles produce the most accurate predictive mean for every component, with relativeℓ2\\ell^\{2\}error of7\.96×10−47\.96\\times 10^\{\-4\}forqq,1\.10×10−21\.10\\times 10^\{\-2\}forpp,6\.43×10−46\.43\\times 10^\{\-4\}forϵ\\epsilon,9\.33×10−49\.33\\times 10^\{\-4\}forn1n\_\{1\},2\.80×10−42\.80\\times 10^\{\-4\}forn2n\_\{2\}, and3\.62×10−43\.62\\times 10^\{\-4\}forn3n\_\{3\}\. S\-PENNs \(half\-normal\) gives corresponding errors of3\.40×10−33\.40\\times 10^\{\-3\},4\.55×10−24\.55\\times 10^\{\-2\},7\.64×10−37\.64\\times 10^\{\-3\},8\.29×10−38\.29\\times 10^\{\-3\},1\.45×10−31\.45\\times 10^\{\-3\}, and4\.23×10−34\.23\\times 10^\{\-3\}\. For both methods, the momentumpphas the largest error because small phase discrepancies are amplified as the trajectory repeatedly changes sign\. MC dropout shows larger phase and amplitude deviations, with errors of4\.86×10−24\.86\\times 10^\{\-2\}forqq,6\.84×10−16\.84\\times 10^\{\-1\}forpp,1\.19×10−21\.19\\times 10^\{\-2\}forϵ\\epsilon,3\.17×10−23\.17\\times 10^\{\-2\}forn1n\_\{1\},2\.39×10−32\.39\\times 10^\{\-3\}forn2n\_\{2\}, and1\.75×10−21\.75\\times 10^\{\-2\}forn3n\_\{3\}\. Additionally, the presented uncertainty bands provide a further distinction between the methods on this representative testing trajectory\. The bands obtained from deep ensembles remain narrow and successfully cover the errors\. MC dropout gives the widest bands, particularly forqqandpp\. However, this larger spread does not fully cover the discrepancies inn1n\_\{1\}andn3n\_\{3\}\. The uncalibrated S\-PENNs bands are comparable to those of deep ensembles\. The coverage curves in the bottom row of Fig\.[4](https://arxiv.org/html/2608.12624#S3.F4)confirm that these raw bands systematically under\-cover, indicating the model is overconfident\. Split conformal calibration successfully brings every component close to the target without altering the underlying thermodynamically admissible samples\.

![Refer to caption](https://arxiv.org/html/2608.12624v1/figs/results/chemicalMotorMain.png)Figure 4:Predictive trajectories and calibration diagnostics for the idealized chemical motor\. Columns correspond toqq,pp,ϵ\\epsilon,n1n\_\{1\},n2n\_\{2\}, andn3n\_\{3\}\. The first three rows show deep ensembles, MC dropout, and S\-PENNs \(half\-normal\), respectively, on a representative testing trajectory\. The predictive mean is plotted as a solid line, the simulated reference as a dashed line, and the uncalibratedμ^±2​σ^\\hat\{\\mu\}\\pm 2\\hat\{\\sigma\}interval as a blue band\. The S\-PENNs row also includes the calibrated 95% conformal interval in red\. Insets enlarge the late\-time response of the three species\. The final row compares empirical and target coverage for the uncalibrated and calibrated S\-PENNs intervals\. Titles report the componentwise relativeℓ2\\ell^\{2\}error \(RL2E\)\.![Refer to caption](https://arxiv.org/html/2608.12624v1/figs/results/chemicalMotorCombinedDiagnostics.png)Figure 5:Predictive accuracy and structure\-preservation diagnostics for the idealized chemical motor\. The top panel shows the componentwise trajectory\-wise relativeℓ2\\ell^\{2\}error \(RL2E\) over the testing dataset forqq,pp,ϵ\\epsilon,n1n\_\{1\},n2n\_\{2\}, andn3n\_\{3\}, comparing deep ensembles, MC dropout, and S\-PENNs using half\-normal, uniform, and exponential distributions for the epistemic index\. The box plots show the medians as solid lines, while the boxes span the interquartile ranges\. The bottom\-left panel shows the energy rateE˙\\dot\{E\}over time for a representative testing trajectory, with the median as a solid line and the 95% prediction interval as a shaded band across sampled S\-PENNs rollouts\. The bottom\-right panel shows the pointwise minimum entropy\-production ratemin⁡\(S˙\)\\min\(\\dot\{S\}\)on a logarithmic scale, where the minimum is taken across the sampled S\-PENNs rollouts at each timestep\.The comparison over the full test dataset is shown in Fig\.[5](https://arxiv.org/html/2608.12624#S3.F5)\. Deep ensembles have the lowest median relativeℓ2\\ell^\{2\}error for every component\. MC dropout constantly produces the largest errors, while all S\-PENNs variants have a lower median than MC dropout across all six variables\. The lower panels verify the thermodynamic behavior of the S\-PENNs samples\. For all three reference distributions, the energy\-rate median remains at zero and the 95% interval stays on the order of10−710^\{\-7\}\. The minimum sampled entropy\-production rate remains positive throughout the rollout\. Thus, as intended by the framework’s design, the uncertainty realizations retain the first\- and second\-law structure even in the presence of nonlinear chemical kinetics and state\-dependent reversible coupling\.

Table 2:Serial wall time and proper scoring rules for the idealized chemical motor example\. The wall time includes model training and generation of predictive samples for all testing trajectories\. ES denotes the energy score over the full multivariate trajectory, and CRPS is averaged for each state variable\. Scores are computed from uncalibrated predictive samples and averaged over the testing dataset\. Lower ES and CRPS values indicate better probabilistic predictions\.∗For deep ensembles, the reported wall time is the serial cost of the 50\-member ensemble, and the scoring rules are computed from the same ensemble\.

To complete the comparison, Table[2](https://arxiv.org/html/2608.12624#S3.T2)reports probabilistic scores and serial computational time\. Deep ensembles give the lowest ES,8\.2×10−18\.2\\times 10^\{\-1\}, and the lowest componentwise CRPS values\. Measured wall times are1\.5×1041\.5\\times 10^\{4\}s for the 50\-member deep ensemble,4\.4×1024\.4\\times 10^\{2\}s for MC dropout, and2\.7×1022\.7\\times 10^\{2\}–3\.2×1023\.2\\times 10^\{2\}s for S\-PENNs\. Thus, S\-PENNs are about two times less costly compared to deep ensembles and outperform MC dropout in probabilistic scores\. Scaling the deep ensemble to 2000 members gives5\.9×1055\.9\\times 10^\{5\}s, which is about three orders of magnitude higher than the cost needed for S\-PENNs\. These results extend the proposed construction to systems with noncanonical reversible coupling and a non\-quadratic reaction potential\.

### 3\.3Viscoplastic model

The third example is designed to assess S\-PENNs on a field\-valued problem characterized by a nonsmooth nonlinear dissipation law, moving beyond finite\-dimensional ODE dynamics, as considered in the first two examples\. More specifically, we consider a one\-dimensional continuum system governed by a viscoplastic constitutive relation of Perzyna type[46](https://arxiv.org/html/2608.12624#bib.bib78),[64](https://arxiv.org/html/2608.12624#bib.bib79), whose governing equations have been put in GENERIC form in[48](https://arxiv.org/html/2608.12624#bib.bib53)\. The state vector𝐱=\(u,p,εv​p,θ\)\\mathbf\{x\}=\(u,p,\\varepsilon^\{vp\},\\theta\), consists of the displacement fieldu⁡\(x,t\)u\(x,t\), momentum fieldp⁡\(x,t\)p\(x,t\), viscoplastic strain fieldεv​p​\(x,t\)\\varepsilon^\{vp\}\(x,t\), and temperature fieldθ⁡\(x,t\)\\theta\(x,t\)\. The total and elastic strains are

ε=∂u∂x,εe=ε−εv​p\.\\varepsilon=\\frac\{\\partial u\}\{\\partial x\},\\qquad\\varepsilon^\{e\}=\\varepsilon\-\\varepsilon^\{vp\}\.\(38\)Neglecting heat conduction, the energy and the entropy functionals are

E⁡\[𝐱\]=∫\(12​ρ​p2\+c​θ\+12​C​\(ε−εv​p\)2\)​𝑑x,S⁡\[𝐱\]=∫c​log⁡θ​𝑑x,E\[\\mathbf\{x\}\]=\\int\\left\(\\frac\{1\}\{2\\rho\}p^\{2\}\+c\\theta\+\\frac\{1\}\{2\}C\(\\varepsilon\-\\varepsilon^\{vp\}\)^\{2\}\\right\)\\,\\mathrm\{d\}x,\\qquad S\[\\mathbf\{x\}\]=\\int c\\log\\theta\\,\\mathrm\{d\}x,\(39\)whereρ\\rho,CC, andccdenote the mass density, elastic modulus, and specific heat, respectively\. The reversible part is governed by the canonical Poisson operator

L=\(0I00−I00000000000\),L=\\begin\{pmatrix\}0&I&0&0\\\\ \-I&0&0&0\\\\ 0&0&0&0\\\\ 0&0&0&0\\end\{pmatrix\},\(40\)whereas the irreversible response is governed by a non\-quadratic Perzyna\-type dissipation potential

Ξ⁡\[𝐱;𝐱∗\]=∫θ2​η​⟨\|ξεv​p\+ξθc​σ\|−σyθ⟩2​𝑑x,\\Xi\[\\mathbf\{x\};\\mathbf\{x\}^\{\*\}\]=\\int\\frac\{\\theta\}\{2\\eta\}\\left\\langle\\left\|\\xi\_\{\\varepsilon^\{vp\}\}\+\\frac\{\\xi\_\{\\theta\}\}\{c\}\\sigma\\right\|\-\\frac\{\\sigma\_\{y\}\}\{\\theta\}\\right\\rangle^\{2\}\\,\\mathrm\{d\}x,\(41\)where𝐱∗=\(ξu,ξp,ξεv​p,ξθ\)\\mathbf\{x\}^\{\*\}=\(\\xi\_\{u\},\\xi\_\{p\},\\xi\_\{\\varepsilon^\{vp\}\},\\xi\_\{\\theta\}\)denotes the dual variables,η\\etais the viscosity,σ=C⁡\(ε−εv​p\)\\sigma=C\(\\varepsilon\-\\varepsilon^\{vp\}\)is the stress,σy\>0\\sigma\_\{y\}\>0is the constant yield stress, and⟨a⟩=max⁡\(a,0\)\\langle a\\rangle=\\max\(a,0\)is the Macaulay bracket\. The potential is convex with respect to the dual variables and non\-quadratic and non\-smooth due to the overstress threshold\. Substitution into the GENERIC evolution equation gives

u˙=pρ,p˙=∂∂x​\(C⁡\(ε−εv​p\)\),ε˙v​p=1η​⟨\|σ\|−σy⟩​sign​\(σ\),θ˙=σ​ε˙v​pc\.\\dot\{u\}=\\frac\{p\}\{\\rho\},\\qquad\\dot\{p\}=\\frac\{\\partial\}\{\\partial x\}\\left\(C\(\\varepsilon\-\\varepsilon^\{vp\}\)\\right\),\\qquad\\dot\{\\varepsilon\}^\{vp\}=\\frac\{1\}\{\\eta\}\\langle\|\\sigma\|\-\\sigma\_\{y\}\\rangle\\,\\mathrm\{sign\}\(\\sigma\),\\qquad\\dot\{\\theta\}=\\frac\{\\sigma\\dot\{\\varepsilon\}^\{vp\}\}\{c\}\.\(42\)
#### 3\.3\.1Data generation and model validation

The dataset is generated by solving the evolution equations above on a bar of lengthl=0\.2l=0\.2m overt∈\[0,2\.0×10−4\]t\\in\[0,2\.0\\times 10^\{\-4\}\]s\. The material parameters are

ρ=7800​kg⋅m−3,C=210​GPa,η=8​GPa⋅s,σy=250​MPa,c=1000​J⋅kg−1⋅K−1\.\\rho=7800\\;\\mathrm\{kg\\cdot m^\{\-3\}\},\\quad C=210\\;\\mathrm\{GPa\},\\quad\\eta=8\\;\\mathrm\{GPa\\cdot s\},\\quad\\sigma\_\{y\}=250\\;\\mathrm\{MPa\},\\quad c=1000\\;\\mathrm\{J\\cdot kg^\{\-1\}\\cdot K^\{\-1\}\}\.\(43\)All simulations start fromu⁡\(x,0\)=0u\(x,0\)=0,p⁡\(x,0\)=0p\(x,0\)=0,εv​p​\(x,0\)=0\\varepsilon^\{vp\}\(x,0\)=0, andθ⁡\(x,0\)=298\\theta\(x,0\)=298K\. A uniform mesh withNx=29N\_\{x\}=29cells givesΔ​x=6\.90×10−3\\Delta x=6\.90\\times 10^\{\-3\}m, with displacement and momentum stored at nodes, and strain, viscoplastic strain, and temperature stored at cell centers\. A Courant–Friedrichs–Lewy \(CFL\) number of0\.50\.5based on the elastic wave speedC/ρ\\sqrt\{C/\\rho\}givesΔ​t=6\.65×10−7\\Delta t=6\.65\\times 10^\{\-7\}s\. The nodal mechanical variables are advanced with velocity–Verlet, and the cell\-centered viscoplastic strain and temperature are advanced with a forward Euler scheme, as in\[[69](https://arxiv.org/html/2608.12624#bib.bib31)\]\.

The left boundary is fixed,u⁡\(0,t\)=0u\(0,t\)=0, while the right boundary follows a smooth ramp to a target engineering strainεt\\varepsilon\_\{t\}and is then held fixed as

u⁡\(l,t\)=\{εt​l​\(3​\(ttramp\)2−2​\(ttramp\)3\),t≤tramp,εt​l,otherwise,u\(l,t\)=\\begin\{cases\}\\varepsilon\_\{t\}\\,l\\left\(3\\left\(\\frac\{t\}\{t\_\{\\mathrm\{ramp\}\}\}\\right\)^\{2\}\-2\\left\(\\frac\{t\}\{t\_\{\\mathrm\{ramp\}\}\}\\right\)^\{3\}\\right\),&t\\leq t\_\{\\mathrm\{ramp\}\},\\\\ \\varepsilon\_\{t\}\\,l,&\\text\{otherwise\},\\end\{cases\}\(44\)wheretramp=1\.6×10−4t\_\{\\mathrm\{ramp\}\}=1\.6\\times 10^\{\-4\}s\. The dataset contains150150loading cases with uniformly spaced target strainsεt∈\[1%,2%\]\\varepsilon\_\{t\}\\in\[1\\%,2\\%\]\. As in the harmonic oscillator example, the split uses7070loading cases for training,5050for conformal calibration, and3030for testing\.

The deterministic backbone follows the N\-GENNs parameterization reviewed in Section[2\.1](https://arxiv.org/html/2608.12624#S2.SS1)\. Because the state variables are fields, the networks learn the local energy, entropy, and dissipation densities rather than the corresponding functionals, which are obtained by integrating these densities over the domain\. The kinetic contribution to the energy density,p2/\(2​ρ\)p^\{2\}/\(2\\rho\), is imposed analytically, while the remaining energy density, and the total entropy and dissipation densities are learned from data\. The Poisson operator is fixed to the canonical form above, consistent with the GENERIC structure for generalized standard materials[48](https://arxiv.org/html/2608.12624#bib.bib53); hence, S\-PENNs perturb only the density networks in this example\. During rollout,u˙=p/ρ\\dot\{u\}=p/\\rhoand the prescribed boundary histories are imposed directly, while the learned model supplies the interior momentum rate and the cell\-centered rates ofεv​p\\varepsilon^\{vp\}andθ\\theta\.

For evaluation, the learned model is fed with the initial conditions of each testing case and rolled out on the same one\-dimensional mesh using an RK4 scheme with the same timestep used to generate the dataset,Δ​t=6\.65×10−7\\Delta t=6\.65\\times 10^\{\-7\}s\. The rollout state consists of the interior nodal values ofuuandpptogether with all cell\-centered values ofεv​p\\varepsilon^\{vp\}andθ\\theta\. Boundary histories foruuandppare imposed from the corresponding simulated case\. S\-PENNs and MC dropout useNs=2000N\_\{s\}=2000stochastic rollouts per testing case, while the deep\-ensemble baseline uses5050independently initialized N\-GENNs with the width choices reported in Table[4](https://arxiv.org/html/2608.12624#A3.T4)\. The full architecture, optimization, and inference settings are summarized in[C](https://arxiv.org/html/2608.12624#A3)\.

#### 3\.3\.2Results and discussion

We present S\-PENNs \(half\-normal\) on a representative testing case for this example in Fig\.[6](https://arxiv.org/html/2608.12624#S3.F6), with the corresponding deep ensemble and MC dropout baseline results provided in[B](https://arxiv.org/html/2608.12624#A2)\. In these field visualizations, columns correspond touu,pp,εv​p\\varepsilon^\{vp\}, andθ\\theta, and the rows show their predictive mean, simulated reference field, pointwise absolute error, and uncalibrated uncertainty\. Fig\.[6](https://arxiv.org/html/2608.12624#S3.F6)further reports the calibrated uncertainty and the empirical coverage against target coverage for S\-PENNs \(half\-normal\)\. For this testing case, S\-PENNs \(half\-normal\) yield relativeℓ2\\ell^\{2\}errors of1\.39×10−31\.39\\times 10^\{\-3\}foruu,1\.27×10−21\.27\\times 10^\{\-2\}forpp,5\.11×10−35\.11\\times 10^\{\-3\}forεv​p\\varepsilon^\{vp\}, and2\.80×10−42\.80\\times 10^\{\-4\}forθ\\theta\. The corresponding errors for deep ensembles are1\.85×10−31\.85\\times 10^\{\-3\},1\.82×10−21\.82\\times 10^\{\-2\},1\.72×10−31\.72\\times 10^\{\-3\}, and1\.30×10−41\.30\\times 10^\{\-4\}, respectively\. MC dropout produces the largest errors overall with5\.13×10−35\.13\\times 10^\{\-3\}foruu,5\.09×10−25\.09\\times 10^\{\-2\}forpp,9\.99×10−39\.99\\times 10^\{\-3\}forεv​p\\varepsilon^\{vp\}, and1\.63×10−31\.63\\times 10^\{\-3\}forθ\\theta\. Despite these quantitative differences, none of the methods qualitatively reproduces the spatial patterns in the pointwise absolute error fields\. The last row of Fig\.[6](https://arxiv.org/html/2608.12624#S3.F6)compares the empirical and target coverage of S\-PENNs \(half\-normal\)\. Before calibration, the uncertainty estimates under\-cover foruuandppbut over\-cover forεv​p\\varepsilon^\{vp\}andθ\\theta, indicating overconfidence in the former two variables and underconfidence in the latter two\. After split conformal calibration, the four coverage curves follow the target diagonal more closely\. After split conformal calibration, all four empirical coverage curves follow the target diagonal more closely\. The calibrated uncertainty fields also provide a closer representation of both the spatial error patterns and the range of error magnitudes for each state variable\. These findings show that conformal prediction calibration improves coverage reliability in this PDE example\.

![Refer to caption](https://arxiv.org/html/2608.12624v1/figs/results/viscoplasticMain.png)Figure 6:S\-PENNs \(half\-normal\) prediction, uncertainty, and calibration fields for the 1D viscoplastic model on a representative testing case\. Columns correspond touu,pp,εv​p\\varepsilon^\{vp\}, andθ\\theta\. From top to bottom, the rows show the predictive mean, the reference field, the pointwise absolute error, the uncalibrated uncertainty given by the pointwise empirical standard deviation across sampled rollouts, the calibrated uncertainty given by the conformal interval half\-widthq^1−α​σ^\\hat\{q\}\_\{1\-\\alpha\}\\hat\{\\sigma\}, and empirical coverage against target coverage\. Titles report the relativeℓ2\\ell^\{2\}error \(RL2E\) for each state variable\.![Refer to caption](https://arxiv.org/html/2608.12624v1/figs/results/viscoplasticRL2Eboxplot.png)Figure 7:Componentwise testing\-case relativeℓ2\\ell^\{2\}error \(RL2E\) over the space–time grid for the 1D viscoplastic testing dataset\. The grouped boxplots compare deep ensembles, MC dropout, and the half\-normal, uniform, and exponential S\-PENNs variants foruu,pp,εv​p\\varepsilon^\{vp\}, andθ\\theta\. Each box shows the median as a solid line and boxes span the interquartile ranges\.Fig\.[7](https://arxiv.org/html/2608.12624#S3.F7)moves from the representative case to the full testing dataset by plotting the statistics of the componentwise relativeℓ2\\ell^\{2\}error for each testing case\. MC dropout has the largest median RL2E for every field\. The deep ensemble has the lowest median errors forεv​p\\varepsilon^\{vp\}andθ\\theta, while S\-PENNs \(exponential\) have the lowest median errors foruuandpp\. Among the S\-PENNs variants, half\-normal gives the lowest median errors forεv​p\\varepsilon^\{vp\}andθ\\theta, whereas uniform and exponential reduce theuuandpperrors relative to half\-normal\. Thus, no S\-PENNs reference distribution dominates across all four variables, but all three variants have lower median errors than MC dropout for every component\.

To compare the quality of the quantified uncertainty and the computational cost, Table[3](https://arxiv.org/html/2608.12624#S3.T3)reports ES and componentwise CRPS from the uncalibrated samples together with serial wall time\. All three S\-PENNs variants achieve lower ES values than the deep ensemble, ranging from4\.8×1044\.8\\times 10^\{4\}to6\.6×1046\.6\\times 10^\{4\}, compared with1\.0×1051\.0\\times 10^\{5\}\. MC dropout performs worst, with an ES of3\.0×1053\.0\\times 10^\{5\}\. The componentwise CRPS results show that S\-PENNs \(exponential\) achieves the lowest values foruuandpp, whereas the deep ensemble performs best forεv​p\\varepsilon^\{vp\}andθ\\theta\. Each S\-PENNs variant also improves all four componentwise CRPS values relative to MC dropout while requiring less than half its wall time\. Specifically, the S\-PENNs runs take1\.0×1031\.0\\times 10^\{3\}s, compared with2\.3×1032\.3\\times 10^\{3\}s for MC dropout\. Measured wall times is7\.5×1047\.5\\times 10^\{4\}s for the 50\-member deep ensemble, which is about one order of magnitude higher than the cost for S\-PENNs\. Scaling the deep ensemble to 2000 members gives3\.0×1063\.0\\times 10^\{6\}s, or approximately3\.0×1033\.0\\times 10^\{3\}times the S\-PENNs cost\.

Table 3:Serial wall time and proper scoring rules for the 1D viscoplastic model\. The wall time includes model training and generation of predictive samples for all testing cases\. For deep ensembles, the reported time is extrapolated from the measured 50\-member ensemble as described in the table note\. ES denotes the energy score over the full\-state rollout, and CRPS is reported for each physical field\. Scores are computed from uncalibrated predictive samples and averaged over the testing dataset\. Lower ES and CRPS values indicate better probabilistic predictions\.∗For deep ensembles, the reported wall time is the serial cost of the 50\-member ensemble, and the scoring rules are computed from the same ensemble\.

Taken together, the 1D viscoplastic results extend the ODE findings from the previous examples to a PDE setting where uncertainty must be quantified over spatiotemporal fields rather than finite\-dimensional trajectories\. All three S\-PENNs variants achieve consistently strong predictive accuracy and comparable overall performance, while outperforming MC dropout in both predictive accuracy and probabilistic scores\. They also provide competitive UQ performance relative to the deep ensemble, with a much lower computational cost\.These results show that the advantages observed in the ODE example carry over to field\-valued PDE dynamics\.

## 4Conclusion

This work introduces Structure\-Preserving Epistemic Neural Networks \(S\-PENNs\) for uncertainty quantification in scientific machine learning models with architecturally enforced physical constraints\. In the GENERIC setting, S\-PENNs attach block\-specific epinets to the thermodynamic building components of N\-GENNs, which serve as base networks\. The epinets share a common epistemic index\. The resulting augmented blocks are then assembled through the same structure\-preserving reparameterizations as in N\-GENNs\. As a result, each sampled vector field remains thermodynamically consistent by construction, preserving total energy and nonnegative entropy production at the level of individual realizations\. Split conformal prediction is then used to calibrate interval widths obtained from the ensemble mean and standard deviation, providing finite\-sample marginal coverage under trajectory\-level exchangeability\. The numerical examples show that S\-PENNs provide a practical balance between structural fidelity, predictive quality, and computational cost across both finite\-dimensional and field\-valued dynamical systems\. S\-PENNs required typically1−31\-3orders of magnitude less serial wall time than deep ensembles, improved substantially over MC dropout in probabilistic scores, and showed only moderate sensitivity to the reference distribution for the epistemic index\.

The present study is limited in two main respects\. First, all three numerical examples use noise\-free state measurements\. With noisy states, the problem becomes an errors\-in\-variables \(EiV\) problem[10](https://arxiv.org/html/2608.12624#bib.bib6),[5](https://arxiv.org/html/2608.12624#bib.bib5),[65](https://arxiv.org/html/2608.12624#bib.bib1),[81](https://arxiv.org/html/2608.12624#bib.bib13): the clean trajectory should satisfy energy conservation and nonnegative entropy production, but the noise\-corrupted trajectory need not\. This differs from intrinsic physical stochasticity, for which individual realizations are governed by the appropriate balance and dissipation laws\. Training the same constrained dynamics directly on noisy measurements would therefore impose the thermodynamic constraints on a noise\-contaminated trajectory, which can bias the learned dynamics\. Recent work in scientific machine learning has begun to address such EiV problems, for example in[81](https://arxiv.org/html/2608.12624#bib.bib13)\. Second, epinets are designed as lightweight auxiliary networks for UQ[53](https://arxiv.org/html/2608.12624#bib.bib44)and have been shown to match large ensembles at orders\-of\-magnitude lower computational cost in benchmarks[52](https://arxiv.org/html/2608.12624#bib.bib45)\. However, the present study evaluates S\-PENNs only on two low\-dimensional ODE systems and a one\-dimensional PDE system, and the computational cost advantage of S\-PENNs relative to deep ensembles is expected to be further amplified in higher\-dimensional settings\.

Future work will first address noisy state measurements by separating the thermodynamically admissible trajectory from the measurement process\. One route is to denoise the measurements before learning the constrained dynamics, while another is to use latent\-variable formulations that represent noisy inputs before dynamics identification[8](https://arxiv.org/html/2608.12624#bib.bib7)\. A second direction is to evaluate scalability on higher\-dimensional thermomechanical systems, where the relative cost of epinet sampling, rollout generation, and deep ensembles may differ from the present benchmarks\.

More broadly, the S\-PENNs principle is not tied to GENERIC dynamics\. It is relevant to other structure\-preserving or physics\-constrained machine learning models in computational mechanics and may also be useful for pretrained foundation models in science and engineering whose architectures or adaptation procedures encode physical priors\. In these settings, uncertainty estimates should remain compatible with the encoded structure\. Because epinets can be attached to frozen or partially frozen building blocks, they provide a lightweight route to physically admissible uncertainty samples for active learning, material design, and reliability assessment\.

## Appendix AEvaluation metrics for the numerical examples

This appendix defines the evaluation metrics used in Section[3](https://arxiv.org/html/2608.12624#S3)\. To distinguish the testing data from the training trajectories in Section[2\.3\.3](https://arxiv.org/html/2608.12624#S2.SS3.SSS3), leti=1,…,Ntesti=1,\\dots,N\_\{\\mathrm\{test\}\}index a trajectory in𝒟test\\mathcal\{D\}\_\{\\mathrm\{test\}\}, let𝐗i\\mathbf\{X\}\_\{i\}denote its reference rollout, and let𝐗^i\(r\)\\hat\{\\mathbf\{X\}\}\_\{i\}^\{\(r\)\},r=1,…,Nsr=1,\\dots,N\_\{s\}, denote itsrr\-th predictive rollout sample\. The scalar indexjjidentifies a state variable at a discrete time point and, for a field\-valued variable, at a spatio\-temporal grid point\. For each reported variableν\\nu,𝒥ν\\mathcal\{J\}\_\{\\nu\}is the set of retained forecast entries used by the componentwise metric: forecast time points for an ODE variable and retained space–time grid points for a field\-valued variable\. The predictive sample mean at entryjjisμ^i​\(j\)=Ns−1​∑r=1Ns𝐗^i\(r\)​\(j\)\\hat\{\\mu\}\_\{i\}\(j\)=N\_\{s\}^\{\-1\}\\sum\_\{r=1\}^\{N\_\{s\}\}\\hat\{\\mathbf\{X\}\}\_\{i\}^\{\(r\)\}\(j\)\. All metrics use the forecast portion of each rollout and therefore exclude the prescribed initial state att=0t=0, and prescribed boundary conditions, if applicable\. Thus, the harmonic\-oscillator and chemical\-motor metrics uset=1,…,Ntt=1,\\dots,N\_\{t\}\. In the viscoplastic example, the imposed boundary conditions ofuuandppare also excluded\.

##### Relativeℓ2\\ell^\{2\}error

Predictive accuracy is measured by the relativeℓ2\\ell^\{2\}error \(RL2E\) of the predictive mean\. For testing trajectoryiiand state variableν\\nu,

RL2Ei\(ν\)=\(∑j∈𝒥ν\(μ^i​\(j\)−𝐗i​\(j\)\)2∑j∈𝒥ν\(𝐗i​\(j\)\)2\)1/2\.\\mathrm\{RL2E\}\_\{i\}^\{\(\\nu\)\}=\\left\(\\frac\{\\sum\_\{j\\in\\mathcal\{J\}\_\{\\nu\}\}\\left\(\\hat\{\\mu\}\_\{i\}\(j\)\-\\mathbf\{X\}\_\{i\}\(j\)\\right\)^\{2\}\}\{\\sum\_\{j\\in\\mathcal\{J\}\_\{\\nu\}\}\\left\(\\mathbf\{X\}\_\{i\}\(j\)\\right\)^\{2\}\}\\right\)^\{1/2\}\.\(45\)Figure titles reportRL2Ei\(ν\)\\mathrm\{RL2E\}\_\{i\}^\{\(\\nu\)\}for the displayed trajectory, and the box plots summarize its distribution over the testing dataset\. For field\-valued variables, the figures also show the pointwise absolute error\|μ^i​\(j\)−𝐗i​\(j\)\|\|\\hat\{\\mu\}\_\{i\}\(j\)\-\\mathbf\{X\}\_\{i\}\(j\)\|\.

##### Empirical coverage

For a nominal coverage level1−α1\-\\alpha, letC^1−α,i​\(j\)\\widehat\{C\}\_\{1\-\\alpha,i\}\(j\)denote the prediction interval at entryjjof testing trajectoryii\. For the uncalibrated rollout samples, its endpoints are the empiricalα/2\\alpha/2and1−α/21\-\\alpha/2quantiles; for calibrated results,C^1−α,i​\(j\)=C^1−α​\(𝐱i\(0\),j\)\\widehat\{C\}\_\{1\-\\alpha,i\}\(j\)=\\widehat\{C\}\_\{1\-\\alpha\}\(\\mathbf\{x\}\_\{i\}^\{\(0\)\};j\)is the split\-conformal interval in Eq\. \([24](https://arxiv.org/html/2608.12624#S2.E24)\)\. Empirical coverage is the fraction of reference rollout entries contained in these intervals[2](https://arxiv.org/html/2608.12624#bib.bib58):

EC\(ν\)\(1−α\)=1Ntest​\|𝒥ν\|∑i=1Ntest∑j∈𝒥ν\{𝐗i\(j\)∈C^1−α,i\(j\)\},\\mathrm\{EC\}^\{\(\\nu\)\}\(1\-\\alpha\)=\\frac\{1\}\{N\_\{\\mathrm\{test\}\}\|\\mathcal\{J\}\_\{\\nu\}\|\}\\sum\_\{i=1\}^\{N\_\{\\mathrm\{test\}\}\}\\sum\_\{j\\in\\mathcal\{J\}\_\{\\nu\}\}\\mathbf\{1\}\\\!\\left\\\{\\mathbf\{X\}\_\{i\}\(j\)\\in\\widehat\{C\}\_\{1\-\\alpha,i\}\(j\)\\right\\\},\(46\)where𝟏​\{⋅\}\\mathbf\{1\}\\\{\\cdot\\\}is the indicator function\. The set𝒥ν\\mathcal\{J\}\_\{\\nu\}contains only forecast entries, namely, it excludes the prescribed initial conditions and boundary conditions\. Ideal empirical calibration corresponds toEC\(ν\)​\(1−α\)=1−α\\mathrm\{EC\}^\{\(\\nu\)\}\(1\-\\alpha\)=1\-\\alpha\. Curves below this diagonal indicate under\-coverage associated with overconfident models, whereas curves above it indicate over\-coverage associated with underconfident models[15](https://arxiv.org/html/2608.12624#bib.bib9),[2](https://arxiv.org/html/2608.12624#bib.bib58),[71](https://arxiv.org/html/2608.12624#bib.bib64)\. We remark that finite calibration and testing datasets can produce deviations from the diagonal even when the conformal procedure satisfies its marginal coverage guarantee[2](https://arxiv.org/html/2608.12624#bib.bib58)\.

##### Proper scoring rules

Probabilistic forecast quality is evaluated on the uncalibrated rollout samples using the continuous ranked probability score \(CRPS\) and the energy score \(ES\)\. These proper scoring rules assess calibration and sharpness jointly, with sharpness considered subject to calibration[16](https://arxiv.org/html/2608.12624#bib.bib8),[15](https://arxiv.org/html/2608.12624#bib.bib9)\. Lower values indicate better probabilistic forecasts\. The CRPS is computed on scalar marginals and averaged separately for each state variable[47](https://arxiv.org/html/2608.12624#bib.bib10),[28](https://arxiv.org/html/2608.12624#bib.bib11):

CRPS\(ν\)=1Ntest​\|𝒥ν\|​∑i=1Ntest∑j∈𝒥ν\[1Ns​∑r=1Ns\|𝐗^i\(r\)​\(j\)−𝐗i​\(j\)\|−12​Ns2​∑r,r′=1Ns\|𝐗^i\(r\)​\(j\)−𝐗^i\(r′\)​\(j\)\|\]\.\\mathrm\{CRPS\}^\{\(\\nu\)\}=\\frac\{1\}\{N\_\{\\mathrm\{test\}\}\|\\mathcal\{J\}\_\{\\nu\}\|\}\\sum\_\{i=1\}^\{N\_\{\\mathrm\{test\}\}\}\\sum\_\{j\\in\\mathcal\{J\}\_\{\\nu\}\}\\left\[\\frac\{1\}\{N\_\{s\}\}\\sum\_\{r=1\}^\{N\_\{s\}\}\\left\|\\hat\{\\mathbf\{X\}\}\_\{i\}^\{\(r\)\}\(j\)\-\\mathbf\{X\}\_\{i\}\(j\)\\right\|\-\\frac\{1\}\{2N\_\{s\}^\{2\}\}\\sum\_\{r,r^\{\\prime\}=1\}^\{N\_\{s\}\}\\left\|\\hat\{\\mathbf\{X\}\}\_\{i\}^\{\(r\)\}\(j\)\-\\hat\{\\mathbf\{X\}\}\_\{i\}^\{\(r^\{\\prime\}\)\}\(j\)\\right\|\\right\]\.\(47\)The ES extends the same assessment to the joint predictive distribution of the full rollout[17](https://arxiv.org/html/2608.12624#bib.bib12)\. Viewing all forecast entries that remain after the exclusions above as vectors𝐗i\\mathbf\{X\}\_\{i\}and𝐗^i\(r\)\\hat\{\\mathbf\{X\}\}\_\{i\}^\{\(r\)\}, respectively, we estimate

ES=1Ntest​∑i=1Ntest\[1Ns​∑r=1Ns‖𝐗^i\(r\)−𝐗i‖2−12​Ns2​∑r,r′=1Ns‖𝐗^i\(r\)−𝐗^i\(r′\)‖2\]\.\\mathrm\{ES\}=\\frac\{1\}\{N\_\{\\mathrm\{test\}\}\}\\sum\_\{i=1\}^\{N\_\{\\mathrm\{test\}\}\}\\left\[\\frac\{1\}\{N\_\{s\}\}\\sum\_\{r=1\}^\{N\_\{s\}\}\\left\\\|\\hat\{\\mathbf\{X\}\}\_\{i\}^\{\(r\)\}\-\\mathbf\{X\}\_\{i\}\\right\\\|\_\{2\}\-\\frac\{1\}\{2N\_\{s\}^\{2\}\}\\sum\_\{r,r^\{\\prime\}=1\}^\{N\_\{s\}\}\\left\\\|\\hat\{\\mathbf\{X\}\}\_\{i\}^\{\(r\)\}\-\\hat\{\\mathbf\{X\}\}\_\{i\}^\{\(r^\{\\prime\}\)\}\\right\\\|\_\{2\}\\right\]\.\(48\)

## Appendix BAdditional baseline results

This appendix reports the representative field visualizations on the same testing case for the deep ensembles and MC dropout baselines referenced in Section[3\.3](https://arxiv.org/html/2608.12624#S3.SS3)\. Figures[8](https://arxiv.org/html/2608.12624#A2.F8)and[9](https://arxiv.org/html/2608.12624#A2.F9)use the same layout as the S\-PENNs results in Fig\.[6](https://arxiv.org/html/2608.12624#S3.F6)\. The four columns correspond touu,pp,εv​p\\varepsilon^\{vp\}, andθ\\theta, and the rows display, from top to bottom, the predictive mean, simulated reference field, pointwise absolute error, and uncalibrated uncertainty\. The uncertainty is the pointwise empirical standard deviation across the 50 independently trained models for the deep ensemble and across 2000 stochastic rollouts for MC dropout\.

![Refer to caption](https://arxiv.org/html/2608.12624v1/figs/appendix/viscoplasticEnsembles.png)Figure 8:Deep ensemble prediction and uncertainty fields for the 1D viscoplastic model on a representative testing instance\. Columns correspond touu,pp,εv​p\\varepsilon^\{vp\}, andθ\\theta\. From top to bottom, the rows show the predictive mean, the simulated solution, the pointwise absolute error, and the uncalibrated uncertainty given by the pointwise empirical standard deviation across ensemble members\. Titles report the relativeℓ2\\ell^\{2\}error \(RL2E\) for each state variable\.![Refer to caption](https://arxiv.org/html/2608.12624v1/figs/appendix/viscoplasticDropout.png)Figure 9:MC dropout prediction and uncertainty fields for the 1D viscoplastic model on a representative testing instance\. Columns correspond touu,pp,εv​p\\varepsilon^\{vp\}, andθ\\theta\. From top to bottom, the rows show the predictive mean, the simulated solution, the pointwise absolute error, and the uncalibrated uncertainty given by the pointwise empirical standard deviation across stochastic forward passes\. Titles report the relativeℓ2\\ell^\{2\}error \(RL2E\) for each state variable\.
## Appendix CNeural\-network setting, training, and inference details

All models are implemented in Python with JAX and Flax\. Training and inference are performed in single precision on a single NVIDIA RTX A6000 GPU\. Unless otherwise stated, the N\-GENNs backbone architectures and base\-network training settings used by each UQ method follow the settings for the corresponding numerical examples in the original N\-GENNs paper[69](https://arxiv.org/html/2608.12624#bib.bib31)\. Table[4](https://arxiv.org/html/2608.12624#A3.T4)summarizes the architecture, optimization, and inference settings used in the harmonic oscillator, idealized chemical motor, and 1D viscoplastic model\. The following paragraphs detail method\-specific implementation choices not fully captured by the table\.

##### S\-PENNs

S\-PENNs’ training follows the two\-stage procedure described in Section[2\.3\.3](https://arxiv.org/html/2608.12624#S2.SS3.SSS3): the deterministic N\-GENNs backbone is trained first, and then frozen while only the epinet parameters are optimized\. In the harmonic\-oscillator and chemical\-motor examples, the epinets augment all thermodynamic blocks\. In the 1D viscoplastic example, they augment only the non\-kinetic contribution to the energy density and the entropy and dissipation potential densities, while the canonical Poisson operator remains fixed\. The same architecture is used for all reference distributions considered in each example\. Each sampled epistemic index produces one rollout sample, and the resulting samples are used to compute the predictive mean, empirical standard deviation, and proper scoring rules\.

##### Deep ensembles

For deep ensembles, the hidden width cycles through the values in Table[4](https://arxiv.org/html/2608.12624#A3.T4)and remains fixed across the hidden layers of each member\. Every member uses the full training dataset, with diversity introduced through the random seed and hidden width\. The wall times in Tables[1](https://arxiv.org/html/2608.12624#S3.T1),[2](https://arxiv.org/html/2608.12624#S3.T2), and[3](https://arxiv.org/html/2608.12624#S3.T3)report the measured computational cost without parallelization\. The predictive statistics and proper scoring rules are computed from the 50 trained members\.

##### MC dropout

MC dropout uses the same N\-GENNs backbone as the deterministic model, with dropout active during training and inference\. A dropout realization is held fixed over all timesteps of a rollout, so each sample follows one sampled vector field\. We note that the 1D viscoplastic example uses a lower dropout rate of 0\.1, rather than the rate of 0\.2 used in the two ODE examples, because the higher rate produced non\-convergent training runs and unstable rollout predictions\.

Table 4:Architectural and training settings for the harmonic oscillator, idealized chemical motor, and 1D viscoplastic examples\. Each method block includes the optimizer, learning\-rate settings, network architecture, and method\-specific training and inference settings\. All neural networks are trained with Softplus activation functions\.MethodSettingHarmonic oscillatorChemical motor1D viscoplastic modelS\-PENNsOptimizerSOAPSOAPSOAPLearning rate \(base networks\)3×10−33\\times 10^\{\-3\}3×10−33\\times 10^\{\-3\}3×10−33\\times 10^\{\-3\}Learning rate \(epinets\)3×10−33\\times 10^\{\-3\}3×10−33\\times 10^\{\-3\}3×10−33\\times 10^\{\-3\}Num\. of hidden layers per block \(base networks\)222Num\. of hidden layers per block \(epinets\)222Hidden\-layer width \(base networks\)303030Hidden\-layer width \(epinets\)101010Training epochs \(base networks\)100001000010000Training epochs \(epinets\)20025001000Epistemic\-index dimension555Prior scale0\.10\.10\.1Deep ensemblesOptimizerSOAPSOAPSOAPLearning rate3×10−33\\times 10^\{\-3\}3×10−33\\times 10^\{\-3\}3×10−33\\times 10^\{\-3\}Ensemble members505050Num\. of hidden layers per block222Hidden\-layer width choices\(20, 30, 40\)\(20, 30, 40\)\(20, 25, 30\)Training epochs per member100001000010000MC dropoutOptimizerSOAPSOAPSOAPLearning rate3×10−33\\times 10^\{\-3\}3×10−33\\times 10^\{\-3\}3×10−33\\times 10^\{\-3\}Dropout rate0\.20\.20\.1Num\. of hidden layers per block222Hidden\-layer width303030Training epochs100001000010000

## Appendix DEpinet construction for inverse problems

The input\-independent epinet branch introduced in Section[2\.3\.1](https://arxiv.org/html/2608.12624#S2.SS3.SSS1)can also be used for inverse problems with unknown global physical parameters\. These benchmarks are reported in the appendix because the main S\-PENNs construction is developed for GENERIC dynamics, whereas the inverse problems isolate the same global\-parameter mechanism in a PINN setting\.

Let𝐬\\mathbf\{s\}denote the independent variables, with𝐬=t\\mathbf\{s\}=tfor an ODE trajectory and𝐬=\(x,t\)\\mathbf\{s\}=\(x,t\)for a spatiotemporal PDE field, and letΩ\\Omegabe the domain on which the governing residual is evaluated\. We consider inverse problems in which the stateu⁡\(𝐬\)u\(\\mathbf\{s\}\)and the unknown global physical parameters𝜷∈ℝdβ\\bm\{\\beta\}\\in\\mathbb\{R\}^\{d\_\{\\beta\}\}are inferred from measurements while violations of the governing equations are penalized through the residual[59](https://arxiv.org/html/2608.12624#bib.bib56),

ℛ⁡\(u,𝜷\)​\(𝐬\)=0,𝐬∈Ω\.\\mathcal\{R\}\(u;\\bm\{\\beta\}\)\(\\mathbf\{s\}\)=0,\\qquad\\mathbf\{s\}\\in\\Omega\.\(49\)Here,ℛ\\mathcal\{R\}denotes the governing residual\. In the deterministic PINN baseline, the state is parameterized by a neural networkμ𝝍u​\(𝐬\)\\mu\_\{\\bm\{\\psi\}\}^\{u\}\(\\mathbf\{s\}\), and the global parameters are represented by the trainable vector𝜷𝝍\\bm\{\\beta\}\_\{\\bm\{\\psi\}\}; both are included in the deterministic parameter set𝝍\\bm\{\\psi\}\. The corresponding loss combines data\-misfit and residual terms, with additional terms included when initial or boundary conditions are applicable\.

Following the ENN notation in Eq\. \([8](https://arxiv.org/html/2608.12624#S2.E8)\), the deterministic PINN is first pretrained and then held fixed while the epinet parametersϕ\\bm\{\\phi\}are optimized\. The state\-branch correctionσϕu\\sigma\_\{\\bm\{\\phi\}\}^\{u\}takes as input the stop\-gradient feature vector𝒉¯𝝍u​\(𝐬\)=\[sg⁡\(𝒉𝝍u​\(𝐬\)\),𝐬\]\\bar\{\\bm\{h\}\}\_\{\\bm\{\\psi\}\}^\{u\}\(\\mathbf\{s\}\)=\[\\mathrm\{sg\}\(\\bm\{h\}\_\{\\bm\{\\psi\}\}^\{u\}\(\\mathbf\{s\}\)\),\\mathbf\{s\}\], which concatenates the frozen base features with the input coordinates\. The parameter\-branch correction𝝈ϕβ\\bm\{\\sigma\}\_\{\\bm\{\\phi\}\}^\{\\beta\}is deliberately input independent and depends only on the shared epistemic index𝐳\\mathbf\{z\}\. Each realization therefore assigns one global perturbation of the parameter vector overΩ\\Omega, rather than a coordinate\-dependent parameter field\. The augmented model is

uϑ​\(𝐬,𝐳\)=μ𝝍u​\(𝐬\)\+σϕu​\(𝒉¯𝝍u​\(𝐬\),𝐳\),𝜷ϑ​\(𝐳\)=𝜷𝝍\+𝝈ϕβ​\(𝐳\),u\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{s\},\\mathbf\{z\}\)=\\mu\_\{\\bm\{\\psi\}\}^\{u\}\(\\mathbf\{s\}\)\+\\sigma\_\{\\bm\{\\phi\}\}^\{u\}\(\\bar\{\\bm\{h\}\}\_\{\\bm\{\\psi\}\}^\{u\}\(\\mathbf\{s\}\),\\mathbf\{z\}\),\\qquad\\bm\{\\beta\}\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{z\}\)=\\bm\{\\beta\}\_\{\\bm\{\\psi\}\}\+\\bm\{\\sigma\}\_\{\\bm\{\\phi\}\}^\{\\beta\}\(\\mathbf\{z\}\),\(50\)whereϑ=\(𝝍,ϕ\)\\bm\{\\vartheta\}=\(\\bm\{\\psi\},\\bm\{\\phi\}\)\. The state\-branch correction follows the input\-dependent epinet construction in Eqs\. \([9](https://arxiv.org/html/2608.12624#S2.E9)\)–\([11](https://arxiv.org/html/2608.12624#S2.E11)\) and the parameter\-branch correction follows the same linear input\-independent form as in Eq\. \([14](https://arxiv.org/html/2608.12624#S2.E14)\),

𝝈ϕβ​\(𝐳\)\\displaystyle\\bm\{\\sigma\}\_\{\\bm\{\\phi\}\}^\{\\beta\}\(\\mathbf\{z\}\)=𝝈ϕβ,learn​\(𝐳\)\+w​𝝈β,prior​\(𝐳\),\\displaystyle=\\bm\{\\sigma\}\_\{\\bm\{\\phi\}\}^\{\\beta,\\mathrm\{learn\}\}\(\\mathbf\{z\}\)\+w\\,\\bm\{\\sigma\}^\{\\beta,\\mathrm\{prior\}\}\(\\mathbf\{z\}\),\(51\)𝝈ϕβ,learn​\(𝐳\)\\displaystyle\\bm\{\\sigma\}\_\{\\bm\{\\phi\}\}^\{\\beta,\\mathrm\{learn\}\}\(\\mathbf\{z\}\)=∑n=1dzznϕnβ,𝝈β,prior\(𝐳\)=∑n=1dzzn𝜻nβ\.\\displaystyle=\\sum\_\{n=1\}^\{d\_\{z\}\}z\_\{n\}\\,\\bm\{\\phi\}\_\{n\}^\{\\beta\},\\qquad\\bm\{\\sigma\}^\{\\beta,\\mathrm\{prior\}\}\(\\mathbf\{z\}\)=\\sum\_\{n=1\}^\{d\_\{z\}\}z\_\{n\}\\,\\bm\{\\zeta\}\_\{n\}^\{\\beta\}\.Here,ϕnβ∈ℝdβ\\bm\{\\phi\}\_\{n\}^\{\\beta\}\\in\\mathbb\{R\}^\{d\_\{\\beta\}\}are learnable coefficient vectors and𝜻nβ∈ℝdβ\\bm\{\\zeta\}\_\{n\}^\{\\beta\}\\in\\mathbb\{R\}^\{d\_\{\\beta\}\}are fixed prior coefficient vectors\. For each sampled𝐳\\mathbf\{z\}, the residualℛ\\mathcal\{R\}is evaluated using the corresponding state realizationuϑ​\(⋅,𝐳\)u\_\{\\bm\{\\vartheta\}\}\(\\cdot,\\mathbf\{z\}\)and global parameter vector𝜷ϑ​\(𝐳\)\\bm\{\\beta\}\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{z\}\)\. The shared epistemic index couples state and parameter uncertainty, while𝜷ϑ​\(𝐳\)\\bm\{\\beta\}\_\{\\bm\{\\vartheta\}\}\(\\mathbf\{z\}\)remains a single global parameter vector overΩ\\Omega\.

We evaluate this construction on two benchmark inverse problems from[82](https://arxiv.org/html/2608.12624#bib.bib66), using the datasets distributed with the open\-source NeuralUQ repository \([https://github\.com/Crunch\-UQ4MI/neuraluq](https://github.com/Crunch-UQ4MI/neuraluq)\)\. The B\-PINNs\-HMC baseline[72](https://arxiv.org/html/2608.12624#bib.bib14)is reproduced from the same NeuralUQ implementation used in[82](https://arxiv.org/html/2608.12624#bib.bib66)\. The proposed construction uses the two\-branch epinet design described above, with the state\-branch associated withuuand the parameter\-branch associated with𝜷\\bm\{\\beta\}\. Table[5](https://arxiv.org/html/2608.12624#A4.T5)summarizes the architecture, training, and inference settings for the two benchmarks shown in the following sections\. At inference, uncertainty is summarized by the empirical predictive mean and standard deviation with the same number of samples generated in the B\-PINNs\-HMC baselines for fair comparison\.

Table 5:Architecture, training, and inference settings for the proposed two\-branch epinet construction in the Kraichnan–Orszag \(KO\) and Korteweg–de Vries \(KdV\) inverse problems\. Tanh is used as the activation function and AdamW is the optimizer used for both training stages\.### D\.1Inverse Kraichnan–Orszag system

As a first benchmark, we consider the inverse Kraichnan–Orszag system, a three\-variable nonlinear dynamical model arising from interactions among inviscid shear waves[70](https://arxiv.org/html/2608.12624#bib.bib72)\. The governing dynamics and initial conditions considered are

x˙1\\displaystyle\\dot\{x\}\_\{1\}=a​x2​x3,\\displaystyle=ax\_\{2\}x\_\{3\},\(52\)x˙2\\displaystyle\\dot\{x\}\_\{2\}=b​x1​x3,\\displaystyle=bx\_\{1\}x\_\{3\},x˙3\\displaystyle\\dot\{x\}\_\{3\}=−\(a\+b\)​x1​x2,\\displaystyle=\-\(a\+b\)x\_\{1\}x\_\{2\},\(x1​\(0\),x2​\(0\),x3​\(0\)\)=\(1\.0,0\.8,0\.5\)\.\\bigl\(x\_\{1\}\(0\),x\_\{2\}\(0\),x\_\{3\}\(0\)\\bigr\)=\(1\.0,0\.8,0\.5\)\.\(53\)The reference coefficients area=1a=1andb=1b=1\. In the inverse setting, both coefficients and the initial conditions are treated as unknowns\. The task is to inferaaandbbtogether with the full trajectoriesx1x\_\{1\},x2x\_\{2\}, andx3x\_\{3\}overt∈\[0,10\]t\\in\[0,10\]from sparse noisy measurements\. Following[82](https://arxiv.org/html/2608.12624#bib.bib66), we use 11 measurements forx1x\_\{1\}andx3x\_\{3\}and 7 measurements forx2x\_\{2\}, each corrupted by zero\-mean Gaussian noise with known standard deviation0\.050\.05\.

Fig\.[10](https://arxiv.org/html/2608.12624#A4.F10)compares the reconstructed trajectories\. Both the reproduced B\-PINNs\-HMC baseline and the proposed two\-branch epinet construction recover the oscillatory dynamics from sparse noisy measurements, and their predictive means remain close to the reference trajectories over the full time interval\. The componentwise relativeℓ2\\ell^\{2\}errors show comparable state\-reconstruction performance\. The proposed construction reduces the errors forx1x\_\{1\}andx3x\_\{3\}, from7\.29×10−27\.29\\times 10^\{\-2\}to3\.93×10−23\.93\\times 10^\{\-2\}and from8\.10×10−28\.10\\times 10^\{\-2\}to6\.69×10−26\.69\\times 10^\{\-2\}, respectively, while B\-PINNs\-HMC is more accurate forx2x\_\{2\}, with an error of6\.30×10−26\.30\\times 10^\{\-2\}compared to7\.39×10−27\.39\\times 10^\{\-2\}\. The estimated parameter distributions in Fig\.[11](https://arxiv.org/html/2608.12624#A4.F11)show that B\-PINNs\-HMC gives sample means and standard deviations ofa=0\.86±0\.06a=0\.86\\pm 0\.06andb=1\.06±0\.03b=1\.06\\pm 0\.03, whereas the proposed two\-branch epinet construction gives more accurate results witha=0\.98±0\.07a=0\.98\\pm 0\.07andb=1\.01±0\.06b=1\.01\\pm 0\.06\.

![Refer to caption](https://arxiv.org/html/2608.12624v1/figs/appendix/inverse_problems/ko_results_comparison.png)Figure 10:Predictive trajectories for the inverse Kraichnan–Orszag problem\. Columns correspond tox1x\_\{1\},x2x\_\{2\}, andx3x\_\{3\}; rows compare B\-PINNs\-HMC and the proposed two\-branch epinet construction\. In each panel, the predictive mean is shown as a solid line, the reference solution as a dashed line, training measurements as orange markers, and the±2\\pm 2standard deviation band as a blue shaded region\. Titles report the relativeℓ2\\ell^\{2\}error \(RL2E\)\.![Refer to caption](https://arxiv.org/html/2608.12624v1/figs/appendix/inverse_problems/ko_parameter_comparison.png)Figure 11:Estimated distributions of the unknown parameters for the inverse Kraichnan–Orszag problem\. Columns correspond toaaandbb; rows compare B\-PINNs\-HMC and the proposed two\-branch epinet construction\. Each panel shows the sample histogram, the sample mean as a solid vertical line, and the reference value as a dashed vertical line\. The sample meanμ\\muand standard deviationσ\\sigmaare annotated above each panel\.
### D\.2Inverse Korteweg–de Vries problem

The second benchmark considers an inverse Korteweg–de Vries \(KdV\) problem\. Compared with the Kraichnan–Orszag system, this example tests the proposed construction on a field\-valued inverse problem with localized wave profiles and soliton interaction\. The KdV equation is a classical dispersive model for nonlinear wave propagation[4](https://arxiv.org/html/2608.12624#bib.bib71), and we write it as

ut=κ1​u​ux\+κ2​ux​x​x\+f,u\_\{t\}=\\kappa\_\{1\}uu\_\{x\}\+\\kappa\_\{2\}u\_\{xxx\}\+f,\(54\)where the reference problem usesκ1=1\.5\\kappa\_\{1\}=1\.5,κ2=0\.25\\kappa\_\{2\}=0\.25andf=0f=0\. Following[4](https://arxiv.org/html/2608.12624#bib.bib71), the exact two\-soliton solution is

u⁡\(x,t\)=2​∂x2log⁡\[exp⁡\(−ω1−ω2\)\+exp⁡\(ω1−ω2\)\+exp⁡\(ω2−ω1\)\+\(a1−a2\)2\(a1\+a2\)2​exp⁡\(ω1\+ω2\)\],u\(x,t\)=2\\partial\_\{x\}^\{2\}\\log\\left\[\\exp\(\-\\omega\_\{1\}\-\\omega\_\{2\}\)\+\\exp\(\\omega\_\{1\}\-\\omega\_\{2\}\)\+\\exp\(\\omega\_\{2\}\-\\omega\_\{1\}\)\+\\frac\{\(a\_\{1\}\-a\_\{2\}\)^\{2\}\}\{\(a\_\{1\}\+a\_\{2\}\)^\{2\}\}\\exp\(\\omega\_\{1\}\+\\omega\_\{2\}\)\\right\],\(55\)where the phase variables are

ωi=aix\+ai3χt\+bi,i=1,2,\\displaystyle\\omega\_\{i\}=a\_\{i\}x\+a\_\{i\}^\{3\}\\chi t\+b\_\{i\},\\qquad i=1,2,\(56\)witha1=1a\_\{1\}=1,a2=2a\_\{2\}=2,b1=b2=log⁡\(3\)/2b\_\{1\}=b\_\{2\}=\\log\(3\)/2, andχ=1\\chi=1\.

The inverse task is to inferκ1\\kappa\_\{1\}andκ2\\kappa\_\{2\}while reconstructingu⁡\(x,t\)u\(x,t\)over the spatiotemporal domain from sparse noisy measurements\. We again use the data from[82](https://arxiv.org/html/2608.12624#bib.bib66), with200200random measurements ofuuand100100measurements offf\. Both measurement types are perturbed by zero\-mean Gaussian noise with known standard deviation0\.10\.1\. Fig\.[12](https://arxiv.org/html/2608.12624#A4.F12)reports the predicted profiles at two representative time slices\. Both B\-PINNs\-HMC and the proposed two\-branch epinet construction recover the dominant two\-soliton structure, including the narrow high\-amplitude peak and the smaller secondary wave\. Their errors differ by regime\. Att=−1\.5t=\-1\.5, where the soliton peaks remain well separated, the proposed construction reduces the relativeℓ2\\ell^\{2\}error from2\.41×10−12\.41\\times 10^\{\-1\}to1\.23×10−11\.23\\times 10^\{\-1\}\. Att=0\.4t=0\.4, closer to the interaction regime, B\-PINNs\-HMC gives the smaller error,4\.82×10−24\.82\\times 10^\{\-2\}compared with5\.37×10−25\.37\\times 10^\{\-2\}for the proposed two\-branch epinet construction\. The parameter distributions in Fig\.[13](https://arxiv.org/html/2608.12624#A4.F13)show that B\-PINNs\-HMC givesκ1=1\.39±0\.03\\kappa\_\{1\}=1\.39\\pm 0\.03andκ2=0\.26±0\.01\\kappa\_\{2\}=0\.26\\pm 0\.01, whereas the proposed two\-branch epinet construction givesκ1=1\.50±0\.05\\kappa\_\{1\}=1\.50\\pm 0\.05andκ2=0\.37±0\.02\\kappa\_\{2\}=0\.37\\pm 0\.02\. Thus, the input\-independent parameter branch identifiesκ1\\kappa\_\{1\}nearly exactly, but it overestimatesκ2\\kappa\_\{2\}\. In contrast, B\-PINNs\-HMC is closer forκ2\\kappa\_\{2\}but biased low forκ1\\kappa\_\{1\}\.

![Refer to caption](https://arxiv.org/html/2608.12624v1/figs/appendix/inverse_problems/kdv_results_comparison.png)Figure 12:Predictive solution profiles for the inverse KdV problem at two representative time slices\. Columns correspond tot=−1\.5t=\-1\.5andt=0\.4t=0\.4; rows compare B\-PINNs\-HMC and the proposed two\-branch epinet construction\. In each panel, the predictive mean is shown as a solid line, the reference solution as a dashed line, training measurements as orange markers, and the±2\\pm 2standard deviation band as a blue shaded region\. Titles report the relativeℓ2\\ell^\{2\}error \(RL2E\)\. The scaling forttfollows[4](https://arxiv.org/html/2608.12624#bib.bib71)\.![Refer to caption](https://arxiv.org/html/2608.12624v1/figs/appendix/inverse_problems/kdv_parameter_comparison.png)Figure 13:Estimated distributions of the unknown parameters for the inverse KdV problem\. Columns correspond toκ1\\kappa\_\{1\}andκ2\\kappa\_\{2\}; rows compare B\-PINNs\-HMC and the proposed two\-branch epinet construction\. Each panel shows the sample histogram, the sample mean as a solid vertical line, and the reference value as a dashed vertical line\. The sample meanμ\\muand standard deviationσ\\sigmaare annotated above each panel\.

## Code availability

The code will be made available upon publication\.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper\.

## Acknowledgments

The authors acknowledge support from the US Department of the Army W911NF2310230\.

## Declaration of generative AI and AI\-assisted technologies in the writing process

During the preparation of this work, the authors used ChatGPT to check grammar and improve sentence clarity\. After using this tool, the authors reviewed and edited the content as needed and took full responsibility for the final published article\.

## References

- B\. Amos, L\. Xu, and J\. Z\. KolterInput convex neural networks\.InInternational conference on machine learning,pp\. 146–155\.Cited by:[§2\.1](https://arxiv.org/html/2608.12624#S2.SS1.p3.2)\.
- Angelopoulos and Bates \(2023\)A\. N\. Angelopoulos and S\. BatesConformal prediction: a gentle introduction\.Foundations and Trends in Machine Learning16\(4\),pp\. 494–591\.Cited by:[Appendix A](https://arxiv.org/html/2608.12624#A1.SS0.SSS0.Px2.p1.1),[Appendix A](https://arxiv.org/html/2608.12624#A1.SS0.SSS0.Px2.p1.2),[§1](https://arxiv.org/html/2608.12624#S1.p6.1),[§2\.4](https://arxiv.org/html/2608.12624#S2.SS4.p3.1),[§2\.4](https://arxiv.org/html/2608.12624#S2.SS4.p4.2)\.
- Bahmani \(2026\)B\. BahmaniConformal quantile regression for neural probabilistic constitutive modeling\.Computer Methods in Applied Mechanics and Engineering457,pp\. 118981\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p4.1)\.
- Beneset al\.\(2006\)N\. Benes, A\. Kasman, and K\. YoungOn decompositions of the KdV 2\-soliton\.Journal of Nonlinear Science16\(2\),pp\. 179–200\.Cited by:[Figure 12](https://arxiv.org/html/2608.12624#A4.F12),[§D\.2](https://arxiv.org/html/2608.12624#A4.SS2.p1.1),[§D\.2](https://arxiv.org/html/2608.12624#A4.SS2.p1.2)\.
- Carrollet al\.\(2006\)R\. J\. Carroll, D\. Ruppert, L\. A\. Stefanski, and C\. M\. CrainiceanuMeasurement error in nonlinear models: a modern perspective\.Chapman and Hall/CRC\.Cited by:[§4](https://arxiv.org/html/2608.12624#S4.p2.1)\.
- Celledoniet al\.\(2021\)E\. Celledoni, M\. J\. Ehrhardt, C\. Etmann, R\. I\. McLachlan, B\. Owren, C\. Schonlieb, and F\. SherryStructure\-preserving deep learning\.European journal of applied mathematics32\(5\),pp\. 888–936\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p1.1)\.
- Chenet al\.\(2018\)R\. T\. Chen, Y\. Rubanova, J\. Bettencourt, and D\. K\. DuvenaudNeural ordinary differential equations\.Advances in neural information processing systems31\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p1.1)\.
- Contiet al\.\(2026\)P\. Conti, J\. Kneifl, A\. Manzoni, A\. Frangi, J\. Fehr, S\. L\. Brunton, and J\. N\. KutzVENI, VINDy, VICI: a generative reduced\-order modeling framework with uncertainty quantification\.Neural Networks,pp\. 108543\.Cited by:[§4](https://arxiv.org/html/2608.12624#S4.p3.1)\.
- Cranmeret al\.\(2019\)M\. Cranmer, S\. Greydanus, S\. Hoyer, P\. Battaglia, D\. Spergel, and S\. HoLagrangian Neural Networks\.InICLR 2020 Workshop on Integration of Deep Neural Models and Differential Equations,External Links:[Link](https://openreview.net/forum?id=iE8tFa4Nq)Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p2.1)\.
- Fuller \(2009\)W\. A\. FullerMeasurement error models\.John Wiley & Sons\.Cited by:[§4](https://arxiv.org/html/2608.12624#S4.p2.1)\.
- Gal and Ghahramani \(2016\)Y\. Gal and Z\. GhahramaniDropout as a Bayesian approximation: representing model uncertainty in deep learning\.Ininternational conference on machine learning,pp\. 1050–1059\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p3.1),[§1](https://arxiv.org/html/2608.12624#S1.p5.1)\.
- Galiotoet al\.\(2024\)N\. Galioto, H\. Sharma, B\. Kramer, and A\. A\. GorodetskyBayesian identification of nonseparable Hamiltonians with multiplicative noise using deep learning and reduced\-order modeling\.Computer Methods in Applied Mechanics and Engineering430,pp\. 117194\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p4.1)\.
- Gao and Ng \(2022\)Y\. Gao and M\. K\. NgWasserstein generative adversarial uncertainty quantification in physics\-informed neural networks\.Journal of Computational Physics463,pp\. 111270\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p3.1)\.
- Gawlikowskiet al\.\(2023\)J\. Gawlikowski, C\. R\. N\. Tassi, M\. Ali, J\. Lee, M\. Humt, J\. Feng, A\. Kruspe, R\. Triebel, P\. Jung, R\. Roscher,et al\.A survey of uncertainty in deep neural networks\.Artificial intelligence review56\(Suppl 1\),pp\. 1513–1589\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p5.1)\.
- Gneitinget al\.\(2007\)T\. Gneiting, F\. Balabdaoui, and A\. E\. RafteryProbabilistic forecasts, calibration and sharpness\.Journal of the Royal Statistical Society Series B: Statistical Methodology69\(2\),pp\. 243–268\.Cited by:[Appendix A](https://arxiv.org/html/2608.12624#A1.SS0.SSS0.Px2.p1.2),[Appendix A](https://arxiv.org/html/2608.12624#A1.SS0.SSS0.Px3.p1.1)\.
- Gneiting and Raftery \(2007\)T\. Gneiting and A\. E\. RafteryStrictly proper scoring rules, prediction, and estimation\.Journal of the American statistical Association102\(477\),pp\. 359–378\.Cited by:[Appendix A](https://arxiv.org/html/2608.12624#A1.SS0.SSS0.Px3.p1.1),[§3](https://arxiv.org/html/2608.12624#S3.p1.1)\.
- Gneitinget al\.\(2008\)T\. Gneiting, L\. I\. Stanberry, E\. P\. Grimit, L\. Held, and N\. A\. JohnsonAssessing probabilistic forecasts of multivariate quantities, with an application to ensemble predictions of surface winds\.Test17\(2\),pp\. 211–235\.Cited by:[Appendix A](https://arxiv.org/html/2608.12624#A1.SS0.SSS0.Px3.p1.2),[§3](https://arxiv.org/html/2608.12624#S3.p1.1)\.
- Goan and Fookes \(2020\)E\. Goan and C\. FookesBayesian Neural Networks: an introduction and survey\.InCase Studies in Applied Bayesian Data Science: CIRM Jean\-Morlet Chair, Fall 2018,pp\. 45–87\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p3.1)\.
- Gopakumaret al\.\(2026\)V\. Gopakumar, A\. Gray, J\. Oskarsson, L\. Zanisi, D\. Giles, M\. J\. Kusner, S\. Pamela, and M\. Peter DeisenrothUncertainty quantification of surrogate models using conformal prediction\.Machine Learning: Science and Technology7\(1\),pp\. 015025\.Cited by:[§2\.4](https://arxiv.org/html/2608.12624#S2.SS4.p2.1)\.
- Greydanuset al\.\(2019\)S\. Greydanus, M\. Dzamba, and J\. YosinskiHamiltonian Neural Networks\.Advances in neural information processing systems32\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p2.1)\.
- Grmela and Öttinger \(1997\)M\. Grmela and H\. C\. ÖttingerDynamics and thermodynamics of complex fluids\. I\. development of a general formalism\.Physical Review E56\(6\),pp\. 6620\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p2.1)\.
- Grmela \(2018\)M\. GrmelaGENERIC guide to the multiscale dynamics and thermodynamics\.Journal of Physics Communications2\(3\),pp\. 032001\.Cited by:[§2\.1](https://arxiv.org/html/2608.12624#S2.SS1.p1.4)\.
- Gruberet al\.\(2025\)A\. Gruber, K\. Lee, H\. Lim, N\. Park, and N\. TraskEfficiently parameterized neural metriplectic systems\.InInternational Conference on Learning Representations,Vol\.2025,pp\. 24835–24859\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p2.1)\.
- Guilhoto and Perdikaris \(2024\)L\. F\. Guilhoto and P\. PerdikarisComposite Bayesian optimization in function spaces using NEON—Neural Epistemic Operator Networks\.Scientific Reports14\(1\),pp\. 29199\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p6.1)\.
- He and Reina \(2025\)Z\. He and C\. ReinaSPIEDiff: robust learning of long\-time macroscopic dynamics from short\-time particle simulations with quantified epistemic uncertainty\.arXiv preprint arXiv:2505\.13501\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p6.1)\.
- He and Reina \(2026\)Z\. He and C\. ReinaEVODMs: variational learning of PDEs for stochastic systems via diffusion models with quantified epistemic uncertainty\.Journal of Computational Physics,pp\. 114722\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p6.1)\.
- Hernándezet al\.\(2021\)Q\. Hernández, A\. Badías, D\. González, F\. Chinesta, and E\. CuetoStructure\-preserving neural networks\.Journal of Computational Physics426,pp\. 109950\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p2.1)\.
- Hersbach \(2000\)H\. HersbachDecomposition of the continuous ranked probability score for ensemble prediction systems\.Weather and Forecasting15\(5\),pp\. 559–570\.Cited by:[Appendix A](https://arxiv.org/html/2608.12624#A1.SS0.SSS0.Px3.p1.1),[§3](https://arxiv.org/html/2608.12624#S3.p1.1)\.
- Huanget al\.\(2022\)S\. Huang, Z\. He, and C\. ReinaVariational Onsager Neural Networks \(VONNs\): a thermodynamics\-based variational learning strategy for non\-equilibrium PDEs\.Journal of the Mechanics and Physics of Solids163,pp\. 104856\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p2.1),[§2\.1](https://arxiv.org/html/2608.12624#S2.SS1.p3.2)\.
- Jacobet al\.\(2025\)B\. Jacob, A\. S\. Nair, A\. A\. Howard, J\. Drgona, and P\. StinisE\-PINNs: Epistemic Physics\-Informed Neural Networks\.arXiv preprint arXiv:2503\.19333\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p6.1),[§2\.3\.1](https://arxiv.org/html/2608.12624#S2.SS3.SSS1.p1.2)\.
- Jinet al\.\(2022\)P\. Jin, Z\. Zhang, I\. G\. Kevrekidis, and G\. E\. KarniadakisLearning Poisson systems and trajectories of autonomous systems via Poisson neural networks\.IEEE Transactions on Neural Networks and Learning Systems34\(11\),pp\. 8271–8283\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p2.1)\.
- Jinet al\.\(2020\)P\. Jin, Z\. Zhang, A\. Zhu, Y\. Tang, and G\. E\. KarniadakisSympNets: intrinsic structure\-preserving symplectic networks for identifying Hamiltonian systems\.Neural Networks132,pp\. 166–179\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p1.1),[§1](https://arxiv.org/html/2608.12624#S1.p2.1)\.
- Joshiet al\.\(2022\)A\. Joshi, P\. Thakolkaran, Y\. Zheng, M\. Escande, M\. Flaschel, L\. De Lorenzis, and S\. KumarBayesian\-EUCLID: discovering hyperelastic material laws with uncertainties\.Computer Methods in Applied Mechanics and Engineering398,pp\. 115225\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p4.1)\.
- Jospinet al\.\(2022\)L\. V\. Jospin, H\. Laga, F\. Boussaid, W\. Buntine, and M\. BennamounHands\-on Bayesian neural networks—a tutorial for deep learning users\.IEEE Computational Intelligence Magazine17\(2\),pp\. 29–48\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p5.1)\.
- Karniadakiset al\.\(2021\)G\. E\. Karniadakis, I\. G\. Kevrekidis, L\. Lu, P\. Perdikaris, S\. Wang, and L\. YangPhysics\-informed machine learning\.Nature Reviews Physics3\(6\),pp\. 422–440\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p1.1)\.
- Kovachkiet al\.\(2023\)N\. Kovachki, Z\. Li, B\. Liu, K\. Azizzadenesheli, K\. Bhattacharya, A\. Stuart, and A\. AnandkumarNeural operator: learning maps between function spaces with applications to pdes\.Journal of Machine Learning Research24\(89\),pp\. 1–97\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p1.1)\.
- Kraaijet al\.\(2020\)R\. C\. Kraaij, A\. Lazarescu, C\. Maes, and M\. PeletierFluctuation symmetry leads to GENERIC equations with non\-quadratic dissipation\.Stochastic Processes and their Applications130\(1\),pp\. 139–170\.Cited by:[§2\.1](https://arxiv.org/html/2608.12624#S2.SS1.p1.4)\.
- Kuleshovet al\.\(2018\)V\. Kuleshov, N\. Fenner, and S\. ErmonAccurate uncertainties for deep learning using calibrated regression\.InInternational conference on machine learning,pp\. 2796–2804\.Cited by:[§2\.4](https://arxiv.org/html/2608.12624#S2.SS4.p1.1)\.
- Lakshminarayananet al\.\(2017\)B\. Lakshminarayanan, A\. Pritzel, and C\. BlundellSimple and scalable predictive uncertainty estimation using deep ensembles\.Advances in neural information processing systems30\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p3.1),[§1](https://arxiv.org/html/2608.12624#S1.p5.1)\.
- Leeet al\.\(2021\)K\. Lee, N\. Trask, and P\. StinisMachine learning structure preserving brackets for forecasting irreversible processes\.Advances in Neural Information Processing Systems34,pp\. 5696–5707\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p2.1)\.
- Leiet al\.\(2018\)J\. Lei, M\. G’Sell, A\. Rinaldo, R\. J\. Tibshirani, and L\. WassermanDistribution\-free predictive inference for regression\.Journal of the American Statistical Association113\(523\),pp\. 1094–1111\.Cited by:[§2\.4](https://arxiv.org/html/2608.12624#S2.SS4.p1.1),[§2\.4](https://arxiv.org/html/2608.12624#S2.SS4.p2.1)\.
- Leviet al\.\(2022\)D\. Levi, L\. Gispan, N\. Giladi, and E\. FetayaEvaluating and calibrating uncertainty prediction in regression tasks\.Sensors22\(15\),pp\. 5540\.Cited by:[§2\.4](https://arxiv.org/html/2608.12624#S2.SS4.p1.1)\.
- Liet al\.\(2020\)Z\. Li, N\. Kovachki, K\. Azizzadenesheli, B\. Liu, K\. Bhattacharya, A\. Stuart, and A\. AnandkumarFourier neural operator for parametric partial differential equations\.arXiv preprint arXiv:2010\.08895\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p1.1)\.
- Linkaet al\.\(2025\)K\. Linka, G\. A\. Holzapfel, and E\. KuhlDiscovering uncertainty: Bayesian constitutive artificial neural networks\.Computer Methods in Applied Mechanics and Engineering433,pp\. 117517\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p4.1)\.
- Luet al\.\(2021\)L\. Lu, P\. Jin, G\. Pang, Z\. Zhang, and G\. E\. KarniadakisLearning nonlinear operators via DeepONet based on the universal approximation theorem of operators\.Nature machine intelligence3\(3\),pp\. 218–229\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p1.1)\.
- Lubliner \(2008\)J\. LublinerPlasticity theory\.Courier Corporation\.Cited by:[§3\.3](https://arxiv.org/html/2608.12624#S3.SS3.p1.1)\.
- Matheson and Winkler \(1976\)J\. E\. Matheson and R\. L\. WinklerScoring rules for continuous probability distributions\.Management science22\(10\),pp\. 1087–1096\.Cited by:[Appendix A](https://arxiv.org/html/2608.12624#A1.SS0.SSS0.Px3.p1.1),[§3](https://arxiv.org/html/2608.12624#S3.p1.1)\.
- Mielke \(2011\)A\. MielkeFormulation of thermoelastic dissipative material behavior using GENERIC\.Continuum Mechanics and Thermodynamics23\(3\),pp\. 233–256\.Cited by:[§3\.3\.1](https://arxiv.org/html/2608.12624#S3.SS3.SSS1.p3.1),[§3\.3](https://arxiv.org/html/2608.12624#S3.SS3.p1.1)\.
- Morrison \(1986\)P\. J\. MorrisonA paradigm for joined Hamiltonian and dissipative systems\.Physica D: Nonlinear Phenomena18\(1\-3\),pp\. 410–419\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p2.1)\.
- Neal \(2012\)R\. M\. NealBayesian learning for neural networks\.Vol\.118,Springer Science & Business Media\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p3.1)\.
- Nováket al\.\(2024\)L\. Novák, H\. Sharma, and M\. D\. ShieldsPhysics\-informed polynomial chaos expansions\.Journal of Computational Physics506,pp\. 112926\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p3.1)\.
- Osbandet al\.\(2023a\)I\. Osband, Z\. Wen, S\. M\. Asghari, V\. Dwaracherla, M\. Ibrahimi, X\. Lu, and B\. Van RoyApproximate thompson sampling via epistemic neural networks\.InUncertainty in Artificial Intelligence,pp\. 1586–1595\.Cited by:[§4](https://arxiv.org/html/2608.12624#S4.p2.1)\.
- Osbandet al\.\(2023b\)I\. Osband, Z\. Wen, S\. M\. Asghari, V\. Dwaracherla, M\. Ibrahimi, X\. Lu, and B\. Van RoyEpistemic neural networks\.Advances in Neural Information Processing Systems36,pp\. 2795–2823\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p6.1),[§2\.2](https://arxiv.org/html/2608.12624#S2.SS2.p1.1),[§2\.2](https://arxiv.org/html/2608.12624#S2.SS2.p1.2),[§2\.3\.2](https://arxiv.org/html/2608.12624#S2.SS3.SSS2.p1.1),[§2\.3\.2](https://arxiv.org/html/2608.12624#S2.SS3.SSS2.p4.1),[§4](https://arxiv.org/html/2608.12624#S4.p2.1)\.
- Öttinger and Grmela \(1997\)H\. C\. Öttinger and M\. GrmelaDynamics and thermodynamics of complex fluids\. II\. illustrations of a general formalism\.Physical Review E56\(6\),pp\. 6633\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p2.1)\.
- Öttinger \(2005\)H\. C\. ÖttingerBeyond equilibrium thermodynamics\.John Wiley & Sons\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p2.1)\.
- Papadopouloset al\.\(2002\)H\. Papadopoulos, K\. Proedrou, V\. Vovk, and A\. GammermanInductive confidence machines for regression\.InEuropean conference on machine learning,pp\. 345–356\.Cited by:[§2\.4](https://arxiv.org/html/2608.12624#S2.SS4.p1.1),[§2\.4](https://arxiv.org/html/2608.12624#S2.SS4.p2.1),[§2\.4](https://arxiv.org/html/2608.12624#S2.SS4.p4.2)\.
- Pavelkaet al\.\(2018\)M\. Pavelka, V\. Klika, and M\. GrmelaMultiscale thermo\-dynamics: introduction to GENERIC\.Walter de Gruyter GmbH & Co KG\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p2.1)\.
- Psaroset al\.\(2023\)A\. F\. Psaros, X\. Meng, Z\. Zou, L\. Guo, and G\. E\. KarniadakisUncertainty quantification in scientific machine learning: methods, metrics, and comparisons\.Journal of Computational Physics477,pp\. 111902\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p3.1)\.
- Raissiet al\.\(2019\)M\. Raissi, P\. Perdikaris, and G\. E\. KarniadakisPhysics\-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations\.Journal of Computational physics378,pp\. 686–707\.Cited by:[Appendix D](https://arxiv.org/html/2608.12624#A4.p2.1),[§1](https://arxiv.org/html/2608.12624#S1.p1.1)\.
- Ross and Heinonen \(2023\)M\. Ross and M\. HeinonenLearning energy conserving dynamics efficiently with Hamiltonian Gaussian processes\.Transactions on Machine Learning Research\.External Links:ISSN 2835\-8856,[Link](https://openreview.net/forum?id=DHEZuKStzH)Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p4.1)\.
- Shafer and Vovk \(2008\)G\. Shafer and V\. VovkA tutorial on conformal prediction\.\.Journal of machine learning research9\(3\)\.Cited by:[§2\.4](https://arxiv.org/html/2608.12624#S2.SS4.p4.2)\.
- Sharmaet al\.\(2024\)H\. Sharma, J\. A\. Gaffney, D\. Tsapetis, and M\. D\. ShieldsLearning thermodynamically constrained equations of state with uncertainty\.APL Machine Learning2\(1\)\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p4.1)\.
- Shin and Choi \(2023\)H\. Shin and M\. ChoiPhysics\-informed variational inference for uncertainty quantification of stochastic differential equations\.Journal of Computational Physics487,pp\. 112183\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p3.1)\.
- Simo and Hughes \(1998\)J\. C\. Simo and T\. J\. HughesComputational inelasticity\.Springer\.Cited by:[§3\.3](https://arxiv.org/html/2608.12624#S3.SS3.p1.1)\.
- Söderström \(2007\)T\. SöderströmErrors\-in\-variables methods in system identification\.Automatica43\(6\),pp\. 939–958\.Cited by:[§4](https://arxiv.org/html/2608.12624#S4.p2.1)\.
- Stankeviciuteet al\.\(2021\)K\. Stankeviciute, A\. M Alaa, and M\. Van der SchaarConformal time\-series forecasting\.Advances in neural information processing systems34,pp\. 6216–6228\.Cited by:[§2\.4](https://arxiv.org/html/2608.12624#S2.SS4.p2.1)\.
- Sun and Yu \(2024\)S\. H\. Sun and R\. YuCopula Conformal prediction for multi\-step time series prediction\.InThe Twelfth International Conference on Learning Representations,External Links:[Link](https://openreview.net/forum?id=ojIJZDNIBj)Cited by:[§2\.4](https://arxiv.org/html/2608.12624#S2.SS4.p2.1)\.
- Tanakaet al\.\(2022\)Y\. Tanaka T\. Iwataet al\.Symplectic Spectrum Gaussian Processes: Learning Hamiltonians from noisy and sparse data\.Advances in Neural Information Processing Systems35,pp\. 20795–20808\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p4.1)\.
- Votrubaet al\.\(2026\)V\. Votruba, Z\. He, W\. Qiu, C\. Reina, and M\. PavelkaNonlinear GENERIC Informed Neural Networks \(N\-GINNs\): learning GENERIC dynamics with non\-quadratic dissipation potentials\.arXiv preprint arXiv:2605\.09058\.Cited by:[Appendix C](https://arxiv.org/html/2608.12624#A3.p1.1),[§1](https://arxiv.org/html/2608.12624#S1.p2.1),[§1](https://arxiv.org/html/2608.12624#S1.p6.1),[§1](https://arxiv.org/html/2608.12624#S1.p7.1),[§2\.1](https://arxiv.org/html/2608.12624#S2.SS1.p2.1),[§3\.1](https://arxiv.org/html/2608.12624#S3.SS1.p1.1),[§3\.2](https://arxiv.org/html/2608.12624#S3.SS2.p1.1),[§3\.2](https://arxiv.org/html/2608.12624#S3.SS2.p2.6),[§3\.3\.1](https://arxiv.org/html/2608.12624#S3.SS3.SSS1.p1.2)\.
- Wan and Karniadakis \(2006\)X\. Wan and G\. E\. KarniadakisMulti\-element generalized polynomial chaos for arbitrary probability measures\.SIAM Journal on Scientific Computing28\(3\),pp\. 901–928\.Cited by:[§D\.1](https://arxiv.org/html/2608.12624#A4.SS1.p1.1)\.
- Xu and Xie \(2021\)C\. Xu and Y\. XieConformal prediction interval for dynamic time\-series\.InInternational Conference on Machine Learning,pp\. 11559–11569\.Cited by:[Appendix A](https://arxiv.org/html/2608.12624#A1.SS0.SSS0.Px2.p1.2)\.
- Yanget al\.\(2021\)L\. Yang, X\. Meng, and G\. E\. KarniadakisB\-PINNs: Bayesian physics\-informed neural networks for forward and inverse PDE problems with noisy data\.Journal of Computational Physics425,pp\. 109913\.Cited by:[Appendix D](https://arxiv.org/html/2608.12624#A4.p4.1),[§1](https://arxiv.org/html/2608.12624#S1.p3.1),[§2\.3\.1](https://arxiv.org/html/2608.12624#S2.SS3.SSS1.p2.1)\.
- Yang and Perdikaris \(2019\)Y\. Yang and P\. PerdikarisAdversarial uncertainty quantification in physics\-informed neural networks\.Journal of Computational Physics394,pp\. 136–152\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p3.1)\.
- Yuet al\.\(2021\)H\. Yu, X\. Tian, W\. E, and Q\. LiOnsagerNet: learning stable and interpretable dynamics using a generalized Onsager principle\.Physical Review Fluids6\(11\),pp\. 114402\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p1.1)\.
- Yuet al\.\(2026\)Y\. Yu, C\. H\. Ho, and Y\. WangA conformal prediction framework for uncertainty quantification in physics\-informed neural networks\.Journal of Computational Physics561,pp\. 114979\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p3.1)\.
- Zelikmanet al\.\(2020\)E\. Zelikman, C\. Healy, S\. Zhou, and A\. AvatiCRUDE: calibrating regression uncertainty distributions empirically\.arXiv preprint arXiv:2005\.12496\.Cited by:[§2\.4](https://arxiv.org/html/2608.12624#S2.SS4.p1.1)\.
- Zhanget al\.\(2019\)D\. Zhang, L\. Lu, L\. Guo, and G\. E\. KarniadakisQuantifying total uncertainty in physics\-informed neural networks for solving forward and inverse stochastic problems\.Journal of Computational Physics397,pp\. 108850\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p3.1)\.
- Zhanget al\.\(2022\)Z\. Zhang, Y\. Shin, and G\. Em KarniadakisGFINNs: GENERIC formalism informed neural networks for deterministic and stochastic dynamical systems\.Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences380\(2229\)\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p1.1),[§1](https://arxiv.org/html/2608.12624#S1.p2.1)\.
- Zhong and Meidani \(2023\)W\. Zhong and H\. MeidaniPI\-VAE: Physics\-Informed Variational Auto\-encoder for stochastic differential equations\.Computer Methods in Applied Mechanics and Engineering403,pp\. 115664\.Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p3.1)\.
- Zhonget al\.\(2020\)Y\. D\. Zhong, B\. Dey, and A\. ChakrabortySymplectic ODE\-Net: learning Hamiltonian dynamics with control\.InInternational Conference on Learning Representations,External Links:[Link](https://openreview.net/forum?id=ryxmb1rKDS)Cited by:[§1](https://arxiv.org/html/2608.12624#S1.p2.1)\.
- Zouet al\.\(2025\)Z\. Zou, X\. Meng, and G\. E\. KarniadakisUncertainty quantification for noisy inputs–outputs in physics\-informed neural networks and neural operators\.Computer Methods in Applied Mechanics and Engineering433,pp\. 117479\.Cited by:[§4](https://arxiv.org/html/2608.12624#S4.p2.1)\.
- Zouet al\.\(2024\)Z\. Zou, X\. Meng, A\. F\. Psaros, and G\. E\. KarniadakisNeuralUQ: a comprehensive library for uncertainty quantification in neural differential equations and operators\.SIAM Review66\(1\),pp\. 161–190\.Cited by:[§D\.1](https://arxiv.org/html/2608.12624#A4.SS1.p1.2),[§D\.2](https://arxiv.org/html/2608.12624#A4.SS2.p2.1),[Appendix D](https://arxiv.org/html/2608.12624#A4.p4.1),[§2\.3\.1](https://arxiv.org/html/2608.12624#S2.SS3.SSS1.p2.1)\.