Multi-step EDMD for Nonlinear Systems

• Multi-step EDMD is a data-driven framework that leverages Koopman operator theory to approximate nonlinear dynamics via linear representations in a lifted function space.
• It minimizes long-horizon errors by directly learning multi-step state outputs, avoiding error compounding seen in recursive single-step methods.
• The framework employs advanced dictionary learning and sparsity techniques, yielding robust surrogates applicable to density forecasting, system identification, and model predictive control.

The multi-step Extended Dynamic Mode Decomposition (EDMD) framework is a family of data-driven numerical methods for forecasting, modeling, and controlling nonlinear dynamical systems via linear approximations in function space. By leveraging action of the Koopman operator on a finite set of observables, multi-step EDMD enables high-fidelity multi-horizon prediction and robust surrogate model construction for stochastic, deterministic, and controlled systems. Key advantages over classical single-step EDMD include direct minimization of long-horizon errors, avoidance of error compounding under recursive prediction, and improved tractability for real-time and control applications.

1. Mathematical Foundations and Koopman Lifting

The multi-step EDMD framework generalizes the spectral approach of Koopman operator theory to nonlinear systems, including both deterministic maps and stochastic differential equations (SDEs) with diffusion. Consider a continuous-time SDE

$dX_t = b(X_t)\,dt + \sigma(X_t)\,dW_t$

with initial distribution $X_0 \sim p_0(x)$, drift $b:\mathbb{R}^d\to\mathbb{R}^d$, and noise $\sigma:\mathbb{R}^d\to\mathbb{R}^{d\times s}$, or a discrete-time system $x_{t+1}=f(x_t, u_t)$. The associated stochastic Koopman semigroup $\{K^t\}$ acts on observables $\psi$ by $K^t\psi(x) = \mathbb{E}[\psi(X_t^x)]$, defining a linear semigroup generated by the infinitesimal generator

$$\mathcal{L}v(x) = b(x)\cdot\nabla v(x) + \frac{1}{2}\Sigma(x):\nabla^2 v(x),\quad \Sigma(x)=\sigma(x)\sigma(x)^\top,$$

with Fokker–Planck adjoint $$\mathcal{L}^*p = -\nabla\cdot[b\,p] + \frac{1}{2}\nabla\cdot\nabla\cdot[\Sigma p]$$ (Zhao et al., 2022). Under ergodicity, there exists a unique invariant density $p_*(x)$, and the $L^2(\mathbb{R}^d; p_*)$ Hilbert space inner product governs the spectral structure of $K^t$ and $\mathcal{L}$.

Deterministic systems use analogous constructions, with observables $\{\psi_j\}_{j=1}^N$ spanning a finite-dimensional subspace, and evolution in this lifted space governed by approximations of the infinite-dimensional Koopman operator (Schurig et al., 10 Nov 2025, Meda et al., 4 Apr 2025).

2. Finite-Dimensional EDMD Approximations and Multi-Step Forecasting

The core idea of EDMD is to approximate $K^{\Delta t}$ by a finite-dimensional matrix acting on the lifted state, using empirical data:

• Snapshot pairs $\{(x_k, y_k)\}_{k=1}^M$ are generated via dynamics and possible control inputs.
• The observables define the lift $\Psi(x)=[\psi_1(x),\dots,\psi_m(x)]^\top$; the data matrices are $\Psi(X), \Psi(Y)\in\mathbb{R}^{m\times M}$.
• Gram matrices are $G = \frac{1}{M}\sum \Psi(x_k)\Psi(x_k)^\top$, $A = \frac{1}{M}\sum \Psi(y_k)\Psi(x_k)^\top$.
• The Koopman approximation is $K = A\,\mathrm{pinv}(G)$, with eigenpairs $(\mu_j, \xi_j)$ giving approximate eigenvalues and eigenfunctions (Zhao et al., 2022, Meda et al., 4 Apr 2025).

Multi-step prediction proceeds by recursive iteration:

$$z_{k+1} = Kz_k + Bu_k,\quad \hat{x}_{k+1} = C z_{k+1},\quad \hat{y}_{k+1} = E z_{k+1},$$

where the lifted state $z_k = \Psi(x_k)$, $C$ selects physical states from the lifted coordinates, and $E$ recovers outputs or measurements (Meda et al., 4 Apr 2025). In stochastic density forecasting, spectral decompositions using $K$ generate multi-step probability evolution via $K^n$ or eigenmode propagation (Zhao et al., 2022).

3. Dictionary Design, Manifold Optimization, and Sparsity

The expressiveness and generalization of EDMD surrogates depend critically on the choice of observables. Standard bases include monomials, radial basis functions (RBFs), and neural networks (Meda et al., 4 Apr 2025). Dictionary learning is further refined by geometric optimization on the Grassmann manifold: selecting a subspace $V=\mathrm{span}\{\psi_1,\dots,\psi_M\}$ that minimizes multi-step forecast error via Riemannian optimization yields approximately invariant, low-dimensional subspaces that are robust to out-of-domain states (Schurig et al., 10 Nov 2025). The error metric integrates projection error over finite horizons and test initial conditions, and the optimization descends to Grassmannian geometry with Riemannian gradients and QR-based retraction.

Explicit dictionary pruning and structure discovery are enabled by including $\ell_1$-type regularization in least-squares identification, which removes irrelevant observables and yields parsimonious, efficient surrogates (Wu et al., 17 Jan 2026). Parallel decomposition across states, steps, and dictionary rows supports scalable computation.

4. Multi-Step Least-Squares Identification and Error Control

