Koopman-Based Learning of Infinitesimal Generators
without Operator Logarithm

Given a pre-compact state space $\mathcal{X}\subseteq\mathbb{R}^{n}$, we consider a continuous-time nonlinear dynamical system of the form

where $x$ denotes the initial condition, and the vector field $f:\mathcal{X}\rightarrow\mathcal{X}$ is assumed to be locally Lipschitz continuous. We denote by $\phi:\mathcal{I}\times\mathcal{X}\rightarrow\mathcal{X}$ the forward flow, also known as the solution map.

The evolution of observable functions of system (1) restricted to $C(\mathcal{X})$ is governed by the family of Koopman operators, defined by

for each $t\geq 0$, where $\circ$ is the composition operator. By the properties of the flow map, it can be easily verified that $\left\{\mathcal{K}_{t}\right\}_{t\geq 0}$ forms a semigroup, with its (infinitesimal) generator defined by

Within the space of continuous observable functions, the limit in Eq. (3) exists on the subspace of all continuously differentiable functions, indicating that $\operatorname{dom}(\mathcal{L})=C^{1}(\mathcal{X})$. In this case, we can verify that the generator of Koopman operators is such that

Note that the finite-difference method [22, 23, 24] indirectly approximates the generator by setting a small terminal time without taking the limit on the r.h.s. of (3). Through this approximation scheme of the time derivative, it can be anticipated that the precision heavily depends on the size of that terminal time.

If $\mathcal{L}$ is a bounded linear operator that generates $\{\mathcal{K}_{t}\}$, then $\mathcal{K}_{t}=e^{t\mathcal{L}}$ for each $t$ in the uniform (operator) topology. However, many of the differential operators that occur in mathematical physics are unbounded. We revisit some concepts to show how to interpret $\mathcal{K}_{t}$ in terms of the possibly unbounded generator on an exponential scale.

We further define the Yosida approximation of $\mathcal{L}$ as

Note that $\left\{\mathcal{L}_{\lambda}\right\}_{\lambda\in\rho(\mathcal{L})}$ is a family of bounded linear operators, and $e^{t\mathcal{L}_{\lambda}}$ is well-defined for each $\lambda\in\rho(\mathcal{L})$.

Let $C(\mathcal{X})$ be endowed with the uniform norm $|\cdot|_{\infty}:=\sup_{x\in\mathcal{X}}|\cdot|$. Suppose that $\{\mathcal{K}_{t}\}_{t\geq 0}$ is a contraction semigroup where

for all $t\geq 0$. Then, $\rho(\mathcal{L})\supseteq(0,\infty)$, and by [25, Lemma 1.3.3], we have that $\left\{\mathcal{L}_{\lambda}\right\}_{\lambda>0}$ converges strongly to $\mathcal{L}$. Specifically, when working on the $C^{2}(\mathcal{X})$ subspace, the convergence rate is given by $\mathcal{O}(\lambda^{-1})$.

To address more general cases, we first present the following facts.

Intuitively, $M$ represents the uniform scaling of the magnitude of the Koopman operator, while $\omega$ indicates the dominant exponential growth or decay rate of the flow on $C(\mathcal{X})$. The above properties provide us with a tool to convert the dominant exponential growth or decay rate of the flow on $C(\mathcal{X})$ by introducing an extra term $\lambda\operatorname{I}$ to the original generator. Specifically, the transformation $e^{-\omega t}\mathcal{K}_{t}$ allows us to work in a topology (by shifting the direction of flows and uniformly compressing by $M$) that is equivalent to the uniform topology of continuous functions, where the semigroup $\{e^{-\omega t}\mathcal{K}_{t}\}$ is a contraction. The convergence of $\mathop{\mathrm{s\text{-}\lim}}_{\lambda\rightarrow\infty}\mathcal{L}_{\lambda}=\mathcal{L}$ for $\left\{\mathcal{L}_{\lambda}\right\}_{\lambda>\omega}$ and its reciprocal convergence rate can similarly be demonstrated in this new topology.

Motivated by representing $\mathcal{L}$ by $\left\{\mathcal{K}_{t}\right\}_{t\geq 0}$ and the Yosida approximation for $\mathcal{L}$ on $\{\lambda>\omega\}$, we establish a connection between $R(\lambda;\mathcal{L})$ and $\left\{\mathcal{K}_{t}\right\}_{t\geq 0}$.

