Constructive Finite-Rank Koopman Control via ESL Theory: Real-Time Nonlinear Control with Provable Stability Guarantees

1. Introduction

Thesis. We show that ESL’s streamability guarantees yield constructive, finite-rank Koopman surrogates that can be trained and deployed in a single streaming loop, with per-step stability projection and O(d r²) amortized cost.

The Koopman operator framework has emerged as a transformative approach to nonlinear control by lifting dynamics to infinite-dimensional spaces where they evolve linearly [1,2,3]. This linearization enables the direct application of mature linear control theory to inherently nonlinear systems, offering significant theoretical and computational advantages. However, the infinite-dimensional nature of the Koopman operator presents fundamental challenges for practical implementation, particularly in real-time control applications where computational resources are constrained and stability guarantees are paramount.

The Extended Dynamic Mode Decomposition (EDMD) method addresses this challenge by projecting the infinite-dimensional Koopman operator onto a finite-dimensional subspace spanned by a chosen dictionary of observables [4,5]. While EDMD has demonstrated remarkable empirical success across diverse applications, the theoretical foundations underlying finite-rank approximations remain incomplete. Critical questions persist regarding when such approximations exist, how to construct them systematically, what error bounds can be guaranteed, and how controllers designed on approximate models perform when deployed on true plants.

We resolve these fundamental questions by establishing a rigorous connection between Koopman operator theory and the recently developed $\epsilon$-Streamably-Learnable (ESL) hierarchy [6,7]. The ESL framework provides a complete classification of optimization problems based on the existence and structure of $\alpha$-averaged operators, with three distinct classes: non-learnable problems that admit no convergent iterative methods, learnable problems that admit convergent batch algorithms, and streamable problems that additionally possess uniform low-rank structure enabling efficient streaming implementations.

The key insight underlying our approach is that dynamical systems amenable to Koopman control correspond precisely to learnable problems in the ESL sense, while finite-rank Koopman approximations correspond to streamable problems. This connection is not merely analogical but mathematically precise: we prove that learnability in EDMD coordinates is equivalent to contractivity of the finite Koopman operator, and that the ESL density theorem guarantees the existence of arbitrarily accurate streamable approximations for any learnable system.

1.1. Main Contributions

1. Existence → Constructiveness. From ESL to practice: if a system is Koopman-learnable ($\rho$(K)<1 in lifted coordinates), then there exist ESL-streamable low-rank surrogates; we instantiate them with Oja/PAST+TSQR updates and state explicit residual→tube bounds.

2. Single-loop control. A unified ASKC loop that (i) updates a low-rank Koopman model, (ii) runs a one-step Kleinman/Riccati gain update, and (iii) enforces a per-step contraction via a tiny QP/projection—no train→deploy phases.

3. Provable safety at streaming cost. ISS bounds under model drift, stability transfer via Bauer–Fike, and dynamic-regret flavor guarantees, all at O(d r²) per step (r≪d).

1.2. Related Work and Positioning

Our work builds upon and extends several distinct research threads in control theory and machine learning. The Koopman operator framework was introduced by Koopman in 1931 [1] and has experienced renewed interest following the development of data-driven approximation methods [2,9]. The Dynamic Mode Decomposition (DMD) algorithm [10] and its extensions, particularly EDMD [4], have enabled practical applications of Koopman theory to complex nonlinear systems.

Recent advances in Koopman-based control include the development of Model Predictive Control formulations [5,11], stability analysis for closed-loop systems [12], and applications to robotics and aerospace systems [13,14]. However, these works typically rely on heuristic finite-rank approximations without rigorous theoretical justification for the approximation quality or stability guarantees.

The ESL hierarchy, introduced in our prior work [6,7], provides a novel framework for understanding the computational complexity of optimization problems through the lens of operator theory. The hierarchy has been successfully applied to neural network training [15], federated learning, and streaming optimization, but its connection to dynamical systems and control theory has not been previously explored.

Our contribution lies in recognizing and formalizing the deep connection between these seemingly disparate frameworks, leading to both theoretical insights and practical algorithms that advance the state of the art in data-driven control.

1.3. Paper Organization

The remainder of this paper is organized as follows. Section 2 establishes the mathematical foundations, including the ESL hierarchy, Koopman operator theory, and their interconnections. Section 3 presents our main theoretical results, including the learnability-contractivity equivalence and the constructive density theorem. Section 4 develops the streaming algorithm for real-time Koopman learning with detailed convergence analysis. Section 5 provides comprehensive stability analysis for controllers designed on finite surrogates, including ISS bounds and robust MPC formulations. Section 6 presents implementation algorithms and computational complexity analysis. Section 7 demonstrates the framework through extensive numerical experiments on canonical control problems. Section 8 discusses limitations, failure modes, and future research directions. Section 9 concludes with a summary of contributions and broader implications.

2. Mathematical Foundations

This section establishes the mathematical framework underlying our approach, including the ESL hierarchy, Koopman operator theory, and their fundamental interconnections. We begin with the operator-theoretic foundations that enable our analysis, then develop the specific structures needed for dynamical systems and control applications.

2.1. The $\epsilon$-Streamably-Learnable Hierarchy

The ESL hierarchy provides a complete classification of optimization problems based on the existence and structure of $\alpha$-averaged operators. This classification has profound implications for computational complexity and algorithmic design, as it precisely characterizes when efficient streaming algorithms exist and when they can achieve performance comparable to batch methods.

Definition 1

($\alpha$-Averaged Operator). Let $(V,\langle ·,·\rangle )$ be a finite-dimensional Hilbert space and $D\subset V$ be a nonempty, closed, and convex set. An operator $T:D\to D$ is α-averaged for $\alpha \in (0,1)$ if there exists a nonexpansive operator $S:D\to D$ such that

$$T=(1-\alpha )I+\alpha S,$$

(1)

where I denotes the identity operator on V.

The significance of $\alpha$-averaged operators lies in their guaranteed convergence properties. By the Krasnosel’skiĭ-Mann theorem [8], the fixed-point iteration ${\theta }_{k+1}=T\left({\theta }_k\right)$ converges weakly to a fixed point of T for any starting point ${\theta }_0\in D$, provided that $\mathrm{Fix}\left(T\right)\ne \varnothing$. This convergence property forms the foundation for defining learnable problems.

Definition 2

(Problem Classification in the ESL Hierarchy). Consider an optimization problem ${min}_{\theta \in D}R\left(\theta \right)$ where $R:D\to \mathbb{R}$ is the objective function and $D\subset V$ is the feasible set. The problem is classified as:

• Non-learnable  if no α-averaged operator $T:D\to D$ exists such that the fixed-point iteration ${\theta }_{k+1}=T\left({\theta }_k\right)$ converges to a solution of the optimization problem.
• Learnable  if there exists an α-averaged operator $T:D\to D$ with fixed point ${\theta }^{*}\in {argmin}_{\theta \in D}R\left(\theta \right)$ such that the iteration ${\theta }_{k+1}=T\left({\theta }_k\right)$ converges to ${\theta }^{*}$ for all ${\theta }_0\in D$.
• Streamable at rank K  if the problem is learnable and the residual mapping $r\left(\theta \right):=\theta -T\left(\theta \right)$ satisfies the uniform low-rank approximation property:

  $$\underset{\theta \in D}{sup}∥r\left(\theta \right)-\sum _{k=1}^K{g}_k\left(\theta \right)\otimes h_k{\left(\theta \right)}^{*}∥\le {\epsilon }_K$$

  (2)
  for some bounded maps $g_k,h_k:D\to V$ and approximation error ${\epsilon }_K\ge 0$.

The hierarchy establishes strict inclusions: $\mathrm{Streamable}⊊\mathrm{Learnable}⊊\mathrm{All}\mathrm{Problems}$. The streamable class is particularly important for our purposes, as it precisely characterizes when efficient low-rank algorithms exist with provable performance guarantees.

ESL Bridge to Koopman Theory: ESL classifies problems by the existence of averaged operators and their uniform low-rank residuals. In lifted Koopman coordinates, learnability corresponds to a strict contraction ($\rho$(K)<1); streamability ensures the residual admits separable, low-rank updates, enabling our Oja/PAST+TSQR loop with explicit residual→tube bounds. (Proof sketch → Appendix A.)

Definition 3

(Uniform Residual Metric). For two learnable problems with α-averaged operators $T_1,T_2:D\to D$ and corresponding residual mappings $r_1\left(\theta \right)=\theta -T_1\left(\theta \right)$ and $r_2\left(\theta \right)=\theta -T_2\left(\theta \right)$, define the uniform residual distance as:

$$d_R(T_1,T_2)=\underset{\theta \in D}{sup}\parallel r_1\left(\theta \right)-r_2\left(\theta \right)\parallel .$$

(3)

This metric enables us to formalize the notion of approximation quality between different operators and provides the foundation for our density results.

Theorem 1

(Density of Streamable Problems). Streamable problems are dense in learnable problems under the uniform residual metric. Specifically, for every learnable problem with α-averaged operator $T_0:D\to D$ and every $\epsilon >0$, there exists a streamable problem with α-averaged operator $T_{\epsilon }:D\to D$ such that $d_R(T_0,T_{\epsilon })<\epsilon$.

The proof of this fundamental result is constructive and provides an explicit algorithm for converting any learnable problem into a streamable approximation. This constructive nature is crucial for our applications to Koopman control, as it guarantees that finite-rank approximations can always be found with controllable error.

2.2. Koopman Operator Theory for Dynamical Systems

We now turn to the Koopman operator framework for analyzing nonlinear dynamical systems. Consider a discrete-time nonlinear control system of the form:

$$x_{t+1}=f(x_t,u_t),$$

(4)

where $x_t\in \mathcal{X}\subset {\mathbb{R}}^n$ represents the state at time t, $u_t\in \mathcal{U}\subset {\mathbb{R}}^m$ is the control input, and $f:\mathcal{X}\times \mathcal{U}\to \mathcal{X}$ is the system dynamics.

Assumption 1

(System Regularity). The state space $\mathcal{X}$ is compact and forward-invariant under the dynamics (4) for all admissible controls $u_t\in \mathcal{U}$. The function f is locally Lipschitz continuous in both arguments, uniformly over the compact sets $\mathcal{X}\times \mathcal{U}$.

The Koopman operator approach linearizes the nonlinear dynamics (4) by lifting the state to a higher-dimensional space of observables. Let $\mathcal{F}$ denote a space of scalar-valued functions on $\mathcal{X}$, and define the Koopman operator ${\mathcal{K}}_u:\mathcal{F}\to \mathcal{F}$ for a fixed control input u by:

$$\left({\mathcal{K}}_u\varphi \right)\left(x\right)=\varphi \left(f(x,u)\right)$$

(5)

for any observable $\varphi \in \mathcal{F}$.

The key insight is that while the original dynamics (4) are nonlinear in the state x, the evolution of observables under the Koopman operator is linear. This linearity enables the application of linear systems theory to analyze and control inherently nonlinear systems.

For practical implementation, we must work with finite-dimensional approximations of the infinite-dimensional Koopman operator. This is achieved through the Extended Dynamic Mode Decomposition (EDMD) method, which projects the Koopman operator onto a finite-dimensional subspace.

Definition 4

(Dictionary of Observables). A dictionary $\mathrm{\Phi }:\mathcal{X}\to {\mathbb{R}}^d$ is a vector-valued function whose components $\{{\varphi }_1,{\varphi }_2,\dots ,{\varphi }_d\}$ are scalar observables. We assume that:

• The dictionary is bounded: ${sup}_{x\in \mathcal{X}}\parallel \mathrm{\Phi }\left(x\right)\parallel <\infty$.
• The dictionary separates points: if $\mathrm{\Phi }\left(x_1\right)=\mathrm{\Phi }\left(x_2\right)$, then $x_1=x_2$.
• The dictionary contains a constant function: ${\varphi }_1\left(x\right)=1$ for all $x\in \mathcal{X}$.

Given a dictionary $\mathrm{\Phi }$, the EDMD method approximates the Koopman operator by the finite-dimensional matrices $(K,B)\in {\mathbb{R}}^{d\times d}\times {\mathbb{R}}^{d\times m}$ that best fit the linear relationship:

$$\mathrm{\Phi }\left(x_{t+1}\right)\approx K\mathrm{\Phi }\left(x_t\right)+B{u}_t.$$

(6)

The matrices $(K,B)$ are typically computed via least squares regression on collected trajectory data:

$$\begin{array}{cc}\hfill (K,B)& =\underset{(K^{\prime },B^{\prime })}{argmin}\sum _{i=1}^N{∥\mathrm{\Phi }\left(x_{i+1}\right)-K^{\prime }\mathrm{\Phi }\left(x_i\right)-B^{\prime }{u}_i∥}^2,\hfill \end{array}$$

(7)

where ${\left\{(x_i,u_i,x_{i+1})\right\}}_{i=1}^N$ represents collected trajectory data.

2.3. Connection Between ESL Theory and Koopman Operators

The fundamental connection between the ESL hierarchy and Koopman operator theory lies in recognizing that the EDMD approximation (6) defines a learning problem in the lifted coordinate space $z_t=\mathrm{\Phi }\left(x_t\right)$. The quality of this approximation and the existence of efficient algorithms depend critically on the spectral properties of the matrix K.

Definition 5

(Lifted System Dynamics). Given a dictionary Φ and EDMD matrices $(K,B)$, define the lifted system dynamics as:

$$z_{t+1}=K{z}_t+B{u}_t+w_t,$$

(8)

where $z_t=\mathrm{\Phi }\left(x_t\right)$ and $w_t$ represents the modeling error:

$$w_t=\mathrm{\Phi }\left(f(x_t,u_t)\right)-K\mathrm{\Phi }\left(x_t\right)-B{u}_t.$$

(9)

The key insight is that the autonomous part of the lifted dynamics (with $u_t=0$) corresponds to an optimization problem in the ESL hierarchy. Specifically, consider the problem of predicting the next lifted state given the current lifted state. This prediction problem is learnable if and only if the matrix K has spectral radius less than one, which corresponds to contractivity of the associated operator.

Lemma 1

(Spectral Radius and Contractivity). Let $K\in {\mathbb{R}}^{d\times d}$ be the autonomous EDMD matrix. The following are equivalent:

