Operator Inference for Model Reduction

• Operator Inference is a non-intrusive, data-driven method that infers reduced operators from snapshot data while retaining projection-based structure.
• It utilizes regression on reduced coordinates to capture the dynamics of high-dimensional PDE/ODE systems, with variants enforcing energy preservation, stability, or Hamiltonian structure.
• Recent extensions address parametric dependence, noisy data, scalability via tensor methods, and the integration of neural operators for flexible modeling.

Searching arXiv for recent Operator Inference papers to ground the article in current literature. Operator Inference is a non-intrusive, data-driven methodology for learning low-dimensional dynamical models from snapshot data while retaining a projection-based reduced-order modeling architecture. In its classical form, it begins from a high-dimensional system obtained, for example, by spatial discretization of a PDE, constructs a reduced basis—typically by Proper Orthogonal Decomposition (POD)—projects state data into reduced coordinates, and then infers reduced operators by regression rather than by explicit projection of inaccessible full-order operators. Across the recent literature, this basic idea appears in continuous-time and discrete-time forms, in first-order and second-order settings, and in variants that enforce energy preservation, Hamiltonian structure, stability, sparsity, or parameter dependence (Gosea et al., 2024, Koike et al., 2024, Filanova et al., 2022, Benner et al., 7 Jul 2025).

1. Definition and conceptual scope

Within the literature represented here, Peherstorfer and Willcox (2016) are identified as introducing foundational Operator Inference for nonintrusive projection-based model reduction, with the central idea of fitting reduced operators of a prescribed polynomial structure directly from data (Gosea et al., 2024). The method occupies a middle position between intrusive projection-based reduced-order modeling and fully black-box learning: it keeps the reduced-state, projection-based viewpoint of POD–Galerkin approaches, but replaces direct access to full-order operators by regression on trajectories.

A standard point of departure is a semi-discrete PDE or ODE system in a high-dimensional state $x(t)\in\mathbb{R}^n$, together with snapshot data $x(t_i)$ and, in many formulations, time-derivative data $\dot{x}(t_i)$. The reduced basis $V\in\mathbb{R}^{n\times r}$ defines reduced coordinates $\hat{x}(t)=V^\top x(t)$, and Operator Inference learns the reduced operators governing $\hat{x}$. In this sense it is both non-intrusive and projection-based: the reduced state evolves in a low-dimensional subspace, but the reduced operators are identified from data rather than obtained by projecting full-order matrices.

The term is used most often in this reduced-operator sense, but the literature also contains a broader operator-learning usage. For elliptic eigenvalue problems, “operator inference” is used for learning domain-to-spectrum and domain-to-eigenfunction maps from pixelized geometries, without constructing reduced matrix operators in POD coordinates (Li et al., 22 Apr 2025). This broader usage differs materially from the classical OpInf program for reduced-order dynamical systems.

2. Core mathematical formulation

A canonical continuous-time formulation assumes that the reduced dynamics obey a structured polynomial ansatz. In the methanation reactor study, the starting point is a quadratic full-order model

$\dot{x}(t)=Ax(t)+H(x(t)\otimes x(t))+f,$

and the inferred reduced model is

$$\dot{\hat{x}}(t)=\hat{A}\hat{x}(t)+\hat{H}(\hat{x}(t)\otimes\hat{x}(t))+\hat{f},$$

with $\hat{x}(t)\in\mathbb{R}^r$, $\hat{A}\in\mathbb{R}^{r\times r}$, $x(t_i)$0, and $x(t_i)$1 (Gosea et al., 2024). Reduced snapshots $x(t_i)$2 and reduced derivatives $x(t_i)$3 are collected, together with the quadratic snapshot matrix built from $x(t_i)$4, and the operators are obtained from a linear least-squares problem in Frobenius norm.

The energy-preserving literature writes the same idea in a slightly different matrix notation. With reduced snapshots $x(t_i)$5, derivative data $x(t_i)$6, and quadratic data $x(t_i)$7, standard OpInf solves

