[orcid=0009-0000-9165-9844]

[orcid = 0000-0003-4564-5999]

Beyond linear subspaces: Nonlinear moment matching meets quadratic manifolds

Reetish Padhi reetishp@vt.edu https://rewtus.github.io/reetish/ organization=Department of Mathematics, Virginia Tech, city=Blacksburg, postcode=24061, state=Virginia, country=USA Serkan Gugercin gugercin@vt.edu https://gugercin.math.vt.edu organization=Department of Mathematics, Virginia Tech, city=Blacksburg, postcode=24061, state=Virginia, country=USA

Abstract

Quadratic manifold-based model order reduction offers a viable pathway to circumvent the limitations of linear subspaces for linear control systems characterized by slow Kolmogorov $n$ -width decay. However, a system-theoretic framework for constructing such quadratic approximations remains absent from the literature. This paper presents a system-agnostic, optimization-free framework for the direct construction of quadratic projection matrices. We prove that the synthesized reduced-order model matches the nonlinear moments of the full-order system and preserves its exact center manifold mapping, thereby ensuring asymptotic tracking of steady-state outputs under specific input classes. Numerical results on transport-dominated benchmark problems, namely, the one-dimensional damped wave and advection equations, show that the proposed framework achieves high-fidelity trajectory reconstruction within a significantly reduced-dimensional state space, yielding substantial online computational savings.

keywords

Quadratic manifolds,Model order reduction ,Nonlinear moment matching ,Kolmogorov

n

-width ,LTI systems ,Transport dominated problems

^†^†credit: Conceptualization, Investigation, Methodology, Software, Formal analysis, Validation, Writing – original draft, Writing – review and editing^†^†credit: Conceptualization, Investigation, Methodology, Supervision, Formal analysis, Validation, Writing – review and editing^†^†corresponding: Corresponding author

1 Introduction

Model Order Reduction (MOR) has emerged as a critical tool for the simulation and control of high-fidelity, large-scale dynamical systems. By approximating the state of a high-dimensional system within a low-dimensional subspace/manifold, MOR enables significant computational savings while maintaining key mathematical properties of the original system. Traditional reduction techniques focus on constructing a reduced model by projecting the state onto a linear subspace. For linear time-invariant (LTI) systems, (linear) projection-based methods such as interpolatory methods morAntBG20, balancing-based approaches BenB17, and Proper Orthogonal Decomposition (POD) volkwein2011pod are well-established. However, the efficacy of linear MOR is fundamentally governed by the Kolmogorov $n$ -width pinkus2012n, which represents the worst-case error arising from the projection of the solution manifold onto the best-possible linear subspace of dimension $r\ll n$ . For many elliptic or parabolic PDEs, the Kolmogorov $n$ -width decays exponentially with $r$ , allowing for low dimensional reduced order models (ROM). Conversely, for hyperbolic or transport-dominated problems, the decay of these widths is significantly slower kol_wave_decay; peherstorfer2022breaking. This “slow decay” constitutes a barrier to reducibility for linear methods. In unger2019kolmogorov, the authors show that for LTI systems, the Kolmogorov $n$ -widths coincide with the Hankel singular values ACA05, thus connecting the concept of the Kolmogorov $n$ -widths to system-theoretic objects.

Nonlinear dimensionality reduction techniques have emerged as a means to circumvent the Kolmogorov $n$ -width barrier using the theory of nonlinear manifolds QM_framework. Among these, machine learning approaches like the autoencoder-based frameworks, utilize deep learning architectures to map high-dimensional data to a latent space via a nonlinear mapping buchfink2023symplectic; otto2023learning; fresca2021comprehensive; fresca2022pod; kim2022fast; kadeethum2022non. Alternatively, dictionary and localized basis methods, e.g., amsallem2012nonlinear; daniel2022physics; geelen2022localized; rewienski2003trajectory, construct a library of localized, typically time-invariant bases that partition the state space, adaptively selecting the optimal local subspace online during system evolution.

Shift based methods consider explicit, time-dependent spatial shifts or transport maps to align moving features along solution trajectories reiss2018shifted; black2020projection; reiss2021optimization; papapicco2022neural; burela2023parametric; krah2025robust. Freezing techniques and symmetry-based reduction frameworks exploit symmetry/equivariance of the original system. They involve Lie group actions and algebraic conditions to transform the governing equations into a moving frame of reference, effectively “freezing” traveling fronts into stationary profiles before performing dimensionality reduction rowley2000reconstruction; rowley2003reduction; beyn2004freezing; ohlberger2013nonlinear. A comprehensive review and classification of these nonlinear approaches, particularly for transport-dominated problems, can be found in hesthaven2026nonlinear.

Existing literature on quadratic manifolds BARNETT_QM1; benner2023quadratic; geelen2023operator; sharma2023symplectic; schwerdtner2024greedy; schwerdtner2025empirical; schwerdtner2025online; paxton2026fast; glas2026structure; rutzmoser2017generalization; schwerdtner2024online has considered different methods of picking the linear and quadratic projection matrices. For instance, rutzmoser2017generalization utilizes structural properties of the governing equations, while geelen2023operator proposes an Operator Inference (OpInf) approach, where the projection subspaces are learned in a data-driven fashion by solving an optimization problem over state snapshots. In benner2023quadratic, the authors propose a quadratic decoder approach for the approximation of nonlinear systems. A greedy approach of selecting quadratic projection bases has been presented in schwerdtner2024greedy which was adapted to an online streaming setting in schwerdtner2024online. On the other hand, in paxton2026fast, the quadratic projection matrices are constructed by solving an optimization problem over the Stiefel manifold. Structure preserving model reduction using quadratic manifolds for Port-Hamiltonian systems has been studied in sharma2023symplectic; glas2026structure. In QM_framework, the authors present a unifying differential-geometric framework for nonlinear manifold model reduction. Higher order polynomial and rational manifold approximations have been explored in geelen2023learning; geelen2024learning; klein2025entropy; buchfink2024approximation.

While quadratic manifold approaches offer a powerful alternative to linear subspaces, they yield reduced systems with higher-order state (nonlinear) dependencies. This structural complexity often makes it difficult to establish direct theoretical comparisons or rigorous statements relating the full- and reduced-order models. Furthermore, existing methods typically rely on optimization formulations to construct these manifolds, which require empirical hyperparameter tuning. To overcome these limitations, the main goal of this work is to provide a purely system-theoretic construction of the quadratic manifold.

A natural choice for establishing system-theoretic guarantees is interpolatory (also known as moment-matching based) model reduction, which traditionally focuses on rational interpolation of transfer functions in the frequency domain; see morAntBG20 for a comprehensive survey on linear time-invariant systems. For structured nonlinear systems, such as bilinear and quadratic-bilinear systems, the input-output behaviour is characterized by Volterra kernels rugh1981nonlinear and the interpolatory methods (moment matching) extends to interpolating these kernels or multivariate transfer functions; see, e.g., flagg2015multipoint; benner2024structured; werner2021structure; gosea2018data; benner2015two; breiten2010krylov; gu2011qlmor; breiten2012interpolation; bai2006projection and the references therein. Alternatively, in astolfi2010, the notion of moment matching was reformulated as matching the steady-state output response of a dynamical system driven by an autonomous signal generator. In this framework, “nonlinear moments” of a dynamical system is characterized by an invariant center manifold mapping satisfying a generalized partial differential equation (invariance equation). This formulation has been expanded to quadratic-bilinear and general nonlinear systems bai2022model; scarciotti2017nonlinear; simard2024parameterization, and extended to non-intrusive, data-driven settings scarciotti2017data; moreschini2025moment. We refer the reader to astolfi_10year_survey; scarciotti2024interconnection for a detailed survey of these methods.

While center-manifold-based moment matching provides an elegant framework for analyzing steady-state output responses, it does not directly yield explicit projection bases for state-space model reduction. In that setting, the invariant manifold mapping remains an abstract time-domain evaluation operator rather than an explicit tool for constructing trial spaces. In this paper, we bridge the gap between nonlinear moment-matching theory and quadratic manifold MOR by reformulating moment matching within a quadratic projection framework. Specifically, Theorem 4.1 establishes that a quadratic manifold structure emerges naturally for a special class of signal generators, directly yielding explicit projection matrices. This formulation provides a practical, computationally efficient approach to construct quadratic reduced-order models with system-theoretic guarantees. By shifting the paradigm from empirical snapshot optimization to system-theoretic interpolation, we prove that the quadratic reduced system matches the nonlinear moments generated by a user-specified signal space (Figure 1). The primary contributions of this work are summarized below:

Figure 1: Schematic representing the nonlinear moment matching framework using quadratic manifolds. The center manifold of the reduced system (dimension

r

) under the quadratic approximation maps to the center manifold of the original system (dimension

n

) and hence has the same asymptotic steady state output as the original system for inputs generated by signal space. The linear and quadratic matrices

{\mathbf{V}}_{1}

and

{\mathbf{V}}_{2}

serve the role of defining the projection matrices for quadratic state reduction.

•

We introduce an optimization-free algorithm for constructing quadratic projection matrices by solving a sequence of decoupled linear Sylvester equations. In this approach, the choice of the signal space acts as a high-level design specification that can be freely adapted to match the characteristic frequencies or transient behaviors of a target application. Once this operating signal space is selected, the projection bases are uniquely determined with a closed-form expression. Unlike prevailing snapshot-based quadratic manifold techniques, this method requires no subsequent empirical hyperparameter tuning, optimization parameter sweeps, or regularization.
•

We establish the core theoretical foundation of this framework in Theorem 4.5. Specifically, we prove that the resulting reduced-order model preserves the nonlinear moments of the full-order system, using the notion of nonlinear moments introduced in astolfi2010. Additionally, we show that the center manifold of the reduced system maps to the quadratic center manifold of the original full-order model (FOM), guaranteeing identical steady-state outputs for signals generated by the specified signal space. This leads to Algorithm 1, which provides a projection-based reformulation of the nonlinear-moment matching framework in the context of quadratic manifolds.
•

We demonstrate the efficacy of the proposed framework on the transport-dominated, one-dimensional damped wave and advection equations. The numerical experiments verify that our quadratic ROM accurately recovers the high-fidelity, full-order steady-state trajectories on the center manifold and thereby, achieves “nonlinear moment matching”. By projecting the high-dimensional state space down to the low dimensional subspace, the framework yields a substantial reduction in simulation time during the ODE solver integration phase, proving that the online savings vastly outweigh the overhead of handling the state-dependent reduced mass matrix. However, state-independent reduced mass matrix can be enforced as in other quadratic manifold approaches.

The remainder of this paper is structured as follows. In Section 2, we recall relevant matrix and tensor definitions and establish the foundational theory behind nonlinear moment matching, invariance equations, and interpolatory model reduction. Section 3 introduces the formulation of the quadratic reduced-order model (23)–(24) resulting from projecting linear full-order dynamics onto quadratic state approximation. In Section 4, we establish the core theoretical foundation of the framework. First, we prove that the center manifold of the full-order system under the driving signal generator is a quadratic manifold and derive closed-form expressions for the projection matrices as solutions of decoupled Sylvester equations. Finally, we prove that the quadratic reduced system constructed with these matrices achieves nonlinear moment matching. The proposed computational framework based on this theoretical analysis is summarized in Algorithm 1. Section 5 provides practical guidelines for parameterizing the signal generator to enforce application-specific interpolation conditions. Section 6 validates the theoretical guarantees on transport-dominated benchmarks, specifically the 1D advection and damped wave equations, and compares performance against standard linear rational interpolation and existing quadratic manifold techniques.

2 Background

In this section, we establish the mathematical notation, definitions, and prerequisite theory for linear and nonlinear moment matching. We begin in Section 2.1 by defining the matrix, vector, and tensor operations used throughout the paper. Section 2.2 reviews center manifold theory and the notion of nonlinear moments for a dynamical system along with connections with linear moment matching framework in Section 2.2.1. Finally, Section 2.2.2 reviews the formal definitions of moment matching and interpolating reduced-order model.

2.1 Notation

In this section, we establish the mathematical notation and tensor operations utilized throughout this manuscript. For a review of these standard definitions and properties, the reader is referred to brewer1978kronecker; graham2018kronecker. We let ${\mathbf{I}}_{n}$ denote the identity matrix of dimension $n\times n$ (where the subscript may be omitted if the dimension is clear from the context) and $\mathbf{0}$ denote a zero matrix of appropriate dimensions. The Kronecker product of two matrices ${\mathbf{A}}\in\mathbb{R}^{m\times n}$ and ${\mathbf{B}}\in\mathbb{R}^{p\times q}$ is denoted by ${\mathbf{A}}\otimes{\mathbf{B}}\in\mathbb{R}^{mp\times nq}$ . The Kronecker sum of two square matrices ${\mathbf{A}}\in\mathbb{R}^{n\times n}$ and ${\mathbf{B}}\in\mathbb{R}^{m\times m}$ is denoted by ${\mathbf{A}}\oplus{\mathbf{B}}\in\mathbb{R}^{nm\times nm}$ and defined as

{\mathbf{A}}\oplus{\mathbf{B}}=({\mathbf{A}}\otimes{\mathbf{I}}_{m})+({\mathbf{I}}_{n}\otimes{\mathbf{B}}).

To facilitate a compact representation of high-order polynomial and multinomial terms, we utilize a shorthand notation for repeated Kronecker operations. For a square matrix ${\mathbf{X}}\in\mathbb{R}^{n\times n}$ , its $k$ -th Kronecker power ${\mathbf{X}}^{(k)}$ and its $k$ -th Kronecker sum ${\mathbf{X}}\oplus_{k}{\mathbf{X}}$ are defined inductively, for $k\geq 2$ , as

	$\displaystyle{\mathbf{X}}^{(k)}$	$\displaystyle={\mathbf{X}}^{(k-1)}\otimes{\mathbf{X}},\quad\text{with }{\mathbf{X}}^{(1)}={\mathbf{X}},$
	$\displaystyle{\mathbf{X}}\oplus_{k}{\mathbf{X}}$	$\displaystyle=({\mathbf{X}}\oplus_{k-1}{\mathbf{X}})\otimes{\mathbf{I}}_{n}+{\mathbf{I}}_{n}^{(k-1)}\otimes{\mathbf{X}},\quad\text{with }{\mathbf{X}}\oplus_{1}{\mathbf{X}}={\mathbf{X}}.$

Equivalently, the $k$ -th Kronecker sum can be expressed explicitly as

{\mathbf{X}}\oplus_{k}{\mathbf{X}}=\sum_{j=1}^{k}{\mathbf{I}}_{n}^{(j-1)}\otimes{\mathbf{X}}\otimes{\mathbf{I}}_{n}^{(k-j)},

(1)

where ${\mathbf{I}}_{n}^{(0)}=1$ . Under this notation, a quadratic state interaction term for a vector $\mathbf{w}$ simplifies directly to $\mathbf{w}^{(2)}=\mathbf{w}\otimes\mathbf{w}$ . Additionally, for matrices of compatible dimensions, the well-known mixed-product property holds

({\mathbf{A}}{\mathbf{B}})\otimes({\mathbf{C}}{\mathbf{D}})=({\mathbf{A}}\otimes{\mathbf{C}})({\mathbf{B}}\otimes{\mathbf{D}}).

(2)

These algebraic properties are used frequently to simplify expressions throughout this paper.

2.2 Nonlinear moments and the invariance equation

In the LTI setting, traditional projection-based interpolatory model reduction methods deal with linear moments, which are discussed in Sec. 2.2.1. The notion of moment of a dynamical system was reformulated using the theory of steady-state responses and invariant manifolds in astolfi2010 leading to a unified framework to both redefine the notion of linear moments for linear and nonlinear systems and introduce the notion of nonlinear moments for both linear and nonlinear systems. To formalize this approach, consider a nonlinear full-order system (FOM) described by

	$\displaystyle\dot{\mathbf{x}}(t)$	$\displaystyle=f(\mathbf{x}(t),\mathbf{u}(t)),\quad\mathbf{x}(0)=\mathbf{x}_{0}\in\mathbb{R}^{n}$		(3)
	$\displaystyle\mathbf{y}(t)$	$\displaystyle=h(\mathbf{x}(t)),$		(3)

where $\mathbf{x}(t)\in\mathbb{R}^{n}$ represents the state, $\mathbf{u}(t)\in\mathbb{R}^{m}$ represents the input vector, $\mathbf{y}(t)\in\mathbb{R}^{p}$ represents the output. The map $f:\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$ describes the nonlinear state evolution and $h:\mathbb{R}^{n}\rightarrow\mathbb{R}^{p}$ represents the output map for the system. Also consider an autonomous (nonlinear) signal generator defined by

	$\displaystyle\dot{\omega}(t)$	$\displaystyle=\mathscr{s}(\omega(t)),\quad\omega(0)=\omega_{0}\in\mathbb{R}^{r}$		(4)
	$\displaystyle\mathbf{y}_{d}(t)$	$\displaystyle=\mathscr{l}(\omega(t)),$		(4)