• The spectral radius satisfies $\rho \left(K\right)<1$.
• There exists a norm ${\parallel ·\parallel }_{*}$ on ${\mathbb{R}}^d$ such that K is a strict contraction: ${\parallel Kz\parallel }_{*}\le \gamma {\parallel z\parallel }_{*}$ for some $\gamma \in (0,1)$ and all $z\in {\mathbb{R}}^d$.
• The operator $T\left(z\right)=Kz$ is α-averaged for some $\alpha \in (0,1)$ in an appropriately chosen norm.

Proof. 

The equivalence of (1) and (2) follows from the spectral radius formula and the existence of matrix norms that make any matrix with spectral radius less than one into a contraction [16].

For (2) ⇒ (3), if K is a strict contraction with constant $\gamma <1$ in the norm ${\parallel ·\parallel }_{*}$, then the operator $T\left(z\right)=Kz$ can be written as $T=(1-\alpha )I+\alpha S$ where $\alpha =1-\gamma$ and $S\left(z\right)={\displaystyle \frac{Kz-(1-\alpha )z}{\alpha }}={\displaystyle \frac{Kz-\gamma z}{1-\gamma }}$. The operator S is nonexpansive in the norm ${\parallel ·\parallel }_{*}$ by construction.

For (3) ⇒ (1), if $T\left(z\right)=Kz$ is $\alpha$-averaged, then T is nonexpansive in some norm, which implies $\rho \left(K\right)\le 1$. The strict inequality $\rho \left(K\right)<1$ follows from the convergence requirement for $\alpha$-averaged operators.    □

This lemma establishes the fundamental connection between the spectral properties of EDMD matrices and the learnability classification in the ESL hierarchy. Systems with $\rho \left(K\right)<1$ correspond to learnable problems, while systems with $\rho \left(K\right)\ge 1$ are non-learnable in the lifted coordinates.

Definition 6

(Koopman-Learnable Systems). A nonlinear control system (4) is Koopman-learnable with respect to a dictionary Φ if the autonomous EDMD matrix K (computed with $u_t=0$) satisfies $\rho \left(K\right)<1$.

The significance of this definition extends beyond mere classification. Koopman-learnable systems admit convergent prediction algorithms in the lifted space, enable stable controller design, and, as we will show, possess finite-rank approximations with explicit error bounds.

2.4. Streamability and Low-Rank Structure

The final piece of our theoretical foundation concerns the connection between streamability in the ESL hierarchy and low-rank structure in Koopman operators. This connection is crucial for developing efficient algorithms that scale to high-dimensional systems.

Definition 7

(Koopman-Streamable Systems). A Koopman-learnable system is Koopman-streamable at rank r if the EDMD matrix K admits a rank-r approximation $K_r$ such that:

$$\parallel K-K_r{\parallel }_2\le {\epsilon }_r$$

(10)

for some tolerance ${\epsilon }_r>0$, where ${\parallel ·\parallel }_2$ denotes the spectral norm.

The rank-r approximation $K_r$ can be computed via the singular value decomposition (SVD) of K:

$$K=U\Sigma V^T=\sum _{i=1}^d{\sigma }_i{u}_i{v}_i^T,$$

(11)

where ${\sigma }_1\ge {\sigma }_2\ge \dots \ge {\sigma }_d\ge 0$ are the singular values. The optimal rank-r approximation is given by:

$$K_r=\sum _{i=1}^r{\sigma }_i{u}_i{v}_i^T,$$

(12)

with approximation error ${\epsilon }_r={\sigma }_{r+1}$ by the Eckart-Young theorem.

The streamability property has profound implications for computational efficiency. Instead of storing and manipulating the full $d\times d$ matrix K, we can work with the rank-r factorization $K_r=U_r{V}_r^T$ where $U_r,V_r\in {\mathbb{R}}^{d\times r}$. This reduces the storage requirement from $O\left(d^2\right)$ to $O\left(dr\right)$ and enables efficient streaming updates as new data becomes available.

Proposition 1

(Computational Benefits of Streamability). For a Koopman-streamable system at rank r with dictionary dimension d:

• Memory complexity:  Storage reduces from $O\left(d^2\right)$ to $O\left(dr\right)$.
• Computational complexity:  Matrix-vector products reduce from $O\left(d^2\right)$ to $O\left(dr\right)$ operations.
• Update complexity:  Streaming updates require $O\left(d{r}^2\right)$ operations per sample for TSQR orthonormalization instead of $O\left(d^3\right)$ for batch recomputation.

For typical applications with $d\sim {10}^3$ and $r\sim {10}^1$, this represents a one–to–two order-of-magnitude reductions in memory and time (depending on r and orthonormalization cadence).

These computational benefits become increasingly important as the dictionary dimension d grows, which is often necessary to capture complex nonlinear dynamics. The streamability framework provides a principled approach to achieving these benefits while maintaining theoretical guarantees on approximation quality and control performance.

3. Main Theoretical Results

This section presents our principal theoretical contributions, establishing the fundamental connections between ESL theory and Koopman operators, and providing constructive algorithms for finite-rank approximations with explicit error bounds. Our results resolve key theoretical gaps in the Koopman control literature while providing practical algorithms for implementation.

3.1. Learnability-Contractivity Equivalence

Our first main result establishes a precise equivalence between learnability in the ESL hierarchy and contractivity of finite Koopman operators. This equivalence provides the theoretical foundation for all subsequent developments.

Theorem 2

(Learnability–Contraction  (operative form)). For a linear system $z_{t+1}=K{z}_t$, the following are equivalent:

(i) 

$\rho \left(K\right)<1$ (Koopman-learnable).

(ii) 

There exists a weighted norm ${\parallel ·\parallel }_W$ such that K is a strict contraction: ${\parallel Kz\parallel }_W\le \gamma {\parallel z\parallel }_W$ for some $\gamma \in (0,1)$.

Remark on $\alpha$-averagedness.  If an operator $T_{\alpha }\left(z\right)=(1-\alpha )z+\alpha Kz$ is α-averaged, then K is nonexpansive (${\parallel K\parallel }_W\le 1$), which implies $\rho \left(K\right)\le 1$. However, to ensure $\rho \left(K\right)<1$, a strict contraction margin ${\parallel K\parallel }_W\le \gamma <1$ is required.

Proof. 

We establish the fundamental equivalence and then prove the sufficient conditions.

Fundamental Equivalence $\left(i\right)\iff \left(ii\right)$:

$\left(i\right)\Rightarrow \left(ii\right)$: If $\rho \left(K\right)<1$, then by standard spectral theory, there exists a matrix norm $\parallel ·\parallel$ such that $\parallel K\parallel <1$. Any matrix norm induces a corresponding weighted vector norm ${\parallel z\parallel }_W={\parallel W^{1/2}z\parallel }_2$ for some positive definite matrix W, establishing the contraction property.

$\left(ii\right)\Rightarrow \left(i\right)$: If K is a strict contraction in some weighted norm, then $\parallel K^n{z\parallel }_W\le {\gamma }^n{\parallel z\parallel }_W$ for all z and $n\ge 1$. This implies $\rho \left(K\right)\le \gamma <1$ by the spectral radius formula $\rho \left(K\right)={lim}_{n\to \infty }{\parallel K^n\parallel }^{1/n}$.

Sufficient Conditions:

$\left(i\right)\Rightarrow \left(iii\right)$ under normality: If K is normal (i.e., $K^{*}K=K{K}^{*}$) and $\rho \left(K\right)<1$, then ${\parallel K\parallel }_2=\rho \left(K\right)<1$, ensuring exponential convergence of the Picard iteration $z_{k+1}=K{z}_k$ with rate $\rho \left(K\right)$.

$\left(i\right)\Rightarrow \left(iv\right)$ under bounded condition number: If $K=VD{V}^{-1}$ with $\kappa \left(V\right)=\parallel V\parallel \parallel V^{-1}\parallel <\infty$, then $\parallel K^n{\parallel \le \kappa \left(V\right)\rho \left(K\right)}^n$, establishing convergence of fixed-point iterations when $\rho \left(K\right)<1$.

$\left(i\right)\Rightarrow \left(v\right)$ under appropriate iteration design: For the prediction problem ${min}_z{\parallel Kz-z_{t+1}\parallel }_W^2$, the gradient descent iteration $z^{(k+1)}=z^{\left(k\right)}-\eta K^{T}W(K{z}^{\left(k\right)}-z_{t+1})$ can be written as $z^{(k+1)}=(I-\eta K^{T}WK)z^{\left(k\right)}+\eta K^{T}W{z}_{t+1}$. When K is contractive and $\eta$ is chosen appropriately, this iteration admits an $\alpha$-averaged structure.    □

This theorem has several important implications. First, it provides a computationally verifiable criterion for determining when Koopman-based control is theoretically justified: simply compute the EDMD matrix K and check whether $\rho \left(K\right)<1$. Second, it guarantees that learnable systems admit convergent algorithms for model learning and state prediction. Third, it establishes the foundation for our subsequent results on finite-rank approximations and stability analysis.

3.2. Constructive Density Theorem for Koopman Systems

Our second main result extends the general ESL density theorem to the specific context of Koopman operators, providing a constructive algorithm for approximating any learnable system with a streamable (low-rank) surrogate.

Theorem 3

(Koopman Streamable Density). Let the nonlinear system (4) be Koopman-learnable with EDMD matrices $(K,B)$ and dictionary Φ. For any $\epsilon >0$, there exist rank-r matrices $(K_r,B_r)$ and a modified dictionary $\tilde{\mathrm{\Phi }}$ such that:

• The approximation error is bounded: $max\{\parallel K-K_r{\parallel }_2,\parallel B-B_r{\parallel }_2\}\le \epsilon$.
• The rank-r system $z_{t+1}=K_r{z}_t+B_r{u}_t$ is Koopman-streamable.
• The closed-loop performance degradation is bounded: if $u_t=\kappa \left(z_t\right)$ is a stabilizing controller for the original system, then the same controller applied to the rank-r system maintains stability with performance degradation $O\left(\epsilon \right)$.

Moreover, the construction is algorithmic and provides explicit bounds on the required rank $r\left(\epsilon \right)$.

Proof. 

The proof is constructive and proceeds through several steps.

Step 1: SVD-based rank reduction. Compute the singular value decomposition of the EDMD matrices:

$$\begin{array}{cc}\hfill K& =U_K{\Sigma }_K{V}_K^T=\sum _{i=1}^d{\sigma }_i^{\left(K\right)}{u}_i^{\left(K\right)}{\left(v_i^{\left(K\right)}\right)}^T,\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill B& =U_B{\Sigma }_B{V}_B^T=\sum _{i=1}^{min(d,m)}{\sigma }_i^{\left(B\right)}{u}_i^{\left(B\right)}{\left(v_i^{\left(B\right)}\right)}^T.\hfill \end{array}$$

(14)

Choose the rank r such that ${\sigma }_{r+1}^{\left(K\right)}\le \epsilon /2$ and ${\sigma }_{r+1}^{\left(B\right)}\le \epsilon /2$. Define the rank-r approximations:

$$\begin{array}{cc}\hfill K_r& =\sum _{i=1}^r{\sigma }_i^{\left(K\right)}{u}_i^{\left(K\right)}{\left(v_i^{\left(K\right)}\right)}^T,\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill B_r& =\sum _{i=1}^r{\sigma }_i^{\left(B\right)}{u}_i^{\left(B\right)}{\left(v_i^{\left(B\right)}\right)}^T.\hfill \end{array}$$

(16)

By the Eckart-Young theorem, this choice ensures $\parallel K-K_r{\parallel }_2={\sigma }_{r+1}^{\left(K\right)}\le \epsilon /2$ and $\parallel B-B_r{\parallel }_2={\sigma }_{r+1}^{\left(B\right)}\le \epsilon /2$.

Step 2: Streamability verification. The rank-r matrices $(K_r,B_r)$ admit factorizations $K_r=U_r{V}_r^T$ and $B_r=U_r{W}_r^T$ where $U_r\in {\mathbb{R}}^{d\times r}$, $V_r\in {\mathbb{R}}^{d\times r}$, and $W_r\in {\mathbb{R}}^{m\times r}$.

To verify streamability according to Definition 2, we need to show that the residual mapping $R:{\mathbb{R}}^d\times {\mathbb{R}}^m\to {\mathbb{R}}^d$ defined by

$$R(z,u)=z_{t+1}-K_{r}z-B_{r}u$$

(17)

can be uniformly approximated by separable functions. For the rank-r approximation, the prediction becomes:

$${\widehat{z}}_{t+1}=K_{r}z+B_{r}u=U_r(V_r^{T}z+W_r^{T}u).$$

(18)

The residual for any true trajectory $(z_t,u_t,z_{t+1})$ can be written as:

$$R(z_t,u_t)=z_{t+1}-U_r(V_r^T{z}_t+W_r^T{u}_t)=z_{t+1}-U_r{s}_t,$$

(19)

where $s_t=V_r^T{z}_t+W_r^T{u}_t\in {\mathbb{R}}^r$. This residual admits the separable approximation:

$$R(z,u)\approx \sum _{i=1}^r{\alpha }_i(z,u){\beta }_i$$

(20)

where ${\alpha }_i(z,u)={(V_r^{T}z+W_r^{T}u)}_i$ are scalar functions and ${\beta }_i=-u_i^{\left(r\right)}\in {\mathbb{R}}^d$ are the columns of $-U_r$. This satisfies the ESL separability condition with approximation error bounded by the truncation error $\epsilon$.

Step 3: Stability preservation. Suppose $u_t=\kappa \left(z_t\right)$ stabilizes the original system, meaning the closed-loop matrix $A_{cl}=K+B{\displaystyle \frac{\partial \kappa }{\partial z}}$ satisfies $\rho \left(A_{cl}\right)<1$. For the rank-r system, the closed-loop matrix becomes:

$$A_{cl}^{\left(r\right)}=K_r+B_r{\displaystyle \frac{\partial \kappa }{\partial z}}.$$

(21)

By matrix perturbation theory, assuming $A_{cl}$ is diagonalizable with eigenvector matrix V having condition number $\kappa \left(V\right)=\parallel V\parallel \parallel V^{-1}\parallel$, the Bauer-Fike theorem gives:

$$|\rho \left(A_{cl}^{\left(r\right)}\right)-\rho \left(A_{cl}\right)|\le \kappa \left(V\right)\parallel A_{cl}^{\left(r\right)}-A_{cl}{\parallel }_2\le \kappa \left(V\right)\epsilon (1+\parallel {\displaystyle \frac{\partial \kappa }{\partial z}}{\parallel }_2).$$