Restricted to the $C^{1}(\mathcal{X})$ subspace, to use the approximation in Section III-A, we can replace $R(\lambda;\mathcal{L})$ with $R(\lambda)$.

The current form of the improper integral in (7) cannot be directly used for a data-driven approximation. We need to further derive an approximation approach based on finite-horizon observable data. Below, we present a direct truncation modification based on (7).

It can be shown that for any $\lambda$, the truncation in Definition III.4 results in an error term $e^{-\lambda\tau}R(\lambda)h(\phi(\tau,x))$ that is uniformly bounded for each observable function $h$, with the dominant term being $e^{-\lambda\tau}$. This indicates that for any arbitrarily large $\lambda$, this truncation does not significantly affect the accuracy compared to the procedure from $\mathcal{L}_{\lambda}$ to $\mathcal{L}$.

For any fixed $\tau>0$, we can therefore use

to approximate $\mathcal{L}$ within a small time-horizon. We illustrate this approximation with the following simple example.

For any fixed $\lambda>0$, given a dictionary $\mathcal{Z}_{N}$ of the form (10), for each $\mathfrak{z}_{i}\in\mathcal{Z}_{N}$ and each $x\in\mathcal{X}$, we consider $\mathfrak{z}_{i}(x)$ as the features and $\lambda^{2}\mathcal{T}_{\tau}\mathfrak{z}_{i}(x)-\lambda\mathfrak{z}_{i}(x)=\lambda^{2}\int_{0}^{\tau}e^{-\lambda s}\mathfrak{z}_{i}(\phi(s,x))ds-\lambda\mathfrak{z}_{i}(x)$ as the labels. To compute the integral, we employ numerical quadrature techniques for approximation. This approach inevitably requires discrete-time observations (snapshots), the number of which is denoted by $\mathbf{N}$, within the interval $[0,\tau]$ of the flow map $\phi(\tau,x)$.

To streamline the evaluation process for numerical examples, drawing inspiration from [27], for any $\tau$ and $i$, we can assess both the trajectory and the integral, i.e. the pair $\left(\phi(\tau,x),\lambda^{2}\int_{0}^{\tau}e^{-\lambda s}\mathfrak{z}_{i}(\phi(s,x))ds\right)$, by numerically solving the following augmented ODE system

We summarize the algorithm for generating training data for one time period as in Algorithm 1.

After obtaining the training data $(\mathfrak{X},\mathfrak{Y})$ using Algorithm 1, we can find $\mathbf{L}$ by $\mathbf{L}=\operatorname{argmin}_{A\in\mathbb{C}^{N\times N}}\|\mathfrak{Y}-\mathfrak{X}A\|_{F}$. The $\mathbf{L}$ is given in closed-form as $\mathbf{L}=\left(\mathfrak{X}^{T}\mathfrak{X}\right)^{\dagger}\mathfrak{X}^{T}\mathfrak{Y}$, where $\dagger$ is the pseudo inverse. Similar to Extended Dynamic Mode Decomposition (EDMD) [26] for learning Koopman operators, the approximations (11) and (12) can be guaranteed. In addition, by the universal approximation theorem, all of the function approximations from above should have uniform convergence to the actual quantities.

It is worth noting that operator learning frameworks, such as autoencoders, which utilize neural networks (NN) as dictionary functions, can reduce human bias in the selection of these functions [5]. Incorporating the proposed learning scheme with NN falls outside the scope of this paper but will be pursued in future work.

The actual vector field is

As we can analytically establish that $\mathcal{L}\mathfrak{z}_{1,0}(x)=f_{1}(x)$ and $\mathcal{L}\mathfrak{z}_{0,1}(x)=f_{2}(x)$, we use the approximation $[\widetilde{\mathcal{L}}_{\lambda}\mathfrak{z}_{1,0},\widetilde{\mathcal{L}}_{\lambda}\mathfrak{z}_{0,1}]$ to conversely obtain $f$.

Note that $\mathfrak{z}_{1,0}(x)=\mathcal{Z}_{N}(x)\mathbf{e}_{2}$ and $\mathfrak{z}_{0,1}(x)=\mathcal{Z}_{N}(x)\mathbf{e}_{5}$, where each $\mathbf{e}_{i}$ for $i\in\left\{0,1,\cdots,N\right\}$ is a column unit vector with all components being $0$ except for the $i$-th component, which is $1$. To apply $\mathbf{L}$, we have that $\widetilde{\mathcal{L}}_{\lambda}\mathfrak{z}_{1,0}(x)\approx\mathcal{Z}_{N}(x)(\mathbf{L}\mathbf{e}_{2})$ and $\widetilde{\mathcal{L}}_{\lambda}\mathfrak{z}_{0,1}(x)\approx\mathcal{Z}_{N}(x)(\mathbf{L}\mathbf{e}_{5})$.