where $\omega(t)\in\mathbb{R}^{r}$ represents the state of the signal generator space, $\mathscr{s}:\mathbb{R}^{r}\rightarrow\mathbb{R}^{r}$ is a nonlinear function that describes the time-evolution of the signal space and $\mathscr{l}:\mathbb{R}^{r}\rightarrow\mathbb{R}^{m}$ is a nonlinear map from the signal generator state, $\omega(t)$ , to the signal generator output, $\mathbf{y}_{d}(t)\in\mathbb{R}^{m}$ . Suppose the nonlinear system in (3) is driven by an input from the autonomous, nonlinear signal generator in (4), i.e., when $\mathbf{u}(t)=\mathbf{y}_{d}(t)$ . Then the combined dynamics of the two systems (3) and (4) can be represented by the interconnected system

$\displaystyle\dot{\omega}(t)$	$\displaystyle=\mathscr{s}(\omega(t)),~$	$\displaystyle\omega(0)=\omega_{0}\in\mathbb{R}^{r}$	(5)
$\displaystyle\dot{\mathbf{x}}(t)$	$\displaystyle=f(\mathbf{x}(t),\mathscr{l}(\omega(t))),~$	$\displaystyle\mathbf{x}(0)=\mathbf{x}_{0}\in\mathbb{R}^{n}$
$\displaystyle\ \mathbf{y}(t)$	$\displaystyle=h(\mathbf{x}(t)).$

Under appropriate assumptions, specifically that the unforced system $\dot{\mathbf{x}}=f(\mathbf{x},0)$ has a locally exponentially stable equilibrium at the origin and the signal generator possesses neutrally stable dynamics (see astolfi2010), the interconnected system in (5) has a locally invariant (center) manifold described by $\mathbf{x}=\mathbf{\Pi}(\omega(t))$ . For more details on center manifold theory and its applications see carr2012applications.

Definition 2.1 (Nonlinear Invariance Equation astolfi2010).

The mapping $\mathbf{\Pi}(\cdot)$ , which parameterizes the invariant manifold associated with $(\mathscr{s},\mathscr{l})$ is the unique local solution to the partial differential equation (PDE)

\frac{\partial\mathbf{\Pi}(\omega)}{\partial\omega}\mathscr{s}(\omega)=f(\mathbf{\Pi}(\omega),l(\omega)),\quad\mathbf{\Pi}(0)=0.

(6)

The PDE (6) is commonly referred to as the nonlinear invariance equation.

Definition 2.2 (Nonlinear Moments astolfi2010).

The nonlinear moment of the system associated with $(\mathscr{s},\mathscr{l})$ is defined as the composite mapping $\mathcal{M}_{\mathscr{s},\mathscr{l}}:\mathbb{R}^{r}\to\mathbb{R}^{p}$ given by the output of the system restricted to the invariant manifold:

\mathcal{M}_{\mathscr{s},\mathscr{l}}(\omega)=h\circ\mathbf{\Pi}(\omega).

(7)

We note that the nonlinear moment and center manifold are associated with a signal generator $(\mathscr{s},\mathscr{l})$ , i.e., the center manifold and nonlinear moment, are different for a different choice of signal generator.

The nonlinear moment matching framework is illustrated in Fig. 2, where the nonlinear moment is a mapping from the signal generator to the steady-state output of the system. The left panel defines the signal space, where the signal generator is governed by the dynamics $\dot{\omega}=\mathscr{s}(\omega)$ and generates the input $\mathbf{u}=\mathscr{l}(\omega)$ . The central panel depicts the center manifold $\mathbf{\Pi}(\omega)$ , which is an invariant geometric surface defined as the solution of the nonlinear invariance equation (6). When the system is driven by this signal generator, an arbitrary state trajectory $\mathbf{x}(t)$ (shown in blue) exhibits initial transient dynamics before eventually converging onto this manifold. The trajectory strictly restricted to this manifold represents the exact steady-state behavior of the system, denoted by $\mathbf{x}_{ss}(t)=\mathbf{\Pi}(\omega(t))$ (shown in black). Finally, the right panel illustrates the output converging to its steady state as the transient decays, matching the steady-state output generated by passing the manifold dynamics through the output function $h(\cdot)$ . The overarching nonlinear moment, defined as $\mathcal{M}_{s,l}(\omega)=h\circ\mathbf{\Pi}(\omega)$ , fundamentally captures this relationship, directly linking the signal space to the steady-state output while entirely bypassing the transient phase.

Figure 2: Schematic representation of the nonlinear moment framework.

2.2.1 Formulations of moments for linear systems

In this section, we consider the special case when $f(\mathbf{x},\mathbf{u})={\mathbf{A}}\mathbf{x}+{\mathbf{B}}\mathbf{u}$ and $g(x)={\mathbf{C}}\mathbf{x}$ where the state-space matrices ${\mathbf{A}}\in\mathbb{R}^{n\times n},{\mathbf{B}}\in\mathbb{R}^{n\times m}$ and ${\mathbf{C}}\in\mathbb{R}^{p\times n}$ are constant. Under this choice the full-order model (FOM) is described by

	$\displaystyle\dot{\mathbf{x}}(t)$	$\displaystyle={\mathbf{A}}\mathbf{x}(t)+{\mathbf{B}}\mathbf{u}(t),\quad\mathbf{x}(0)=\mathbf{x}_{0}$		(8)
	$\displaystyle\mathbf{y}(t)$	$\displaystyle={\mathbf{C}}\mathbf{x}(t),$		(8)

where $\mathbf{x}(t)\in\mathbb{R}^{n}$ represents the state, $\mathbf{u}(t)\in\mathbb{R}^{m}$ represents the input, and $\mathbf{y}(t)\in\mathbb{R}^{p}$ represents the output. Assuming zero initial conditions ( $\mathbf{x}(0)=\mathbf{0}$ ), taking the Laplace transform converts the differential equations (8) into an algebraic equation that reads ${\mathbf{Y}}(s)={\mathbf{C}}(s\mathbf{I}-{\mathbf{A}})^{-1}{\mathbf{B}}{\mathbf{U}}(s)$ where ${\mathbf{Y}}(s)$ and ${\mathbf{U}}(s)$ denote the Laplace transforms of the output $\mathbf{y}(t)$ and input $\mathbf{u}(t)$ , respectively. The function ${\mathbf{H}}(s)={\mathbf{C}}(s\mathbf{I}-{\mathbf{A}})^{-1}{\mathbf{B}}\in\mathbb{R}^{p\times m}$ denotes the associated transfer function of the LTI system and represents the input-output behaviour of the system in the frequency domain. In this LTI setting, the transfer function is a rational function of $s$ .

Definition 2.3 (Linear Moments ACA05).

Let ${\mathbf{H}}(s)$ denote the transfer function of the linear time-invariant (LTI) system in (8). The $0$ -th linear moment at a frequency $s^{*}\in\mathbb{C}$ is defined as the transfer function evaluated at that point, denoted by $\eta_{0}(s^{*})={\mathbf{H}}(s^{*})$ . For any integer $k\geq 1$ , the $k$ -th linear moment at $s^{*}$ is defined as the $k$ -th derivative of the transfer function evaluated at $s^{*}$ , given by

\eta_{k}(s^{*})=(-1)^{k}\left.\frac{d^{k}}{ds^{k}}{\mathbf{H}}(s)\right|_{s=s^{*}}.

Linear moment matching methods construct reduced systems, i.e., systems of the form (8) with a reduced state dimension, whose now lower-degree rational transfer function interpolates the moments of the original linear system as defined in Definition 2.3. This frequency-domain formulation is standard in interpolatory projection-based model reduction framework morAntBG20. Astolfi astolfi2010 demonstrated that these linear moments can be equivalently characterized in the time domain as the steady-state output response generated by an autonomous, linear signal generator of dimension $r$ ,

	$\displaystyle\dot{\omega}(t)$	$\displaystyle={\mathbf{S}}\omega(t),\quad\omega(0)=\omega_{0}\in\mathbb{R}^{r},$		(9)
	$\displaystyle\mathbf{u}(t)$	$\displaystyle={\mathbf{L}}\omega(t).$		(9)

where ${\mathbf{S}}\in\mathbb{R}^{r\times r}$ , ${\mathbf{L}}\in\mathbb{R}^{m\times r}$ are matrices chosen such that the pair ${\mathbf{S}}$ and ${\mathbf{L}}$ are locally observable, i.e., when

\displaystyle\operatorname{rank}\left(\left[\begin{array}[]{c}{\mathbf{L}}\\ {\mathbf{L}}{\mathbf{S}}\\ \vdots\\ {\mathbf{L}}{\mathbf{S}}^{r-1}\end{array}\right]\right)=r.

Intuitively, this condition ensures that the input $\mathbf{u}(t)$ depends on every component of $\omega(t)$ with no hidden or redundant components in the signal state. Similar to (5), we consider the interconnection of the LTI system in (8) and linear signal generator in (9), which leads to the interconnected system

$\displaystyle\dot{\omega}(t)$	$\displaystyle={\mathbf{S}}\omega(t),~$	$\displaystyle\omega(0)=\omega_{0}\in\mathbb{R}^{r}$	(10)
$\displaystyle\dot{\mathbf{x}}(t)$	$\displaystyle={\mathbf{A}}\mathbf{x}(t)+{\mathbf{B}}\mathscr{{\mathbf{L}}}\omega(t),~$	$\displaystyle\mathbf{x}(0)=\mathbf{x}_{0}\in\mathbb{R}^{n}$
$\displaystyle\ \mathbf{y}(t)$	$\displaystyle={\mathbf{C}}\mathbf{x}(t).$

In this linear setting, the invariance equation (6) reduces to a linear matrix equation (11). The center manifold $\mathbf{\Pi}(\omega)=\mathbf{\Pi}\omega$ where $\mathbf{\Pi}\in\mathbb{R}^{n\times r}$ is a linear subspace obtained as the solution of the Sylvester equation

{\mathbf{A}}\mathbf{\Pi}+{\mathbf{B}}{\mathbf{L}}=\mathbf{\Pi}{\mathbf{S}}.

(11)

Moreover, the moment $\mathcal{M}_{{\mathbf{S}},{\mathbf{L}}}(\cdot)$ associated with $({\mathbf{S}},{\mathbf{L}})$ in the sense of astolfi2010 is given as

\displaystyle\mathcal{M}_{{\mathbf{S}},{\mathbf{L}}}(\omega)={\mathbf{C}}\mathbf{\Pi}(\omega)={\mathbf{C}}\mathbf{\Pi}\omega.

To see how this fits into the classical projection-based framework in morAntBG20, consider the case where we want to interpolate the transfer function ${\mathbf{H}}(s)$ at $r$ distinct complex points $\{\mu_{1},\mu_{2},\dots,\mu_{r}\}$ along the tangential directions $\{\ell_{1},\dots,\ell_{r}\}$ . In the interpolatory projection-based framework, this is referred to as one-sided tangential interpolation where we construct a projection matrix $\mathbf{V}\in\mathbb{R}^{n\times r}$ whose columns span the Krylov subspace,

\text{Im}(\mathbf{V})=\text{span}\left\{(\mu_{1}\mathbf{I}-{\mathbf{A}})^{-1}{\mathbf{B}}\ell_{1},\dots,(\mu_{r}\mathbf{I}-{\mathbf{A}})^{-1}{\mathbf{B}}\ell_{r}\right\}.

(12)

If we choose the signal generator matrix ${\mathbf{S}}$ to be diagonal, ${\mathbf{S}}=\text{diag}(\mu_{1},\dots,\mu_{r})$ , and set ${\mathbf{L}}=[\ell_{1},\dots,\ell_{r}]$ such that the columns correspond to tangential directions in $\mathbb{R}^{m}$ , the Sylvester equation (11) can be solved column-by-column. For the $i$ -th column $\mathbf{\Pi}_{i}$ , the equation yields:

{\mathbf{A}}\mathbf{\Pi}_{i}+{\mathbf{B}}\ell_{i}=\mu_{i}\mathbf{\Pi}_{i}\implies\mathbf{\Pi}_{i}=(\mu_{i}\mathbf{I}-{\mathbf{A}})^{-1}{\mathbf{B}}\ell_{i}.

(13)

Thus, the columns of the steady-state mapping matrix $\mathbf{\Pi}$ are precisely the exact vectors that form the interpolating subspace in the projection-interpolatory framework, $\boldsymbol{\Pi}={\mathbf{V}}$ . Thus, the moment $\mathcal{M}_{{\mathbf{S}},{\mathbf{L}}}(\cdot)$ reduces to

\mathcal{M}_{{\mathbf{S}},{\mathbf{L}}}(\omega)={\mathbf{C}}\mathbf{\Pi}\omega=[{\mathbf{H}}(\mu_{1})\ell_{1}\quad{\mathbf{H}}(\mu_{2})\ell_{2}\quad\dots\quad{\mathbf{H}}(\mu_{r})\ell_{r}]\omega

(14)

From (14), we see that the linear moment (and the steady state response) for an LTI system is determined by transfer function evaluations at the frequencies $\mu_{1},\dots,\mu_{r}$ and the nonlinear moment matching framework in astolfi2010 boils down to the classical interpolatory projection framework morAntBG20.

Figure 3: Schematic representation of linear moment for LTI systems.

When both the system and the signal generator are linear, the nonlinear moment-matching framework in Fig. 2 boils down to the linear framework depicted in Fig. 3. In this setting, the nonlinear invariance equation (2.1) reduces to the linear Sylvester equation (11), causing the center manifold to simplify to an $r$ -dimensional linear subspace. Due to the invariance property of this subspace, trajectories originating on it remain confined to it for all time. State trajectories starting off on the subspace (such as $x(t)$ , shown in blue) undergo an initial transient before decaying onto the subspace to yield the steady-state state trajectory $\mathbf{x}_{ss}(t)={\mathbf{C}}\mathbf{\Pi}\omega(t)$ . Consequently, the system output exhibits a corresponding transient phase before converging to its steady-state response. The overarching moment mapping $\mathcal{M}_{{\mathbf{S}},{\mathbf{L}}}(\omega)={\mathbf{C}}\mathbf{\Pi}\omega$ directly links the signal space to this steady-state output, which, for appropriate choices of $({\mathbf{S}},{\mathbf{L}})$ , simplifies to classical transfer function interpolation as shown in (14).

Remark 2.4.

Thus, in the linear setting with the above choice of signal generator, the matrices ${\mathbf{S}}$ and ${\mathbf{L}}$ in the nonlinear moment matching framework act as the direct analogs to the interpolation points $\mu_{i}$ and tangential directions $\ell_{i}$ in the classical projection-interpolatory framework in morAntBG20. Because the steady-state dynamics restricted to the manifold $\mathbf{\Pi}$ are driven purely by the signal generator (i.e., $\dot{\mathbf{x}}=\mathbf{\Pi}{\mathbf{S}}\omega$ ), a reduced-order model constructed with these matrices intrinsically achieves the desired one-sided tangential interpolation. Consequently, throughout this paper, we treat the signal generator not merely as a fixed input structure, but as tunable design parameters that can be explicitly chosen to enforce targeted interpolation requirements and steady-state fidelity in the reduced system.

2.2.2 Moment-matching (interpolatory) reduced systems

This section reviews the formal definition of moment matching for nonlinear systems and establishes the general structure of an interpolatory reduced-order model (ROM). Consider the full-order system described in (3), and let a candidate reduced system be given by

	$\displaystyle\dot{\hat{\mathbf{x}}}(t)$	$\displaystyle=\hat{f}(\hat{\mathbf{x}}(t),\mathbf{u}(t)),\quad\hat{\mathbf{x}}(0)=\hat{\mathbf{x}}_{0},$		(15)
	$\displaystyle\mathbf{y}_{r}(t)$	$\displaystyle=\hat{h}(\hat{\mathbf{x}}(t)),$		(15)

where $\hat{\mathbf{x}}(t)\in\mathbb{R}^{r}$ represents the reduced state vector with $r<n$ , $\mathbf{u}(t)\in\mathbb{R}^{m}$ is the input, and $\mathbf{y}_{r}(t)\in\mathbb{R}^{p}$ is the ROM output. The map $\hat{f}:\mathbb{R}^{r}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{r}$ governs the reduced state evolution, and $\hat{h}:\mathbb{R}^{r}\rightarrow\mathbb{R}^{p}$ defines the output map. Suppose the system in (15) is driven by the signal generator $(\mathscr{s},\mathscr{l})$ , yielding the interconnected system:

$\displaystyle\dot{\omega}(t)$	$\displaystyle=\mathscr{s}(\omega(t)),\quad$	$\displaystyle\omega(0)$	$\displaystyle=\omega_{0},$	(16)
$\displaystyle\dot{\hat{\mathbf{x}}}(t)$	$\displaystyle=\hat{f}(\hat{\mathbf{x}}(t),\mathscr{l}(\omega(t))),\quad$	$\displaystyle\hat{\mathbf{x}}(0)$	$\displaystyle=\hat{\mathbf{x}}_{0},$
$\displaystyle\mathbf{y}_{r}(t)$	$\displaystyle=\hat{h}(\hat{\mathbf{x}}(t)).$

This interconnected system possesses its own center manifold, denoted by $\boldsymbol{\pi}_{r}(\omega)$ , which is obtained by solving the corresponding invariance equation in Definition 2.1. The nonlinear moment of the reduced system at $(\mathscr{s},\mathscr{l})$ is then given as

\widehat{\mathcal{M}}_{\mathscr{s},\mathscr{l}}(\omega)=\hat{h}\circ\boldsymbol{\pi}_{r}(\omega).

(17)

Definition 2.5 (Nonlinear moment matching astolfi2010).