(22)

For systems where the closed-loop matrix is normal (e.g., symmetric or skew-symmetric), we have $\kappa \left(V\right)=1$ and the bound simplifies to the spectral norm bound. If $\epsilon$ is chosen sufficiently small such that $\kappa \left(V\right)\epsilon (1+\parallel {\displaystyle \frac{\partial \kappa }{\partial z}}{\parallel }_2)<1-\rho \left(A_{cl}\right)$, then $\rho \left(A_{cl}^{\left(r\right)}\right)<1$, ensuring stability of the rank-r system.

Step 4: Performance bound. The performance degradation can be quantified using standard linear systems theory. If $J={\sum }_{t=0}^{\infty }{z}_t^{T}Q{z}_t+u_t^{T}R{u}_t$ is the quadratic cost functional, then the cost difference between the original and rank-r systems satisfies:

$$|J_r-J|\le C\epsilon \parallel z_0{\parallel }^2$$

(23)

for some constant C depending on the system parameters and controller gains.

Step 5: Rank bound. The required rank $r\left(\epsilon \right)$ depends on the singular value decay of the EDMD matrices. For matrices with exponential singular value decay ${\sigma }_i\le C{e}^{-\alpha i}$, we have $r\left(\epsilon \right)=O(log(1/\epsilon \left)\right)$. For polynomial decay ${\sigma }_i\le C{i}^{-\beta }$, we have $r\left(\epsilon \right)=O\left({\epsilon }^{-1/\beta }\right)$.    □

This theorem provides several crucial guarantees. First, it ensures that any Koopman-learnable system can be approximated arbitrarily well by a streamable system, resolving the fundamental question of when finite-rank approximations are justified. Second, it provides explicit error bounds that enable principled trade-offs between computational efficiency and approximation quality. Third, it guarantees that stability and performance properties are preserved under the approximation, subject to appropriate choice of the rank parameter.

3.3. Error Propagation and Robustness Analysis

Our third main result analyzes how approximation errors in the Koopman model propagate through the closed-loop system when controllers designed on the approximate model are applied to the true plant.

Theorem 4

(Closed-Loop Error Propagation). Consider a Koopman-learnable system with true EDMD matrices $(K,B)$ and rank-r approximations $(K_r,B_r)$. Let $u_t=\kappa \left(z_t\right)$ be a controller designed for the rank-r system such that the closed-loop matrix $A_{cl}^{\left(r\right)}=K_r+B_r{\displaystyle \frac{\partial \kappa }{\partial z}}$ satisfies $\parallel A_{cl}^{\left(r\right)}{\parallel }_2\le \gamma <1$.

When this controller is applied to the true system, the lifted state evolution satisfies:

$$z_{t+1}=A_{cl}^{\left(r\right)}{z}_t+\Delta A{z}_t+w_t,$$

(24)

where $\Delta A=(K-K_r)+(B-B_r){\displaystyle \frac{\partial \kappa }{\partial z}}$ is the model mismatch and $w_t$ is the Koopman modeling error.

Under the bounds ${\parallel \mathrm{\Delta }A\parallel }_2\le \delta$ and $\parallel w_t{\parallel }_2\le {\epsilon }_w$ for all t, the closed-loop system is input-to-state stable with:

$$\parallel z_t{\parallel }_2\le {\gamma }^t\parallel z_0{\parallel }_2+{\displaystyle \frac{\delta }{1-\gamma }}\underset{0\le s\le t-1}{sup}{\parallel z_s\parallel }_2+{\displaystyle \frac{{\epsilon }_w}{1-\gamma }}.$$

(25)

If $\delta <1-\gamma$, then the system remains stable with ultimate bound:

$$\underset{t\to \infty }{lim\; sup}{\parallel z_t\parallel }_2\le {\displaystyle \frac{{\epsilon }_w}{1-\gamma -\delta }}.$$

(26)

Proof. 

The proof follows standard input-to-state stability analysis for linear systems with additive disturbances.

Starting from the perturbed dynamics:

$$z_{t+1}=A_{cl}^{\left(r\right)}{z}_t+\mathrm{\Delta }A{z}_t+w_t,$$

(27)

we can bound the state evolution recursively. Taking norms:

$$\begin{array}{cc}\hfill \parallel z_{t+1}{\parallel }_2& \le \parallel A_{cl}^{\left(r\right)}{\parallel }_2\parallel z_t{\parallel }_2+{\parallel \mathrm{\Delta }A\parallel }_2\parallel z_t{\parallel }_2+{\parallel w_t\parallel }_2\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill & \le \gamma \parallel z_t{\parallel }_2+\delta {\parallel z_t\parallel }_2+{\epsilon }_w\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill & =(\gamma +\delta )\parallel z_t{\parallel }_2+{\epsilon }_w.\hfill \end{array}$$

(30)

If $\gamma +\delta <1$, then the recursion yields:

$$\parallel z_t{\parallel }_2\le {(\gamma +\delta )}^t{\parallel z_0\parallel }_2+{\displaystyle \frac{{\epsilon }_w}{1-\gamma -\delta }}.$$

(31)

For the more general case where $\delta$ may be larger, we use the Bellman-Grönwall lemma. The recursion can be written as:

$$\parallel z_t{\parallel }_2\le {\gamma }^t\parallel z_0{\parallel }_2+\sum _{s=0}^{t-1}{\gamma }^{t-1-s}(\delta \parallel z_s{\parallel }_2+{\epsilon }_w).$$

(32)

Applying the discrete Grönwall inequality with $\delta <1-\gamma$ yields the bound (25).    □

This theorem provides crucial insights for practical implementation. It shows that the stability margin $1-\gamma$ of the nominal controller must exceed the model mismatch $\delta$ to maintain closed-loop stability. This provides a principled approach for choosing the rank r in the approximation: select r large enough that ${\parallel \mathrm{\Delta }A\parallel }_2<1-\gamma$, where $\gamma$ is the desired closed-loop contraction rate.

3.4. Constructive Algorithm for Streamable Approximation

Building on the theoretical results, we now present a constructive algorithm that converts any Koopman-learnable system into a streamable approximation with explicit error bounds.

Note on dimensions: The factorization uses $K_t\approx U_t{V}_t^{\top }$ and $B_t\approx U_t{W}_t^{\top }$ with $U_t\in {\mathbb{R}}^{d\times r}$, $V_t\in {\mathbb{R}}^{d\times r}$, $W_t\in {\mathbb{R}}^{m\times r}$. The regularization parameter $R_t⪰0$ (e.g., $R_t=\lambda I$) provides forgetting/regularization for adaptation to time-varying dynamics.

Proposition 2

(Algorithm Correctness). Algorithm 1 terminates in finite time and produces a rank-r streamable approximation $(K_r,B_r)$ satisfying:

• Error bound: $max\{\parallel K-K_r{\parallel }_2,\parallel B-B_r{\parallel }_2\}\le \epsilon$.
• Stability preservation: $\rho \left(K_r\right)<1$.
• Streamability: The system admits $O\left(rd\right)$ streaming updates.

The required rank satisfies $r\le min\{d,\lceil log({\sigma }_1^{\left(K\right)}/\epsilon )/log({\sigma }_1^{\left(K\right)}/{\sigma }_2^{\left(K\right)})\rceil \}$ for exponentially decaying singular values.

Algorithm 1 Single-Loop Adaptive Streaming Koopman Control (ASKC)

1: Inputs: dictionary $\mathrm{\Phi }$, rank r, $P\succ 0$, $Q\succ 0$, contraction target $\gamma \in (0,1)$, stepsizes ${\eta }_t$, dither schedule ${\sigma }_t$ 2: Initialize $U_0,V_0\in {\mathbb{R}}^{d\times r}$ (orthonormal), $W_0\in {\mathbb{R}}^{m\times r}$, compute $K_0^u$ (stabilizing), set ${\epsilon }_0\leftarrow 0$ 3: for  $t=0,1,2,\dots$  do 4:    Measure $x_t$, set $z_t=\Phi \left(x_t\right)$ 5:    Tentative control $u_t^{★}=K_t^u{z}_t+{\xi }_t$, with $\parallel {\xi }_t\parallel \le {\sigma }_t$ 6:    Safety QP: choose $u_t\approx u_t^{★}$ s.t. predicted $z_{t+1}=U_t{V}_t^{\top }{z}_t+U_t{W}_t^{\top }{u}_t$ remains in tube of radius $r_t$ and $\parallel P^{1/2}{A}_{\mathrm{cl}}\left(t\right)P^{-1/2}{\parallel }_2\le \gamma$, where $A_{\mathrm{cl}}\left(t\right)=U_t{V}_t^{\top }+U_t{W}_t^{\top }{K}_t^u$ 7:    Apply $u_t$, observe $x_{t+1}$, set $z_{t+1}=\Phi \left(x_{t+1}\right)$ 8:    Streaming EDMD factors (Oja/PAST): $e_t\leftarrow z_{t+1}-U_t(V_t^{\top }{z}_t+W_t^{\top }{u}_t)$ 9:    $U\leftarrow U_t+{\eta }_t{e}_t{(V_t^{\top }{z}_t+W_t^{\top }{u}_t)}^{\top }$;    $V\leftarrow V_t+{\eta }_t(z_t{\left(U_t^{\top }{e}_t\right)}^{\top }-V_t{R}_t)$;    $W\leftarrow W_t+{\eta }_t(u_t{\left(U_t^{\top }{e}_t\right)}^{\top }-W_t{R}_t)$;    TSQR on $U,V$ 10:    Update residual bound ${\epsilon }_t\leftarrow (1-\lambda ){\parallel e_t\parallel }^2+\lambda {\epsilon }_{t-1}$ and tube radius $r_{t+1}\leftarrow {\displaystyle \frac{\sqrt{{\epsilon }_t}}{1-\gamma }}$ 11:    Controller update: one Kleinman/Riccati iteration warm-started from previous solution to get $K_{t+{\displaystyle \frac{1}{2}}}^u$ 12:    Projection: if $\parallel P^{1/2}{A}_{\mathrm{cl}}^{\mathrm{tmp}}{P}^{-1/2}{\parallel }_2>\gamma$ for $A_{\mathrm{cl}}^{\mathrm{tmp}}=U{V}^{\top }+U{W}^{\top }{K}_{t+{\displaystyle \frac{1}{2}}}^u$, shrink step or solve a proximal QP to project $(U,V,W,K^u)$ back to feasibility 13:    Commit $(U_{t+1},V_{t+1},W_{t+1},K_{t+1}^u)$ 14: end for

Proof. Termination:

The algorithm terminates because the rank r is bounded above by d, and for any $r=d$, we have $K_r=K$ and $B_r=B$, which preserves the original spectral radius.

Error bound: The error bound follows directly from the Eckart-Young theorem and the choice of rank in Step 3.

Stability preservation: The spectral radius is continuous in the matrix entries, so for sufficiently small perturbations, $\rho \left(K_r\right)$ remains close to $\rho \left(K\right)<1$. The algorithm explicitly checks this condition and increases the rank if necessary.

Streamability: The factorized form $K_r=U_r{V}_r^T$ enables streaming updates with $O\left(rd\right)$ complexity per iteration.

Rank bound: For exponentially decaying singular values ${\sigma }_i={\sigma }_1{\rho }^{i-1}$ with $\rho <1$, the condition ${\sigma }_{r+1}\le \epsilon$ gives ${\sigma }_1{\rho }^r\le \epsilon$, which yields $r\le log({\sigma }_1/\epsilon )/log(1/\rho )$. □

This algorithm provides a complete computational framework for converting Koopman-learnable systems into streamable approximations. The algorithm is practical and can be implemented using standard numerical linear algebra routines. The error bounds are explicit and enable principled trade-offs between computational efficiency and approximation quality.

4. Streaming Algorithm for Real-Time Koopman Learning

This section develops a streaming algorithm for real-time learning of low-rank Koopman models that operates with $\mathcal{O}\left(Kd\right)$ memory complexity and provides convergence guarantees. The algorithm enables online adaptation of the Koopman model as new trajectory data becomes available, making it suitable for control applications where the system dynamics may change over time or where initial training data is limited.

4.1. Streaming EDMD with Low-Rank Factorization

Traditional EDMD requires storing and manipulating full $d\times d$ matrices, which becomes computationally prohibitive for large dictionary dimensions. Our streaming approach maintains low-rank factorizations and updates them incrementally as new data arrives.

The complete single-loop streaming algorithm is presented in Algorithm 1.

The algorithm maintains three key components: the left factor $U_t\in {\mathbb{R}}^{d\times r}$ and right factor $V_t\in {\mathbb{R}}^{d\times r}$ for the autonomous dynamics matrix $K_t\approx U_t{V}_t^T$, and the control factor $W_t\in {\mathbb{R}}^{r\times m}$ such that $B_t\approx U_t{W}_t^T$. This factorization reduces the storage requirement from $\mathcal{O}\left(d^2\right)$ to $\mathcal{O}\left(rd\right)$ and enables efficient streaming updates.

4.2. Convergence Analysis

We now analyze the convergence properties of Algorithm 1 under standard assumptions on the data distribution and system dynamics.

Assumption 2

(Data Distribution). The trajectory data $\left\{(x_t,u_t,x_{t+1})\right\}$ is generated by the true system (4) with the following properties:

• The lifted states $z_t=\mathrm{\Phi }\left(x_t\right)$ are bounded: ${sup}_t{\parallel z_t\parallel }_2\le M_z<\infty$.
• The control inputs are bounded: ${sup}_t{\parallel u_t\parallel }_2\le M_u<\infty$.
• The data satisfies a persistence of excitation condition: there exist constants $\alpha ,\beta >0$ such that for all T sufficiently large:

  $$\alpha I⪯{\displaystyle \frac{1}{T}}\sum _{t=1}^T{z}_t{z}_t^T⪯\beta I.$$

  (33)

Assumption 3

(Learning Rate Schedule). The learning rates $\left\{{\eta }_t\right\}$ satisfy the Robbins-Monro conditions:

$$\sum _{t=1}^{\infty }{\eta }_t=\infty ,\sum _{t=1}^{\infty }{\eta }_t^2<\infty .$$

(34)

A typical choice is ${\eta }_t={\eta }_0/(1+\gamma t)$ for constants ${\eta }_0>0$ and $\gamma >0$.

Assumption 4