$x(t_i)$8

where $x(t_i)$9 and $\dot{x}(t_i)$0 (Koike et al., 2024). In this form, each row is decoupled, so the unconstrained problem reduces to $\dot{x}(t_i)$1 independent least-squares solves.

Discrete-time variants also exist. In noisy-data settings, one writes a polynomial reduced model

$\dot{x}(t_i)$2

and solves a least-squares problem in the reduced coordinates with a regression matrix built from monomial features of $\dot{x}(t_i)$3 and the inputs $\dot{x}(t_i)$4 (Uy et al., 2021). This formulation is central to the analysis of unbiasedness and mean-squared-error bounds under re-projection-based sampling.

A distinct reformulation replaces derivative fitting by trajectory fitting. In the differentiable-solver approach, the reduced dynamics are embedded as equality constraints in an ODE-constrained optimization problem, and the reduced operators are calibrated by differentiating through the time-stepping solver rather than by fitting the right-hand side pointwise (Hartmann et al., 2021). The objective becomes the mismatch between observed and simulated trajectories over the entire time horizon.

3. Reduced spaces, coordinates, and system classes

The reduced basis is most commonly constructed by POD from an SVD of the snapshot matrix. In the methanation reactor example, the snapshot matrix $\dot{x}(t_i)$5 is decomposed as $\dot{x}(t_i)$6, the reduced basis is $\dot{x}(t_i)$7, and the retained dimension $\dot{x}(t_i)$8 is chosen by an energy criterion; for the main model, $\dot{x}(t_i)$9 captures $V\in\mathbb{R}^{n\times r}$0 of the energy, with $V\in\mathbb{R}^{n\times r}$1 modes for conversion and $V\in\mathbb{R}^{n\times r}$2 modes for temperature (Gosea et al., 2024). In energy-preserving OpInf, a related criterion

$V\in\mathbb{R}^{n\times r}$3

is used to select $V\in\mathbb{R}^{n\times r}$4 from the POD spectrum (Koike et al., 2024).

Second-order mechanical formulations preserve the mass–damping–stiffness structure instead of rewriting the system as first order. For linear elastodynamics, the full model

$V\in\mathbb{R}^{n\times r}$5

is reduced directly in second-order form, and the inferred normalized reduced model is

$V\in\mathbb{R}^{n\times r}$6

with regression performed on projected displacement, velocity, and acceleration data (Filanova et al., 2022). When the complete external force is available, the same paper infers $V\in\mathbb{R}^{n\times r}$7 directly under positive-definiteness constraints.

Operator Inference has also been adapted to constrained mechanical systems governed by index-2 and index-3 DAEs. There, the key observation is that proper DAEs admit an underlying ODE on a hidden constraint manifold. By building the POD basis from constraint-consistent displacement or velocity snapshots, one obtains a reduced coordinate system already contained in the manifold, and OpInf then identifies the hidden second-order ODE directly in reduced space (Benner et al., 7 Jul 2025). In this setting, the inferred reduced operators $V\in\mathbb{R}^{n\times r}$8 are obtained from reduced displacement, velocity, acceleration, input, and output data.

The incompressible-flow literature adds a further structural specialization. For semi-discrete Navier–Stokes equations written as a quadratic DAE, pressure can be eliminated through a discrete pressure Poisson equation, producing an underlying quadratic ODE for the velocity alone and a separate algebraic relation for pressure. Operator Inference is then applied to the reduced velocity ODE, while pressure is inferred through a separate algebraic map from reduced velocity and input (Benner et al., 2020).

4. Structure preservation, stability, and physical constraints

A central theme in recent OpInf research is that least-squares fitting alone does not generally preserve structure. In the energy-preserving setting, the quadratic term is required to satisfy

$V\in\mathbb{R}^{n\times r}$9

or, in tensor components under symmetry in the last two indices,