The system in (15) achieves nonlinear moment matching of the full-order system (3) at $(\mathscr{s},\mathscr{l})$ if, for all $\omega$ in a neighborhood of the origin, the moments coincide, i.e.,

\widehat{\mathcal{M}}_{\mathscr{s},\mathscr{l}}(\omega)=\hat{h}\circ\boldsymbol{\pi}_{r}(\omega)=h\circ\boldsymbol{\pi}(\omega)=\mathcal{M}_{\mathscr{s},\mathscr{l}}(\omega).

(18)

Additionally, the system in (15) is formally considered a reduced-order model of (3) if $r<n$ .

In astolfi2010, the author introduce a parameterized family of reduced models that inherently match the moments of (3) at $(\mathscr{s},\mathscr{l})$ . This family is defined by

	$\displaystyle\dot{\mathbf{x}}_{r}(t)$	$\displaystyle=\mathscr{s}(\mathbf{x}_{r}(t))-\delta(\mathbf{x}_{r}(t))\mathscr{l}(\mathbf{x}_{r}(t))+\delta(\mathbf{x}_{r}(t))\mathbf{u}(t),$		(19)
	$\displaystyle\mathbf{y}_{r}(t)$	$\displaystyle=h(\boldsymbol{\pi}(\mathbf{x}_{r}(t))),$		(19)

where $\mathbf{x}_{r}\in\mathbb{R}^{r}$ is the reduced state vector, and $\delta(\cdot)$ is an arbitrary mapping chosen such that the associated invariance equation,

\frac{\partial p(\omega)}{\partial\omega}\mathscr{s}(\omega)=\mathscr{s}(p(\omega))-\delta(p(\omega))\mathscr{l}(p(\omega))+\delta(p(\omega))\mathscr{l}(\omega),

(20)

admits the unique solution $p(\omega)=\omega$ . Satisfying this condition guarantees that the center manifold of the reduced system collapses to the identity map ( $\boldsymbol{\pi}_{r}(\omega)=\omega$ ). Consequently, the reduced model evaluates the exact full-order output map along the FOM’s center manifold, ensuring that the ROM perfectly replicates the steady-state output response of the high-dimensional plant for the prescribed class of inputs.

Remark 2.6.

For a given FOM and chosen signal generator, there exist an infinite number of reduced systems that satisfy the moment-matching condition. Furthermore, the FOM and ROM are not required to share the same structural form; a linear FOM may be reduced to a nonlinear ROM, and vice versa. The fundamental requirement is simply that the ROM dimension is strictly smaller than the FOM dimension ( $r\ll n$ ) while fulfilling the exact moment-matching property in (18). In this paper, we present a method of constructing a quadratic ROM that interpolates the nonlinear moments of a linear (LTI) full order system.

3 Defining the quadratic reduced system

We consider the problem of constructing reduced systems for linear time-invariant (LTI) systems of the form

	$\displaystyle\dot{\mathbf{x}}(t)$	$\displaystyle={\mathbf{A}}\mathbf{x}(t)+{\mathbf{B}}\mathbf{u}(t),\quad\mathbf{x}(0)=\mathbf{x}_{0}\in\mathbb{R}^{n}$		(21)
	$\displaystyle\mathbf{y}(t)$	$\displaystyle={\mathbf{C}}\mathbf{x}(t),$		(21)

where $\mathbf{x}(t)\in\mathbb{R}^{n}$ represents the state, $\mathbf{u}(t)\in\mathbb{R}^{m}$ represents the input and $\mathbf{y}(t)\in\mathbb{R}^{p}$ represents the output and the state matrices are ${\mathbf{A}}\in\mathbb{R}^{n\times n},{\mathbf{B}}\in\mathbb{R}^{n\times m}$ and ${\mathbf{C}}\in\mathbb{R}^{p\times n}$ . In this paper, we focus on quadratic approximations, where the state $\mathbf{x}(t)\in\mathbb{R}^{n}$ is approximated as

\mathbf{x}(t)\approx\boldsymbol{\Phi}(\mathbf{x}_{r}(t))=\mathbf{V}_{1}\mathbf{x}_{r}(t)+\mathbf{V}_{2}(\mathbf{x}_{r}(t)\otimes\mathbf{x}_{r}(t)),

(22)

where $\mathbf{x}_{r}(t)\in\mathbb{R}^{r}$ is the reduced state with $r\ll n$ . The matrices $\mathbf{V}_{1}\in\mathbb{R}^{n\times r}$ and $\mathbf{V}_{2}\in\mathbb{R}^{n\times r^{2}}$ define the linear and quadratic components of the manifold, respectively. Substituting the quadratic approximation (22) into the linear FOM (21) yields the residual equation

\displaystyle\mathbf{V}_{1}\frac{d\mathbf{x}_{r}}{dt}+\mathbf{V}_{2}\frac{d}{dt}(\mathbf{x}_{r}(t)\otimes\mathbf{x}_{r}(t))=\mathbf{A}\mathbf{V}_{1}\mathbf{x}_{r}(t)+\mathbf{A}\mathbf{V}_{2}(\mathbf{x}_{r}(t)\otimes\mathbf{x}_{r}(t))+\mathbf{B}\mathbf{u}(t)+\mathbf{r}(t),

where $\mathbf{r}(t)\in\mathbb{R}^{n}$ is the residual vector. Within the framework of classical interpolatory MOR, we introduce a test projection matrix ${\mathbf{W}}\in\mathbb{R}^{n\times r}$ such that ${\mathbf{W}}^{\top}{\mathbf{V}}_{1}$ is invertible and $\mathbf{W}^{\top}\mathbf{r}(t)=\mathbf{0}$ . In other words, we solve the FOM along the columns of ${\mathbf{W}}$ . Thus, projecting the residual equation with ${\mathbf{W}}^{\top}$ from the left yields

\displaystyle\mathbf{W}^{\top}\mathbf{V}_{1}\frac{d\mathbf{x}_{r}}{dt}+\mathbf{W}^{\top}\mathbf{V}_{2}\frac{d}{dt}(\mathbf{x}_{r}\otimes\mathbf{x}_{r})

\displaystyle=\mathbf{W}^{\top}\mathbf{A}\mathbf{V}_{1}\mathbf{x}_{r}+\mathbf{W}^{\top}\mathbf{A}\mathbf{V}_{2}(\mathbf{x}_{r}\otimes\mathbf{x}_{r})+\mathbf{W}^{\top}\mathbf{B}\mathbf{u}.

This formulation results in a reduced system given by


$\displaystyle({\mathbf{I}}_{r}+{\mathbf{E}}_{r}(\mathbf{x}_{r}\oplus\mathbf{x}_{r}))\frac{d\mathbf{x}_{r}(t)}{dt}$	$\displaystyle=\mathbf{A}_{r}\mathbf{x}_{r}(t)+\mathbf{H}_{r}(\mathbf{x}_{r}(t)\otimes\mathbf{x}_{r}(t))+\mathbf{B}_{r}\mathbf{u}(t)$	(23a)
$\displaystyle\mathbf{y}_{r}(t)$	$\displaystyle=\widehat{h}(\mathbf{x}_{r})=\mathbf{C}_{r}\mathbf{x}_{r}(t)+\mathbf{K}_{r}(\mathbf{x}_{r}(t)\otimes\mathbf{x}_{r}(t)),$	(23b)

where the reduced matrices and the state-dependent reduced mass matrix term denoted by $\mathbf{M}_{r}(\mathbf{x}_{r})$ are defined as


$\displaystyle\mathbf{A}_{r}$	$\displaystyle=(\mathbf{W}^{\top}\mathbf{V}_{1})^{-1}\mathbf{W}^{\top}\mathbf{A}\mathbf{V}_{1},\quad\mathbf{B}_{r}=(\mathbf{W}^{\top}\mathbf{V}_{1})^{-1}\mathbf{W}^{\top}\mathbf{B},\quad\mathbf{C}_{r}=\mathbf{C}\mathbf{V}_{1}\quad\mathbf{H}_{r}=(\mathbf{W}^{\top}\mathbf{V}_{1})^{-1}\mathbf{W}^{\top}\mathbf{A}\mathbf{V}_{2},$	(24a)
$\displaystyle\mathbf{E}_{r}$	$\displaystyle=(\mathbf{W}^{\top}\mathbf{V}_{1})^{-1}\mathbf{W}^{\top}\mathbf{V}_{2},\qquad\mathbf{M}_{r}(\mathbf{x}_{r})={\mathbf{I}}_{r}+{\mathbf{E}}_{r}(\mathbf{I}_{r}\otimes\mathbf{x}_{r}+\mathbf{x}_{r}\otimes\mathbf{I}_{r})\quad\mathbf{K}_{r}=\mathbf{C}\mathbf{V}_{2}.$	(24b)

In the proposed framework of this paper, the choice of the left projection basis ${\mathbf{W}}$ is uncoupled and independent of the choice of ${\mathbf{V}}_{1},{\mathbf{V}}_{2}$ . This flexibility allows one to tailor the algebraic properties of the reduced-order model to specific computational requirements. For example, by enforcing the orthogonality condition, $\mathbf{W}^{\top}\mathbf{V}_{2}=\mathbf{0}$ (like in geelen2023operator; QM_framework for example), $\mathbf{M}_{r}$ simplifies to the identity matrix ${\mathbf{I}}_{r}$ . This is particularly advantageous as it simplifies the reduced dynamics, ensuring that the resulting ROM avoids state dependent terms multiplying the derivative of the reduced state. In our framework, we consider the more general case of where ${\mathbf{W}}^{\top}{\mathbf{V}}_{2}\neq 0$ , while naturally retaining both the orthogonalized formulation and the classical Galerkin projection ${\mathbf{W}}={\mathbf{V}}_{1}$ as special cases. We note that the mass-matrix term like in (23) also appears in the context of nonlinear decoder approximation in benner2023quadratic.

In the next section, we provide a system-theoretic way of choosing these projection matrices ${\mathbf{V}}_{1},{\mathbf{V}}_{2}$ and ${\mathbf{W}}$ such that the ROM defined in (23)–(24) preserves the nonlinear moments of the FOM defined in (21).

4 Constructing interpolatory ROM using quadratic approximations

This section establishes the theoretical foundations for constructing an interpolating reduced-order model (ROM) that achieves moment matching over a prescribed signal space. We begin by characterizing the center manifold obtained as the solution to the invariance equation under a specialized class of signal generators. Specifically, Theorem 4.1 demonstrates that this center manifold is quadratic, providing theoretical motivation for the choice of the quadratic projection matrices ${\mathbf{V}}_{1}$ and ${\mathbf{V}}_{2}$ . These developments culminate in Theorem 4.5, the primary theoretical result of this work, which states that the ROM defined in (23) and (24) matches the moments of the full-order system in the sense of astolfi2010, given an appropriate selection of ${\mathbf{V}}_{1}$ , ${\mathbf{V}}_{2}$ , and ${\mathbf{W}}$ . To facilitate readability, the proof of this main result is structured into two sequential steps. First, Theorem 4.3 establishes that the ROM constructed with these specific projection matrices guarantees that its own center manifold reduces to the identity map. Second, the overarching moment-matching result is achieved by combining the quadratic center manifold property from Theorem 4.1 with the identity mapping property from Theorem 4.3.

4.1 Quadratic Center Manifolds for Model Reduction

We start by introducing an autonomous signal generator whose state $\omega(t)\in\mathbb{R}^{r}$ evolves linearly, while its output $\mathbf{u}(t)$ incorporates both linear and quadratic terms. Let $\omega_{0}=\omega(0)\in\mathbb{R}^{r}$ denote the initial state of the generator. Suppose the signal generator dynamics are chosen as

	$\displaystyle\dot{\omega}(t)$	$\displaystyle=\mathbf{S}\omega(t),\quad\quad\omega(0)=\omega_{0},$		(25)
	$\displaystyle\mathbf{u}(t)$	$\displaystyle=\mathbf{L}(\omega(t))=\mathbf{L}_{1}\omega(t)+\mathbf{L}_{2}\omega^{(2)}(t),$		(25)

where $\mathbf{S}\in\mathbb{R}^{r\times r}$ is a full-rank matrix, $\mathbf{L}_{1}\in\mathbb{R}^{m\times r}$ and $\mathbf{L}_{2}\in\mathbb{R}^{m\times r^{2}}$ . Consequently, the explicit, time-dependent input that will be fed into the FOM system (21) is given by

\mathbf{u}(t)=\mathbf{L}_{1}e^{\mathbf{S}t}\omega_{0}+\mathbf{L}_{2}\left(e^{\mathbf{S}t}\otimes e^{\mathbf{S}t}\right)(\omega_{0}\otimes\omega_{0}).

(26)

To provide intuition for the upcoming formulation, we note that selecting the signal generator parameters ${\mathbf{S}}$ , ${\mathbf{L}}_{1}$ , and ${\mathbf{L}}_{2}$ in (25) is functionally equivalent to choosing the interpolation frequencies and tangential directions, much like discussed in Section 2.2.1. Consequently, the quadratic ROM constructed later in this section explicitly interpolates the exact nonlinear moments of the full-order system at these targeted frequencies and directions. We rigorously formalize this connection between choice of signal generators and nonlinear moment interpolation and provide canonical choices for selecting these parameters in Section 5. Assuming an interconnection of this signal space (25) with the LTI system from (21), we obtain

$\displaystyle\dot{}\omega(t)$	$\displaystyle={\mathbf{S}}\omega(t),\quad\omega(0)=\omega_{0},$	(27)
$\displaystyle\dot{\mathbf{x}}(t)$	$\displaystyle={\mathbf{A}}\mathbf{x}(t)+{\mathbf{B}}({\mathbf{L}}_{1}\omega(t)+{\mathbf{L}}_{2}(\omega(t)\otimes\omega(t))),\quad\mathbf{x}(0)=\mathbf{x}_{0}\in\mathbb{R}^{n},$
$\displaystyle\mathbf{y}(t)$	$\displaystyle={\mathbf{C}}\mathbf{x}(t).$

This interconnected system admits a center invariant manifold which, as discussed in Section 2, is given as the solution of the invariance equation. For this interconnected system, the invariance equation has the form

\frac{\partial\Pi(\omega)}{\partial\omega}\mathbf{S}\omega=\mathbf{A}\Pi(\omega)+\mathbf{B}\left(\mathbf{L}_{1}\omega+\mathbf{L}_{2}\omega^{(2)}\right).

(28)

In the following theorem, we show that for a signal generator of the form (25), the solution of the invariance equation (28) is a quadratic manifold. The matrices that define this quadratic center manifold mapping are a natural choice for the projection matrices $\mathbf{V}_{1}$ and $\mathbf{V}_{2}$ used to construct the ROM matrices in (24).

Theorem 4.1.

Suppose ${\mathbf{S}}\in\mathbb{R}^{r\times r}$ is chosen such that the eigenvalues of ${\mathbf{A}}$ are disjoint from those of ${\mathbf{S}}$ and its Kronecker sums ( ${\mathbf{S}}\oplus_{k}{\mathbf{S}}$ ) for $k\geq 2$ . Then the unique, analytic solution to the invariance equation (28) is given by

\displaystyle\mathbf{\Pi}(\omega)=\mathbf{\Pi}_{1}\omega+\mathbf{\Pi}_{2}\omega^{(2)},

where $\mathbf{\Pi}_{1}\in\mathbb{R}^{n\times r}$ and $\mathbf{\Pi}_{2}\in\mathbb{R}^{n\times r^{2}}$ are constant matrices given as the solution of the Sylvester equations

	$\displaystyle\boldsymbol{\Pi}_{1}\mathbf{S}-\mathbf{A}\boldsymbol{\Pi}_{1}$	$\displaystyle=\mathbf{B}\mathbf{L}_{1},$		(29)
	$\displaystyle\boldsymbol{\Pi}_{2}(\mathbf{S}\oplus\mathbf{S})-\mathbf{A}\boldsymbol{\Pi}_{2}$	$\displaystyle=\mathbf{B}\mathbf{L}_{2}.$		(30)

Proof.

The center manifold under the signal generator is given as the solution of the following invariance equation

\displaystyle\frac{\partial\mathbf{\Pi}(\omega)}{\partial\omega}={\mathbf{A}}\Pi(\omega)+{\mathbf{B}}({\mathbf{L}}_{1}\omega+{\mathbf{L}}_{2}\omega^{(2)}).

We follow an idea similar to the derivation of Volterra series expansion for structured dynamical systems in rugh1981nonlinear and the power series approach to solving (regulator) invariance equations in huang2004nonlinear; bai2022model. Suppose the center manifold $\boldsymbol{\mathbf{\Pi}}(\omega)$ of the FOM (21) under signal generated by (25) is given by the power series expansion (around $0$ ), i.e.,

\boldsymbol{\Pi}(\omega)=\sum_{i=1}^{\infty}\boldsymbol{\Pi}_{i}\omega^{(i)}=\boldsymbol{\Pi}_{1}\omega+\boldsymbol{\Pi}_{2}\omega^{(2)}+\boldsymbol{\Pi}_{3}\omega^{(3)}+\dots.

(31)

Substituting the power series expansion (31) into the invariance equation (28) yields

\sum_{i=1}^{\infty}\frac{\partial(\boldsymbol{\Pi}_{i}\omega^{(i)})}{\partial\omega}\mathbf{S}\omega=\mathbf{A}\sum_{i=1}^{\infty}\boldsymbol{\Pi}_{i}\omega^{(i)}+\mathbf{B}\mathbf{L}_{1}\omega+\mathbf{B}\mathbf{L}_{2}\omega^{(2)}.