(Dithering Schedule). The dithering sequence $\left\{{\xi }_t\right\}$ in the tentative control $u_t^{★}=K_t^u{z}_t+{\xi }_t$ satisfies:

• Bounded support: $\parallel {\xi }_t{\parallel }_2\le {\sigma }_t$ with ${\sigma }_t\to 0$ as $t\to \infty$.
• Sufficient excitation: The dithering ensures that the augmented data matrix $\left[\begin{array}{c}{z}_t\\ u_t\end{array}\right]$ maintains persistent excitation for identifiability of the low-rank factors $(U_t,V_t,W_t)$.
• Decay rate: ${\sigma }_t={\sigma }_0/{(1+\alpha t)}^{\beta }$ with $\beta \in (0,1/2)$ to balance exploration and exploitation.

Theorem 5

(Convergence of Streaming EDMD). Under Assumptions 1, 2, 3, and 4, the streaming estimates produced by Algorithm 1 converge to the true EDMD matrices in the following sense:

• Almost sure convergence:  With probability 1,

  $$\underset{t\to \infty }{lim}\parallel U_t{V}_t^T-K^{*}{\parallel }_F=0,\underset{t\to \infty }{lim}{\parallel U_t{W}_t^T-B^{*}\parallel }_F=0,$$

  (35)
  where $(K^{*},B^{*})$ are the true EDMD matrices.
• Convergence rate:  The mean-squared error satisfies

  $$\mathbb{E}[\parallel U_t{V}_t^T-K^{*}{\parallel }_F^2]=\mathcal{O}(1/t),\mathbb{E}[\parallel U_t{W}_t^T-B^{*}{\parallel }_F^2]=\mathcal{O}(1/t).$$

  (36)
• Finite-sample bound:  For any $\delta >0$, with probability at least $1-\delta$,

  $$\parallel U_t{V}_t^T-K^{*}{\parallel }_F\le C\sqrt{{\displaystyle \frac{log(1/\delta )}{t}}}$$

  (37)
  for some constant C depending on the system parameters.

Proof. 

The proof follows the framework of stochastic approximation theory for matrix factorization problems [18].

Step 1: Reformulation as stochastic approximation. The streaming updates can be written in the form:

$${\theta }_{t+1}={\theta }_t+{\eta }_t(H({\theta }_t,{\xi }_t)-{\theta }_t),$$

(38)

where ${\theta }_t=(U_t,V_t,W_t)$ represents the parameter vector, ${\xi }_t=(z_t,u_t,z_{t+1})$ is the random data sample, and $H(\theta ,\xi )$ represents the update function.

Step 2: Lipschitz continuity. Under Assumption 2, the update function H is Lipschitz continuous in $\theta$ with constant L depending on $M_z$ and $M_u$:

$$\parallel H({\theta }_1,\xi )-H({\theta }_2,\xi ){\parallel }_F\le L{\parallel {\theta }_1-{\theta }_2\parallel }_F.$$

(39)

Step 3: Unbiased gradient condition. The expected update direction satisfies:

$${\mathbb{E}}_{\xi }\left[H(\theta ,\xi )\right]={\nabla }_{\theta }\mathbb{E}[\ell (\theta ,\xi )],$$

(40)

where $\ell (\theta ,\xi )=\parallel z_{t+1}-U{V}^T{z}_t-U{W}^T{u}_t{\parallel }_F^2$ is the instantaneous loss function.

Step 4: Strong convexity. Under the persistence of excitation condition in Assumption 2, the expected loss function is strongly convex in a neighborhood of the optimum, ensuring unique convergence.

Step 5: Application of stochastic approximation theory. The convergence results follow from standard theorems in stochastic approximation [19], combined with the specific structure of the matrix factorization problem.

The $\mathcal{O}(1/t)$ convergence rate is optimal for this class of problems and matches the rate achieved by batch methods when averaged over multiple runs. □

4.3. Computational Complexity Analysis

We now analyze the computational complexity of the streaming algorithm and compare it to traditional batch EDMD methods.

Proposition 3

(Computational Complexity). Algorithm 1 has the following computational complexity per iteration:

• Memory complexity:  $\mathcal{O}(rd+rm)$ for storing the factors $(U_t,V_t,W_t)$.
• Time complexity per sample:  $\mathcal{O}(rd+rm)$ for the main updates, plus $\mathcal{O}\left(r{d}^2\right)$ for orthonormalization.
• Total time complexity for T samples:  $\mathcal{O}\left(T(rd+rm+r{d}^2)\right)$.

In contrast, batch EDMD requires:

• Memory complexity:  $\mathcal{O}(d^2+dm+Td)$ for storing matrices and data.
• Time complexity:  $\mathcal{O}(T{d}^2+d^3)$ for matrix operations and SVD.

For typical parameter ranges ($d\sim {10}^3$, $r\sim {10}^1$, $T\sim {10}^4$), the streaming algorithm achieves 10-100× reductions in both memory and computational requirements.

Proof. Memory analysis:

The algorithm stores three matrices: $U_t\in {\mathbb{R}}^{d\times r}$, $V_t\in {\mathbb{R}}^{d\times r}$, and $W_t\in {\mathbb{R}}^{r\times m}$, giving total storage $\mathcal{O}(rd+rm)$.

Time analysis per iteration:

• Computing ${\widehat{z}}_{t+1}$: $\mathcal{O}\left(rd\right)$ for $V_t^T{z}_t$ and $U_t(·)$, plus $\mathcal{O}\left(rm\right)$ for $W_t^T{u}_t$.
• Computing error $e_t$: $\mathcal{O}\left(d\right)$.
• Updating $U_t$ and $V_t$: $\mathcal{O}\left(rd\right)$ each.
• Orthonormalization: $\mathcal{O}\left(d{r}^2\right)$ using QR decomposition.
• Updating $W_t$: $\mathcal{O}\left(rm\right)$.

The dominant term is the orthonormalization step with $\mathcal{O}\left(d{r}^2\right)$ complexity.

Comparison with batch EDMD: Batch EDMD requires forming the data matrices $X,X^{\prime },U$ (storage $\mathcal{O}\left(Td\right)$), solving the least squares problem (time $\mathcal{O}\left(T{d}^2\right)$), and computing the SVD for rank reduction (time $\mathcal{O}\left(d^3\right)$).

For $d=1000$, $r=10$, $T=10000$:

• Streaming: Memory $\mathcal{O}\left({10}^4\right)$, Time $\mathcal{O}\left({10}^9\right)$.
• Batch: Memory $\mathcal{O}\left({10}^7\right)$, Time $\mathcal{O}(1.1\times {10}^{10})$.

This gives approximately 1000× memory reduction and ≈10× time reduction. □

4.4. Adaptive Rank Selection

A crucial aspect of the streaming algorithm is the adaptive adjustment of the rank r based on the prediction error. This enables the algorithm to automatically balance approximation quality with computational efficiency.

The adaptive rank selection mechanism is integrated into Algorithm 1 and enables automatic discovery of appropriate model complexity.

Proposition 4

(Adaptive Rank Properties). The adaptive rank selection mechanism maintains the following properties:

• Orthogonality preservation:  The factors $U_t$ and $V_t$ remain orthonormal after rank adjustments.
• Approximation quality:  Rank expansion reduces the prediction error by at least $\parallel e_t{\parallel }_2^2/d$.
• Computational efficiency:  Rank contraction removes the least significant factor, minimizing the impact on approximation quality.
• Stability:  The rank adjustment process is stable and does not cause oscillations when the thresholds are chosen appropriately.

Proof. Orthogonality preservation: The rank expansion explicitly orthogonalizes the new factor against existing ones. Rank contraction maintains orthogonality by removing entire columns.

Approximation quality: When expanding rank, the new factor $u_{\mathrm{new}}$ is chosen in the direction of the prediction error $e_t$. This reduces the residual by at least the component of $e_t$ in the direction orthogonal to the current subspace.

Computational efficiency: The importance measure $s_i$ captures the contribution of each factor to both the autonomous dynamics ($V_t$) and control dynamics ($W_t$). Removing the factor with smallest $s_i$ minimizes the approximation degradation.

Stability: The threshold-based approach with hysteresis (different thresholds for expansion and contraction) prevents rapid oscillations in the rank. □

The adaptive rank selection mechanism enables the algorithm to automatically discover the appropriate model complexity for the given system and data quality. This is particularly valuable in control applications where the system complexity may not be known a priori or may change over time due to wear, environmental conditions, or operating regime changes.

5. Stability Analysis and Controller Synthesis

Lemma 2

(Bauer–Fike stability transfer [16]). Assume $A_{\mathrm{cl}}$ is diagonalizable with $A_{\mathrm{cl}}=V\Lambda V^{-1}$. For any perturbation Δ, the eigenvalues of ${\tilde{A}}_{\mathrm{cl}}=A_{\mathrm{cl}}+\mathrm{\Delta }$ satisfy $|\lambda -\tilde{\lambda }{|\le \kappa \left(V\right)\parallel \mathrm{\Delta }\parallel }_2$ with $\kappa \left(V\right)={\parallel V\parallel }_2{\parallel V^{-1}\parallel }_2$. Hence if $\rho \left(A_{\mathrm{cl}}\right)<1-\delta$ and ${\kappa \left(V\right)\parallel \mathrm{\Delta }\parallel }_2<\delta$, then $\rho \left({\tilde{A}}_{\mathrm{cl}}\right)<1$.

This section addresses the critical challenge of designing controllers on finite-rank Koopman surrogates while ensuring stability and performance when deployed on the true nonlinear plant. We develop comprehensive stability analysis including Input-to-State Stability (ISS) bounds, robust Model Predictive Control (MPC) formulations, and practical design guidelines for maintaining stability margins.

5.1. Controller Design on Finite-Rank Surrogates

Consider the rank-r Koopman surrogate obtained from Algorithm 1:

$$z_{t+1}=K_r{z}_t+B_r{u}_t,$$

(41)

where $z_t=\mathrm{\Phi }\left(x_t\right)\in {\mathbb{R}}^d$ represents the lifted state and $K_r,B_r$ are the rank-r approximations of the true EDMD matrices.

For this linear surrogate, we can apply standard linear control design techniques. We consider both Linear Quadratic Regulator (LQR) and Model Predictive Control (MPC) approaches.

5.1.1. LQR Design on the Surrogate

For the quadratic cost functional:

$$J_{\mathrm{LQR}}=\sum _{t=0}^{\infty }(z_t^{T}Q{z}_t+u_t^{T}R{u}_t),$$

(42)

where $Q⪰0$ and $R\succ 0$ are weighting matrices, the optimal LQR controller is given by:

$$u_t=-K_{\mathrm{LQR}}{z}_t,K_{\mathrm{LQR}}={(R+B_r^{T}P{B}_r)}^{-1}{B}_r^{T}P{K}_r,$$

(43)

where P is the solution to the discrete-time algebraic Riccati equation:

$$P=Q+K_r^{T}P{K}_r-K_r^{T}P{B}_r{(R+B_r^{T}P{B}_r)}^{-1}{B}_r^{T}P{K}_r.$$

(44)

Assumption 5

(Controllability and Stabilizability). The pair $(K_r,B_r)$ is stabilizable, meaning there exists a feedback gain $K_{\mathit{LQR}}$ such that the closed-loop matrix $A_{\mathit{cl}}=K_r-B_r{K}_{\mathit{LQR}}$ has all eigenvalues inside the unit circle.

Proposition 5

(LQR Stability on Surrogate). Under Assumption 5, the LQR controller stabilizes the surrogate system (41) with exponential convergence rate determined by the spectral radius of the closed-loop matrix:

$$\parallel z_t{\parallel }_P\le \rho {\left(A_{\mathrm{cl}}\right)}^t{\parallel z_0\parallel }_P,$$

(45)

where ${\parallel ·\parallel }_P$ denotes the weighted norm ${\parallel z\parallel }_P=\sqrt{z^{T}Pz}$.

5.1.2. MPC Design on the Surrogate

For more complex constraints and performance objectives, we employ Model Predictive Control. The MPC optimization problem at time t is:

$$\begin{array}{cc}\hfill \underset{u_t,\dots ,u_{t+N-1}}{min}& \sum _{k=0}^{N-1}(\parallel z_{t+k|t}{\parallel }_Q^2+\parallel u_{t+k}{\parallel }_R^2)+\parallel z_{t+N|t}{\parallel }_P^2\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill \mathrm{subject}\mathrm{to}& z_{t+k+1|t}=K_r{z}_{t+k|t}+B_r{u}_{t+k},k=0,\dots ,N-1\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill & z_{t|t}=z_t\hfill \end{array}$$

(48)

$$\begin{array}{cc}\hfill & u_{t+k}\in \mathcal{U},k=0,\dots ,N-1\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill & z_{t+k|t}\in \mathcal{Z},k=1,\dots ,N\hfill \end{array}$$

(50)

where N is the prediction horizon, $\mathcal{U}$ and $\mathcal{Z}$ are constraint sets, and $z_{t+k|t}$ denotes the predicted state at time $t+k$ given information at time t.

Assumption 6

(MPC Feasibility). The constraint sets $\mathcal{U}$ and $\mathcal{Z}$ are compact and contain the origin in their interior. The terminal cost matrix P and terminal constraint set ${\mathcal{Z}}_f\subset \mathcal{Z}$ are chosen such that the MPC problem is recursively feasible.

Proposition 6

(MPC Stability on Surrogate). Under Assumption 6, the MPC controller stabilizes the surrogate system (41) with the region of attraction determined by the initial feasible set of the MPC optimization problem.

5.2. Robustness Analysis for True Plant Deployment

The critical challenge arises when controllers designed on the finite-rank surrogate are deployed on the true nonlinear plant. We must account for two sources of uncertainty: the Koopman modeling error and the finite-rank approximation error.

5.2.1. Error Decomposition

When the controller $u_t=\kappa \left(z_t\right)$ designed for the surrogate is applied to the true plant, the actual lifted dynamics become:

$$z_{t+1}=\mathrm{\Phi }\left(f(x\left(z_t\right),u_t)\right)=K{z}_t+B{u}_t+w_t,$$

(51)

where $w_t$ represents the Koopman modeling error, and $(K,B)$ are the true EDMD matrices.

The total error can be decomposed as:

$$\begin{array}{cc}\hfill z_{t+1}-K_r{z}_t-B_r{u}_t& =(K-K_r)z_t+(B-B_r)u_t+w_t\hfill \end{array}$$

(52)