In other words, we approximate $\widetilde{\mathcal{L}}_{\lambda}\mathfrak{z}_{1,0}$ and $\widetilde{\mathcal{L}}_{\lambda}\mathfrak{z}_{0,1}$ (and hence $f_{1}$ and $f_{2}$) using a linear combination of functions within $\mathcal{Z}_{N}$. The weights for these approximations are given by $\mathbf{L}\mathbf{e}_{2}$ and $\mathbf{L}\mathbf{e}_{5}$, respectively. We report the corresponding results in Table I and II.

We compare the aforementioned results with those obtained using the benchmark approach as described in [16] with the same $M$ and $\mathcal{Z}_{N}$. As described in Remark IV.1, $\mathbf{K}\in\mathbb{C}^{N\times N}$ can be obtained such that $\mathcal{K}_{s}\mathfrak{z}_{1,0}(x)\approx\mathcal{Z}_{N}(x)(\mathbf{K}\mathbf{e}_{2})$ and $\mathcal{K}_{s}\mathfrak{z}_{0,1}(x)\approx\mathcal{Z}_{N}(x)(\mathbf{K}\mathbf{e}_{5})$. Then, we have the approximation, as claimed by [16], $f_{1}(x)\approx\mathcal{Z}_{N}(x)(\log(\mathbf{K})\mathbf{e}_{2}/s)$ and $f_{2}(x)\approx\mathcal{Z}_{N}(x)(\log(\mathbf{K})\mathbf{e}_{5}/s)$ w.r.t. $|\cdot|_{\infty}$. According to [16, Section VI.A], we set $s=0.5$ to avoid the multi-valued matrix logarithm. We report the weights obtained by taking the real parts of $\log(\mathbf{K})\mathbf{e}_{2}/s$ and $\log(\mathbf{K})\mathbf{e}_{5}/s$ in Table III and IV, respectively. Multiple orders of magnitude in accuracy have been established using the proposed method.

It is worth noting that, unlike the experiment in [16] where $i,j\in\left\{0,1,2,3\right\}$, we deliberately choose different numbers for $i$ and $j$, leading to non-negligible imaginary parts after taking the matrix logarithm of $\mathbf{K}$. The imaginary parts of the learned weights using the Koopman-logarithm approach are reported in Table V and VI. The presence of non-negligible imaginary parts indicates that the underlying system is sensitive to the selection of dictionary functions when employing the Koopman-logarithm approach, an effect that is significant and cannot be overlooked when aiming to minimize human intervention in identifying unknown systems in practice. Furthermore, the original experiment in [16] used $M=300^{2}$ samples for a data fitting, and when $M$ is reduced to $100^{2}$, the quality of the operator learning deteriorates.

Observing the identified system dynamics, we anticipate that the time-reversed system is (locally) asymptotically stable w.r.t. the origin. We use the learned $\mathbf{L}$ to construct polynomial Lyapunov functions based on the Lyapunov equation $\mathcal{L}V(x)=|x|^{2}$, where the sign on the r.h.s. is reversed due to the dynamics being reversed. To use $\mathbf{L}$ and the library functions from $\mathcal{Z}_{N}$, we define the weights as $\theta:=[\theta_{0},\theta_{1},\cdots,\theta_{N-1}]^{T}$ and seek a $\theta$ such that $V(x;\theta):=\mathcal{Z}_{N}(x)\theta$ , with the objective of minimizing $|\mathcal{Z}_{N}(x)\mathbf{L}\theta-x_{1}^{2}-x_{2}^{2}|$. Ignoring the small terms of the magnitude $10^{-9}$, the constructed Lyapunov function is given by

Additionally, $\mathcal{L}V(x;\theta)$ is approximated by $\mathcal{Z}_{N}(x)\mathbf{L}\theta$, with the maximal value observed at $\mathcal{Z}_{N}(0)\mathbf{L}\theta=-0.068$, which indicates that stability predition using the data-driven Lyapunov function is verified to be valid. This indicates that the stability prediction, as inferred using the data-driven Lyapunov function, is verified to be valid. The visualization can be found in Fig 1.