(32)

By collecting the terms corresponding to each degree of $\omega^{(i)}$ , we derive a sequence of Sylvester equations

$\displaystyle\boldsymbol{\Pi}_{1}\mathbf{S}-\mathbf{A}\boldsymbol{\Pi}_{1}$	$\displaystyle=\mathbf{B}\mathbf{L}_{1},$	(33)
$\displaystyle\boldsymbol{\Pi}_{2}(\mathbf{S}\oplus\mathbf{S})-\mathbf{A}\boldsymbol{\Pi}_{2}$	$\displaystyle=\mathbf{B}\mathbf{L}_{2},$	(34)
$\displaystyle\boldsymbol{\Pi}_{k}(\mathbf{S}\oplus_{k}{\mathbf{S}})-\mathbf{A}\boldsymbol{\Pi}_{k}$	$\displaystyle=\mathbf{0},\quad\text{for }k\geq 3,$	(35)

where $\oplus_{k}$ denotes the $k$ Kronecker sum. From the theory of Sylvester equations (see, for example, ACA05), for $k\geq 3$ , $\boldsymbol{\Pi}_{k}=\mathbf{0}$ is the unique solution of the Sylvester equation if the eigenvalues of $\mathbf{A}$ are disjoint from the eigenvalues of the $k$ -fold Kronecker sum of $\mathbf{S}$ . Thus, by the hypothesis of the theorem, we have $\boldsymbol{\Pi}_{k}=\mathbf{0}$ as the unique solution to (35) for $k\geq 3$ . Thus, the center manifold for the FOM under the signal generator in (25) is given by

\mathbf{\Pi}(\omega)=\mathbf{\Pi}_{1}\omega+\mathbf{\Pi}_{2}\omega^{(2)},

where $\mathbf{\Pi}_{1}\in\mathbb{R}^{n\times r}$ and $\mathbf{\Pi}_{2}\in\mathbb{R}^{n\times r^{2}}$ are constant matrices given as the solution of the Sylvester equations (29) and (30) respectively. ∎

The statement of Theorem 4.1 provides a strong motivation for the use of quadratic manifolds to approximate the FOM state since the center manifold, obtained as the solution of the nonlinear-invariance equation in (28), is quadratic (see Fig. 4) in $\omega$ and hence can be reduced to a $r$ dimensional reduced space using a linear and quadratic projection matrix. Furthermore, the theorem also provides a natural choice for the linear and quadratic projection bases in the quadratic approximation:

\mathbf{V}_{1}:=\boldsymbol{\Pi}_{1}\in\mathbb{R}^{n\times r},\quad\mathbf{V}_{2}:=\boldsymbol{\Pi}_{2}\in\mathbb{R}^{n\times r^{2}}.

(36)

Because these matrices are obtained by directly solving linear Sylvester equations, our approach provides an optimization-free, closed form expressions for ${\mathbf{V}}_{1}$ and ${\mathbf{V}}_{2}$ . The QM ROM (23)–(24) constructed using ${\mathbf{V}}_{1}$ and ${\mathbf{V}}_{2}$ achieves nonlinear moment matching as shown in the following subsections (see Theorem 4.5). This construction also provides a significant practical advantage: it isolates the left projection matrix ${\mathbf{W}}$ as a completely free design parameter, constrained only by the mild requirement that ${\mathbf{W}}^{\top}{\mathbf{V}}_{1}$ remains invertible. By completely decoupling the moment-matching conditions from the choice of ${\mathbf{W}}$ , our framework allows flexibility to enforce secondary numerical or physical properties (such as stability or structure preservation) without affecting the moment matching property. This can be considered analogous to the one-sided vs two-sided interpolatory projections in the linear projection-based interpolatory MOR; see, e.g., (morAntBG20, Thm. 3.3.1). In the subsequent theorems, we make no additional assumptions on ${\mathbf{W}}$ to preserve this generality. While existing literature typically relies on standard choices, such as enforcing a Galerkin-like projection ( ${\mathbf{W}}={\mathbf{V}}_{1}$ ) or orthogonalizing against the quadratic basis ( ${\mathbf{W}}^{\top}{\mathbf{V}}_{2}=\mathbf{0}$ ), our framework automatically accommodates these choices of ${\mathbf{W}}$ .

Remark 4.2.

The assumption in Theorem 4.1 (which carries over to Theorems 4.3 and 4.5) that the eigenvalues of ${\mathbf{A}}$ are disjoint from those of ${\mathbf{S}}$ and its Kronecker sums ${\mathbf{S}}\oplus_{k}{\mathbf{S}}$ has two key implications. Mathematically, it guarantees unique solutions to the Sylvester equations, ensuring that the invariant quadratic center manifold is uniquely defined. From a system-theoretic perspective, it prevents placing interpolation frequencies at the poles of the transfer function of the original LTI system, thereby ensuring well-posed steady-state dynamics.

Figure 4: Schematic showing the quadratic center manifold of the LTI system (21). An arbitrary trajecteory

\mathbf{x}(t)

converges to its corresponding steady state

\mathbf{x}_{ss}(t)

trajectory on the center manifold.

4.2 Theoretical interpolation guarantees

Based on Theorem 4.1, we construct the reduced system (23)-(24) using projection matrices $\mathbf{V}_{1}$ and $\mathbf{V}_{2}$ defined in (36) and ${\mathbf{W}}$ chosen such that ${\mathbf{W}}^{\top}{\mathbf{V}}_{1}$ is invertible. In this section we provide the theoretical guarantees that the resulting reduced-order model (ROM) preserves the nonlinear moments of the full-order model (FOM). We first demonstrate that the center manifold of the ROM constructed using ${\mathbf{V}}_{1}$ and ${\mathbf{V}}_{2}$ as in (36) maps to the center manifold of the FOM, thereby ensuring nonlinear moment matching along ${\mathbf{W}}^{\top}$ . We denote the nonlinear moments of the system as $\boldsymbol{\mathcal{M}}(\omega)$ and suppress the dependence on ${\mathbf{S}},{\mathbf{L}}_{1},{\mathbf{L}}_{2}$ to avoid notational clutter. For the reader’s convenience, we recall the equations of the ROM (23) interconnected with the signal generator from (25) given by


	$\displaystyle\dot{\omega}(t)$	$\displaystyle={\mathbf{S}}\omega(t),\quad\omega(0)=\omega_{0}\in\mathbb{R}^{r}$
	$\displaystyle{\mathbf{M}}(\mathbf{x}_{r})\frac{d\mathbf{x}_{r}(t)}{dt}$	$\displaystyle=\mathbf{A}_{r}\mathbf{x}_{r}(t)+\mathbf{H}_{r}(\mathbf{x}_{r}(t)\otimes\mathbf{x}_{r}(t))+\mathbf{B}_{r}\left({\mathbf{L}}_{1}\omega(t)+{\mathbf{L}}_{2}(\omega(t)\otimes\omega(t))\right),\quad\mathbf{x}_{r}(0)\in\mathbb{R}^{r}$
	$\displaystyle\mathbf{y}_{r}(t)$	$\displaystyle=\widehat{h}(\mathbf{x}_{r})=\mathbf{C}_{r}\mathbf{x}_{r}(t)+\mathbf{K}_{r}(\mathbf{x}_{r}(t)\otimes\mathbf{x}_{r}(t)),$

with the matrices of the ROM constructed as in (24) with the choice ${\mathbf{V}}_{1}=\boldsymbol{\Pi}_{1}$ and ${\mathbf{V}}_{2}=\boldsymbol{\Pi}_{2}$ as shown in (36). The reduced-order model (ROM) naturally has a mass-matrix term ${\mathbf{M}}_{r}(\mathbf{x}_{r})$ multiplying the time derivative. To incorporate this into the standard nonlinear moment matching framework defined in (2.1), we recast the system into a form without the mass-matrix term. By our specific choice of the projection matrix $\mathbf{W}$ , the product $\mathbf{W}^{\top}\mathbf{V}_{1}$ is nonsingular. Consequently, the reduced mass matrix $\mathbf{M}_{r}(\mathbf{x}_{r})$ is invertible at the equilibrium point $\mathbf{x}_{r}=0$ . By continuity, there exists a neighborhood around the origin where $\mathbf{M}_{r}(\mathbf{x}_{r})$ remains invertible. This local invertibility allows us to formulate the invariance equation for the ROM as

	$\displaystyle\frac{\partial\boldsymbol{\pi}_{r}}{\partial\omega}\mathbf{S}\omega$	$\displaystyle=(\mathbf{M}_{r}(\boldsymbol{\pi}_{r}(\omega)))^{-1}\left(\mathbf{A}_{r}\boldsymbol{\pi}_{r}(\omega)+\mathbf{H}_{r}(\boldsymbol{\pi}_{r}(\omega)\otimes\boldsymbol{\pi}_{r}(\omega))+\mathbf{B}_{r}(\mathbf{L}_{1}\omega+\mathbf{L}_{2}\omega^{(2)})\right)$
	$\displaystyle\implies\mathbf{M}_{r}(\boldsymbol{\pi}_{r}(\omega))\frac{\partial\boldsymbol{\pi}_{r}}{\partial\omega}\mathbf{S}\omega$	$\displaystyle=\mathbf{A}_{r}\boldsymbol{\pi}_{r}(\omega)+\mathbf{H}_{r}(\boldsymbol{\pi}_{r}(\omega)\otimes\boldsymbol{\pi}_{r}(\omega))+\mathbf{B}_{r}(\mathbf{L}_{1}\omega+\mathbf{L}_{2}\omega^{(2)}).$		(38)

Thus, we adopt the implicit formulation in (38) as the defining invariance equation for the ROM in the subsequent results. Additionally using (7) and the definition of the reduced system in (23), we have that the moment of the ROM for ${\mathbf{S}},{\mathbf{L}}_{1}$ and ${\mathbf{L}}_{2}$ is given by

\widehat{\mathcal{M}}(\omega)=\widehat{h}\circ\boldsymbol{\pi}_{r}(\omega)={\mathbf{C}}_{r}\boldsymbol{\pi}_{r}(\omega)+{\mathbf{K}}_{r}(\boldsymbol{\pi}_{r}(\omega)\otimes\boldsymbol{\pi}_{r}(\omega)).

(39)

Theorem 4.3.

Consider the ROM defined in (23)–(24) with $\mathbf{V}_{1},\mathbf{V}_{2}$ chosen as the solution to (29)–(30) and let ${\mathbf{W}}$ be chosen such that $\mathbf{W}^{\top}\mathbf{V}_{1}$ is invertible. Let $\boldsymbol{\pi}_{r}(\omega)$ be the solution to the ROM invariance equation (38), i.e.,

{\mathbf{M}}(\boldsymbol{\pi}_{r}(\omega))\frac{\partial\boldsymbol{\pi}_{r}}{\partial\omega}\mathbf{S}\omega=\mathbf{A}_{r}\boldsymbol{\pi}_{r}(\omega)+\mathbf{H}_{r}(\boldsymbol{\pi}_{r}(\omega)\otimes\boldsymbol{\pi}_{r}(\omega))+\mathbf{B}_{r}(\mathbf{L}_{1}\omega+\mathbf{L}_{2}\omega^{(2)})

in a neighbourhood of $0$ . If the eigenvalues of $\mathbf{A}_{r}$ are disjoint from the eigenvalues of $\mathbf{S}$ and its Kronecker sums ( ${\mathbf{S}}\oplus_{k}{\mathbf{S}}$ ) for $k\geq 2$ , then the identity mapping

\boldsymbol{\pi}_{r}(\omega)=\omega

(40)

is the unique solution to the invariance equation (38) in a neighbourhood of $0$ .

Proof.

We follow a similar power series approach used in Theorem 4.1. Consider the invariance equation for the ROM given by

\displaystyle{\mathbf{M}}(\boldsymbol{\pi}_{r})\frac{\partial\boldsymbol{\pi}_{r}}{\partial\omega}\mathbf{S}\omega=\mathbf{A}_{r}\boldsymbol{\pi}_{r}(\omega)+\mathbf{H}_{r}(\boldsymbol{\pi}_{r}(\omega)\otimes\boldsymbol{\pi}_{r}(\omega))+\mathbf{B}_{r}(\mathbf{L}_{1}\omega+\mathbf{L}_{2}\omega^{(2)}).

Consider the power series expansion for $\boldsymbol{\pi_{r}}$ in a neighbourhood around 0 given as

\displaystyle\boldsymbol{\pi_{r}}(\omega)

\displaystyle={\mathbf{P}}_{1}\omega+{\mathbf{P}}_{2}\omega^{(2)}+{\mathbf{P}}_{3}\omega^{(3)}+\dots=\sum_{i=1}^{\infty}{\mathbf{P}}_{i}\omega^{(i)}.

We will show that under the conditions stated in the theorem, $\boldsymbol{\pi_{r}}(\omega)=\omega$ is the unique solution to the invariance equation, i.e., ${\mathbf{P}}_{1}={\mathbf{I}}_{r}$ and ${\mathbf{P}}_{k}=0$ for $k>1$ . Using the power series expansion in the invariance equation for the ROM (38) and collecting the terms corresponding to $\omega$ , we have the Sylvester equation

	$\displaystyle{\mathbf{P}}_{1}{\mathbf{S}}$	$\displaystyle={\mathbf{A}}_{r}{\mathbf{P}}_{1}+{\mathbf{B}}_{r}{\mathbf{L}}_{1}$
	$\displaystyle\implies{\mathbf{P}}_{1}{\mathbf{S}}$	$\displaystyle=({\mathbf{W}}^{\top}{\mathbf{V}}_{1})^{-1}{\mathbf{W}}^{\top}{\mathbf{A}}{\mathbf{V}}_{1}{\mathbf{P}}_{1}+({\mathbf{W}}^{\top}{\mathbf{V}}_{1})^{-1}{\mathbf{W}}^{\top}{\mathbf{B}}{\mathbf{L}}_{1}.$

We have that ${\mathbf{P}}_{1}={\mathbf{I}}_{r}$ is a solution to this equation since

\displaystyle{\mathbf{W}}^{\top}({\mathbf{V}}_{1}{\mathbf{I}}_{r}{\mathbf{S}})

\displaystyle={\mathbf{W}}^{\top}({\mathbf{A}}{\mathbf{V}}_{1}+{\mathbf{B}}{\mathbf{L}}_{1})

from the definition of ${\mathbf{V}}_{1}$ in (29). By the hypothesis of the theorem, we have that eigenvalues of ${\mathbf{A}}_{r}$ are disjoint from the eigenvalues of ${\mathbf{S}}$ and hence this is the unique solution of the Sylvester equation. Similarly, collecting the terms corresponding to $\omega^{(2)}$ in the invariance equation, we obtain

\displaystyle{\mathbf{P}}_{2}({\mathbf{S}}\oplus{\mathbf{S}})+({\mathbf{W}}^{\top}{\mathbf{V}}_{1})^{-1}{\mathbf{W}}^{\top}{\mathbf{V}}_{2}\left(({\mathbf{P}}_{1}{\mathbf{S}}\otimes{\mathbf{P}}_{1})+({\mathbf{P}}_{1}\otimes{\mathbf{P}}_{1}{\mathbf{S}})\right)={\mathbf{A}}_{r}{\mathbf{P}}_{2}+{\mathbf{H}}_{r}({\mathbf{P}}_{1}\otimes{\mathbf{P}}_{1})+{\mathbf{B}}_{r}{\mathbf{L}}_{2}.

Since we just established ${\mathbf{P}}_{1}={\mathbf{I}}_{r}$ , this simplifies to

\displaystyle{\mathbf{P}}_{2}({\mathbf{S}}\oplus{\mathbf{S}})+({\mathbf{W}}^{\top}{\mathbf{V}}_{1})^{-1}{\mathbf{W}}^{\top}{\mathbf{V}}_{2}\left({\mathbf{S}}\oplus{\mathbf{S}}\right)={\mathbf{A}}_{r}{\mathbf{P}}_{2}+{\mathbf{H}}_{r}+{\mathbf{B}}_{r}{\mathbf{L}}_{2}.

Rearranging terms and substituting the expressions for the ROM matrices, we get

\displaystyle{\mathbf{P}}_{2}({\mathbf{S}}\oplus{\mathbf{S}})-({\mathbf{W}}^{\top}{\mathbf{V}}_{1})^{-1}{\mathbf{W}}^{\top}{\mathbf{A}}{\mathbf{V}}_{1}{\mathbf{P}}_{2}

\displaystyle=({\mathbf{W}}^{\top}{\mathbf{V}}_{1})^{-1}{\mathbf{W}}^{\top}\left(-{\mathbf{V}}_{2}({\mathbf{S}}\oplus{\mathbf{S}})+{\mathbf{A}}{\mathbf{V}}_{2}+{\mathbf{B}}{\mathbf{L}}_{2}\right)=0

since ${\mathbf{V}}_{2}$ satisfies the Sylvester equation in (30). Thus, ${\mathbf{P}}_{2}=0$ is the only solution of this Sylvester equation since the eigenvalues of $\mathbf{A}_{r}$ are disjoint from the eigenvalues of $\mathbf{S}\oplus{\mathbf{S}}$ . Similarly, collecting terms corresponding to $\omega^{(k)}$ for $k\geq 3$ and substituting the expressions for ${\mathbf{P}}_{1},\dots,{\mathbf{P}}_{k-1}$ , gives us