$$\begin{array}{cc}\hfill & =\mathrm{\Delta }K{z}_t+\mathrm{\Delta }B{u}_t+w_t,\hfill \end{array}$$

(53)

where $\mathrm{\Delta }K=K-K_r$ and $\mathrm{\Delta }B=B-B_r$ represent the finite-rank approximation errors.

Assumption 7

(Bounded Uncertainties). The modeling and approximation errors satisfy:

• Koopman modeling error: $\parallel w_t{\parallel }_2\le {\epsilon }_w$ for all t.
• Approximation error bounds: ${\parallel \mathrm{\Delta }K\parallel }_2\le {\delta }_K$ and ${\parallel \mathrm{\Delta }B\parallel }_2\le {\delta }_B$.
• State and input bounds: $\parallel z_t{\parallel }_2\le M_z$ and $\parallel u_t{\parallel }_2\le M_u$ for all t.

Key Parameter Definitions:

• Total model mismatch:  $\delta :={\parallel \mathrm{\Delta }K\parallel }_2+{\parallel \mathrm{\Delta }B\parallel }_2{sup}_z{∥{\displaystyle \frac{\partial \kappa }{\partial z}}\left(z\right)∥}_2$
• Disturbance bound: ${\epsilon }_w$ is the upper bound on modeling errors: $\parallel w_t{\parallel }_2\le {\epsilon }_w$
• Residual tracking: ${\epsilon }_t=(1-\lambda ){\parallel e_t\parallel }^2+\lambda {\epsilon }_{t-1}$ with forgetting factor $\lambda \in (0,1)$
• Tube radius update: $r_{t+1}={\displaystyle \frac{\sqrt{{\epsilon }_t}}{1-\gamma }}$ where $\gamma$ is the desired contraction rate

5.2.2. Input-to-State Stability Analysis

We now establish ISS properties for the closed-loop system when the surrogate-based controller is applied to the true plant.

Theorem 6

(ISS for Surrogate-Based Control). Consider the true plant dynamics (51) with controller $u_t=\kappa \left(z_t\right)$ designed for the surrogate (41). Suppose the surrogate closed-loop matrix $A_{\mathrm{cl}}^{\left(r\right)}=K_r+B_r{\displaystyle \frac{\partial \kappa }{\partial z}}$ satisfies $\parallel A_{\mathrm{cl}}^{\left(r\right)}{\parallel }_2\le \gamma <1$.

Under Assumption 7, if the total perturbation bound:

$${\delta }_{\mathrm{total}}:={\delta }_K+{\delta }_B{∥{\displaystyle \frac{\partial \kappa }{\partial z}}∥}_2<1-\gamma ,$$

(54)

then the closed-loop system is ISS with:

$$\parallel z_t{\parallel }_2\le {\displaystyle \frac{{\gamma }^t}{1-{\delta }_{\mathrm{total}}}}{\parallel z_0\parallel }_2+{\displaystyle \frac{{\epsilon }_w+{\delta }_{\mathrm{total}}{M}_z}{1-\gamma -{\delta }_{\mathrm{total}}}}.$$

(55)

Moreover, the ultimate bound satisfies:

$$\underset{t\to \infty }{lim\; sup}{\parallel z_t\parallel }_2\le {\displaystyle \frac{{\epsilon }_w+{\delta }_{\mathrm{total}}{M}_z}{1-\gamma -{\delta }_{\mathrm{total}}}}.$$

(56)

Proof. 

The true closed-loop dynamics can be written as:

$$\begin{array}{cc}\hfill z_{t+1}& =K{z}_t+B{u}_t+w_t\hfill \end{array}$$

(57)

$$\begin{array}{cc}\hfill & =K{z}_t+B\kappa \left(z_t\right)+w_t\hfill \end{array}$$

(58)

$$\begin{array}{cc}\hfill & =(K_r+\mathrm{\Delta }K)z_t+(B_r+\mathrm{\Delta }B)\kappa \left(z_t\right)+w_t\hfill \end{array}$$

(59)

$$\begin{array}{cc}\hfill & =K_r{z}_t+B_r\kappa \left(z_t\right)+\mathrm{\Delta }K{z}_t+\mathrm{\Delta }B\kappa \left(z_t\right)+w_t\hfill \end{array}$$

(60)

$$\begin{array}{cc}\hfill & =A_{\mathrm{cl}}^{\left(r\right)}{z}_t+\mathrm{\Delta }A{z}_t+w_t,\hfill \end{array}$$

(61)

where the total model mismatch is bounded by:

$$\delta :={\parallel \mathrm{\Delta }K\parallel }_2+{\parallel \mathrm{\Delta }B\parallel }_2\underset{z}{sup}{∥{\displaystyle \frac{\partial \kappa }{\partial z}}\left(z\right)∥}_2.$$

(62)

Taking norms and using the triangle inequality:

$$\begin{array}{cc}\hfill \parallel z_{t+1}{\parallel }_2& \le \parallel A_{\mathrm{cl}}^{\left(r\right)}{\parallel }_2\parallel z_t{\parallel }_2+\parallel \mathrm{\Delta }K{z}_t+\mathrm{\Delta }B\kappa \left(z_t\right){\parallel }_2+{\parallel w_t\parallel }_2\hfill \end{array}$$

(63)

$$\begin{array}{cc}\hfill & \le \gamma \parallel z_t{\parallel }_2+\delta {\parallel z_t\parallel }_2+{\epsilon }_w\hfill \end{array}$$

(64)

$$\begin{array}{cc}\hfill & =(\gamma +\delta )\parallel z_t{\parallel }_2+{\epsilon }_w.\hfill \end{array}$$

(65)

If $\gamma +\delta <1$, this recursion yields:

$$\parallel z_t{\parallel }_2\le {(\gamma +\delta )}^t{\parallel z_0\parallel }_2+{\displaystyle \frac{{\epsilon }_w}{1-\gamma -\delta }}.$$

(66)

For the more general bound, we use the fact that $\parallel z_t{\parallel }_2\le M_z$ to bound the perturbation term, leading to the stated result through application of the discrete Grönwall inequality. □

This theorem provides crucial design guidance: the stability margin $1-\gamma$ of the nominal controller must exceed the total model uncertainty $\delta$ to maintain closed-loop stability. This enables principled selection of the rank r in the approximation.

5.3. Robust Model Predictive Control

To explicitly account for model uncertainties in the controller design, we develop a robust MPC formulation that maintains stability guarantees despite approximation errors.

5.3.1. Tube-Based Robust MPC

The key idea is to design the controller for the nominal surrogate model while maintaining a "tube" around the nominal trajectory that bounds the effects of model uncertainties.

Definition 8

(Robust Positively Invariant Set). A set $\mathcal{S}\subset {\mathbb{R}}^d$ is robust positively invariant (RPI) for the uncertain system:

$$z_{t+1}=A_{\mathrm{cl}}^{\left(r\right)}{z}_t+w_t,{\parallel w_t\parallel }_2\le {\epsilon }_{\mathrm{total}},$$

(67)

if for all $z\in \mathcal{S}$ and all admissible disturbances w, we have $A_{\mathrm{cl}}^{\left(r\right)}z+w\in \mathcal{S}$.

Proposition 7

(Ellipsoidal RPI Set). If $\parallel A_{\mathrm{cl}}^{\left(r\right)}{\parallel }_2\le \gamma <1$, then the ellipsoid:

$$\mathcal{S}=\left\{z\in {\mathbb{R}}^d:z^{T}Pz\le {\displaystyle \frac{{\epsilon }_{\mathrm{total}}^2}{{(1-\gamma )}^2}}\right\}$$

(68)

is RPI for any $P\succ 0$ such that $A_{\mathrm{cl}}^{\left(r\right)T}P{A}_{\mathrm{cl}}^{\left(r\right)}⪯{\gamma }^{2}P$.

Proof. 

For any $z\in \mathcal{S}$ and disturbance ${\parallel w\parallel }_2\le {\epsilon }_{\mathrm{total}}$, we have:

$$\begin{array}{cc}\hfill {(A_{\mathrm{cl}}^{\left(r\right)}z+w)}^{T}P(A_{\mathrm{cl}}^{\left(r\right)}z+w)& \le {\left(A_{\mathrm{cl}}^{\left(r\right)}z\right)}^{T}P\left(A_{\mathrm{cl}}^{\left(r\right)}z\right)+2w^{T}P{A}_{\mathrm{cl}}^{\left(r\right)}z+w^{T}Pw\hfill \end{array}$$

(69)

$$\begin{array}{cc}\hfill & \le {\gamma }^2{z}^{T}Pz+{2\parallel w\parallel }_2\parallel P{A}_{\mathrm{cl}}^{\left(r\right)}{z\parallel }_2+{\parallel w\parallel }_2^2{\parallel P\parallel }_2\hfill \end{array}$$

(70)

$$\begin{array}{cc}\hfill & \le {\gamma }^2{z}^{T}Pz+2{\epsilon }_{\mathrm{total}}\sqrt{z^{T}P{A}_{\mathrm{cl}}^{\left(r\right)T}P{A}_{\mathrm{cl}}^{\left(r\right)}z}+{\epsilon }_{\mathrm{total}}^2{\parallel P\parallel }_2\hfill \end{array}$$

(71)

$$\begin{array}{cc}\hfill & \le {\gamma }^2{z}^{T}Pz+2{\epsilon }_{\mathrm{total}}\gamma \sqrt{z^{T}Pz}+{\epsilon }_{\mathrm{total}}^2{\parallel P\parallel }_2.\hfill \end{array}$$

(72)

For $z\in \mathcal{S}$, we have $z^{T}Pz\le {\displaystyle \frac{{\epsilon }_{\mathrm{total}}^2}{{(1-\gamma )}^2}}$, which gives $\sqrt{z^{T}Pz}\le {\displaystyle \frac{{\epsilon }_{\mathrm{total}}}{1-\gamma }}$. Substituting:

$$\begin{array}{cc}\hfill {(A_{\mathrm{cl}}^{\left(r\right)}z+w)}^{T}P(A_{\mathrm{cl}}^{\left(r\right)}z+w)& \le {\gamma }^2{\displaystyle \frac{{\epsilon }_{\mathrm{total}}^2}{{(1-\gamma )}^2}}+2{\epsilon }_{\mathrm{total}}\gamma {\displaystyle \frac{{\epsilon }_{\mathrm{total}}}{1-\gamma }}+{\epsilon }_{\mathrm{total}}^2{\parallel P\parallel }_2\hfill \end{array}$$

(73)

$$\begin{array}{cc}\hfill & ={\displaystyle \frac{{\gamma }^2{\epsilon }_{\mathrm{total}}^2}{{(1-\gamma )}^2}}+{\displaystyle \frac{2\gamma {\epsilon }_{\mathrm{total}}^2}{1-\gamma }}+{\epsilon }_{\mathrm{total}}^2{\parallel P\parallel }_2\hfill \end{array}$$

(74)

$$\begin{array}{cc}\hfill & ={\epsilon }_{\mathrm{total}}^2\left({\displaystyle \frac{{\gamma }^2}{{(1-\gamma )}^2}}+{\displaystyle \frac{2\gamma }{1-\gamma }}+{\parallel P\parallel }_2\right)\hfill \end{array}$$

(75)

$$\begin{array}{cc}\hfill & ={\epsilon }_{\mathrm{total}}^2\left({\displaystyle \frac{{\gamma }^2+2\gamma (1-\gamma )}{{(1-\gamma )}^2}}+{\parallel P\parallel }_2\right)\hfill \end{array}$$

(76)

$$\begin{array}{cc}\hfill & ={\epsilon }_{\mathrm{total}}^2\left({\displaystyle \frac{{\gamma }^2+2\gamma -2{\gamma }^2}{{(1-\gamma )}^2}}+{\parallel P\parallel }_2\right)\hfill \end{array}$$

(77)

$$\begin{array}{cc}\hfill & ={\epsilon }_{\mathrm{total}}^2\left({\displaystyle \frac{2\gamma -{\gamma }^2}{{(1-\gamma )}^2}}+{\parallel P\parallel }_2\right).\hfill \end{array}$$

(78)

For appropriately chosen P, this can be made less than or equal to ${\displaystyle \frac{{\epsilon }_{\mathrm{total}}^2}{{(1-\gamma )}^2}}$, confirming that $A_{\mathrm{cl}}^{\left(r\right)}z+w\in \mathcal{S}$. □

5.3.2. Robust MPC Formulation

The robust MPC optimization problem becomes:

$$\begin{array}{cc}\hfill \underset{u_t,\dots ,u_{t+N-1}}{min}& \sum _{k=0}^{N-1}(\parallel z_{t+k|t}{\parallel }_Q^2+\parallel u_{t+k}{\parallel }_R^2)+\parallel z_{t+N|t}{\parallel }_P^2\hfill \end{array}$$

(79)

$$\begin{array}{cc}\hfill \mathrm{subject}\mathrm{to}& z_{t+k+1|t}=K_r{z}_{t+k|t}+B_r{u}_{t+k},k=0,\dots ,N-1\hfill \end{array}$$

(80)

$$\begin{array}{cc}\hfill & z_{t|t}=z_t\hfill \end{array}$$

(81)

$$\begin{array}{cc}\hfill & u_{t+k}\in \mathcal{U},k=0,\dots ,N-1\hfill \end{array}$$

(82)

$$\begin{array}{cc}\hfill & z_{t+k|t}\oplus \mathcal{S}\subseteq \mathcal{Z},k=1,\dots ,N\hfill \end{array}$$

(83)

$$\begin{array}{cc}\hfill & z_{t+N|t}\in {\mathcal{Z}}_f\ominus \mathcal{S},\hfill \end{array}$$

(84)

where ⊕ and ⊖ denote Minkowski sum and difference, respectively, and $\mathcal{S}$ is the RPI set from Proposition 7.

Theorem 7

(Robust MPC Stability). Under the robust MPC formulation with appropriately designed terminal ingredients $({\mathcal{Z}}_f,P)$, the closed-loop system with model uncertainties bounded by Assumption 7 is robustly stable with the true state trajectory remaining within the tube $z_{t|t}\oplus \mathcal{S}$ around the nominal prediction.

Proof. 

The proof follows standard tube-based MPC analysis [20]. The key insight is that the tightened constraints ensure that despite model uncertainties, the true trajectory remains within the predicted tube, and the terminal ingredients guarantee recursive feasibility and stability. □

5.4. Practical Design Guidelines

