4 Moment Matching Method for Model Reduction

• 4 Moment Matching Method is a technique that aligns the steady-state output moments of complex systems with those of simplified models using polynomial basis functions.
• It employs a Galerkin residual approach to solve invariance equations, ensuring the reduced-order model accurately reproduces the original system's response.
• Validated on systems up to 1000 dimensions, the method offers efficient computation and high-fidelity steady-state matching for nonlinear dynamic systems.

The 4 Moment Matching Method encompasses a spectrum of mathematical and algorithmic techniques for constructing, estimating, or reducing models by enforcing the agreement of moments (expectations of various orders) between modeled (or empirical) distributions and the target or observed distributions. In the context of nonlinear model reduction, moment matching refers to matching the steady-state or asymptotic output response of the reduced-order model with that of the original high-dimensional system when driven by a specified class of signal generators. Recent developments leverage polynomial approximation schemes and Galerkin residual techniques to approximate solutions to the invariance equations that characterize these moments, enabling efficient reduced-order modeling even at large system dimensions (Doebeli et al., 17 Dec 2024).

1. Moment Matching Technique: Concept and Relevance to Model Reduction

The foundational idea of moment matching in model reduction is to parametrize the steady-state output response (the “moment”) of a full-order system in terms of its behavior when excited by an exogenous nonlinear or linear signal generator. Formally, if the full-order system is given as $\dot{x} = f(x,u)$, $y = h(x)$, and a signal generator as $\dot{\omega} = s(\omega)$, $v = \ell(\omega)$, with the interconnection $u = v$, then moments are defined via a mapping $\pi(\omega)$ such that $x_\infty = \pi(\omega)$ describes the invariant manifold on which the system output tracks the signal generator trajectory in stationarity.

For high-dimensional systems ($n$ up to 1000), the goal is to identify a polynomial mapping $\pi: \mathbb{R}^d \to \mathbb{R}^n$ that reduces the complex dynamics to an invariant subspace while reproducing the same steady-state output. By ensuring that the reduced-order model produces outputs $y_r = h(\pi(r))$ matching those of the original system for the same excitation, the essential dynamics and response characteristics are preserved even with a dramatically lower-dimensional state representation (Doebeli et al., 17 Dec 2024).

2. Polynomial Approximation and Numerical Scheme

The method approximates the invariant mapping $\pi(\omega)$ by expanding each of its $n$ components in a polynomial basis over the generator state domain $\Omega \subset \mathbb{R}^d$. Specifically, the approximation takes the form

$$\pi_i^N(\omega) = \sum_{k=1}^N c_{i,k} \phi_k(\omega),$$

where the $\{\phi_k\}$ are a (tensor) basis of degree $M$ polynomials or other complete function families. The choice of a global polynomial basis—especially in low $d$—permits the exact or quadrature-based evaluation of high-order moment integrals that arise in the Galerkin projection.

To solve the coupling introduced by inserting this expansion into the invariance equations, a Galerkin residual method is employed. The residuals—differences between the derivative and the vector field—are projected onto the basis: $$R_i(\omega) = \frac{\partial \pi_i^N}{\partial \omega} s(\omega) - f_i(\pi^N(\omega), \ell(\omega));$$

$$\langle R_i, \phi_j \rangle_{\Omega} = \int_\Omega R_i(\omega) \phi_j(\omega) d\omega = 0 \quad \forall j.$$

This results in a system of nonlinear algebraic equations in the coefficients $c_{i, k}$, which is solved using Newton's method, with the Jacobian evaluated as needed (using a Moore–Penrose pseudoinverse when singular). This approach is robust even for large $n$ due to the block structure and tensorization available for polynomial representations.

3. Invariance Equations and Their Role in Moment Matching

The invariance equations define the mapping $\pi$ characterizing the locally attractive center manifold connecting the full system and the signal generator. For the interconnected system: $$\dot{x} = f(x, \ell(\omega)), \qquad \dot{\omega} = s(\omega),$$ the invariance equation is

$$\frac{\partial \pi}{\partial \omega} s(\omega) = f(\pi(\omega), \ell(\omega)), \qquad \pi(0) = 0.$$