\displaystyle{\mathbf{P}}_{k}({\mathbf{S}}\oplus_{k}{\mathbf{S}})-{\mathbf{A}}_{r}{\mathbf{P}}_{k}=0,

which has the unique solution ${\mathbf{P}}_{k}=0$ under the hypothesis of the theorem. Thus, we conclude that $\boldsymbol{\pi_{r}}(\omega)=\omega$ is the unique solution to (38) in a neighbourhood of 0. ∎

Beyond its role as a step toward proving the moment-matching result in Theorem 4.5, the preceding theorem provides geometric insights into the ROM center manifold. Eq. (40) demonstrates that the ROM state directly and perfectly tracks the coordinates of the driving signal generator. This identity map simplifies the approximation of the full-order center manifold, leading directly to the following corollary.

Corollary 4.4.

Let $\boldsymbol{\pi_{r}}(\omega)$ denote the solution of the ROM invariance equation (38). Then the map $\mathbf{p}(\omega)$ defined as

\mathbf{p}(\omega)={\mathbf{V}}_{1}\boldsymbol{\pi_{r}}(\omega)+{\mathbf{V}}_{2}(\boldsymbol{\pi_{r}}(\omega)\otimes\boldsymbol{\pi_{r}}(\omega))

solves the FOM invariance equation (28) along direction ${\mathbf{W}}^{\top}$ .

The significance of Corollary 4.4 lies in its direct connection to projection-based model order reduction. The statement of the corollary implies that the residual of the high-dimensional invariance equation is forced to zero when projected onto the subspace defined by the left projection matrix ${\mathbf{W}}$ . This is equivalent to a generalized Petrov-Galerkin projection condition applied directly to the nonlinear manifold equations. Here, the trial space for the full-order state is defined by a quadratic manifold expansion spanned by ${\mathbf{V}}_{1}$ and ${\mathbf{V}}_{2}$ , while the test subspace spanned by ${\mathbf{W}}$ , allows one to enforce additional desirable properties in the ROM. This provides an important link between nonlinear center manifold theory and classical projection-based frameworks. Theorems 4.1 and 4.3 allow us to prove that the ROM achieves nonlinear moment matching for the signal generator $({\mathbf{S}},{\mathbf{L}}_{1},{\mathbf{L}}_{2})$ .

Theorem 4.5.

Suppose the assumptions of Theorems 4.1 and 4.3 hold and the ROM is constructed in (23)–(24) with the choice of ${\mathbf{V}}_{1},{\mathbf{V}}_{2}$ as in (36) and ${\mathbf{W}}$ chosen so that ${\mathbf{W}}^{\top}{\mathbf{V}}_{1}$ is invertible. Then, the center manifold $\boldsymbol{\pi_{r}}(\omega)$ of the ROM maps to the center manifold $\mathbf{\Pi}(\omega)$ of the LTI system in (21) under the map

\mathbf{\Pi}(\omega)={\mathbf{V}}_{1}\boldsymbol{\pi_{r}}(\omega)+{\mathbf{V}}_{2}(\boldsymbol{\pi_{r}}(\omega)\otimes\boldsymbol{\pi_{r}}(\omega)).

Moreover, the nonlinear moment, $\widehat{\mathcal{M}}(\cdot)$ , of the ROM defined in (23)–(24) matches the nonlinear moment, $\mathcal{M}(\cdot)$ , of the LTI system in (21), i.e.,

\displaystyle\mathcal{M}(\omega)=\widehat{\mathcal{M}}(\omega),

for all $\omega$ in a neighborhood of $\omega=0$ .

Proof.

The nonlinear moment for the ROM is given as,

\displaystyle\widehat{\mathcal{M}}(\omega)={\mathbf{C}}_{r}\boldsymbol{\pi}_{r}(\omega)+{\mathbf{K}}_{r}(\boldsymbol{\pi}_{r}(\omega)\otimes\boldsymbol{\pi}_{r}(\omega)).

Since $\mathbf{\Pi}(\omega)={\mathbf{V}}_{1}\omega+{\mathbf{V}}_{2}\omega^{(2)}$ from Theorem 4.1 and $\boldsymbol{\pi_{r}}(\omega)=\omega$ in a neighbourhood around 0 from Theorem 4.3, we have that

\displaystyle\mathcal{M}(\omega)={\mathbf{C}}\left({\mathbf{V}}_{1}\omega+{\mathbf{V}}_{2}\omega^{(2)}\right)={\mathbf{C}}_{r}\boldsymbol{\pi_{r}}(\omega)+{\mathbf{K}}_{r}\boldsymbol{\pi_{r}}^{(2)}(\omega)=\widehat{\mathcal{M}}(\omega)

for all $\omega$ in a neighborhood of $0$ . ∎

This theorem shows that the choice of ${\mathbf{V}}_{1}$ and ${\mathbf{V}}_{2}$ made in (36) achieves nonlinear moment matching in the sense of astolfi2010. We also note that this result hold for any choice of ${\mathbf{W}}$ such that ${\mathbf{W}}^{\top}{\mathbf{V}}_{1}$ is invertible. Corollary 4.4 provides some insight into the role of ${\mathbf{W}}$ in the moment matching framework since we can interpret ${\mathbf{W}}$ to be the directions along which we solve the FOM invariance equations. Because this choice of ${\mathbf{W}}$ is arbitrary, we have freedom in enforcing additional properties in our ROM. One such condition is forcing the ROM to have no mass matrix term, i.e., pick ${\mathbf{W}}$ such that ${\mathbf{W}}^{\top}{\mathbf{V}}_{2}=0$ so that ${\mathbf{M}}_{r}(\mathbf{x}_{r})={\mathbf{I}}_{r}$ . Combining the theoretical results established in this section and the reduced system we derived in Section 3, we propose an algorithm to compute the quadratic reduced system (23)–(24).

Input: Linear Time-Invariant (LTI) system matrices

({\mathbf{A}},{\mathbf{B}},{\mathbf{C}})

as given in (21).

Output: Reduced Order Model (ROM) that achieves moment matching as shown in Theorem 4.5

1.

Choose matrices ${\mathbf{S}}\in\mathbb{R}^{r\times r}$ , ${\mathbf{L}}_{1}\in\mathbb{R}^{m\times r}$ , and ${\mathbf{L}}_{2}\in\mathbb{R}^{m\times r^{2}}$ to interpolate desired input signals generated by the signal generator. Some canonical choices for the signal generator are discussed in Section 5.

Solve Sylvester equations for ${\mathbf{V}}_{1}\in\mathbb{R}^{n\times r}$ and ${\mathbf{V}}_{2}\in\mathbb{R}^{n\times r^{2}}$ :

	$\displaystyle{\mathbf{A}}{\mathbf{V}}_{1}+{\mathbf{B}}{\mathbf{L}}_{1}$	$\displaystyle={\mathbf{V}}_{1}{\mathbf{S}}$
	$\displaystyle{\mathbf{A}}{\mathbf{V}}_{2}+{\mathbf{B}}{\mathbf{L}}_{2}$	$\displaystyle={\mathbf{V}}_{2}({\mathbf{S}}\oplus{\mathbf{S}})$

These Sylvester equations define the center manifold of the full order system as shown in Theorem 4.1.

3.
Choose a projection matrix ${\mathbf{W}}\in\mathbb{R}^{n\times r}$ such that ${\mathbf{W}}^{\top}{\mathbf{V}}_{1}$ is invertible. Some choices include:
- •
  
  Enforce ${\mathbf{W}}^{\top}{\mathbf{V}}_{2}=\mathbf{0}$ (Eliminates state-dependent mass-matrix term)
- •
  
  Set ${\mathbf{W}}={\mathbf{V}}_{1}$ (Galerkin projection)

Compute the reduced-order model matrices:

	$\displaystyle{\mathbf{A}}_{r}$	$\displaystyle=({\mathbf{W}}^{\top}{\mathbf{V}}_{1})^{-1}{\mathbf{W}}^{\top}{\mathbf{A}}{\mathbf{V}}_{1},\quad{\mathbf{B}}_{r}=({\mathbf{W}}^{\top}{\mathbf{V}}_{1})^{-1}{\mathbf{W}}^{\top}{\mathbf{B}},\quad{\mathbf{C}}_{r}={\mathbf{C}}{\mathbf{V}}_{1},\quad{\mathbf{K}}_{r}={\mathbf{C}}{\mathbf{V}}_{2},$
	$\displaystyle{\mathbf{H}}_{r}$	$\displaystyle=({\mathbf{W}}^{\top}{\mathbf{V}}_{1})^{-1}{\mathbf{W}}^{\top}{\mathbf{A}}{\mathbf{V}}_{2},\quad{\mathbf{E}}_{r}=({\mathbf{W}}^{\top}{\mathbf{V}}_{1})^{-1}{\mathbf{W}}^{\top}{\mathbf{V}}_{2},\quad{\mathbf{M}}_{r}(\mathbf{x}_{r})={\mathbf{I}}_{r}+{\mathbf{E}}_{r}(\mathbf{x}_{r}\otimes\mathbf{x}_{r})$

Construct the interpolating QM ROM given by:

	$\displaystyle{\mathbf{M}}_{r}(\mathbf{x}_{r})\,\dot{\mathbf{x}}_{r}(t)$	$\displaystyle={\mathbf{A}}_{r}\mathbf{x}_{r}(t)+{\mathbf{B}}_{r}\mathbf{u}(t)+{\mathbf{H}}_{r}\big(\mathbf{x}_{r}(t)\otimes\mathbf{x}_{r}(t)\big)$
	$\displaystyle\mathbf{y}_{r}(t)$	$\displaystyle={\mathbf{C}}_{r}\mathbf{x}_{r}(t)+{\mathbf{K}}_{r}(\mathbf{x}_{r}(t)\otimes\mathbf{x}_{r}(t))$

Theoretical Guarantee: By construction, the resulting QM ROM matches the nonlinear moment of the full-order model, i.e., $\hat{\mathcal{M}}(\omega)=\mathcal{M}(\omega)$ . This ensures exact steady-state dynamic reproduction for any input signal generated by $({\mathbf{S}},{\mathbf{L}}_{1},{\mathbf{L}}_{2})$ (Theorem 4.5).

Algorithm 1 Moment matching quadratic manifold approach

Algorithm 1 summarizes the step-by-step procedure for constructing the optimization-free quadratic reduced-order model. We begin, in Step 1, by defining the signal generator matrices ${\mathbf{S}},{\mathbf{L}}_{1}$ , and ${\mathbf{L}}_{2}$ , which encode the targeted interpolation frequencies and tangential directions. Then, in Step 2, these parameters are used to determine the linear and quadratic basis matrices, ${\mathbf{V}}_{1}$ and ${\mathbf{V}}_{2}$ , by solving the two decoupled linear Sylvester equations (29) and (30). Step 3 chooses the left projection matrix ${\mathbf{W}}$ , which can be tailored to enforce specific geometric properties (such as a standard Galerkin projection where ${\mathbf{W}}={\mathbf{V}}_{1}$ or an orthogonality condition ${\mathbf{W}}^{\top}{\mathbf{V}}_{2}=0$ ), provided that the invertibility condition on ${\mathbf{W}}^{\top}{\mathbf{V}}_{1}$ is satisfied. Using these bases, in Step 4, one projects the full-order system to efficiently compute the reduced matrices, which notably includes assembling the terms for the state-dependent mass matrix ${\mathbf{M}}(\mathbf{x}_{r})$ . Finally, in Step 5, the resulting nonlinear reduced-order model is computed, which is theoretically guaranteed to match the targeted nonlinear moments of the original full-order system as per Theorem 4.5.

Remark 4.6.

The online storage and evaluation of the quadratic ROM constructed via Algorithm 1 are independent of the full-order dimension $n$ . Although the projection matrices ${\mathbf{V}}_{1}\in\mathbb{R}^{n\times r}$ and ${\mathbf{V}}_{2}\in\mathbb{R}^{n\times r^{2}}$ are computed during the offline phase, they are never loaded during online execution. By precomputing all reduced-order operators offline, the online memory complexity scales strictly as $\mathcal{O}(r^{3})$ .

5 Enforcing application specific interpolation conditions

The selection of an appropriate signal generator matrix $\mathbf{S}$ for moment matching is inherently problem-dependent, as it dictates the specific family of trajectories the reduced-order model (ROM) is constrained to replicate. In this section, we detail canonical choices for $\mathbf{S}$ designed to capture distinct dynamic behaviors, ranging from persistent steady-state oscillations to quasi-polynomial transient signals, and provide explicit analytical solutions for the resulting manifold templates.

5.1 Real canonical forms for sinusoidal inputs

To interpolate purely sinusoidal steady-state signals at a given set of driving frequencies, the signal generator matrix must possess purely imaginary eigenvalues $\{\pm j\mu_{1},\pm j\mu_{2},\dots,\pm j\mu_{k}\}$ . To ensure that all subsequent reduction stages, including the derivation of the projection matrices $\mathbf{V}_{1}$ and $\mathbf{V}_{2}$ , and consequently the ROM operators themselves, remain strictly within the real domain, we enforce a real canonical form for $\mathbf{S}$ astolfi2010; morAntBG20:

\mathbf{S}=\operatorname{blkdiag}\left(\begin{bmatrix}0&\mu_{1}\\ -\mu_{1}&0\end{bmatrix},\,\dots,\,\begin{bmatrix}0&\mu_{k}\\ -\mu_{k}&0\end{bmatrix}\right)\in\mathbb{R}^{r\times r},

(41)

where $r=2k$ . By matching the moment with this block-diagonal structure, the ROM accurately replicates the frequency response of the original high-dimensional system at the selected frequencies. This choice of ${\mathbf{S}}$ results in inputs of the form

\displaystyle\mathbf{u}(t)=\left(\sum_{i=1}^{r}c_{i,1}\sin(\mu_{i}t)+c_{i,2}\cos(\mu_{i}t)\right)+\left(\sum_{i,j=1}^{r}c_{ij,1}\sin((\mu_{i}+\mu_{j})t)+c_{ij,2}\cos((\mu_{i}+\mu_{j})t)\right),

(42)

where $c_{i,1},c_{i,2},c_{ij,1},c_{ij,2}$ are constants dependent on ${\mathbf{L}}_{1},{\mathbf{L}}_{2},\omega_{0}$ for $1\leq i,j,\leq r$ . The constructed ROM, thereby, matches the steady-state behaviour of the FOM for all inputs of the form in (42). We note that the ROM not only interpolates frequencies $\mu_{i},~1\leq i\leq r$ , but also the mixed interpolation frequencies $\mu_{i}+\mu_{j},1\leq i,j\leq r$ . Furthermore, this real canonical representation of ${\mathbf{S}}$ eliminates the need for complex arithmetic, directly yielding a real-valued reduced-order system. This construction is conceptually analogous to the basis transformations employed to enforce realness within interpolatory model reduction and Loewner frameworks morAntBG20.

5.2 Closed-form solutions for diagonal signal generators

To illustrate the explicit algebraic bridge between center manifold theory and classical transfer function interpolation, consider a diagonalized signal generator configuration where $\mathbf{S}=\operatorname{diag}(s_{1},\dots,s_{r})\in\mathbb{C}^{r\times r}$ . Let $\mathbf{\ell}_{i}\in\mathbb{R}^{m}$ denote the $i$ -th column of $\mathbf{L}_{1}$ and let $\mathbf{b}_{ij}\in\mathbb{R}^{m}$ represent the column of $\mathbf{L}_{2}$ associated with the interacting state components $\omega_{i}\omega_{j}$ (for $i,j=1,\dots,r$ ). Under this diagonal framework, the Sylvester equations defining the linear and quadratic manifold mappings decouple column-by-column. This yields the following exact, closed-form expressions for the columns of $\boldsymbol{\Pi}_{1}$ and $\boldsymbol{\Pi}_{2}$ :

	$\displaystyle\boldsymbol{\Pi}_{1}$	$\displaystyle=\left[(s_{1}\mathbf{I}-\mathbf{A})^{-1}\mathbf{B}\mathbf{\ell}_{1},\,\dots,\,(s_{r}\mathbf{I}-\mathbf{A})^{-1}\mathbf{B}\mathbf{\ell}_{r}\right]\in\mathbb{R}^{n\times r},$		(43)
	$\displaystyle\boldsymbol{\Pi}_{2}$	$\displaystyle=\left[\dots,\,\left((s_{i}+s_{j})\mathbf{I}-\mathbf{A}\right)^{-1}\mathbf{B}\mathbf{b}_{ij},\,\dots\right]\in\mathbb{R}^{n\times r^{2}}.$		(44)

Recalling the FOM transfer function $\mathbf{H}(s)=\mathbf{C}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B}$ , the steady-state nonlinear system output profile $\mathcal{M}(\omega)=\mathbf{C}\boldsymbol{\Pi}_{1}\omega+\mathbf{C}\boldsymbol{\Pi}_{2}(\omega\otimes\omega)$ can be mapped directly to frequency-domain evaluations of the full-order model as

\mathcal{M}(\omega)=\sum_{i=1}^{r}\mathbf{H}(s_{i})\mathbf{\ell}_{i}\omega_{i}+\sum_{i=1}^{r}\sum_{j=1}^{r}\mathbf{H}(s_{i}+s_{j})\mathbf{b}_{ij}\omega_{i}\omega_{j}.