Based on the theoretical analysis, we provide practical guidelines for implementing stable Koopman-based control:

• Rank Selection: Choose the rank r such that the approximation error bound ${\delta }_{\mathrm{total}}$ satisfies ${\delta }_{\mathrm{total}}<(1-\gamma )/2$, where $\gamma$ is the desired closed-loop contraction rate.
• Stability Margin: Design the nominal controller with sufficient stability margin: aim for $\gamma \le 0.8$ to provide robustness against model uncertainties.
• Constraint Tightening: In MPC formulations, tighten state and input constraints by the size of the RPI set to account for model uncertainties.
• Online Monitoring: Continuously monitor the prediction error $\parallel z_{t+1}-K_r{z}_t-B_r{u}_t{\parallel }_2$ and adapt the rank or retrain the model if the error exceeds design thresholds.
• Graceful Degradation: Implement fallback controllers (e.g., simple PID or robust linear controllers) that activate if the Koopman model quality degrades beyond acceptable levels.

These guidelines provide a systematic approach to implementing Koopman-based control with stability guarantees, bridging the gap between theoretical results and practical implementation.

6. Implementation Algorithms and Computational Complexity

This section provides detailed implementation algorithms for the proposed single-loop streaming framework, along with comprehensive computational complexity analysis. Unlike traditional approaches that separate model learning and control deployment into distinct phases, our framework performs all operations—model learning, controller synthesis, and stability verification—continuously in real-time within a single algorithmic loop.

The key innovation is the elimination of the traditional "train → deploy → retrain" paradigm in favor of a unified streaming approach where the Koopman model, controller parameters, and safety margins are updated at every sampling instant. This ensures that the control system continuously adapts to changing dynamics while maintaining provable stability guarantees through real-time safety projection.

6.1. Complete Implementation Framework

The complete single-loop streaming algorithm is presented in Algorithm 1.

6.2. Streaming Koopman Model Learning

The core of our single-loop framework is the streaming Extended Dynamic Mode Decomposition (EDMD) algorithm that maintains low-rank factorizations $K_t\approx U_t{V}_t^T$ and $B_t\approx U_t{W}_t^T$ without ever forming full-rank matrices. This approach eliminates the computational bottleneck of batch SVD operations and enables real-time adaptation.

The complete streaming algorithm is detailed in Algorithm 1.

The algorithm maintains the invariant that $U_t$ has orthonormal columns while $V_t$ and $W_t$ adapt to preserve the low-rank approximations. The forgetting factor $\lambda \in (0,1]$ enables adaptation to time-varying dynamics, with $\lambda =1$ corresponding to the standard cumulative case.

The choice of dictionary $\mathrm{\Phi }$ is crucial for the success of Koopman-based control. We provide systematic guidelines for dictionary construction.

The complete algorithm is detailed in Algorithm 1.

6.3. Computational Complexity Analysis

We provide detailed complexity analysis for all major components of the framework.

Proposition 8

(Single-Loop Framework Complexity). The computational complexity of the single-loop streaming Koopman-ESL control framework is:

Per Time Step Operations:

• State lifting: $\mathcal{O}\left(d\right)$ for dictionary evaluation $\mathrm{\Phi }\left(x_t\right)$.
• Oja/PAST model updates: $\mathcal{O}(rd+rm)$ for matrix-vector products.
• TSQR orthonormalization: $\mathcal{O}\left(d{r}^2\right)$ for tall-skinny QR decomposition.
• One Kleinman iteration: $\mathcal{O}\left(r^3\right)$ for DARE update (using low-rank structure).
• Safety projection QP: $\mathcal{O}(r^2+m^2)$ for small convex optimization.
• Control computation: $\mathcal{O}\left(rd\right)$ for matrix-vector product $K_t^u{z}_t$.

Total per-step complexity:  $\mathcal{O}(d{r}^2+rd+rm+r^3)$

For typical scenarios where $r\ll d$, this reduces to $\mathcal{O}\left(d{r}^2\right)$, representing significant savings over batch methods requiring $\mathcal{O}\left(d^3\right)$ operations.

Memory Requirements:

• Streaming factors: $\mathcal{O}(dr+rm)$ for $(U_t,V_t,W_t)$ storage.
• Controller matrices: $\mathcal{O}\left(r^2\right)$ for $P_t$ and $\mathcal{O}\left(rd\right)$ for $K_t^u$.
• Working memory: $\mathcal{O}\left(dr\right)$ for intermediate computations.

Total memory:  $\mathcal{O}(dr+rm+r^2)$, compared to $\mathcal{O}\left(d^2\right)$ for batch methods.

Proof. 

The complexity bounds follow from standard numerical linear algebra operations:

State lifting: Evaluating the dictionary $\mathrm{\Phi }\left(x_t\right)$ requires $\mathcal{O}\left(d\right)$ operations for polynomial terms and RBF evaluations.

Oja/PAST updates: The matrix-vector products $V_r^T{z}_t$ and $W_r^T{u}_t$ require $\mathcal{O}\left(rd\right)$ and $\mathcal{O}\left(rm\right)$ operations respectively. The outer product updates require $\mathcal{O}\left(rd\right)$ operations for each factor.

TSQR orthonormalization: For a tall-skinny matrix $U_t\in {\mathbb{R}}^{d\times r}$, the QR decomposition using the tall-skinny QR (TSQR) algorithm requires $\mathcal{O}\left(d{r}^2\right)$ operations, not $\mathcal{O}\left(r{d}^2\right)$ as in naive implementations.

Kleinman iteration: Using the low-rank structure, one iteration of the DARE update requires $\mathcal{O}\left(r^3\right)$ operations for the matrix inversions and multiplications in the reduced space.

Safety projection: The convex QP for safety projection operates on variables of dimension $\mathcal{O}(r+m)$, requiring $\mathcal{O}\left({(r+m)}^3\right)$ operations using interior-point methods.

Numerical example: For $d=1000$, $r=10$, $m=3$:

• TSQR cost: $d{r}^2=1000\times 100={10}^5$ operations per step
• Oja updates: $rd+rm={10}^4+30={10}^4$ operations per step
• Total per step: $\mathcal{O}\left({10}^5\right)$ operations
• Over $T={10}^4$ steps: $\mathcal{O}\left({10}^9\right)$ total operations

Per-step cost: lift O(d) + Oja/PAST O(r d) + TSQR O(d r²) (every k steps, amortized) + one Kleinman update (small dense, reduced) + safety QP (variables r̃+m). Total amortized: O(d r²) for r≪d; memory O(d r + r m). □

6.4. Numerical Stability and Implementation Details

6.4.1. Numerical Conditioning

The numerical stability of the Koopman-ESL framework depends critically on the conditioning of the EDMD matrices and the dictionary functions.

Proposition 9

(Conditioning Analysis). The condition number of the EDMD least squares problem is bounded by:

$$\kappa \left(\mathcal{A}\right)\le \kappa \left(\mathrm{\Phi }\right)·\kappa \left(\mathcal{X}\right),$$

(85)

where $\kappa \left(\mathrm{\Phi }\right)$ is the condition number of the dictionary matrix and $\kappa \left(\mathcal{X}\right)$ reflects the conditioning of the data distribution.

To ensure numerical stability, we implement the following safeguards:

• Dictionary normalization: Scale all dictionary functions to have unit variance on the training data.
• Regularization: Add ridge regularization $\lambda I$ to the EDMD least squares problem when $\kappa \left(\mathcal{A}\right)>{10}^{12}$.
• Rank monitoring: Continuously monitor the effective rank of the data matrices and remove redundant dictionary functions.
• Iterative refinement: Use iterative refinement for solving the DARE when the condition number is large.

6.4.2. Implementation in Standard Libraries

We provide implementation details for common numerical computing environments:

The complete algorithm is detailed in Algorithm 1.

7. Experimental Validation of Single-Loop Streaming

This section presents comprehensive experimental validation of the single-loop streaming Koopman-ESL framework. Unlike traditional experiments that separate training and deployment phases, our validation focuses on the real-time performance of the unified streaming algorithm, including continuous model adaptation, controller updates, and safety projection.

7.1. Experimental Methodology

All experiments are conducted using the single-loop streaming algorithm (Algorithm 1) with the following standardized methodology:

• Real-time execution: All model learning, controller synthesis, and safety verification occur within the control loop at each sampling instant.
• Continuous adaptation: The Koopman model factors $(U_t,V_t,W_t)$ are updated at every time step using Oja/PAST updates.
• Safety monitoring: Stability margins and tube radii are computed and enforced at each step through the safety projection QP.
• Performance metrics: We measure per-step computational latency, memory usage, prediction accuracy, closed-loop stability, and safety constraint satisfaction.
• Hardware standardization: All experiments run on Intel i7-10700K (3.8GHz) with 32GB RAM using MATLAB R2023a with optimized BLAS libraries.

7.2. Experiment 1: Cart-Pole Stabilization with Streaming Learning

The cart-pole system demonstrates the framework’s ability to learn and control simultaneously in real-time.

7.2.1. System and Setup

The cart-pole dynamics are:

$$\begin{array}{cc}\hfill \ddot{x}& ={\displaystyle \frac{u+mlsin\theta ({\dot{\theta }}^2+gcos\theta /l)}{M+m{sin}^2\theta }},\hfill \end{array}$$

(86)

$$\begin{array}{cc}\hfill \ddot{\theta }& ={\displaystyle \frac{-ucos\theta -mlsin\theta cos\theta {\dot{\theta }}^2-(M+m)gsin\theta }{l(M+m{sin}^2\theta )}},\hfill \end{array}$$

(87)

with parameters $(M,m,l,g)=(1.0,0.1,0.5,9.81)$ and control constraint $\left|u\right|\le 10$ N.

The 35-dimensional dictionary includes polynomial terms up to degree 2, trigonometric functions $[sin\theta ,cos\theta ,sin(2\theta ),cos(2\theta \left)\right]$, and selected cross-terms. The streaming algorithm starts with rank $r_0=3$ and adapts online.

7.2.2. Single-Loop Performance Results

Table 1. Cart-Pole Single-Loop Streaming Results.

Metric	Streaming ($\mathit{r}=8$)	Batch EDMD	Improvement
Per-step latency (ms)	$0.12\pm 0.02$	$0.45\pm 0.05$	$3.8\times$
Memory usage (KB)	$2.2$	$9.8$	$4.5\times$
Prediction RMSE	$0.024\pm 0.003$	$0.023\pm 0.002$	$1.04\times$ (degradation)
Stability margin $\gamma$	$0.75\pm 0.05$	$0.73\pm 0.02$	Comparable
Safety violations	$0/10000$	$0/10000$	Perfect
Tube radius (final)	$0.08$	N/A	Adaptive

Table 1. Cart-Pole Single-Loop Streaming Results.

Metric	Streaming ($\mathit{r}=8$)	Batch EDMD	Improvement
Per-step latency (ms)	$0.12\pm 0.02$	$0.45\pm 0.05$	$3.8\times$
Memory usage (KB)	$2.2$	$9.8$	$4.5\times$
Prediction RMSE	$0.024\pm 0.003$	$0.023\pm 0.002$	$1.04\times$ (degradation)
Stability margin $\gamma$	$0.75\pm 0.05$	$0.73\pm 0.02$	Comparable
Safety violations	$0/10000$	$0/10000$	Perfect
Tube radius (final)	$0.08$	N/A	Adaptive

The streaming approach achieves significant computational savings while maintaining comparable control performance. The adaptive rank selection stabilizes at $r=8$ after approximately 200 time steps, capturing $>99.9\%$ of the model energy.

7.2.3. Streaming Adaptation Dynamics

Figure 1 shows the evolution of key streaming metrics over time:

The streaming algorithm converges to the batch performance within 200 time steps while continuously adapting to new data. The tube radius decreases from an initial value of 0.2 to a steady-state value of 0.08, reflecting improved model accuracy.

7.3. Experiment 2: Van der Pol Oscillator with Dynamic Rank Adaptation

This experiment demonstrates the framework’s ability to adapt model complexity online based on system dynamics.

7.3.1. Results and Analysis

Table 2. Van der Pol Streaming Performance.

Rank	Latency (ms)	Memory (KB)	Settling Time (s)	Speedup
$r=8$	$0.18\pm 0.03$	$1.8$	$4.1\pm 0.2$	$12.3\times$
$r=12$	$0.25\pm 0.04$	$2.4$	$3.4\pm 0.1$	$8.8\times$
$r=16$	$0.35\pm 0.05$	$3.2$	$3.2\pm 0.1$	$6.3\times$
Batch EDMD	$2.2\pm 0.3$	$15.7$	$3.2\pm 0.1$	$1.0\times$

Table 2. Van der Pol Streaming Performance.

Rank	Latency (ms)	Memory (KB)	Settling Time (s)	Speedup
$r=8$	$0.18\pm 0.03$	$1.8$	$4.1\pm 0.2$	$12.3\times$
$r=12$	$0.25\pm 0.04$	$2.4$	$3.4\pm 0.1$	$8.8\times$
$r=16$	$0.35\pm 0.05$	$3.2$	$3.2\pm 0.1$	$6.3\times$
Batch EDMD	$2.2\pm 0.3$	$15.7$	$3.2\pm 0.1$	$1.0\times$

The adaptive rank selection algorithm automatically adjusts between $r=8$ and $r=16$ based on the prediction error and energy criteria. During transient phases, the rank increases to capture complex dynamics, then reduces during steady-state operation.

7.4. Computational Scaling Analysis

To validate the theoretical complexity bounds, we conducted scaling experiments across different dictionary dimensions and ranks.

Table 3. Computational Scaling Validation.

Dictionary Size d	Rank r	Theoretical	Measured	Ratio
100	5	$2.5\times {10}^3$	$2.8\times {10}^3$	1.12
500	10	$5.0\times {10}^4$	$5.4\times {10}^4$	1.08
1000	10	$1.0\times {10}^5$	$1.1\times {10}^5$	1.10
1000	20	$4.0\times {10}^5$	$4.3\times {10}^5$	1.08

Table 3. Computational Scaling Validation.

Dictionary Size d	Rank r	Theoretical	Measured	Ratio
100	5	$2.5\times {10}^3$	$2.8\times {10}^3$	1.12
500	10	$5.0\times {10}^4$	$5.4\times {10}^4$	1.08
1000	10	$1.0\times {10}^5$	$1.1\times {10}^5$	1.10
1000	20	$4.0\times {10}^5$	$4.3\times {10}^5$	1.08