$\hat{x}(t)=V^\top x(t)$0

Standard OpInf does not enforce this and can therefore introduce artificial energy production or dissipation; EP-OpInf imposes the condition through constrained least squares and a compressed representation of the quadratic operator (Koike et al., 2024).

Hamiltonian operator inference enforces a different structure. For canonical Hamiltonian systems $\hat{x}(t)=V^\top x(t)$1, the inferred reduced operator is constrained so that the reduced system remains Hamiltonian, and for noncanonical systems the inferred reduced Poisson operator is constrained to be skew-symmetric. The resulting C-H-OpInf and NC-H-OpInf formulations are linear solves with symmetry or skew-symmetry constraints, and the reduced models preserve Hamiltonian structure by construction (Gruber et al., 2023).

Stability-preserving formulations are especially prominent in process engineering and structural dynamics. In the methanation-reactor study, the reduced operators are learned with a stability-oriented parameterization following Goyal–Duff–Benner (2023): $\hat{x}(t)=V^\top x(t)$2 is parameterized using symmetric positive definite and skew-symmetric components so that its eigenvalues lie in the left half-plane, and $\hat{x}(t)=V^\top x(t)$3 is parameterized so that the quadratic form defines a global Lyapunov function, yielding an asymptotically stable reduced model (Gosea et al., 2024). The optimization is performed in PyTorch with Adam, a CyclicLR scheduler in triangular2 mode, early stopping, and a regularization factor $\hat{x}(t)=V^\top x(t)$4 on $\hat{x}(t)=V^\top x(t)$5.

For nonlinear structural dynamics, stability is enforced through polynomial Lyapunov structure. The reduced model is written as

$\hat{x}(t)=V^\top x(t)$6

where $\hat{x}(t)=V^\top x(t)$7 is a polynomial potential. Global uniform boundedness or input-to-state stability is enforced by requiring the potential and related coercivity expressions to belong to the sum-of-squares cone, which converts the inference problem into a convex semidefinite program (Boef et al., 2024). In the DAE mechanical setting, symmetric positive (semi)definiteness of the reduced mass, damping, and stiffness operators is imposed by semidefinite programming, again to guarantee stability and interpretability (Benner et al., 7 Jul 2025).

These structure-preserving formulations also clarify a common misconception: OpInf is not inherently structure-preserving. Conservation, Hamiltonian form, symmetry, energy preservation, or stability must be encoded by constrained optimization, operator parameterization, or problem-specific ansätze; otherwise the learned reduced operators are determined only by data fit (Koike et al., 2024).

5. Parametric, noisy, and scalable extensions

A major extension is parametric Operator Inference. In the OpInf-LLM framework, a shared POD basis $\hat{x}(t)=V^\top x(t)$8 is built across parameter values for each PDE family, and for each training parameter $\hat{x}(t)=V^\top x(t)$9 one learns reduced operators $\hat{x}$0. At query time, operator entries are inferred at new parameters by element-wise interpolation or linear regression in parameter space, after which the reduced ODE is integrated and the field is reconstructed (Wang et al., 2 Feb 2026). This is a nonintrusive parametric OpInf strategy rather than an affine intrusive projection.

Noisy-data OpInf motivates a different set of ideas. Under re-projection, the regression matrix can be made deterministic while only the response is noisy, and the inferred operators become unbiased estimators of the projection-based reduced operators. The operator mean-squared error scales with $\hat{x}$1, and the resulting prediction error bounds motivate active operator inference: snapshot selection is designed to maximize the smallest singular value of the regression matrix, thereby reducing the noise-to-signal ratio and the prediction MSE (Uy et al., 2021).

Scalability has been addressed through tensor methods. Tensor-Train Operator Inference introduces full-order TT-OpInf, full-order QTT-OpInf, and reduced-order TT-ROM OpInf. In these variants, snapshot data, quadratic terms, design tensors, and even inferred operators are represented in tensor-train format, and least-squares problems are solved through TT pseudoinverse and TT contractions (Danis et al., 9 Sep 2025). The reduced-order TT-ROM variant preserves the standard reduced OpInf least-squares stage but replaces matrix POD by TT-based compression of the snapshot tensor.