(45)

Equation (45) provides a clear interpretation of the framework: the linear manifold component interpolates the system’s response at the fundamental driving frequencies $s_{i}$ , while the quadratic component captures the ”mixed” frequencies $s_{i}+s_{j}$ . Thus, the reduced-order model incorporates the mathematical effects of $r$ frequencies as well as $(r^{2}+r)/2$ mixed frequencies while maintaining a strictly $r$ -dimensional reduced state-space.

Remark 5.1.

For readers acquainted with the projection-based interpolatory model reduction for matching linear moments, the structure of $\boldsymbol{\Pi}_{1}$ and $\boldsymbol{\Pi}_{2}$ would look familiar. However, the structural formulation in (45) reveals a crucial distinction regarding the interpolation properties of the reduced-order model. Under the classical linear rational interpolation paradigm morAntBG20, matching the full-order transfer function at the fundamental frequencies $s_{i}$ and the mixed harmonic frequencies $s_{i}+s_{j}$ would require a linear ROM of total dimension $r+(r+r^{2})/2=\mathcal{O}(r^{2})$ . A linear subspace treats the mixed-frequency modes $s_{i}+s_{j}$ as independent, uncoupled coordinates. While such a high-dimensional linear system achieves identical one-sided interpolation points, it is not guaranteed to preserve the system’s intrinsic center manifold geometry.

5.3 Matching higher-order derivatives

The proposed framework can be systematically extended to capture complex transient dynamics such as quasipolynomial signals that naturally arise in systems characterized by repeated eigenvalues or damped modes. This is accomplished by constructing a signal space that interpolates the higher-order derivatives of the transfer function. To illustrate this construction, consider a signal generator defined by the matrices

\mathbf{S}=\begin{bmatrix}\mu&1\\ 0&\mu\end{bmatrix}\in\mathbb{C}^{2\times 2},\quad\mathbf{L}_{1}=\begin{bmatrix}\ell_{1}&\ell_{2}\end{bmatrix}\in\mathbb{C}^{m\times 2},\quad\mathbf{L}_{2}=\begin{bmatrix}b_{1}&b_{2}&b_{3}&b_{4}\end{bmatrix}\in\mathbb{C}^{m\times 4}.

(46)

The resulting nonlinear moment $\mathcal{M}(\mathbf{w})$ can then be evaluated as

$\displaystyle\mathcal{M}(\mathbf{w})$	$\displaystyle={\mathbf{C}}\mathbf{V}_{1}\mathbf{w}+{\mathbf{C}}\mathbf{V}_{2}(\mathbf{w}\otimes\mathbf{w})$
	$\displaystyle=\mathbf{H}(\mu)(\ell_{1}w_{1}+\ell_{2}w_{2})+\mathbf{H}(2\mu)(b_{1}w_{1}^{2}+b_{2}w_{1}w_{2}+b_{3}w_{2}w_{1}+b_{4}w_{2}^{2})$
	$\displaystyle\quad+\mathbf{H}^{\prime}(\mu)\ell_{1}w_{2}+\mathbf{H}^{\prime}(2\mu)\left(b_{1}w_{1}w_{2}+b_{1}w_{2}w_{1}+(b_{2}+b_{3})w_{2}^{2}\right)+\mathbf{H}^{\prime\prime}(2\mu)b_{1}w_{2}^{2}$
	$\displaystyle=c_{1}\mathbf{H}(\mu)+c_{2}\mathbf{H}(2\mu)+c_{3}\mathbf{H}^{\prime}(\mu)+c_{4}\mathbf{H}^{\prime}(2\mu)+c_{5}\mathbf{H}^{\prime\prime}(2\mu).$	(47)

This formulation yields a basis of amplitude-modulated, exponentially decaying or oscillating signals. Consequently, any scalar control input $\mathbf{u}(t)$ generated within this signal space takes the explicit form

\mathbf{u}(t)=c_{1}e^{\mu t}+c_{2}te^{\mu t}+c_{3}e^{2\mu t}+c_{4}te^{2\mu t}+c_{5}t^{2}e^{2\mu t},

(48)

where the weighting coefficients $c_{1},\dots,c_{5}$ are uniquely determined by the initial conditions $\mathbf{w}(0)$ and the matrices $\mathbf{L}_{1},\mathbf{L}_{2}$ . Because $s\in\mathbb{C}$ , the inputs and projection matrices are, in general, complex-valued. To enforce real-valued physical quantities, $\mathbf{S}$ can be represented in real canonical form to interpolate $s=i\mu$ ( $\mu\in\mathbb{R}$ ):

\mathbf{S}=\begin{bmatrix}0&\mu&1&0\\ -\mu&0&0&1\\ 0&0&0&\mu\\ 0&0&-\mu&0\end{bmatrix}\in\mathbb{R}^{4\times 4}.

(49)

By parameterizing the signal generator in this real-valued canonical structure with a Jordan block, the nonlinear moment matching condition reduces directly into solving matrix Sylvester equations. This aligns with the foundational Sylvester equation formulations for projection based model reduction for LTI system in gallivan2004sylvester.

6 Numerical results

While the theoretical results and the reduction framework proposed in this work are broadly applicable to any linear time-invariant (LTI) system of the form (21) including regimes where conventional linear subspaces are proven to perform effectively, linear projection methods remain fundamentally ill-suited for efficiently capturing shift-variant dynamics or LTI systems with slow decay of Hankel singular values (i.e., slow Kolmogorov $n$ -width decay). Thus, we consider classical transport-dominated problems, specifically the one-dimensional advection-diffusion equation and the one-dimensional wave equation, with their characteristic slow decay of the (normalized) Hankel singular values shown in Figure 5. Consequently, such systems serve as ideal benchmarks to demonstrate the efficacy and computational gains of the proposed quadratic manifold framework. All numerical simulations were implemented on MATLAB R2023b (version 23.2.0, Update 4) on a laptop equipped with an Apple M3 Pro chip, running macOS 14.6.1. The repository to generate all plots and simulations in this paper can be found in padhi_code.