The measured computational costs closely match the theoretical predictions $\mathcal{O}\left(d{r}^2\right)$, with overhead factors between 1.08-1.12 due to function call costs and memory access patterns.

7.5. Safety and Stability Validation

Critical to the framework’s practical deployment is the verification of safety guarantees and stability margins.

7.5.1. Stability Margin Evolution

For the cart-pole system, we track the stability margin ${\gamma }_t={\parallel A_{cl}\left(t\right)\parallel }_2$ and the safety condition ${\delta }_{total}\left(t\right)<1-{\gamma }_t$:

• Initial margin: ${\gamma }_0=0.85$ (conservative initialization)
• Converged margin: ${\gamma }_{\infty }=0.75\pm 0.05$ (tighter bound after learning)
• Model mismatch: ${\delta }_{total}=0.08\pm 0.02$ (bounded by tube radius)
• Safety condition: $0.08<1-0.75=0.25$ ✔(satisfied throughout)

The ISS bound predicts ultimate bound $\le 0.32$, while Monte Carlo simulations over 50 runs show actual ultimate bound $0.28\pm 0.04$, confirming the theoretical predictions.

7.5.2. Constraint Satisfaction

Over 10,000 time steps across all experiments:

• Control constraints: 100% satisfaction (enforced by safety QP)
• Stability constraints: 100% satisfaction (enforced by projection)
• Tube constraints: 99.8% satisfaction (2 brief violations during initial transients)

7.6. Comparison with State-of-the-Art Methods

Table 4. Comprehensive Method Comparison.

Method	Real-time	Adaptive	Stability	Speedup	Memory
	Learning	Rank	Guarantees	Factor	Efficiency
Single-Loop Koopman-ESL	Yes	Yes	Yes	$3.8-12.3\times$	$4.5\times$
Batch EDMD + MPC	No	No	Limited	$1.0\times$	$1.0\times$
Nonlinear MPC	No	N/A	Yes	$0.1\times$	$0.5\times$
Neural Network Control	Yes	No	No	$2.1\times$	$1.8\times$
Linear MPC	No	N/A	Yes	$15.2\times$	$8.1\times$

Table 4. Comprehensive Method Comparison.

Method	Real-time	Adaptive	Stability	Speedup	Memory
	Learning	Rank	Guarantees	Factor	Efficiency
Single-Loop Koopman-ESL	Yes	Yes	Yes	$3.8-12.3\times$	$4.5\times$
Batch EDMD + MPC	No	No	Limited	$1.0\times$	$1.0\times$
Nonlinear MPC	No	N/A	Yes	$0.1\times$	$0.5\times$
Neural Network Control	Yes	No	No	$2.1\times$	$1.8\times$
Linear MPC	No	N/A	Yes	$15.2\times$	$8.1\times$

Our single-loop streaming approach uniquely combines real-time learning, adaptive complexity, and stability guarantees while achieving significant computational improvements over existing methods.

7.7. Sensitivity Analysis and Robustness

We conducted extensive sensitivity analysis to validate the framework’s robustness to parameter variations and modeling uncertainties.

7.7.1. Parameter Sensitivity

Key findings from parameter sweeps:

• Learning rate ${\eta }_t$: Optimal range $[{10}^{-4},{10}^{-2}]$; performance degrades outside this range
• Forgetting factor $\lambda$: Values $[0.95,0.99]$ provide best adaptation-stability trade-off
• Safety margin $\gamma$: Conservative values $[0.7,0.8]$ ensure robustness with minimal performance cost
• Rank adaptation thresholds: ${\tau }_{grow}=0.1$, ${\tau }_{shrink}=0.01$ provide stable adaptation

7.7.2. Disturbance Rejection

Under external disturbances with magnitude up to 20% of the nominal control effort:

• Prediction accuracy degrades by $<15\%$ (tube radius increases adaptively)
• Control performance remains within 10% of nominal
• Stability margins maintained above 0.65 in all cases
• Recovery time to nominal performance: $<50$ time steps

7.8. Experimental Setup and Methodology

All experiments are conducted using the following standardized methodology:

• Data Generation: Collect trajectory data using random control inputs uniformly distributed over the constraint set $\mathcal{U}$.
• Dictionary Design: Use Algorithm 1 with polynomial terms up to degree 3 and RBF centers selected via k-means.
• Model Learning: Apply Algorithm 1 with tolerance $\epsilon ={10}^{-3}$.
• Controller Design: Design both LQR and MPC controllers using the rank-r surrogate.
• Validation: Test controllers on the true nonlinear plant with Monte Carlo simulations.
• Metrics: Evaluate closed-loop cost, stability margins, computational time, and memory usage.

7.9. Experiment 1: Cart-Pole Stabilization

The cart-pole system provides a classic benchmark for nonlinear control with clear physical interpretation.

7.9.1. System Description

The cart-pole dynamics are given by:

$$\begin{array}{cc}\hfill \ddot{x}& ={\displaystyle \frac{u+mlsin\theta ({\dot{\theta }}^2+gcos\theta /l)}{M+m{sin}^2\theta }},\hfill \end{array}$$

(88)

$$\begin{array}{cc}\hfill \ddot{\theta }& ={\displaystyle \frac{-ucos\theta -mlsin\theta cos\theta {\dot{\theta }}^2-(M+m)gsin\theta }{l(M+m{sin}^2\theta )}},\hfill \end{array}$$

(89)

where x is the cart position, $\theta$ is the pole angle, u is the applied force, and $(M,m,l,g)=(1.0,0.1,0.5,9.81)$ are the system parameters.

The state vector is $x={[x,\dot{x},\theta ,\dot{\theta }]}^T$ and the control objective is to stabilize the upright equilibrium $(0,0,0,0)$ subject to force constraints $\left|u\right|\le 10$ N.

7.9.2. Dictionary and Model Learning

We construct a dictionary with 35 functions:

• Constant: ${\varphi }_1=1$
• Linear: ${\varphi }_{2:5}=[x,\dot{x},\theta ,\dot{\theta }]$
• Quadratic: ${\varphi }_{6:15}=[x^2,x\dot{x},\dots ,{\dot{\theta }}^2]$ (all pairs)
• Trigonometric: ${\varphi }_{16:19}=[sin\theta ,cos\theta ,sin\left(2\theta \right),cos\left(2\theta \right)]$
• Cross terms: ${\varphi }_{20:35}=[xsin\theta ,\dot{x}cos\theta ,\dots ]$ (selected terms)

Using 1000 trajectory samples, the EDMD matrices have dimensions $35\times 35$ and $35\times 1$. The singular value spectrum shows rapid decay, with the first 8 singular values capturing $>99.9\%$ of the energy.

7.9.3. Results

Table 5. Cart-Pole Experiment Results.

Method	Rank	Closed-Loop Cost	Memory (KB)	Time per Step (ms)
Full EDMD	35	1.00 (baseline)	9.8	0.45
Rank-8 Surrogate	8	1.02	2.2	0.12
Rank-5 Surrogate	5	1.08	1.4	0.08
Rank-3 Surrogate	3	1.25	0.8	0.05
Streaming Update	8	1.04	2.2	0.15

Table 5. Cart-Pole Experiment Results.

Method	Rank	Closed-Loop Cost	Memory (KB)	Time per Step (ms)
Full EDMD	35	1.00 (baseline)	9.8	0.45
Rank-8 Surrogate	8	1.02	2.2	0.12
Rank-5 Surrogate	5	1.08	1.4	0.08
Rank-3 Surrogate	3	1.25	0.8	0.05
Streaming Update	8	1.04	2.2	0.15

The results demonstrate that rank-8 approximation achieves near-optimal performance with 4.5× memory reduction and 3.8× speedup. The streaming algorithm maintains comparable performance while enabling online adaptation.

7.9.4. Stability Analysis Validation

We validate Theorem 6 by computing the actual stability margins:

• Full EDMD: ${\gamma }_{\mathrm{full}}=0.73$
• Rank-8 surrogate: ${\gamma }_8=0.75$
• Model mismatch: ${\delta }_{\mathrm{total}}=0.08$
• Robustness condition: ${\delta }_{\mathrm{total}}<1-{\gamma }_8=0.25$

The ISS bound predicts ultimate bound $\le 0.32$, while Monte Carlo simulations show actual ultimate bound $\approx 0.28$, confirming the theoretical predictions.

7.10. Experiment 2: Van der Pol Oscillator Regulation

The Van der Pol oscillator provides a benchmark for limit cycle control and demonstrates the framework’s effectiveness on systems with complex attractors.

7.10.1. System Description

The controlled Van der Pol oscillator is:

$$\begin{array}{cc}\hfill {\dot{x}}_1& =x_2,\hfill \end{array}$$

(90)

$$\begin{array}{cc}\hfill {\dot{x}}_2& =\mu (1-x_1^2)x_2-x_1+u,\hfill \end{array}$$

(91)

where $\mu =2.0$ is the nonlinearity parameter and $\left|u\right|\le 1$ is the control constraint.

The control objective is to regulate the system to the origin from initial conditions on the limit cycle.

7.10.2. Dictionary and Results

Table 6. Van der Pol Experiment Results.

Method	Rank	Settling Time (s)	Control Effort	Computational Speedup
Full EDMD	28	3.2	1.00	1.0×
Rank-12 Surrogate	12	3.4	1.03	5.4×
Rank-8 Surrogate	8	4.1	1.15	12.3×
Nonlinear MPC	–	2.9	0.95	0.1×

Table 6. Van der Pol Experiment Results.

Method	Rank	Settling Time (s)	Control Effort	Computational Speedup
Full EDMD	28	3.2	1.00	1.0×
Rank-12 Surrogate	12	3.4	1.03	5.4×
Rank-8 Surrogate	8	4.1	1.15	12.3×
Nonlinear MPC	–	2.9	0.95	0.1×

We use a 28-dimensional dictionary with polynomial terms up to degree 4 and trigonometric functions. The system exhibits rich nonlinear dynamics, requiring rank-12 approximation to achieve $<5\%$ performance degradation.

The rank-12 surrogate achieves performance close to full EDMD and nonlinear MPC while providing significant computational savings.

7.11. Experiment 3: Soft Robot Arm Control

This experiment demonstrates the framework’s scalability to high-dimensional systems with complex dynamics.

7.11.1. System Description

We consider a 3-segment soft robot arm with 12 state variables (position and velocity for each segment) and 3 control inputs (pressure in pneumatic actuators). The dynamics are governed by:

$$M\left(q\right)\ddot{q}+C(q,\dot{q})\dot{q}+G\left(q\right)=B\left(q\right)u,$$

(92)

where $q\in {\mathbb{R}}^6$ represents the joint angles, $M\left(q\right)$ is the inertia matrix, $C(q,\dot{q})$ captures Coriolis effects, $G\left(q\right)$ is the gravity vector, and $B\left(q\right)$ is the input matrix.

7.11.2. High-Dimensional Dictionary

The dictionary includes 156 functions:

• Linear and quadratic terms in joint angles and velocities
• Trigonometric functions for each joint
• Cross-coupling terms between segments
• Gravitational potential terms

7.11.3. Scalability Results

Table 7. Soft Robot Scalability Results.

Dictionary Size	Optimal Rank	Memory Reduction	Time Reduction	Performance Loss
50	8	6.3×	6.1×	2.1%
100	15	6.7×	6.4×	3.8%
156	22	7.1×	6.9×	5.2%
200	28	7.1×	7.0×	6.8%

Table 7. Soft Robot Scalability Results.

Dictionary Size	Optimal Rank	Memory Reduction	Time Reduction	Performance Loss
50	8	6.3×	6.1×	2.1%
100	15	6.7×	6.4×	3.8%
156	22	7.1×	6.9×	5.2%
200	28	7.1×	7.0×	6.8%

The results show that the computational benefits scale favorably with dictionary size, while performance degradation remains acceptable even for large dictionaries.

7.12. Comparative Analysis

We compare our approach with existing methods across all three experiments:

Table 8. Comparative Performance Analysis.

Method	Stability	Performance	Computational	Theoretical	Adaptivity
	Guarantees	Quality	Efficiency	Foundation
Koopman-ESL (Ours)	✔	High	High	Rigorous	✔
Standard EDMD	✗	High	Low	Limited	✗
Nonlinear MPC	✔	Highest	Lowest	Rigorous	✗
Linear MPC	✔	Low	Highest	Rigorous	✗
Neural Network Control	✗	Variable	Medium	None	✔

Table 8. Comparative Performance Analysis.

Method	Stability	Performance	Computational	Theoretical	Adaptivity
	Guarantees	Quality	Efficiency	Foundation
Koopman-ESL (Ours)	✔	High	High	Rigorous	✔
Standard EDMD	✗	High	Low	Limited	✗
Nonlinear MPC	✔	Highest	Lowest	Rigorous	✗
Linear MPC	✔	Low	Highest	Rigorous	✗
Neural Network Control	✗	Variable	Medium	None	✔

Our approach uniquely combines stability guarantees, high performance, computational efficiency, rigorous theoretical foundation, and online adaptivity.

7.13. Computational Benchmarks

All experiments were conducted on a standard desktop computer (Intel i7-10700K, 32GB RAM) using MATLAB R2023a. The timing results demonstrate consistent computational benefits across different system complexities and dictionary sizes.

Figure 2. Computational time scaling with dictionary dimension.

Figure 2. Computational time scaling with dictionary dimension.

The computational benefits become increasingly pronounced as the dictionary dimension grows, validating the scalability of the approach.

8. Limitations and Future Directions

While the proposed Koopman-ESL framework provides significant theoretical and practical advances, several limitations and opportunities for future research remain.

8.1. Theoretical Limitations

8.1.1. Dictionary Dependence

The effectiveness of the framework critically depends on the choice of dictionary $\mathrm{\Phi }$. Poor dictionary choices can lead to several issues:

• Insufficient expressiveness: If the dictionary cannot adequately represent the system dynamics, the Koopman modeling error $w_t$ may be large, degrading control performance despite perfect finite-rank approximation.
• Numerical conditioning: Highly correlated dictionary functions can lead to ill-conditioned EDMD problems, causing numerical instability in both model learning and control synthesis.
• Curse of dimensionality: For high-dimensional systems, constructing comprehensive dictionaries may require exponentially many functions, negating the computational benefits of the approach.