In the linear case, this reduces to a Sylvester equation: $\Pi S = A\Pi + BL.$ The solution $\pi(\omega)$ parametrizes the pullback of the generator’s state onto the center manifold of the full-order system, ensuring that $x(t) = \pi(\omega(t))$ and the steady-state output matches for all initial conditions on the manifold.

By numerically solving these invariance PDEs, the polynomial approximation scheme achieves a reduced-order surrogate that exactly matches the system’s moment (steady-state response) for arbitrary generator excitations within $\Omega$.

4. Construction and Validation of Reduced-Order Models

The approximate mapping $\pi^N(\omega)$ is used to construct ROMs by embedding it as the core of the reduced dynamics, e.g.,

$$\dot{r} = s(r) - \bar{g}(r, \ell(r)) + \bar{g}(r, u), \qquad y_r = h(\pi^N(r)),$$

with $r \in \mathbb{R}^d$ directly representing the generator state and $\bar{g}$ selected to ensure reduced model stability. This construction ensures that the reduced-order model matches the output behavior of the original for all signal generator trajectories.

Empirical validation includes the computation of the L2 norm of the residuals of the invariance equation, as well as comparison of steady-state output trajectories between the ROM and the high-dimensional system. For low $n$ and low-degree polynomials, residual errors close to machine precision are observed (e.g., $10^{-16}$ in some 2-state test cases), and for large $n$ (e.g., $n=1000$) residuals remain small and decrease with increase in polynomial degree.

5. Numerical Results and Applicability

Extensive numerical results demonstrate the scalability and accuracy of the method for both linear and nonlinear high-dimensional systems. For instance, in a nonlinear RL ladder circuit with $n=2$ to $n=1000$ and either linear or Van der Pol signal generators, the method efficiently recovers the invariant mapping $\pi$ and associated ROMs that faithfully reproduce the full model's steady-state output. For moderate polynomial degree ($M \sim 5$), computation times are practical (from several seconds for $n=2$ to under 2 hours for $n=1000$ with $d=2$).

Key practical findings include:

• The polynomial Galerkin method enables precision matching of moments even at large scale.
• The use of a global polynomial basis (tensor monomials or orthogonal polynomials) interacts well with quadrature and tensor contraction techniques, bypassing the need for gridding or discretizing the state space.
• The reduced models constructed in this way fail to capture only transient dynamics not on the center manifold, but strongly match the steady-state output signal response, as per the invariance equation predictions.

6. Applications, Scalability, and Advantages

The method addresses prominent challenges in nonlinear model reduction for high-dimensional systems—where full identification or simulation is intractable—by focusing on steady-state, or moment, agreement. Applications include:

• Large-scale electrical circuits, where reduced models are needed for simulation and design, but only input-output steady-state behavior is relevant.
• Nonlinear mechanical systems (e.g., cart-pendulum, RL-ladders) where model order reduction is required for controller synthesis or simulation of nonlinear oscillatory behaviors.
• Problems involving nonlinear or periodic excitation of highly over-parameterized systems, where the moment matching approach directly yields low-rank surrogates for long-term system identification and prediction.

The primary advantage lies in the direct attack on the invariance equation using spectral (polynomial) expansion as opposed to data-driven fitting or state-space identification. This not only provides theoretical guarantees of “moment fidelity” but also enables algorithmic tractability, as exhibited for $n$ in the hundreds or thousands. The method’s framework also allows flexible generalization to settings with alternative generator classes, other basis functions, or approximation methods.

7. Limitations and Further Developments

While highly effective for matching steady-state moments, moment matching cannot guarantee fidelity for transient or strongly nonlinear non-steady-state input responses not spanned by the generator excitation space. Additionally, the method’s efficiency hinges on the dimensionality of the generator ($d$)—for large $d$, the volume of the function space and complexity of the expansion grows rapidly. Further improvements could include:

• Use of adaptive or sparse polynomial bases to mitigate curse-of-dimensionality effects for very high $d$.
• Integration with machine learning approaches for “data-driven” basis selection and residual minimization.
• Extension to operator-valued moments for transfer function matching when the system is excited across a broad range of frequencies.

This polynomial approximation scheme for moment matching thus provides an effective and scalable approach to nonlinear model reduction, supporting accurate ROM construction for large-scale dynamical systems while maintaining strict fidelity to prescribed steady-state responses (Doebeli et al., 17 Dec 2024).