(a) One-dimensional advection equation (

n=1024

(b) One-dimensional damped wave equation (

n=4096

Figure 5: Decay behavior of the normalized Hankel singular values (HSV), illustrating the characteristically slow decay rates for the respective systems with the chosen parameters.

6.1 One-dimensional linear transport equation

As the first benchmark, we consider the dynamics of a one-dimensional linear transport (advection) equation defined on a spatial domain of length 1. This problem models the pure translation of a scalar field with a constant velocity $v=1$ directed toward the right. The system is governed by the following first-order hyperbolic partial differential equation (PDE):

\frac{\partial\xi}{\partial t}+\frac{\partial\xi}{\partial z}=0,\quad z\in(0,1),\quad t>0.

(50)

Here, $\xi(z,t)$ represents the continuous scalar state (for example displacement) at spatial position $z$ and time $t$ . The system is driven by a time-dependent boundary condition on the left ( $z=0$ ), given by

\xi(0,t)=\mathbf{u}(t),

(51)

where $\mathbf{u}(t)$ is a scalar input. The system is observed at the right boundary ( $z=1$ ), with $\mathbf{y}(t)=\xi(1,t)$ representing the scalar output of the system. By explicitly enforcing the Dirichlet boundary condition at the inlet node ( $z=0$ ), the remaining $n=N$ interior and outlet nodes form the active state vector $\mathbf{x}(t)\in\mathbb{R}^{n}$ . The resulting LTI system is in the form of (21) where the state-space matrices $\mathbf{A}\in\mathbb{R}^{n\times n}$ and $\mathbf{B}\in\mathbb{R}^{n\times 1}$ are obtained using Chebyshev pseudospectral method for spatial discretization. This ensures we do not introduce numerical diffusion to the system, thereby making the system diffusion-dominated. The output matrix $\mathbf{C}\in\mathbb{R}^{1\times n}$ acts as a selector vector isolating the node corresponding to the observation node at $z=1$ and we set $n=1024$ to obtain a LTI system of that dimension. As we see in Figure 5(a), the system showcases extreme slow decay in the singular values making it a challenging problem for model reduction.

Constructing the interpolating reduced system

We construct a reduced-order system as defined in Section 3 that achieves moment matching for a specified signal space $({\mathbf{S}},{\mathbf{L}}_{1},{\mathbf{L}}_{2})$ as shown in Theorem. 4.5. The signal generator matrices are chosen as,

	$\displaystyle{\mathbf{S}}$	$\displaystyle=\operatorname{blkdiag}\left(\left[\begin{array}[]{cc}0&\mu_{1}\\ \overline{\mu}_{1}&0\end{array}\right],\dots,\left[\begin{array}[]{cc}0&\mu_{r/2}\\ \overline{\mu}_{r/2}&0\end{array}\right]\right)\in\mathbb{R}^{r\times r},$		(52)
	$\displaystyle{\mathbf{L}}_{1}$	$\displaystyle=[1~1~\dots~1]\in\mathbb{R}^{1\times r},\quad{\mathbf{L}}_{2}=[1~1~\dots~1]\in\mathbb{R}^{1\times r^{2}}.$		(52)

We set the reduction order to $r=20$ . The interpolation frequencies are chosen as complex conjugate pairs, spaced logarithmically along the imaginary axis between $10^{1}$ and $10^{3}$ . This selection also allows us to study dynamics with different timescales using slow varying and fast varying input signals. For simplicity and uniformity, the tangential weights are set to ones. Under this choice of signal generator, we construct a quadratic ROM of the form (23)–(24) using Algorithm 1.

Refer to caption — Figure 6: Steady state reconstruction of the discretized linear transport equation system on the center manifold.

To show that the quadratic ROM (23)– (24) constructed from the algorithm indeed achieves moment matching, we generate a signal from the signal generator in (52) with a random initial condition $\omega_{0}\in\mathbb{R}^{20}$ . This input signal is then used to simulate the QM ROM and FOM from $t\in[0,1]$ . We plot the state evolution and output of both the full and reduced order systems assuming that the two systems start on their respective center manifolds, i.e.,

\displaystyle\mathbf{x}(0)

\displaystyle=\boldsymbol{\pi}(\omega_{0})={\mathbf{V}}_{1}\omega_{0}+{\mathbf{V}}_{2}\omega_{0}^{(2)}\in\mathbb{R}^{1024},\quad\quad\mathbf{x}_{r}(0)=\boldsymbol{\pi}_{r}(\omega_{0})=\omega_{0}\in\mathbb{R}^{20}.

We then lift the reduced states $\mathbf{x}_{r}(t)\in\mathbb{R}^{20}$ (for $t>0$ ) to the full space using the quadratic mapping,

\displaystyle\boldsymbol{\Phi}(\mathbf{x}_{r}(t))=\mathbf{V}_{1}\mathbf{x}_{r}(t)+\mathbf{V}_{2}(\mathbf{x}_{r}(t)\otimes\mathbf{x}_{r}(t))\in\mathbb{R}^{1024},

allowing us to evaluate the state reconstruction error against the true FOM solution. Figure 6 illustrates the space-time evolution of the state for the full-order model (FOM), the proposed quadratic manifold reduced-order model (ROM), and the resulting absolute reconstruction error under the input signal generated by the signal generator in (52). The distinct diagonal bands in the left and middle panels prominently capture the continuous wave propagation of the FOM ( $n=1024$ ) and the ROM ( $r=20$ ), respectively, across the spatial grid over the simulation time $t\in[0,1]$ . Despite a drastic reduction in dimensionality, the quadratic manifold ROM successfully tracks the travelling wavefronts accurately.

Figure 7: Numerical verification of moment matching property for the one-dimensional advection equation

We recall that the steady state output of the FOM under the signal generator is given by the mapping $\mathcal{M}(\omega(t))$ where $\omega(t)$ represents the signal state vector. The top panel in Figure 7 displays the time history of the steady-state system outputs $\mathcal{M}(\omega(t))$ and $\widehat{\mathcal{M}}(\omega(t))$ for the FOM and the quadratic manifold ROM (QM ROM), respectively, over $t\in[0,1]$ . The trajectories of the full and reduced models are virtually indistinguishable, with the red dashed line perfectly tracking the complex, multi-frequency oscillations of the reference full-order system. The bottom panel plots the absolute output error $|\mathcal{M}(\omega(t))-\widehat{\mathcal{M}}(\omega(t))|$ for both system components on a semi-logarithmic scale. Theoretically, the quadratic ROM tracks the steady state outputs and center manifold of the full order system exactly. However, minor discrepancies accumulate in practice due to numerical integration. Nevertheless, the relative error remains strictly bounded and low throughout the entire simulation horizon. These results showcase the capability of the proposed framework to accurately capture transport-dominated dynamics within an exceptionally low-dimensional state space. Table 1(a) reports the empirical online runtimes comparing the QM-ROM against the FOM, with the QM-ROM offering over a $380\times$ online speedup. Further computational profiling and theoretical time-complexity details are provided in Section 6.3.

6.2 One dimensional damped wave equation

As the next benchmark, we consider the dynamics of a one-dimensional damped wave equation on a spatial domain of length $1$ . The system is governed by the PDE given below.

\frac{\partial^{2}\xi}{\partial t^{2}}+\gamma\frac{\partial\xi}{\partial t}=c^{2}\frac{\partial^{2}\xi}{\partial z^{2}}+f(z,t),

(53)

where $\xi(z,t)$ is the displacement at spatial position $z$ and time $t$ , $c$ represents the speed of wave propagation, $\gamma>0$ is the viscous damping coefficient, and $f(z,t)$ represents the external applied force. The PDE is spatially discretized over $n$ internal grid points using a standard central finite difference method. Let $\mathbf{w}(t)\in\mathbb{R}^{n}$ represent the vector of displacements at these grid points. To cast this second-order system into a standard first-order form, we augment the displacement and velocity vectors to define the state vector as $\mathbf{x}(t)=\begin{bmatrix}\mathbf{w}(t)^{T}&\dot{\mathbf{w}}(t)^{T}\end{bmatrix}^{T}\in\mathbb{R}^{2n}$ . This yields a continuous-time linear time-invariant (LTI) system of dimension $2n$ in the form of (21) where $\mathbf{u}(t)$ is the scalar input and $\mathbf{y}(t)$ is the scalar output (the measured displacement). The spatial distribution of the input, originally $f(z,t)$ in the PDE, is captured by the input matrix ${\mathbf{B}}$ . The system matrix ${\mathbf{A}}\in\mathbb{R}^{2n\times 2n}$ , the input matrix ${\mathbf{B}}\in\mathbb{R}^{2n\times 1}$ , and the output matrix ${\mathbf{C}}\in\mathbb{R}^{1\times 2n}$ are defined as

{\mathbf{A}}=\begin{bmatrix}\mathbf{0}_{n\times n}&{\mathbf{I}}_{n\times n}\\ c^{2}{\mathbf{L}}_{n}&-\gamma{\mathbf{I}}_{n\times n}\end{bmatrix},\quad{\mathbf{B}}=\begin{bmatrix}\mathbf{0}_{n\times 1}\\ \tilde{{\mathbf{B}}}\end{bmatrix},\quad{\mathbf{C}}=\begin{bmatrix}\tilde{{\mathbf{C}}}&\mathbf{0}_{1\times n}\end{bmatrix},

(54)

where ${\mathbf{L}}_{n}$ is the discrete Laplacian matrix representing the finite difference approximation of the spatial derivative. In this specific implementation, the input force is applied at the center of the domain. Therefore, $\tilde{{\mathbf{B}}}$ places a value of $1$ at the index corresponding to the middle node and zeros everywhere else, meaning the input directly affects only the velocity of that specific node. Similarly, the measurement matrix $\tilde{{\mathbf{C}}}$ is designed to extract the displacement of this exact same middle node.

For the numerical experiments, the physical parameters are chosen as $c=1$ , and $\gamma=10^{-8}$ . Spatially discretizing the domain with $n=2048$ internal grid points yields an LTI system of dimension $2n=4096$ . The selection of a small, non-zero damping coefficient $\gamma$ ensures that the system dynamics remain transport-dominated, while strictly displacing the eigenvalues of ${\mathbf{A}}$ from the imaginary axis into the open left-half plane, a standard and realistic assumption when numerically simulating such systems. The Hankel singular values (HSVs) for this system are shown in Figure 5(b). As illustrated, the HSVs exhibit an exceptionally slow decay rate. This characteristic is typical for wave propagation problems and indicates that traditional linear model reduction methods would require a high-dimensional state to capture the dynamics accurately.

Constructing the interpolating reduced system

In this section, we construct a reduced-order system that achieves moment matching for a specified signal space defined by the tuple $({\mathbf{S}},{\mathbf{L}}_{1},{\mathbf{L}}_{2})$ . The signal generator matrix ${\mathbf{S}}$ is selected to match the imaginary frequencies close to the eigenvalues of the original system matrix ${\mathbf{A}}$ . For simplicity and uniformity, the tangential weights are set to ones. The resulting signal generator matrices are chosen as,

	$\displaystyle{\mathbf{S}}$	$\displaystyle=\operatorname{blkdiag}\left(\left[\begin{array}[]{cc}0&\mu_{1}\\ \overline{\mu}_{1}&0\end{array}\right],\dots,\left[\begin{array}[]{cc}0&\mu_{r/2}\\ \overline{\mu}_{r/2}&0\end{array}\right]\right)\in\mathbb{R}^{r\times r},$
	$\displaystyle{\mathbf{L}}_{1}$	$\displaystyle=[1~1~\dots~1]\in\mathbb{R}^{1\times r},\quad\quad{\mathbf{L}}_{2}=[1~1~\dots~1]\in\mathbb{R}^{1\times r^{2}}.$

We set the reduction order to $r=16$ . The interpolation frequencies are chosen as complex conjugate pairs, spaced logarithmically along the imaginary axis between $10$ and $10^{2.5}$ . This selection is motivated by the fact that the poles of the original system are located around this region. By spanning this frequency range, the reduced model accurately captures both the high-frequency components, which represent short-scale phenomena, and the low-frequency components, which dictate slow-scale dynamics. Following the procedure in Algorithm 1, we compute the projection matrices ${\mathbf{V}}_{1}$ and ${\mathbf{V}}_{2}$ by solving the Sylvester equations in (29) and (30) and setting ${\mathbf{W}}={\mathbf{V}}_{1}$ .

Figure 9: Numerical verification of the moment matching property in the one-dimesional wave equation

Just like in the advection equation example, we try to show that the quadratic ROM ((23)– (24)) constructed from the algorithm indeed achieves moment matching at $({\mathbf{S}},{\mathbf{L}}_{1},{\mathbf{L}}_{2})$ . We generate an input signal generated by the signal generator with a random initial condition $\omega_{0}\in\mathbb{R}^{r}$ and plot the state evolution and output of both the full and reduced order systems, over $t\in[0,2]$ , assuming that the two systems start on their respective center manifolds, i.e.,

\displaystyle\mathbf{x}(0)=\boldsymbol{\pi}(\omega_{0})={\mathbf{V}}_{1}\omega_{0}+{\mathbf{V}}_{2}\omega_{0}^{(2)}\in\mathbb{R}^{4096}~~\mbox{and}~~\mathbf{x}_{r}(0)=\boldsymbol{\pi}_{r}(\omega_{0})=\omega_{0}\in\mathbb{R}^{16}.

To facilitate direct comparison between the low-dimensional state trajectory $\mathbf{x}_{r}(t)\in\mathbb{R}^{16}$ and the full-order state trajectory $\mathbf{x}(t)\in\mathbb{R}^{4096}$ , the ROM state is lifted to the full-order space via the mapping

\displaystyle\boldsymbol{\Phi}(\mathbf{x}_{r}(t))=\mathbf{V}_{1}\mathbf{x}_{r}(t)+\mathbf{V}_{2}(\mathbf{x}_{r}(t)\otimes\mathbf{x}_{r}(t))\in\mathbb{R}^{4096}.

Figure 8 illustrates the state evolution of the FOM, state reconstruction of the QM-ROM and the absolute error in the state reconstruction of the QM-ROM. Note that the full state vector $\mathbf{x}(t)$ is partitioned such that the first $n$ components correspond to the nodal displacements, while the remaining $n$ components represent the nodal velocities. The left panel displays the evolution of the state displacement and velocity versus the reduced-order approximation $\boldsymbol{\Phi}(\mathbf{x}_{r}(t))$ . It is evident that the quadratic ROM successfully captures the transient dynamics of the wave equation, remaining tightly bound to the full-order trajectory even as the wave propagates. The right panel quantifies this performance through the absolute error. The error remains small and bounded, confirming that the solution to the Sylvester equations (29) and (30) effectively identifies the nonlinear center manifold of the system.

In addition to state-space reconstruction, Figure 9 evaluates the steady-state output response $\mathcal{M}(\boldsymbol{\omega}(t))$ of the full-order system alongside the reduced-order output $\widehat{\mathcal{M}}(\boldsymbol{\omega}(t))$ obtained from the quadratic moment-matching ROM (QM ROM). As depicted in the top panel, the output trajectory predicted by the QM ROM is visually indistinguishable from the full-order response over the entire time interval $t\in[0,2]$ . The bottom panel provides the absolute error $\|\mathcal{M}(\boldsymbol{\omega}(t))-\widehat{\mathcal{M}}(\boldsymbol{\omega}(t))\|_{2}$ , which remains bounded between $10^{-12}$ and $10^{-11}$ . This serves as concrete numerical verification that the proposed ROM satisfies the moment-matching conditions in the steady-state output space for the one-dimensional wave equation. Table 1(b) reports the empirical online runtimes comparing the QM-ROM against the FOM, with the QM-ROM offering over a $330\times$ online speedup.

6.3 Computational profiling and online performance

In this section, we address the computational scaling of the framework with system dimension $n$ . We use the term “offline time” to refer to the time taken to compute the reduced system and the term “online time” to mean the time taken to simulate the system once it is formed. We show that our method not only has a competitive online time scaling but also quick offline time. A potential computational bottleneck in quadratic manifold-based ROMs stems from the presence of the state-dependent reduced mass matrix, given in (24) as

\mathbf{M}_{r}(\mathbf{x}_{r})=\mathbf{I}_{r}+\mathbf{E}_{r}(\mathbf{I}_{r}\otimes\mathbf{x}_{r}+\mathbf{x}_{r}\otimes\mathbf{I}_{r}).

(55)

Since $\mathbf{M}_{r}(\mathbf{x}_{r})$ varies continuously with the reduced state $\mathbf{x}_{r}(t)$ , a naive implementation would explicitly assemble the large, highly redundant Kronecker product matrix of size $r^{2}\times r$ and solve a dense $r\times r$ linear system at every time step or stage of an ODE solver. Such an unoptimized approach would heavily penalize the execution runtime and jeopardize any computational gains. We provide two methods to effectively handle the simulation of this term without losing online speedups.

We exploit the underlying algebraic structure of the Kronecker product, allowing us to evaluate the action of $\mathbf{M}_{r}(\mathbf{x}_{r})$ without ever explicitly forming or storing the constituent high-dimensional matrices. Consider the term ${\mathbf{M}}_{r}(\hat{\mathbf{x}})$

\displaystyle{\mathbf{M}}_{r}(\hat{\mathbf{x}}):=({\mathbf{I}}_{r}+{\mathbf{E}}_{r}(\hat{\mathbf{x}}\otimes{\mathbf{I}}_{r}+{\mathbf{I}}_{r}\otimes\hat{\mathbf{x}}))\quad\text{with }{\mathbf{E}}_{r}=\left[{\mathbf{E}}_{r}^{(1)}\quad{\mathbf{E}}_{r}^{(2)}\quad\dots\quad{\mathbf{E}}_{r}^{(r)}\right],

where ${\mathbf{E}}_{r}\in\mathbb{R}^{r\times r^{2}}$ as defined in (24), $\hat{\mathbf{x}}\in\mathbb{R}^{r}$ is an arbitrary vector, and ${\mathbf{E}}_{r}^{(i)}\in\mathbb{R}^{r\times r}$ represent the matrix blocks of ${\mathbf{E}}_{r}$ . Then using the properties of the Kronecker product, we have that

\displaystyle{\mathbf{M}}_{r}(\hat{\mathbf{x}})=\sum_{i=1}^{r}\mathbf{x}_{i}{\mathbf{E}}_{r}^{(i)}+\left[{\mathbf{E}}_{r}^{(1)}\hat{\mathbf{x}}\quad{\mathbf{E}}_{r}^{(2)}\hat{\mathbf{x}}\quad\dots\quad{\mathbf{E}}_{r}^{(r)}\hat{\mathbf{x}}\right].

(56)

As we see in (56), the computation of ${\mathbf{M}}_{r}(\hat{\mathbf{x}})$ can be done in $\mathcal{O}(r^{3})$ flops where $r\ll n$ , e.g., $r\approx 20$ . Direct LU or Cholesky factorization of this assembled matrix also scales as $\mathcal{O}(r^{3})$ , thereby ensuring the online time complexity is cubic in $r$ and significantly cheaper than solving the full order system. The implementation details can be found in the repository padhi_code.

As an alternative, enforcing the condition ${\mathbf{W}}^{\top}{\mathbf{V}}_{2}=0$ (a standard practice in quadratic manifold approaches) results in ${\mathbf{E}}_{r}=0$ , simplifying the reduced mass matrix to the identity matrix. While our framework retains the generality of a state-dependent mass matrix to ensure strict moment matching, this choice of ${\mathbf{W}}$ offers a viable path for scenarios where online efficiency is the paramount priority.

Table 1 summarizes the online simulation runtimes for the full-order models and the proposed QM ROMs across both benchmark problems. The QM ROM achieves substantial online speedups, reducing execution times from hundreds of seconds down to under one second and thus yielding an approximate $388\times$ speedup for the 1D advection equation and a $330\times$ speedup for the 1D damped wave equation. This reduction in computational cost stems directly from confining the online time integration to the low-dimensional state space $\mathbb{R}^{r}$ , effectively decoupling the online solver cost from the large full-order dimension $n$ .

Table 1: Computational wall-clock online runtime comparison between the full-order model (FOM) and proposed QM ROM across benchmark problems.

(a) 1D Advection Equation

Method	Online Time (s)
FOM ( $n=1024$ )	378.4343
QM ROM ( $r=20$ )	0.9752

(b) 1D Damped Wave Equation

Method	Online Time (s)
FOM ( $n=4096$ )	167.8696
QM ROM ( $r=16$ )	0.5042

6.4 Theoretical and numerical comparisons with linear methods

To contextualize the performance benefits of the proposed Quadratic Manifold (QM) framework, we evaluate its online time complexity and empirical performance against classical linear projection-based reduced-order models (ROMs), such as One-Sided Rational interpolation (OSR) morAntBG20.

For transport-dominated or wave-propagation problems, e.g., the 1D advection equation in Section 6.1, linear projection methods suffer from a fundamental barrier due to the slow decay of the Kolmogorov $n$ -width (the Hankel singular values for the LTI systems). Thus, to reach an acceptable error tolerance, a linear subspace requires a large subspace dimension. As discussed in Section 5.2, the nonlinear moment $\mathcal{M}(\omega)$ of the system is given by

\mathcal{M}(\omega)=\sum_{i=1}^{r}\mathbf{H}(s_{i})\mathbf{\ell}_{i}\omega_{i}+\sum_{i=1}^{r}\sum_{j=1}^{r}\mathbf{H}(s_{i}+s_{j})\mathbf{b}_{ij}\omega_{i}\omega_{j}.

(57)

To match these dynamics, a linear ROM must interpolate the transfer function at $r+(r^{2}+r)/2$ unique points, increasing the required linear subspace dimension to $r_{\text{linear}}=\mathcal{O}(r^{2})$ . The online time-complexity for linear methods scales as $\mathcal{O}(r_{\text{linear}}^{2})$ since their reduced matrices can be pre-computed and there is no time-dependent mass-matrix term. However, due to larger state dimension, the execution cost of an equivalent linear ROM scales as $\mathcal{O}(r_{\text{linear}}^{2})=\mathcal{O}(r^{4})$ . In contrast, the QM framework confines the online ODE solver strictly to a low-dimensional state vector $\mathbf{x}_{r}(t)\in\mathbb{R}^{r}$ . By exploiting tensor-structured evaluations in (56), the QM ROM achieves an online time complexity of $\mathcal{O}(r^{3})$ .

Figure 10: Comparative analysis of the offline and online runtime as a function of the reduced manifold dimension

r

for One-Sided Rational Interpolation (OSR) and the proposed system-theoretic Quadratic Manifold (QM) framework.

To empirically validate these online runtime properties, we benchmark both frameworks on the 1D linear advection problem ( $n=1024$ ) across manifold dimensions $r\in[4,26]$ . We ensure a fair comparison by comparing a quadratic ROM of dimension $r$ with a linear ROM that interpolates the transfer function ${\mathbf{H}}(s)$ at the respective $r+(r^{2}+r)/2$ points. For instance, a QM model with $r=20$ is evaluated against a linear OSR model interpolating the transfer function at the corresponding 230 points. The OSR implementation constructs a linear projection matrix followed by orthogonalization for numerical stability, selecting only directions corresponding to non-zero singular values in rank-deficient cases. This results in an added offline cost when $r_{\text{linear}}$ is large. Empirical online runtimes are measured using the CPU time required by the ODE solver to integrate the ROM solution over $t\in[0,1]$ , as shown in Figure 10.

In the online phase (Figure 10, left), the linear OSR model exhibits rapid (theoretically $\mathcal{O}(r^{4})$ ) runtime growth. Conversely, the QM ROM maintains a milder $\mathcal{O}(r^{3})$ growth, delivering up to an order-of-magnitude wall-clock speedup at higher dimensions $r$ . In the offline phase (Figure 10, right), QM remains highly scalable by solving decoupled Sylvester equations (29)–(30). We note that for some choices of signal generator, ${\mathbf{V}}_{1}$ and ${\mathbf{V}}_{2}$ can be computed using their closed-form solutions as shown in Section 5.2. Because both QM and OSR theoretically guarantee exact interpolation of the targeted nonlinear moments (as per Theorem 4.5 and (57)), the QM framework achieves superior online computational efficiency without compromising on accuracy in this example.

6.5 Comparison with existing quadratic manifold methods

To contextualize the efficacy and system-theoretic guarantees of the proposed interpolatory quadratic manifold approach, we benchmark it against state-of-the-art data-driven quadratic manifold techniques. In current literature, non-linear reduced-order models (ROMs) defined on quadratic manifolds are frequently constructed by optimizing for the linear and quadratic projection matrices, ${\mathbf{V}}_{1}$ and ${\mathbf{V}}_{2}$ , using state snapshot trajectories. As a benchmark, we consider the data-driven greedy snapshot optimization framework introduced in schwerdtner2024greedy. To ensure a consistent baseline within an intrusive model reduction setting, both the greedy method and our proposed interpolatory formulation are provided full access to the full-order model (FOM) system operators ${\mathbf{A}},{\mathbf{B}},$ and ${\mathbf{C}}$ .

We evaluate both methods on the 1D advection equation benchmark detailed in Section 6.1. A training dataset of $1{,}000$ state snapshots is collected over the initial interval $t\in[0,1]$ by driving the FOM with the signal generator

\displaystyle\begin{aligned} {\mathbf{S}}&=\operatorname{blkdiag}\left(\left[\begin{array}[]{cc}0&\mu_{1}\\ \overline{\mu}_{1}&0\end{array}\right],\dots,\left[\begin{array}[]{cc}0&\mu_{r/2}\\ \overline{\mu}_{r/2}&0\end{array}\right]\right)\in\mathbb{R}^{r\times r},\end{aligned}

which results in inputs of the form

\displaystyle\mathbf{u}(t)=\left(\sum_{i=1}^{r}c_{i,1}\sin(\mu_{i}t)+c_{i,2}\cos(\mu_{i}t)\right)+\left(\sum_{i,j=1}^{r}c_{ij,1}\sin((\mu_{i}+\mu_{j})t)+c_{ij,2}\cos((\mu_{i}+\mu_{j})t)\right),

where $c_{i,1},c_{i,2},c_{ij,1},c_{ij,2}$ are constants dependent on ${\mathbf{L}}_{1},{\mathbf{L}}_{2},\omega_{0}$ for $1\leq i,j,\leq r$ as presented in (42). Following the greedy optimization algorithm from schwerdtner2024greedy with our implementation provided in padhi_code, projection matrices ${\mathbf{V}}_{g,1}$ and ${\mathbf{V}}_{g,2}$ are computed to minimize empirical state trajectory reconstruction error over the $t\in[0,1]$ training window. Model accuracy is measured across the same training snapshot dataset. The test basis (left projection matrix) is set to ${\mathbf{W}}_{g}={\mathbf{V}}_{g,1}$ , yielding the greedy quadratic ROM (gQM-ROM) of dimension $r=20$ , which has the form (23)–(24) with an identity mass-matrix term since ${\mathbf{W}}_{g}^{\top}{\mathbf{V}}_{g,2}=0$ . The singular value decay of the snapshot data used to construct the greedy ROM projection matrices is shown in bottom left panel of Figure 11. On the other hand, Algorithm 1 is used to construct the proposed interpolatory QM-ROM ( $r=20$ ).

Table 2: Computational wall-clock time and state reconstruction error comparison between the Full-Order Model (FOM), the proposed interpolatory QM-ROM, and the greedy data-driven QM ROM (gQM-ROM) at reduced order

r=20

Method	Offline Time (s)	Online Time (s)	Absolute Error
Full order model (Full)	—	147.0238	—
Proposed interpolatory (QM-ROM)	3.0298	0.7266	$\mathbf{6.196\times 10^{-5}}$
Greedy QM-ROM (gQM-ROM) schwerdtner2024greedy	23.4082	0.4709	$1.622\times 10^{-3}$

Figure 11 illustrates the spatio-temporal dynamics and corresponding absolute reconstruction errors across the simulation window $t\in[0,1]$ . While both reduced models successfully capture the wave propagation, for this example, the proposed interpolatory QM-ROM achieves more than two orders of magnitude higher accuracy ( $\mathcal{O}(10^{-5})$ compared to $\mathcal{O}(10^{-3})$ for gQM-ROM). This might be due to the fact that the snapshot-based approaches are based on minimizing an empirical least-squares state-reconstruction fitting over a finite ensemble of snapshots. While the use of a quadratic manifold significantly relaxes the dependence on slow singular value decay relative to linear subspace methods, namely an order- $r$ quadratic manifold aims to capture information up to approximately the $(r+r^{2})$ -th singular direction rather than just the $r$ -th, the fitting residual at order $r$ is still governed by how quickly the neglected singular values decay beyond that effective order, as shown in the bottom left panel of Figure 11.

On the other hand, for LTI systems, the proposed framework is based on interpolating the steady-state trajectory for a special class of inputs (see Section 5). Consequently, the reduced order $r$ in our setting parameterizes the state dimension of the target exogenous signal generator (and thus the richness of the interpolated input family), instead of truncation rank for snapshot data.

Table 2 details the computational statistics and reconstruction errors for both the proposed interpolatory QM-ROM and the greedy quadratic ROM (gQM-ROM). While the gQM-ROM provides a robust framework that can be tailored via specific hyperparameter selections, the optimization-free nature of the proposed formulation (Algorithm 1) serves as a complementary approach, naturally yielding consistently low offline construction times. During the online phase, both methods share an identical theoretical time complexity of $\mathcal{O}(r^{3})$ due to their quadratic state operators. In practice, the gQM-ROM achieves excellent online execution speeds (as reflected in Table 2) by enforcing an identity mass-matrix structure ( ${\mathbf{M}}_{r}(\mathbf{x}_{r})=\mathbf{0}$ ). Notably, this structural advantage can be integrated into the proposed interpolatory QM-ROM (by picking ${\mathbf{W}}$ such that ${\mathbf{W}}^{\top}{\mathbf{V}}_{2}=0$ ) if comparable runtime efficiency is desired.

Finally, we emphasize that many parameters appear in the optimization-based QM methods, such as snapshot sample density, time horizons, and regularization parameters. Thus, hyperparameter tuning could potentially yield further improvements in the accuracy of gQM-ROM beyond what we obtained here. However, the primary objective of this comparison is not to conduct an exhaustive optimization benchmark, especially considering that gQM-ROM is applicable to general nonlinear dynamical systems rather than the LTI systems we consider here. Moreover, gQM-ROM and other snapshot based quadratic manifold approaches are not generally employed in a projection setting, but rather in a data-driven modeling framework. Our main goal here, instead, is to demonstrate that the proposed interpolatory QM-ROM provides a compelling, system-theoretic alternative to existing data-driven quadratic manifold approaches. Using the theory of nonlinear moment matching, our framework bypasses snapshot collection and hyperparameter tuning, delivering certified system-theoretical guarantees alongside computational savings for this example. Extending our framework to a data-driven setting and general nonlinear dynamical systems remain an ongoing research direction.

7 Conclusions and future directions

We introduced an optimization-free algebraic framework for constructing quadratic manifold-based reduced-order models (ROMs). By moving beyond the inherent limitations of traditional linear subspaces, our approach provides a principled, system-theoretic pathway to overcome the slow decay of Kolmogorov $n$ -widths that typically affects advection- and transport-dominated problems. We established rigorous theoretical guarantees for the proposed framework, proving that the constructed ROM successfully matches the nonlinear moments of the full-order model. Furthermore, we demonstrated that the exact center manifold mapping is preserved, thereby ensuring the asymptotic tracking of steady-state outputs under specific classes of inputs. The theoretical results were validated through numerical experiments on a transport-dominated one-dimensional damped wave equation and one-dimensional advection equation. The numerical results show that the quadratic manifold ROM captures the complex transient dynamics with high fidelity. This reduction directly translates to substantial online and offline computational savings. A natural extension is to generalize this system-theoretic quadratic manifold framework to broader classes of nonlinear or parameterized partial differential equations. Additionally, exploring schemes to optimally select the signal generator parameters will further expand the practical utility of quadratic manifold-based model order reduction.

Acknowledgements

We thank Prof. Benjamin Peherstorfer and Prof. Steffen Werner for their feedback on the manuscript.

Data availability

The MATLAB code for reproducing the plots and comparisons in this paper are available on Github padhi_code.

References

[1] D. Amsallem, M. J. Zahr, and C. Farhat. Nonlinear model order reduction based on local reduced-order bases. International Journal for Numerical Methods in Engineering, 92(10):891–916, 2012.
[2] A. C. Antoulas. Approximation of Large-Scale Dynamical Systems. Adv. Des. Control 6. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2005.
[3] A. C. Antoulas, C. A. Beattie, and S. Gugercin. Interpolatory Methods for Model Reduction. Computational Science & Engineering. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2020.
[4] A. Astolfi. Model reduction by moment matching for linear and nonlinear systems. IEEE Transactions on Automatic Control, 55(10):2321–2336, 2010.
[5] A. Astolfi, G. Scarciotti, J. Simard, N. Faedo, and J. V. Ringwood. Model reduction by moment matching: Beyond linearity a review of the last 10 years. In 2020 59th IEEE conference on decision and control (CDC), pages 1–16. IEEE, 2020.
[6] H. Bai, T. Mylvaganam, and G. Scarciotti. Model reduction for quadratic-bilinear systems using nonlinear moments. In 2022 European Control Conference (ECC), pages 1702–1707. IEEE, 2022.
[7] Z. Bai and D. Skoogh. A projection method for model reduction of bilinear dynamical systems. Linear algebra and its applications, 415(2-3):406–425, 2006.
[8] J. Barnett and C. Farhat. Quadratic approximation manifold for mitigating the Kolmogorov barrier in nonlinear projection-based model order reduction. Journal of Computational Physics, 464:111348, 2022.
[9] P. Benner and T. Breiten. Two-sided projection methods for nonlinear model order reduction. SIAM Journal on Scientific Computing, 37(2):B239–B260, 2015.
[10] P. Benner and T. Breiten. Chapter 6: Model Order Reduction Based on System Balancing, pages 261–295. SIAM, 2017.
[11] P. Benner, P. Goyal, J. Heiland, and I. Pontes Duff. A quadratic decoder approach to nonintrusive reduced-order modeling of nonlinear dynamical systems. PAMM, 23(1):e202200049, 2023.
[12] P. Benner, S. Gugercin, and S. W. Werner. Structured interpolation for multivariate transfer functions of quadratic-bilinear systems. Advances in Computational Mathematics, 50(2):18, 2024.
[13] W.-J. Beyn and V. Thümmler. Freezing solutions of equivariant evolution equations. SIAM Journal on Applied Dynamical Systems, 3(2):85–116, 2004.
[14] F. Black, P. Schulze, and B. Unger. Projection-based model reduction with dynamically transformed modes. ESAIM: Mathematical Modelling and Numerical Analysis, 54(6):2011–2043, 2020.
[15] T. Breiten and P. Benner. Interpolation-based $\mathcal{H}_{2}$ -model reduction of bilinear control system. SIAM J. Matrix Anal. Appl, 33(3):859–885, 2012.
[16] T. Breiten and T. Damm. Krylov subspace methods for model order reduction of bilinear control systems. Systems & Control Letters, 59(8):443–450, 2010.
[17] J. Brewer. Kronecker products and matrix calculus in system theory. IEEE Transactions on Circuits and Systems, 25:772–781, 1978.
[18] P. Buchfink, S. Glas, and B. Haasdonk. Symplectic model reduction of Hamiltonian systems on nonlinear manifolds and approximation with weakly symplectic autoencoder. SIAM Journal on Scientific Computing, 45(2):A289–A311, 2023.
[19] P. Buchfink, S. Glas, and B. Haasdonk. Approximation bounds for model reduction on polynomially mapped manifolds. Comptes Rendus. Mathématique, 362(G13):1881–1891, 2024.
[20] P. Buchfink, S. Glas, B. Haasdonk, and B. Unger. Model reduction on manifolds: A differential geometric framework. Physica D: Nonlinear Phenomena, 468:134299, 2024.
[21] S. Burela, P. Krah, and J. Reiss. Parametric model order reduction for a wildland fire model via the shifted POD based deep learning method. arXiv preprint arXiv:2304.14872, 2023.
[22] J. Carr. Applications of centre manifold theory. Springer Science & Business Media, 2012.
[23] T. Daniel, F. Casenave, N. Akkari, A. Ketata, and D. Ryckelynck. Physics-informed cluster analysis and a priori efficiency criterion for the construction of local reduced-order bases. Journal of Computational Physics, 458:111120, 2022.
[24] G. Flagg and S. Gugercin. Multipoint volterra series interpolation and $\mathcal{H}_{2}$ optimal model reduction of bilinear systems. SIAM Journal on Matrix Analysis and Applications, 36(2):549–579, 2015.
[25] S. Fresca, L. Dede’, and A. Manzoni. A comprehensive deep learning-based approach to reduced order modeling of nonlinear time-dependent parametrized PDEs. Journal of Scientific Computing, 87(2):61, 2021.
[26] S. Fresca and A. Manzoni. POD-DL-ROM: Enhancing deep learning-based reduced order models for nonlinear parametrized PDEs by proper orthogonal decomposition. Computer Methods in Applied Mechanics and Engineering, 388:114181, 2022.
[27] K. Gallivan, A. Vandendorpe, and P. Van Dooren. Sylvester equations and projection-based model reduction. Journal of Computational and Applied Mathematics, 162(1):213–229, 2004.
[28] R. Geelen, L. Balzano, and K. Willcox. Learning latent representations in high-dimensional state spaces using polynomial manifold constructions. In 2023 62nd IEEE Conference on Decision and Control (CDC), pages 4960–4965. IEEE, 2023.
[29] R. Geelen, L. Balzano, S. Wright, and K. Willcox. Learning physics-based reduced-order models from data using nonlinear manifolds. Chaos: An Interdisciplinary Journal of Nonlinear Science, 34(3), 2024.
[30] R. Geelen and K. Willcox. Localized non-intrusive reduced-order modelling in the operator inference framework. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 380(2229), 2022.
[31] R. Geelen, S. Wright, and K. Willcox. Operator inference for non-intrusive model reduction with quadratic manifolds. Computer Methods in Applied Mechanics and Engineering, 403:115717, 2023.
[32] S. Glas and H. Mu. Structure-preserving model reduction on manifolds of port-Hamiltonian systems. arXiv preprint arXiv:2603.08656, 2026.
[33] I. V. Gosea and A. C. Antoulas. Data-driven model order reduction of quadratic-bilinear systems. Numerical Linear Algebra with Applications, 25(6):e2200, 2018.
[34] A. Graham. Kronecker products and matrix calculus with applications. Courier Dover Publications, 2018.
[35] C. Greif and K. Urban. Decay of the Kolmogorov $n$ -width for wave problems. Applied Mathematics Letters, 96:216–222, 2019.
[36] C. Gu. Qlmor: A projection-based nonlinear model order reduction approach using quadratic-linear representation of nonlinear systems. IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 30(9):1307–1320, 2011.
[37] J. S. Hesthaven, B. Peherstorfer, and B. Unger. Nonlinear model reduction for transport-dominated problems. Acta Numerica, 35:173–272, 2026.
[38] J. Huang. Nonlinear output regulation: theory and applications. SIAM, 2004.
[39] T. Kadeethum, F. Ballarin, Y. Choi, D. O’Malley, H. Yoon, and N. Bouklas. Non-intrusive reduced order modeling of natural convection in porous media using convolutional autoencoders: comparison with linear subspace techniques. Advances in Water Resources, 160:104098, 2022.
[40] Y. Kim, Y. Choi, D. Widemann, and T. Zohdi. A fast and accurate physics-informed neural network reduced order model with shallow masked autoencoder. Journal of Computational Physics, 451:110841, 2022.
[41] R. Klein, B. Sanderse, P. Costa, R. Pecnik, and R. Henkes. Entropy-stable model reduction of one-dimensional hyperbolic systems using rational quadratic manifolds. Journal of Computational Physics, 528:113817, 2025.
[42] P. Krah, A. Marmin, B. Zorawski, J. Reiss, and K. Schneider. A robust shifted proper orthogonal decomposition: Proximal methods for decomposing flows with multiple transports. SIAM Journal on Scientific Computing, 47(2):A633–A656, 2025.
[43] A. Moreschini, M. Scandella, A. Astolfi, and T. Parisini. Moment matching by kernel-based learning. IEEE Transactions on Automatic Control, 2025.
[44] M. Ohlberger and S. Rave. Nonlinear reduced basis approximation of parameterized evolution equations via the method of freezing. Comptes Rendus Mathematique, 351(23-24):901–906, 2013.
[45] S. E. Otto, G. R. Macchio, and C. W. Rowley. Learning nonlinear projections for reduced-order modeling of dynamical systems using constrained autoencoders. Chaos: An Interdisciplinary Journal of Nonlinear Science, 33(11), 2023.
[46] R. Padhi. Beyond linear subspaces: Nonlinear moment matching meets quadratic manifolds. https://github.com/Rewtus/Moment-matching-quadratic-manifolds, 2026.
[47] D. Papapicco, N. Demo, M. Girfoglio, G. Stabile, and G. Rozza. The neural network shifted-proper orthogonal decomposition: a machine learning approach for non-linear reduction of hyperbolic equations. Computer Methods in Applied Mechanics and Engineering, 392:114687, 2022.
[48] G. Paxton, S. Cheon, R. Geelen, and S. A. McQuarrie. Fast quadratic manifold learning for nonlinear dimensionality reduction in large-scale systems using Riemannian optimization. arXiv preprint arXiv:2605.26039, 2026.
[49] B. Peherstorfer. Breaking the Kolmogorov barrier with nonlinear model reduction. Notices of the American Mathematical Society, 69(5):725–733, 2022.
[50] A. Pinkus. $N$ -width in Approximation Theory. Springer Science & Business Media, 2012.
[51] J. Reiss. Optimization-based modal decomposition for systems with multiple transports. SIAM Journal on Scientific Computing, 43(3):A2079–A2101, 2021.
[52] J. Reiss, P. Schulze, J. Sesterhenn, and V. Mehrmann. The shifted proper orthogonal decomposition: A mode decomposition for multiple transport phenomena. SIAM Journal on Scientific Computing, 40(3):A1322–A1344, 2018.
[53] M. Rewienski and J. White. A trajectory piecewise-linear approach to model order reduction and fast simulation of nonlinear circuits and micromachined devices. IEEE Transactions on computer-aided design of integrated circuits and systems, 22(2):155–170, 2003.
[54] C. W. Rowley, I. G. Kevrekidis, J. E. Marsden, and K. Lust. Reduction and reconstruction for self-similar dynamical systems. Nonlinearity, 16(4):1257–1275, 2003.
[55] C. W. Rowley and J. E. Marsden. Reconstruction equations and the Karhunen–Loève expansion for systems with symmetry. Physica D: Nonlinear Phenomena, 142(1-2):1–19, 2000.
[56] W. J. Rugh. Nonlinear system theory. Johns Hopkins University Press Baltimore, 1981.
[57] J. B. Rutzmoser, D. J. Rixen, P. Tiso, and S. Jain. Generalization of quadratic manifolds for reduced order modeling of nonlinear structural dynamics. Computers & Structures, 192:196–209, 2017.
[58] G. Scarciotti and A. Astolfi. Data-driven model reduction by moment matching for linear and nonlinear systems. Automatica, 79:340–351, 2017.
[59] G. Scarciotti and A. Astolfi. Nonlinear model reduction by moment matching. Foundations and Trends in System and Control, 4(3-4):224–409, 2017.
[60] G. Scarciotti and A. Astolfi. Interconnection-based model order reduction-a survey. European Journal of Control, 75:100929, 2024.
[61] P. Schwerdtner, S. Gugercin, and B. Peherstorfer. Empirical sparse regression on quadratic manifolds. SIAM Journal on Scientific Computing, 47(6):A3085–A3107, 2025.
[62] P. Schwerdtner, P. Mohan, A. Pachalieva, J. Bessac, D. O’Malley, and B. Peherstorfer. Online learning of quadratic manifolds from streaming data for nonlinear dimensionality reduction and nonlinear model reduction. arXiv preprint arXiv:2409.02703, 2024.
[63] P. Schwerdtner, P. Mohan, A. Pachalieva, J. Bessac, D. O’Malley, and B. Peherstorfer. Online learning of quadratic manifolds from streaming data for nonlinear dimensionality reduction and nonlinear model reduction. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 481(2314), 2025.
[64] P. Schwerdtner and B. Peherstorfer. Greedy construction of quadratic manifolds for nonlinear dimensionality reduction and nonlinear model reduction. SIAM Journal on Mathematics of Data Science, 8(3):793–819, 2026.
[65] H. Sharma, H. Mu, P. Buchfink, R. Geelen, S. Glas, and B. Kramer. Symplectic model reduction of Hamiltonian systems using data-driven quadratic manifolds. Computer Methods in Applied Mechanics and Engineering, 417:116402, 2023.
[66] J. D. Simard, A. Moreschini, and A. Astolfi. Parameterization of all differential-algebraic moment matching interpolants. IEEE Transactions on Automatic Control, 70(3):1875–1882, 2024.
[67] B. Unger and S. Gugercin. Kolmogorov $n$ -widths for linear dynamical systems. Advances in Computational Mathematics, 45(5):2273–2286, 2019.
[68] S. Volkwein. Model reduction using proper orthogonal decomposition. Lecture notes, Institute of Mathematics and Scientific Computing, University of Graz, 2011.
[69] S. W. Werner. Structure-preserving model reduction for mechanical systems. PhD thesis, Otto-von-Guericke-Universität Magdeburg, Fakultät für Mathematik, 2021.