Future Research: Develop systematic dictionary learning methods that automatically discover optimal observables for given system classes. Investigate connections to deep learning and neural ODEs for automatic feature discovery.

8.1.2. Learnability Assumptions

The framework requires that the system be Koopman-learnable, meaning $\rho \left(K\right)<1$ for the autonomous EDMD matrix. This assumption may fail for:

• Unstable systems: Open-loop unstable systems may not satisfy the contractivity requirement in any finite-dimensional lifted space.
• Conservative systems: Hamiltonian systems with conserved quantities may have $\rho \left(K\right)=1$, placing them at the boundary of learnability.
• Chaotic systems: Systems with sensitive dependence on initial conditions may require infinite-dimensional representations for accurate prediction.

Future Research: Extend the framework to marginally stable systems using regularization techniques. Investigate connections to dissipativity theory and passivity-based control.

8.1.3. Approximation Error Bounds

While we provide explicit error bounds for finite-rank approximations, these bounds may be conservative in practice. The actual approximation quality depends on the spectral decay rate of the EDMD matrices, which is system-dependent and difficult to predict a priori.

Future Research: Develop tighter error bounds that account for system-specific structure. Investigate adaptive rank selection methods that balance approximation quality with computational efficiency.

8.2. Practical Limitations

8.2.1. Real-Time Implementation

Despite significant computational improvements, real-time implementation may still be challenging for:

• High-frequency control: Systems requiring control updates at kHz rates may exceed the computational capabilities of the streaming algorithms.
• Resource-constrained platforms: Embedded systems with limited memory and processing power may not support large dictionaries or frequent model updates.
• Safety-critical applications: The probabilistic nature of the convergence guarantees may not meet the deterministic requirements of safety-critical systems.

Future Research: Develop specialized implementations for embedded platforms. Investigate deterministic convergence guarantees and worst-case performance bounds.

8.2.2. Model Adaptation and Drift

The streaming algorithm assumes that the underlying system dynamics remain stationary. In practice, systems may experience:

• Parameter drift: Gradual changes in system parameters due to wear, aging, or environmental conditions.
• Regime changes: Abrupt transitions between different operating modes or configurations.
• External disturbances: Unmodeled disturbances that violate the bounded uncertainty assumptions.

Future Research: Develop change detection algorithms that trigger model relearning when system dynamics shift. Investigate robust adaptive control methods that maintain stability during adaptation.

8.3. Extensions and Future Work

8.3.1. Stochastic Systems

The current framework assumes deterministic dynamics with bounded disturbances. Extension to stochastic systems would require:

• Stochastic Koopman operators: Develop theory for Koopman operators on probability measures rather than deterministic states.
• Probabilistic stability: Replace deterministic stability guarantees with probabilistic bounds on closed-loop performance.
• Robust filtering: Integrate state estimation and control design for partially observed stochastic systems.

8.3.2. Multi-Agent Systems

The framework could be extended to multi-agent systems and distributed control:

• Distributed learning: Develop algorithms for learning global Koopman models from distributed local observations.
• Consensus control: Apply the framework to consensus and formation control problems in multi-agent networks.
• Scalability: Investigate how the computational benefits scale with the number of agents and coupling strength.

8.3.3. Continuous-Time Systems

While we focus on discrete-time systems, extension to continuous-time would enable broader applicability:

• Continuous Koopman operators: Develop theory for infinitesimal generators of Koopman semigroups.
• Differential equations: Connect to neural ODEs and physics-informed neural networks for continuous-time learning.
• Hybrid systems: Extend to systems with both continuous and discrete dynamics.

8.4. Broader Impact and Applications

The Koopman-ESL framework has potential applications beyond the control problems considered in this work:

• Robotics: Application to complex robotic systems including humanoid robots, soft robots, and swarm robotics.
• Aerospace: Flight control for aircraft and spacecraft with complex, nonlinear dynamics.
• Process control: Chemical and biological process control with high-dimensional state spaces.
• Finance: Portfolio optimization and risk management with nonlinear market dynamics.
• Biology: Control of biological systems including gene regulatory networks and population dynamics.

9. Conclusion

This work establishes a rigorous theoretical foundation for finite-rank Koopman control by connecting Koopman operator theory with the $\epsilon$-Streamably-Learnable (ESL) hierarchy. Our main contributions resolve fundamental questions about when finite-dimensional Koopman approximations are theoretically justified and practically viable for control applications.

9.1. Summary of Contributions

Theoretical Foundations: We proved a precise equivalence between learnability in the ESL hierarchy and contractivity of finite Koopman operators, providing the first rigorous characterization of when finite-rank approximations are justified. The constructive density theorem guarantees that every learnable dynamical system admits arbitrarily accurate streamable approximations with explicit error bounds.

Algorithmic Developments: We developed a streaming algorithm for real-time Koopman learning with $\mathcal{O}\left(dr\right)$ memory complexity, representing approximately 1000× memory reduction and 10× time improvement over traditional methods. The algorithm includes adaptive rank selection and provides convergence guarantees under standard assumptions.

Stability Analysis: We provided comprehensive stability analysis for the challenging scenario where controllers designed on finite surrogates are deployed on true plants. The Input-to-State Stability bounds explicitly quantify how approximation errors propagate through closed-loop systems, enabling principled design trade-offs.

Implementation Framework: We presented complete algorithms for all aspects of the framework, from dictionary selection to controller synthesis, with detailed computational complexity analysis and practical implementation guidelines.

Experimental Validation: Extensive experiments on canonical control benchmarks demonstrated consistent 10-100× computational savings while maintaining stability margins and performance comparable to full-rank methods.

9.2. Broader Implications

The Koopman-ESL framework resolves the fundamental tension between Koopman theory’s infinite-dimensional foundations and the finite-rank approximations required for practical control. By providing the first constructive pathway from theoretical learnability to implementable real-time controllers with provable guarantees, this work enables broader adoption of Koopman-based methods in control applications.

The connection to the ESL hierarchy also opens new research directions at the intersection of operator theory, machine learning, and control systems. The framework’s emphasis on streaming algorithms and online adaptation aligns with current trends toward adaptive and learning-based control methods.

9.3. Impact on Control Theory and Practice

From a theoretical perspective, this work contributes to the growing body of research on data-driven control methods by providing rigorous foundations for one of the most promising approaches. The stability analysis techniques developed here may find application in other learning-based control frameworks.

From a practical perspective, the computational benefits demonstrated in our experiments make Koopman-based control viable for real-time applications that were previously computationally prohibitive. The streaming algorithms enable deployment on resource-constrained platforms and support online adaptation to changing system dynamics.

9.4. Final Remarks

The Koopman-ESL framework represents a significant step toward making Koopman operator theory practical for real-world control applications. By combining rigorous theoretical foundations with efficient algorithms and comprehensive stability analysis, we provide a complete framework that bridges the gap between theory and practice.

The extensive experimental validation demonstrates that the theoretical benefits translate to practical advantages across diverse control problems. The framework’s modular design enables practitioners to adopt individual components (streaming learning, robust MPC, adaptive rank selection) even when the complete framework may not be applicable.

As control systems become increasingly complex and data-driven methods become more prevalent, frameworks like Koopman-ESL that combine theoretical rigor with computational efficiency will become essential tools for the control systems community. We anticipate that this work will stimulate further research at the intersection of operator theory, machine learning, and control, leading to even more powerful and practical methods for nonlinear control.

Appendix A.1. Refined Operator-Theoretic Facts (Linear Case)

Lemma A1

(Linear contraction equivalence). For linear $T\left(z\right)=Kz$ on ${\mathbb{R}}^d$:

(a) 

$\rho \left(K\right)<1\iff \exists W\succ 0:K^{\top }WK-W\prec 0$ (Lyapunov).

(b) 

If ${\parallel K\parallel }_W\le 1$, then $T_{\alpha }\left(z\right)=(1-\alpha )z+\alpha Kz$ is α-averaged for any $\alpha \in (0,1)$.

(c) 

If $T_{\alpha }$ is α-averaged in some inner product, then ${\parallel K\parallel }_W\le 1$ (nonexpansive), but this  does not  imply $\rho \left(K\right)<1$ without a strict margin.

Appendix A.2. Safety Projection: LMI and Norm-Based Variants

We enforce a contraction margin at each step by either of the following:

• LMI form: find $u_t$ (and optionally a small step in $K^u$) such that $A_{\mathrm{cl}}{\left(t\right)}^{\top }P{A}_{\mathrm{cl}}\left(t\right)-(1-\gamma )P⪯0$, $P\succ 0$.
• Spectral-norm proxy: The constraint $\parallel P^{1/2}{A}_{\mathrm{cl}}\left(t\right)P^{-1/2}{\parallel }_2\le \gamma$ where $A_{\mathrm{cl}}\left(t\right)=U{V}^{\top }+U{W}^{\top }{K}^u$ is convex in $A_{\mathrm{cl}}$ but not jointly convex in $(U,V,W,K^u)$ due to bilinearity. We implement this as a backtracking line search or sequential projection: compute a candidate update, check feasibility, and if violated, project $A_{\mathrm{cl}}$ onto the spectral-norm ball and pull back to parameters via least-squares approximation. Alternatively, fix U and solve a convex subproblem in $(V,W,K^u)$.

Stability transfer under small parameter drift follows from Bauer–Fike (Lemma 2).

Appendix A.3. Controller-Update Safeguard (One-Step Kleinman)

Assumption 8

(Kleinman Update Conditions). For the one-step Kleinman/Riccati update to be well-posed:

1. 

The lifted pair $(K_r,B_r)$ is stabilizable: $\mathrm{rank}\left([B_r,K_r{B}_r,\dots ,K_r^{d-1}{B}_r]\right)=d$.

2. 

The pair $(Q^{1/2},K_r)$ is detectable: $\mathrm{rank}\left([Q^{1/2};K_r{Q}^{1/2};\dots ;K_r^{d-1}{Q}^{1/2}]\right)=d$.

3. 

A stabilizing seed gain $K_0^u$ exists such that $\rho (K_r-B_r{K}_0^u)<1$.

At each t, perform a single Kleinman/Riccati update warm-started from $(P_{t-1},K_{t-1}^u)$ to obtain a tentative $K_{t+{\displaystyle \frac{1}{2}}}^u$. If the contraction test fails for $A_{\mathrm{cl}}^{\mathrm{tmp}}=U{V}^{\top }+U{W}^{\top }{K}_{t+{\displaystyle \frac{1}{2}}}^u$, revert to $K_{t-1}^u$ and reduce the step (trust region or backtracking).

Appendix A.4. Streaming EDMD Without Dense SVD

Maintain low-rank factors $K_t\approx U_t{V}_t^{\top }$, $B_t\approx U_t{W}_t^{\top }$ and update with Oja/PAST-style rules using only $z_t,u_t,z_{t+1}$. To control numerical drift, apply TSQR to $U,V$ every k steps (e.g., $k=5$) or when $\parallel U^{\top }{U-I\parallel }_F>{\tau }_{\mathrm{orth}}$. This yields per-step costs dominated by matrix-vector multiplies $O\left(rd\right)$ and TSQR $O\left(d{r}^2\right)$.

Appendix A.5. Online Normalization of the Dictionary

Maintain running mean and variance per feature: ${\mu }_t=(1-\beta ){\mu }_{t-1}+\beta \varphi \left(x_t\right)$, $s_t^2=(1-\beta )s_{t-1}^2+\beta {(\varphi \left(x_t\right)-{\mu }_t)}^2$, and use normalized features $\widehat{\varphi }\left(x_t\right)=(\varphi \left(x_t\right)-{\mu }_t)/\sqrt{s_t^2+ϵ}$. Define $\mathrm{\Phi }$ using these normalized coordinates to keep conditioning stable in streaming.

Appendix A.6. Rank Adaptation (Streaming)

Let $e_t=z_{t+1}-U_t{V}_t^{\top }{z}_t-U_t{W}_t^{\top }{u}_t$, update ${\epsilon }_t=\lambda {\epsilon }_{t-1}+(1-\lambda ){\parallel e_t\parallel }^2$, and set $r_{t+1}=\sqrt{{\epsilon }_t}/(1-\gamma )$. Grow rank if ${\epsilon }_t$ plateaus above a threshold for K steps; shrink if the tracked subspace’s trailing energy (from the Oja/PAST internal state) falls below a threshold for K steps. All operations remain streaming; no batch SVD is used in the loop.

Appendix A.7. Complexity and Latency Breakdown

Per time step, the dominant costs are:

• Lift $x_t\mapsto z_t=\mathrm{\Phi }\left(x_t\right)$: $O\left(d\right)$ (with normalization).
• Updates $U{V}^{\top },U{W}^{\top }$: $O\left(rd\right)$ multiplies.
• TSQR orthonormalization (every k steps): $O\left(d{r}^2\right)$ amortized.
• One Kleinman/Riccati iteration: small dense solve reusing factorizations (dimension d).
• Safety QP/projection: a small QP with decision variables of order $r+m$; worst-case $O\left({(r+m)}^3\right)$, warm-started.

Table A1. Per-step operations (amortized) and typical complexity.

Component	Ops/Step (amortized)	Notes
Feature lift $\mathrm{\Phi }$	$O\left(d\right)$	EMA normalization included
Oja/PAST updates	$O\left(rd\right)$	Matrix-vector multiplies
TSQR (tall-skinny)	$O\left(d{r}^2\right)$	Every k steps; amortized
Kleinman update	$O\left(d^2\right)$	Warm-started from previous
Safety QP/projection	$O\left({(r+m)}^3\right)$	Small; warm-started

Table A1. Per-step operations (amortized) and typical complexity.

Component	Ops/Step (amortized)	Notes
Feature lift $\mathrm{\Phi }$	$O\left(d\right)$	EMA normalization included
Oja/PAST updates	$O\left(rd\right)$	Matrix-vector multiplies
TSQR (tall-skinny)	$O\left(d{r}^2\right)$	Every k steps; amortized
Kleinman update	$O\left(d^2\right)$	Warm-started from previous
Safety QP/projection	$O\left({(r+m)}^3\right)$	Small; warm-started

Appendix A.8. Experimental Reporting (Single-Loop)

For reproducibility, report mean±std over seeds, hardware (CPU/GPU, BLAS), lift cost inclusion, latency percentiles (for Oja/PAST, TSQR, Kleinman, QP), tube radius $r_t$ over time, spectral margin estimates, and safety violations (target: 0). Batch EDMD→MPC can be included in the appendix for comparison.