A key advance of the multi-step EDMD framework is direct learning of the condensed multi-step state-output map. Rather than identifying an operator $A$ for single-step propagation and composing $A^k$ (which leads to error compounding and potential instability for $\Vert A\Vert_2>1$), the framework fits for each prediction step $k$ a linear map from dictionary lifts and past inputs to the $k$-step-ahead state:

$$x_{j,k} = E_k \psi(x_{j,0}) + F_k [u_{j,0};\dots;u_{j,k-1}] + \mathrm{residual}.$$

The identification is performed via convex least-squares:

$$\min_{E_k,F_k} \sum_{j=1}^M \|x_{j,k} - E_k\psi(x_{j,0}) - F_k[u_{j,0};\dots;u_{j,k-1}]\|_2^2 + \beta\|[E_{k};F_{k}]\|_2^2 + \tau\|E_{k}\|_1,$$

with row-wise $\ell_1$ sparsity for dictionary pruning (Wu et al., 17 Jan 2026). This produces stable, direct multi-step surrogates whose error does not grow with prediction horizon $k$. Theorem 3 in (Wu et al., 17 Jan 2026) provides, under Sobolev regularity and bounded-dictionary assumptions, high-probability error bounds independent of $k$:

$$\|f_i^{(k)} - E_{k,i}^\top\psi\|_H \leq C_{ms}(p^{-s}\|f_i^{(k)}\|_{H^s} + \sqrt{N^2\log M/M}),$$

where $N$ is the dictionary size, $M$ the number of data, and $p$ the degree for polynomial lifts. In contrast, one-step EDMD bounds exhibit exponential-in-$k$ error blow-up when $\|K\|_H > 1$.

5. Applications, Numerical Performance, and Practical Considerations

Multi-step EDMD has been demonstrated for prediction, nonlinear system identification, density forecasting, and model predictive control (MPC) across a wide range of systems:

• Density forecast for SDEs: Approximates the semigroup generated by the Fokker–Planck operator, achieving weak convergence of the density representation as data and dictionary sizes increase (Zhao et al., 2022).
• Nonlinear system surrogates: For electric vehicle (EV) cabin climate, multi-step EDMD with RBF dictionaries achieves state RMSE $<5\%$ and power RMSE $<3\%$ for $N=35$ basis, outperforming polynomial and neural-network dictionaries at small $N$ (Meda et al., 4 Apr 2025).
• Koopman MPC: Direct multi-step identification solves the full condensed map for state and control over control horizon $T$, enabling stable closed-loop control where one-step EDMD leads to diverging errors and controller instability. In benchmarks, multi-step EDMD achieves stable MSE under open and closed loop for $k\gg 10$, while one-step MSE diverges for $k>10$ due to spectral radius $>1$ (Wu et al., 17 Jan 2026).
• Dictionary learning: Grassmannian-based shaping yields low-dimensional surrogates with out-of-domain generalization and large speedups in predictor evaluation (Schurig et al., 10 Nov 2025).

Per-step computational costs are governed by $O(N^2)$ matrix-vector products and $O(N\phi)$ dictionary evaluation (with $\phi$ cost per basis). For $N\sim 30$–$40$, real-time execution on embedded platforms is practical (Meda et al., 4 Apr 2025).

6. Algorithmic and Implementation Summary

A generic multi-step EDMD workflow proceeds as:

1. Dictionary selection: Choose functions $\{\psi_j\}$—polynomial, RBF, or learned—appropriate to system structure and data availability.
2. Data acquisition: Collect snapshot (state, next-state) or (trajectory) data under representative state and control distributions.
3. Lifting and matrix assembly: Formulate $\Psi(X)$, $\Psi(Y)$ and associated Gram matrices; assemble multi-step input-output structures for supervised learning.
4. Identification: Solve (regularized) least-squares for operator matrices, possibly under parallelizable decomposition across step $k$ and state $i$; perform dictionary pruning as necessary.
5. Prediction: For given initial condition and input sequence, propagate via

$$z_{k+1} = K z_k + B u_k\,,\quad\text{or}\quad x_{k} = E_k \psi(x_0) + F_k [u_0;\ldots;u_{k-1}]$$

up to desired horizon, reconstruct state and output via $C$ and $E$ matrices.

1. Validation: Assess with trajectory RMSE, Consistency Index (CI), or finite-horizon forecast errors; compare performance with alternative dictionaries or identification schemes (Zhao et al., 2022, Meda et al., 4 Apr 2025, Schurig et al., 10 Nov 2025, Wu et al., 17 Jan 2026).

When relevant (e.g., Grassmannian dictionary learning), implement Riemannian optimization to shape active observables for improved multi-step invariance and forecast accuracy (Schurig et al., 10 Nov 2025). Closed-loop extensions with linear MPC/LQR or online adaptation using recursive least squares are direct, leveraging the lifted linear structure of the surrogate.

7. Convergence Guarantees and Limitations

Convergence properties are established in terms of increasing data ($M\to\infty$), basis set expressiveness ($m,N\to\infty$), and regularity of the lifting functions. For stochastic systems, EDMD approximations converge to the Galerkin projection of the Fokker–Planck semigroup in $L^2(p_\ast)$, and the truncated density forecast $p_m(x,t)$ converges weakly to the true solution (Zhao et al., 2022). In deterministic and control settings, multi-step EDMD avoids error blow-up inherent to one-step recursion, with finite-time sample error scaling optimally in dictionary size and sample count (Wu et al., 17 Jan 2026). Overfitting and lack of generalization can arise with excessive dictionary dimension, motivating use of Grassmannian shaping and sparsity regularization (Schurig et al., 10 Nov 2025, Wu et al., 17 Jan 2026).

A plausible implication is that careful system-specific dictionary design, appropriate regularization, and direct multi-horizon identification are essential for the robust practical application of EDMD in high-dimensional settings and real-time control scenarios.