The operator class itself has also been generalized. NN-OpInf replaces fixed polynomial operators by structure-preserving composable neural operators, including skew-symmetric, symmetric positive semidefinite, and potential-based blocks, while still operating in POD coordinates and learning reduced latent dynamics from projected data (Parish et al., 9 Mar 2026). For parameter-dependent matrix equations, OpInf has been adapted to vectorized algebraic problems such as Lyapunov and Riccati equations by rewriting the full problem as a polynomial algebraic equation in the vectorized state and inferring reduced polynomial operators from solution snapshots (Wen et al., 20 Nov 2025).

6. Applications, performance, limitations, and terminological breadth

The application range is broad. In process engineering, Operator Inference has been applied to a fixed-bed CO$\hat{x}$2 methanation reactor modeled by a 1D pseudo-homogeneous PDE system with energy and mass balances, spatially discretized by finite volumes on 200 control volumes. The resulting quadratic–linear ROM with $\hat{x}$3 achieves a relative Frobenius error of approximately $\hat{x}$4 and requires only $\hat{x}$5 of the full-model runtime, while capturing conversion and temperature profiles through the startup phase (Gosea et al., 2024). The same paper explicitly positions this as “an important milestone towards the implementation of fast and reliable digital twin architectures.”

In design and control for fluid-like PDEs, EP-OpInf preserves the energy-preserving structure of quadratic convection terms. For Burgers and Kuramoto–Sivashinsky, intrusive and EP-OpInf reduced operators satisfy the energy-preserving constraint up to machine precision, whereas standard OpInf exhibits significant violations, even when its short-horizon data fit is strong (Koike et al., 2024). In incompressible flows, velocity-only OpInf ROMs with separate pressure reconstruction provide competitive or superior performance relative to intrusive POD and DMD-based methods on driven-cavity and cylinder-wake benchmarks (Benner et al., 2020).

Mechanical applications include linear elastodynamics, constrained multibody DAEs, and nonlinear structural dynamics. In the constrained DAE setting, reduced dimensions of $\hat{x}$6 and $\hat{x}$7 are reported for the anchored and triple-chain examples, with good performance under both training and test loads (Benner et al., 7 Jul 2025). In nonlinear structural dynamics, quartic stable sparse OpInf ROMs substantially outperform linear ROMs and remain stable under reduced training data, whereas unconstrained variants can become unstable on validation trajectories (Boef et al., 2024).

Limitations recur across the literature. Several papers emphasize that classical continuous-time OpInf assumes access to derivative information or accurate numerical differentiation, that regression matrices may be ill-conditioned, and that extrapolation beyond the training regime is not guaranteed (Gosea et al., 2024, Uy et al., 2021). Structure-preserving variants often incur higher computational cost because constrained least-squares, semidefinite programs, or nonconvex neural optimization replace straightforward row-decoupled regression (Koike et al., 2024, Benner et al., 7 Jul 2025, Parish et al., 9 Mar 2026). Generalization across operating conditions remains a distinct issue: the methanation study, for example, validates a single startup scenario and explicitly leaves variable parameters and variable input loads to future work (Gosea et al., 2024).

A final terminological point is that the phrase “operator inference” is no longer confined to reduced polynomial ODE learning. In the elliptic eigenvalue literature, it denotes direct learning of domain-to-eigenvalue and domain-to-eigenfunction maps using CNNs and Fourier Neural Operators, with geometry preprocessing such as domain scaling, main-axis alignment, and detailed pixelization (Li et al., 22 Apr 2025). That usage shares the idea of learning an operator-valued map from data, but it differs from classical OpInf in its absence of POD coordinates, reduced dynamical operators, or least-squares identification of projected ODE coefficients.