Deep Koopman-layered Model with Universal Property Based on Toeplitz Matrices

## Abstract

Abstract We propose deep Koopman-layered models with learnable parameters in the form of Toeplitz matrices for analyzing the dynamics of time-series data.
The proposed model has both theoretical solidness and flexibility.
By virtue of the universal property of Toeplitz
matrices and the reproducing property underlined in the model, we can show its universality and the generalization property.
In addition, the flexibility of the proposed model enables the model to fit time-series data coming from nonautonomous dynamical systems.
When training the model, we apply Krylov subspace methods for efficient computations.
In addition, the proposed model can be regarded as a neural ODE-based model.
In this sense, the proposed model establishes a new connection among Koopman operators, neural ODEs, and numerical linear algebraic methods.

## 1 Introduction

Koopman operator has been one of the important tools in machine learning (Kawahara, ; Ishikawa et al., ; Lusch et al., ; Brunton & Kutz, ; Hashimoto et al., ).
Koopman operators are linear operators that describe the composition of functions and are applied to analyzing time-series data generated by nonlinear dynamical systems (Koopman, ; Budišić et al., ; Klus et al., ; Giannakis & Das, ; Mezić, ).
By computing the eigenvalues of Koopman operators, we can understand the long-term behavior of the undelined dynamical systems.
An important feature of Applying Koopman operators is that we can estimate them with given time-series data through fundamental linear algebraic tools such as projection.
A typical approach to estimate Koopman operators is extended dynamical mode decomposition (EDMD) (Williams et al., ).
For EDMD, we need to choose the dictionary functions to determine the representation space of the Koopman operator, and what choice of them gives us a better estimation is far from trivial.
In addition, since we construct the estimation in an analytical way, the model is not flexible enough to incorporate additional information about dynamical systems.
With EDMD as a starting point, many DMD-based methods are proposed (Kawahara, ; Colbrook & Townsend, ; Schmid, ).
For autonomous systems, we need to estimate a single Koopman operator.
In this case, Ishikawa et al. () proposed to choose derivatives of kernel functions as dictionary functions based on the theory of Jet spaces.
Several works deal with nonautonomous systems.
Maćešić et al. () applied EDMD to estimate a time-dependent Koopman operator for each time window.
However, as far as we know, no existing works show proper choices of dictionary functions for nonautonomous systems based on theoretical analysis.
In addition, since each Koopman operator for a time window is estimated individually in these frameworks, we cannot take the information of other Koopman operators into account.

To find a proper representation space and gain the flexibility of the model, neural network-based Koopman methods have been proposed (Lusch et al., ; Azencot et al., ; Shi & Meng, ).
These methods set the encoder from the data space to the representation space where the Koopman operator is defined, and the decoder from the representation space to the data space, as deep neural networks.
Then, we train them.
Neural network-based Koopman methods for nonautonomous systems have also been proposed.
Liu et al. () proposed to decompose the Koopman operator into a time-invariant part and a time-variant part.
The time-variant part of the Koopman operator is constructed individually for each time window using EDMD.
Xiong et al. () assumed the ergodicity of the dynamical system and considered time-averaged Koopman operators for nonautonomous dynamical systems.
However, their theoretical properties have not been fully understood, and since the representation space changes as the learning process proceeds, their theoretical analysis is challenging.

In this work, we propose a framework that estimates multiple Koopman operators over time with the Fourier basis representation space and learnable Toeplitz matrices.
Using our framework, we can estimate multiple Koopman operators simultaneously and can capture the transition of properties of data along time via multiple Koopman operators.
We call each Koopman operator the Koopman-layer, and the whole model the deep Koopman-layered model.
The proposed model has both theoretical solidness and flexibility.
We show that the Fourier basis is a proper basis for constructing the representation space even for nonautonomous dynamical systems in the sense that we can show its theoretical properties such as universality and generalization bound.
In addition, the proposed model has learnable parameters, which makes the model more flexible to fit nonautonomous dynamical systems than the analytical methods such as EDMD.
The proposed model resolves the issue of theoretical analysis for the neutral network-based methods and that of the flexibility for the analytical methods simultaneously.

We show that each Koopman operator is represented by the exponential of a matrix constructed with Toeplitz matrices and diagonal matrices.
This allows us to apply Krylov subspace methods (Gallopoulos & Saad, ; Güttel, ; Hashimoto & Nodera, ) to compute the estimation of Koopman operators with low computational costs.
By virtue of the universal property of Toeplitz matrices (Ye & Lim, ), we can show the universality of the proposed model with a linear algebraic approach.
We also show a generalization bound of the proposed model using a reproducing kernel Hilbert space (RKHS) associated with the Fourier functions.
We can analyze both the universality and generalization error with the same framework.

The proposed model can also be regarded as a neural ODE-based model (Chen et al., ; Teshima et al., ; Li et al., ).
While in the existing method, we train the models with numerical analysis approaches, in the proposed method, we train the models with a numerical linear algebraic approach.
The universality and generalization results of the proposed model can also be seen as those for the neural ODE-based models.
Our method sheds light on a new linear algebraic approach to the design of neural ODEs.

Our contributions are summarized as follows:

- •
We propose a model for analyzing nonautonomous dynamical systems that has both theoretical solidness and flexibility.
We show that the Fourier basis provides us with a proper representation space, in the sense that we can show the universality and the generalization bound regarding the model.
As for the flexibility, we can learn multiple Koopman operators simultaneously, which enables us to extract the transition of properties of dynamical systems along time.
- •
We apply Krylov subspace methods to compute the estimation of Koopman operators.
This establishes a new connection between Koopman operator theoretic approaches and Krylov subspace methods, which opens up future directions for extracting further information about dynamical systems using numerical linear algebraic approaches.
- •
We provide a new implementation method for neural ODEs purely with numerical linear algebraic approaches, not with numerical analysis approaches.

## 2 Preliminary

In this section, we provide notations and fundamental notions required in this paper.

### 2.1 Notations

In this paper, we use a generalized concept of matrices.
For a finite index set $N\subset\mathbb{Z}^{d}$ and $a_{j,l}\in\mathbb{C}$ $(j,l\in N)$, we call $A=[a_{j,l}]_{j,l\in N}$ an $N$ by $N$ matrix and denote by $\mathbb{C}^{N\times N}$ the space of all $N$ by $N$ matrices.
Indeed, by constructing a bijection $I:N\to\{1,\ldots,|N|\}$ and setting $\tilde{a}_{I(j),I(l)}=a_{j,l}$, we obtain a standard matrix $[\tilde{a}_{I(j),I(l)}]_{I(j),I(l)}$ corresponding to $A$.

### 2.2 L 2 superscript 𝐿 2 L^{2} italic_L start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT space and Reproducing kernel Hilbert space on the torus

We consider two function spaces, the $L^{2}$ space and RKHS, in this paper.
Let $\mathbb{T}$ be the torus $\mathbb{R}/2\pi\mathbb{Z}$.
We denote by $L^{2}(\mathbb{T}^{d})$ the space of square-integrable complex-valued functions on $\mathbb{T}^{d}$, equipped with the Lebesgue measure.
As for the RKHS, let $\kappa:\mathbb{T}^{d}\times\mathbb{T}^{d}\to\mathbb{C}$ be a positive definite kernel, which satisfies the following two properties:

- 1.
${\kappa}(x,y)=\overline{\kappa(y,x)}$  for $x,y\in\mathbb{T}^{d}$,
- 2.
$\sum_{n,m=1}^{N}\!\overline{c_{n}}c_{m}\kappa(x_{n},x_{m})\!\geq\!0$  for $N\in\mathbb{N}$, $c_{n}\in\mathbb{C}$, $x_{n}\in\mathbb{T}$.

Let ${\phi}$ be a map defined as $\phi(x)=\kappa(\cdot,x)$, which is called the feature map.
The RKHS $\mathcal{H}_{\kappa}$ is the Hilbert space spanned by $\{\phi(x)\,\mid\,x\in\mathbb{T}^{d}\}$.
The inner product $\left\langle\cdot,\cdot\right\rangle:\mathcal{H}_{\kappa}\times\mathcal{H}_{
\kappa}\to\mathbb{C}$ in $\mathcal{H}_{\kappa}$ is defined as

$$ $\bigg{\langle}\sum_{n=1}^{N}c_{n}{\phi}(x_{n}),\sum_{m=1}^{M}d_{m}{\phi}(y_{m} )\bigg{\rangle}=\sum_{n=1}^{N}\sum_{m=1}^{M}\overline{c_{n}}d_{m}\kappa(x_{n}, y_{m})$ $$

for $c_{n},d_{n}\in\mathbb{C}$ and $x_{n},y_{n}\in\mathbb{T}^{d}$.
Note that by the definition of $\kappa$, $\left\langle\cdot,\cdot\right\rangle$ is well-defined and satisfies the axiom of inner products.
An important property for RKHSs is the reproducing property.
For $x\in\mathbb{T}^{d}$ and $v\in\mathcal{H}_{\kappa}$, we have $\left\langle\phi(x),v\right\rangle=v(x)$, which is useful for deriving a generalization bound.

### 2.3 Koopman generator and operator

Consider an ODE $\frac{\mathrm{d}x}{\mathrm{d}t}(t)=f(x(t))$ on $\mathbb{T}^{d}$.
Let $g:\mathbb{R}\times\mathbb{T}^{d}$ be the flow of the ODE, that is, $g$ satisfies $g(0,x)=x$ and $g(s,g(t,x))=g(s+t,x)$ for $x\in\mathbb{T}^{d}$.
We assume $g$ is continuous and invertible.
For $t\in\mathbb{R}$, we define the Koopman operator $K^{t}$ on $L^{2}(\mathbb{T}^{d})$ by the composition with $g(t,\cdot)$ as $K^{t}h(x)=h(g(t,x))$ for $h\in L^{2}(\mathbb{T}^{d})$ and $x\in\mathbb{T}^{d}$.
The Koopman operator is a linear operator that maps a function $h$ to a function $h(g(t,\cdot))$.
Note that the Koopman operator $K^{t}$ is linear even if $g(t,\cdot)$ is nonlinear.
Since $K^{t}$ depends on $t$, we can consider the family of Koopman operators $\{K^{t}\}_{t\in\mathbb{R}}$.
For $h\in C^{1}(\mathbb{T}^{d})$, where $C^{1}(\mathbb{T}^{d})$ is the space of continuous differentiable functions on $\mathbb{T}^{d}$, let $L$ be the linear operator defined as

$$ $\displaystyle Lh=\lim_{t\to\infty}\frac{K^{t}h-h}{t},$ $$

where the limit is by means of $L^{2}(\mathbb{T})$.
We call $L$ the Koopman generator.
We write $K^{t}=\mathrm{e}^{tL}$.
If $L$ is bounded, then it coincides with the standard definition $\mathrm{e}^{tL}=\sum_{i=1}^{\infty}(tL)^{i}/i!$.
If $L$ is unbounded, it can be justified by approximating $L$ by a sequence of bounded operators and considering the strong limit of the sequence of the exponential of the bounded operators (Yosida, ).

## 3 Deep Koopman-layered model

We propose deep Koopman-layered models based on the Koopman operator theory, which have both theoretical solidness and flexibility.

### 3.1 Multiple dynamical systems and Koopman generators

Consider $J$ ODEs $\frac{\mathrm{d}x}{\mathrm{d}t}(t)=f_{j}(x(t))$ on $\mathbb{T}^{d}$ for $j=1,\ldots,J$.
Let $g_{j}:\mathbb{R}\times\mathbb{T}^{d}$ be the flow of the $j$th ODE. For $v\in L^{2}(\mathbb{T}^{d})$, consider the following model:

$$ $\displaystyle G(x)=v\circ g_{J}(t_{J},\cdot)\circ\cdots\circ g_{1}(t_{1},\cdot )(x)=v(g_{J}(t_{J},\cdots g_{1}(t_{1},x))).$ (1) $$

This model describes a switching dynamical system, and also is regarded as a discrete approximation of a nonautonomous dynamical system.

###### Remark 1

Since we are focusing on the complex-valued function space $L^{2}(\mathbb{T})$, $G$ itself is a complex-valued function.
However, we can easily extend the model to the flow $g_{J}(t_{J},\cdot)\circ\cdots\circ g_{1}(t_{1},\cdot)$, which is a map from $\mathbb{T}^{d}$ to $\mathbb{T}^{d}$.
We can obtain a complex-valued function on $\mathbb{T}^{d+1}$ that describes a map from $\mathbb{T}^{d}$ to $\mathbb{T}^{d}$.
Indeed, let $\tilde{g}_{j}(x,y)=[g_{j}(t_{j},x),y]$ for $x\in\mathbb{T}^{d}$ and $y\in\mathbb{T}$.
Let $\tilde{v}$ be a function that satisfies $\tilde{v}(x,k/d)=x_{k}$, and let $G=\tilde{v}\circ\tilde{g}_{J}\circ\cdots\circ\tilde{g}_{1}$.
Then, $G(\cdot,k/d)$ is the $k$th element of $g_{J}(t_{J},\cdot)\circ\cdots\circ g_{1}(t_{1},\cdot)$.

###### Remark 1

### 3.2 Approximation of Koopman generators using Toeplitz matrices

We consider training the model (1) using given time-series data.
For this purpose, we apply the Koopman operator theory.
Let $L_{j}$ be the Koopman generator associated with the flow $g_{j}$.
Since the Koopman operator $K^{t_{j}}_{j}$ of $g_{j}$ is represented as $\mathrm{e}^{t_{j}L_{j}}$, the model (1) is represented as

$$ $\displaystyle G=\mathrm{e}^{t_{1}L_{1}}\cdots\mathrm{e}^{t_{J}L_{J}}v.$ $$

To deal with the Koopman generators defined on the infinite-dimensional space, we approximate them using Fourier functions.
Let $q_{n}(x)=\mathrm{e}^{\mathrm{i}n\cdot x}$ for $n\in\mathbb{Z}^{d}$ and $x\in\mathbb{T}^{d}$, where $\mathrm{i}$ is the imaginary unit.
Let $M^{j}_{r}\subset\mathbb{Z}^{d}$ be a finite index set for $r=1,\ldots,R_{j}$.
We set the $k$th element of the function $f_{j}$ in the ODE as

$$ $\displaystyle\sum_{m_{R_{j}}\in M^{j}_{R_{j}}}a^{j,k}_{m_{R_{j}},R_{j}}q_{m_{R _{j}}}\cdots\sum_{m_{1}\in M^{j}_{1}}a^{j,k}_{m_{1},1}q_{m_{1}}$ (2) $$

with $a^{j,k}_{m_{r},r}\in\mathbb{C}$,
the product of weighted sums of Fourier functions.
Then, we approximate the Koopman generator $L_{j}$ by projecting the input vector onto the space $V_{N}:=\operatorname{Span}\{q_{n}\,\mid\,n\in N\}$, where $N\subset\mathbb{Z}^{d}$ is a finite index set, applying $L_{j}$, and projecting it back to $V_{N}$ as $Q_{N}Q_{N}^{*}L_{j}Q_{N}Q_{N}^{*}$.
Here, $Q_{N}:\mathbb{C}^{N}\to V_{N}$ is the linear operator defined as $Q_{N}c=\sum_{n\in N}c_{n}q_{n}$ for $c=(c_{n})_{n\in N}\in\mathbb{C}^{N}$ and ^∗ is the adjoint.
Note that $Q_{N}Q_{N}^{*}$ is the projection onto $V_{N}$.
Then, the representation matrix $Q_{N}^{*}L_{j}Q_{N}$ of the approximated Koopman generator $Q_{N}Q_{N}^{*}L_{j}Q_{N}Q_{N}^{*}$ is written as follows.
Throughout the paper, all the proofs are documented in Appendix A.

###### Proposition 2

The $(n,l)$-entry of the representation matrix $Q_{N}^{*}L_{j}Q_{N}$ of the approximated operator is

$$ $\displaystyle\sum_{k=1}^{d}\sum_{n_{R_{j}}-l\in M^{j}_{R_{j}}}\sum_{n_{R_{j}-1 }-n_{R_{j}}\in M^{j}_{R_{j}-1}}\cdots\sum_{n_{2}-n_{3}\in M^{j}_{2}}\sum_{n-n_ {2}\in M^{j}_{1}}$ $\displaystyle\qquad\qquad\qquad\qquad a^{j,k}_{n_{R_{j}}-l,R_{j}}a^{j,k}_{n_{R _{j}-1}-n_{R_{j}},R_{j-1}}\cdots a^{j,k}_{n_{2}-n_{3},2}a^{j,k}_{n-n_{2},1} \mathrm{i}l_{k},$ (3) $$

where $l_{k}$ is the $k$th element of the index $l\in\mathbb{Z}^{d}$.
Moreover, we set $n_{r}=m_{R_{j}}+\cdots+m_{r}+l$, thus $n_{1}=n$, $m_{r}=n_{r}-n_{r+1}$ for $r=1,\ldots,R_{j}-1$, and $m_{R_{j}}=n_{R_{j}}-l$.

Note that since the sum involves the differences of indices, it can be written using Toeplitz matrices, whose $(n,l)$-entry depends only on $n-l$.
We approximate the sum appearing in Eq. (3) by restricting the index $n_{r}$ to $N$, combine with the information of time $t_{j}$, and set a matrix $\mathbf{L}_{j}\in\mathbb{C}^{N\times N}$ as

$$ $\displaystyle\mathbf{L}_{j}=t_{j}\sum_{k=1}^{d}A_{1}^{j,k}\cdots A_{R_{j}}^{j, k}D_{k},$ $$

where $A_{r}^{j,k}$ is the Toeplitz matrix defined as $A_{r}^{j,k}=[a_{n-l,r}^{j,k}]_{n,l\in N}$ and $D_{k}$ is the diagonal matrix defined as $(D_{k})_{l,l}=\mathrm{i}l_{k}$.
We finally regard $Q_{N}\mathbf{L}_{j}Q_{N}^{*}$ as an approxitation of the Koopman genertor $L_{j}$.

Then, we construct the approximation $\mathbf{G}$ of $G$, defined in Eq. (1), as

$$ $\displaystyle\mathbf{G}=\mathrm{e}^{Q_{N}\mathbf{L}_{1}Q_{N}^{*}}\cdots\mathrm {e}^{Q_{N}\mathbf{L}_{J}Q_{N}^{*}}v=Q_{N}\mathrm{e}^{\mathbf{L}_{1}}\cdots \mathrm{e}^{\mathbf{L}_{J}}Q_{N}^{*}v.$ (4) $$

We call the model $\mathbf{G}$ deep Koopman-layered model.

To compute the product of the matrix exponential $\mathrm{e}^{\mathbf{L}_{j}}$ and the vector $\mathrm{e}^{\mathbf{L}_{j+1}}\cdots\mathrm{e}^{\mathbf{L}_{J}}Q_{N}^{*}v$, we can use Krylov subspace methods.
If the number of indices for describing $f_{j}$ is smaller than that for describing the whole model, i.e., $|M^{j}_{r}|\ll|N|$, then the Toeplitz matrix $A^{j,k}_{r}$ is sparse.
In this case, the matrix-vector product can be computed with the computational cost of $O(\sum_{r=1}^{R_{j}}|M^{j}_{r}||N|)$.
Thus, one iteration of the Krylov subspace method costs $O(\sum_{r=1}^{R_{j}}|M^{j}_{r}||N|)$, which makes the computation efficient compared to direct methods without taking the structure of the matrix into account, whose computational cost results in $O(|N|^{3})$.

###### Remark 3

To restrict $f_{j}$ to be a real-valued map and reduce the number of parameters $a_{m,r}^{j,k}$, we can set $M^{j}_{r}$ as $\{-m^{j}_{1,r},\ldots,m^{j}_{1,r}\}\times\cdots\times\{-m^{j}_{d,r},\ldots,m^{
j}_{d,r}\}$ for $m^{j}_{k,r}\in\mathbb{N}$ for $k=1,\ldots,d$.
In addition, we set $a_{m,r}^{j,k}=\overline{a_{-m,r}^{j,k}}$ for $m\in M_{r}^{j}$.
Then, we have $a_{m,r}^{j,k}q_{m}=\overline{a_{-m,r}^{j,k}q_{-m}}$, and $f_{j}$ is real-valued.

###### Remark 4

An advantage of applying Koopman operators is that their spectra describe the properties of dynamical systems.
For example, if the dynamical system is measure preserving, then the corresponding Koopman operator is unitary.
Since each Koopman layer is an estimation of the Koopman operator, we can analyze time-series data coming from nonsutonomous dynamical systems by computing the eigenvalues of the Koopman layers.
We will observe the eigenvalues of Koopman layers numerically in Subsection 6.2.

###### Proposition 2

###### Remark 3

###### Remark 4

## 4 Universality

In this section, we show the universal property of the proposed deep Koopman-layered model.
We can interpret the model $\mathbf{G}$ as the approximation of the target function by transforming the function $v$ into the target function using the linear operator $Q_{N}\mathrm{e}^{\mathbf{L}_{1}}\cdots\mathrm{e}^{\mathbf{L}_{J}}Q_{N}^{*}$.
If we can represent any linear operator by $\mathrm{e}^{\mathbf{L}_{1}}\cdots\mathrm{e}^{\mathbf{L}_{J}}$, then we can transform $v$ into any target function in $V_{N}$, which means we can approximate any function as $N$ goes to the whole set $\mathbb{Z}^{d}$.
Thus, this property corresponds to the universality of the model.
In Section 3, by constructing the model with the matrix $\mathrm{e}^{\mathbf{L}_{1}}\cdots\mathrm{e}^{\mathbf{L}_{J}}$ based on the Koopman operators with the Fourier functions, we restrict the number of parameters of the linear operator that transforms $v$ into the target function.
The universality of the model means that this restriction is reasonable in the sense of representing the target functions using the deep Koopman-layered model.

Let $T(N,\mathbb{C})=\{\sum_{k=1}^{d}A_{1}^{k}\cdots A_{R_{k}}^{k}D_{k},\,\mid\,R_{
k}\in\mathbb{N},\ A_{1}^{k}\cdots A_{R_{k}}^{k}\in\mathbb{C}^{N\times N}\mbox{
: Toeplitz}\}$.
Let $L_{0}^{2}(\mathbb{T}^{d})=\overline{\operatorname{Span}\{q_{n}\,\mid\,n\neq 0\}}$ be the space of $L^{2}$ functions whose average is 0 0.
We show the following fundamental result of the universality of the model:

###### Theorem 5

Assume $v\in L_{0}^{2}(\mathbb{T}^{d})$ and $v\neq 0$.
For any $f\in L_{0}^{2}(\mathbb{T}^{d})$ with $f\neq 0$ and for any $\epsilon>0$, there exist a finite set $N\subset\mathbb{Z}\setminus\{0\}$, a positive integer $J$, and matrices $\mathbf{L}_{1},\ldots,\mathbf{L}_{J}\in T(N,\mathbb{C})$ such that $\|f-\mathbf{G}\|\leq\epsilon$ and $\mathbf{G}=Q_{N}\mathrm{e}^{\mathbf{L}_{1}}\cdots\mathrm{e}^{\mathbf{L}_{J}}Q_
{N}^{*}v$.

Theorem 5 is for a single function $f$, but applying Theorem 5 for each component of $\mathbf{G}$, we obtain the following result for the flow $g_{\tilde{J}}(t_{\tilde{J}},\cdot)\circ\cdots\circ g_{1}(t_{1},\cdot)$ with $\tilde{J}\in\mathbb{N}$, which is considered in Eq. (1).

###### Corollary 6

Assume $v\in L_{0}^{2}(\mathbb{T}^{d})$ and $v\neq 0$.
For any sequence $g_{1}(t_{1},\cdot),\ldots,g_{\tilde{J}}(t_{\tilde{J}},\cdot)$ of flows that satisfies $v\circ g_{\tilde{J}}(t_{\tilde{J}},\cdot)\circ\cdots\circ g_{j}(t_{j},\cdot)
\in L^{2}_{0}(\mathbb{T}^{d})$ and $v\circ g_{\tilde{J}}(t_{\tilde{J}},\cdot)\circ\cdots\circ g_{j}(t_{j},\cdot)\neq
0$ for $j=1,\ldots\tilde{J}$, and for any $\epsilon>0$, there exist a finite set $N\subset\mathbb{Z}\setminus\{0\}$, integers $0<J_{1}<\cdots<J_{\tilde{J}}$, and matrices $\mathbf{L}_{1},\ldots,\mathbf{L}_{J_{\tilde{J}}}\in T(N,\mathbb{C})$ such that $\|v\circ g_{\tilde{J}}(t_{\tilde{J}},\cdot)\circ\cdots\circ g_{j}(t_{j},\cdot)
-\mathbf{G}_{j}\|\leq\epsilon$ and $\mathbf{G}_{j}=Q_{N}\mathrm{e}^{\mathbf{L}_{J_{j-1}+1}}\cdots\mathrm{e}^{
\mathbf{L}_{J_{\tilde{J}}}}Q_{N}^{*}v$ for $j=1,\ldots,\tilde{J}$, where $J_{0}=1$.

###### Remark 7

The function space $L_{0}^{2}(\mathbb{T}^{d})$ for the target function is not restrictive.
By adding a constant to the functions in $L_{0}^{2}(\mathbb{T}^{d})$, we can represent any function in $L^{2}(\mathbb{T}^{d})$.
Thus, by adding one additional learnable parameter $c\in\mathbb{C}$ to the model $\mathbf{G}$ in Theorem 5 and consider the model $\mathbf{G}(x)+c$ for an input $x\in\mathbb{T}^{d}$, we can represent any function in $L^{2}(\mathbb{T}^{d})$.

The proof of Theorem 5 is obtained by a linear algebraic approach.
By virtue of setting $f_{j}$ as the product of weighted sums of Fourier functions as explained in Eq. (2), the approximation of the Koopman generator is composed of Toeplitz matrices.
As a result, we can apply the following proposition regarding Toeplitz matrices by Ye & Lim (, Theorem 2).

###### Proposition 8

For any $B\in\mathbb{C}^{N\times N}$, there exists $R=\lfloor|N|\rfloor+1$ Toeplitz matrices $A_{1},\ldots,A_{R}$ such that $B=A_{1}\cdots A_{R}$.

We use Proposition 8 to show the following lemma regarding the representation with $T(N,\mathbb{C})$.

###### Lemma 9

Assume $N\subset\mathbb{Z}^{d}\setminus\{0\}$.
Then, we have $\mathbb{C}^{N\times N}=T(N,\mathbb{C})$.

Since $\mathbb{C}^{N\times N}$ is a Lie algebra and the corresponding Lie group $GL(N,\mathbb{C})$, the group of nonsingular $N$ by $N$ matrices, is connected, we have the following lemma (Hall, , Corollary 3.47).

###### Lemma 10

We have $GL(N,\mathbb{C})=\{\mathrm{e}^{\mathbf{L}_{1}}\cdots\mathrm{e}^{\mathbf{L}_{J}
}\,\mid\,J\in\mathbb{N},\ \mathbf{L}_{1},\ldots,\mathbf{L}_{J}\in\mathbb{C}^{N
\times N}\}$.

We also use the following transitive property of $GL(N,\mathbb{C})$ and finally obtain Theorem 5.

###### Lemma 11

For any $\mathbf{u},\mathbf{v}\in\mathbb{C}^{N}\setminus\{0\}$, there exists $A\in GL(N,\mathbb{C})$ such that $\mathbf{u}=A\mathbf{v}$.

###### Theorem 5

###### Corollary 6

###### Remark 7

###### Proposition 8

###### Lemma 9

###### Lemma 10

###### Lemma 11

##### Training deep Koopman-layered model with time-series data

Based on Corollary 6, we can train the deep Koopman-layered model using time-series data as follows:
We first fix the final nonlinear transform $v$ in the model $\mathbf{G}$ taking Remark 1 into account, the number of layers $\tilde{J}$, and the index sets $N$, $M_{r}^{j}$.
We input a family of time-series data $\{x_{s,0},\ldots,x_{s,\tilde{J}}\}_{s=1}^{S}$ to $\mathbf{G}$.
We learn the parameter $a_{m,r}^{j,k}$ in $\mathbf{G}$ by minimizing $\sum_{s=1}^{S}\ell(v(x_{s,\tilde{J}}),\mathbf{G}_{j}(x_{s,j-1}))$ for $j=1,\ldots,\tilde{J}$ using an optimization method.
For example, we can set an objective function $\sum_{j=1}^{\tilde{J}}\sum_{s=1}^{S}\ell(v(x_{s,\tilde{J}}),\mathbf{G}_{j}(x_{
s,j-1}))$.
Here $\ell:\mathbb{C}\times\mathbb{C}\to\mathbb{R}$ is a loss function.
For example, we can set $\ell$ as the squared error.
We will show examples of the training of deep Koopman-layered models and numerical results in Section 6.

## 5 Generalization bound

We investigate the generalization property of the proposed deep Koopman-layered model in this section.
Our framework with Koopman operators enables us to derive a generalization bound involving the norms of Koopman operators.

Let $\mathcal{G}_{N}=\{Q_{N}\mathrm{e}^{\mathbf{L}_{1}}\cdots\mathrm{e}^{\mathbf{L}
_{J}}Q_{N}^{*}v\,\mid\,\mathbf{L}_{1},\ldots,\mathbf{L}_{J}\in T(N,\mathbb{C})\}$ be the function class of deep Koopman-layered model (4).
Let $\ell(\mathcal{G}_{N})=\{(x,y)\mapsto\ell(f(x),y)\,\mid\,f\in\mathcal{G}_{N}\}$ for a function $\ell$ that is bounded by $C>0$.
Then, we have the following result of a generalization bound for the deep Koopman-layered model.

###### Proposition 12

Let $h\in\ell(\mathcal{G}_{N})$, $x$ and $y$ be random variables, $S\in\mathbb{N}$, and $x_{1},\ldots,x_{S}$ and $y_{1},\ldots,y_{S}$ be i.i.d. samples drawn from the distributions of $x$ and $y$, respectively.
For any $\delta>0$, with probability at least $1-\delta$, we have

$$ $\displaystyle\mathrm{E}[h(x,y)]\leq\frac{1}{S}\sum_{s=1}^{S}h(x_{n},y_{n})+ \frac{\alpha}{\sqrt{S}}\max_{j\in N}\mathrm{e}^{\tau\|j\|_{1}}\!\!\!\!\!\!\!\! \!\sup_{\mathbf{L}_{1},\ldots,\mathbf{L}_{J}\in T(N,\mathbb{C})}\!\!\!\!\!\| \mathrm{e}^{\mathbf{L}_{1}}\|\cdots\|\mathrm{e}^{\mathbf{L}_{J}}\|\,\|v\|+3C \sqrt{\frac{\log(\delta/2)}{S}}.$ $$

We use the Rademacher complexity to derive Proposition 12.
For this purpose, we regard the model (1) as a function in an RKHS.
For $j\in\mathbb{Z}^{d}$ and $x\in\mathbb{T}^{d}$, let $\tilde{q}_{j}(x)=\mathrm{e}^{-\tau\|j\|_{1}}\mathrm{e}^{\mathrm{i}j\cdot x}$, where $\tau>0$ is a fixed parameter and $\|[j_{1},\ldots,j_{d}]\|_{1}=|j_{1}|+\cdots+|j_{d}|$ for $[j_{1},\ldots,j_{d}]\in\mathbb{Z}^{d}$.
Let $\kappa(x,y)=\sum_{j\in\mathbb{Z}^{d}}\overline{\tilde{q}_{j}(x)}\tilde{q}_{j}(y)$, and consider the RKHS $\mathcal{H}_{\kappa}$ associated with the kernel $\kappa$.
Note that $\{\tilde{q}_{j}\,\mid\,j\in\mathbb{Z}^{d}\}$ is an orthonormal basis of $\mathcal{H}_{\kappa}$.
Giannakis et al. () and Das et al. () used this kind of RKHSs for simulating dynamical systems on a quantum computer based on the Koopman operator theory and for approximating Koopman operators by a sequence of compact operators.
Here, we use the RKHS $\mathcal{H}_{\kappa}$ for deriving a generalization bound.
To regard the function $\mathbf{G}\in V_{N}=\operatorname{Span}\{q_{j}\,\mid\,j\in N\}\subset L^{2}(
\mathbb{T}^{d})$ as a function in $\mathcal{H}_{\kappa}$, we define an inclusion map $\iota_{N}:V_{N}\to\mathcal{H}_{\kappa}$ as $\iota_{N}q_{j}=\mathrm{e}^{\tau\|j\|_{1}}\tilde{q}_{j}$ for $j\in N$.
Then, the operator norm of $\iota_{N}$ is $\|\tau_{N}\|=\max_{j\in N}\mathrm{e}^{\tau\|j\|_{1}}$.

Let $S\in\mathbb{N}$, $\sigma_{1},\ldots,\sigma_{S}$ be i.i.d. Rademacher variables, and $x_{1},\ldots,x_{S}$ be given samples.
Then, the empirical Rademacher complexity $\hat{R}_{S}(\mathcal{G}_{N})$ is bounded as follows.

###### Lemma 13

We have

$$ $\displaystyle\hat{R}_{S}(\mathcal{G}_{N})\leq\frac{\alpha}{\sqrt{S}}\max_{j\in N }\mathrm{e}^{\tau\|j\|_{1}}\sup_{\mathbf{L}_{1},\ldots,\mathbf{L}_{J}\in T(N, \mathbb{C})}\|\mathrm{e}^{\mathbf{L}_{1}}\|\cdots\|\mathrm{e}^{\mathbf{L}_{J}} \|\,\|v\|,$ $$

where $\alpha=\sum_{j\in\mathbb{Z}^{d}}\mathrm{e}^{-2\tau\|j\|_{1}}$.

We can see that the complexity of the model depends exponentially on both $N$ and $J$.
Combining Lemma 4.2 in Mohri et al. () and Lemma 13, we can derive Proposition 12.

###### Remark 14

The exponential dependence of the generalization bound on the number of layers is also typical for standard neural networks (Neyshabur et al., ; Bartlett et al., ; Golowich et al., ; Hashimoto et al., ).

###### Remark 15

Based on Proposition 12, we can control the generalization error by adding a regularization term to the loss function to make $\|\mathrm{e}^{\mathbf{L}_{1}}\|\cdots\|\mathrm{e}^{\mathbf{L}_{J}}\|$ smaller.

###### Proposition 12

###### Lemma 13

###### Remark 14

###### Remark 15

## 6 Numerical results

We empirically confirm the fundamental properties of the proposed deep Koopman-layered model.

### 6.1 Representation power and generalization

To confirm the fundamental property of the Koopman layer, we first consider an autonomous system.
Consider the van der Pol oscillator on $\mathbb{T}$

$$ $\frac{\mathrm{d}^{2}x(t)}{\mathrm{d}t^{2}}=-\mu(1-x(t)^{2})\frac{\mathrm{d}x(t )}{\mathrm{d}t}+x(t),$ (5) $$

where $\mu=3$.
By setting $\mathrm{d}x/\mathrm{d}t$ as a new variable, we regard Eq. (5) as a first-ordered system on the two-dimensional space.
We discretized Eq. (5) with the time-interval $\Delta t=0.01$, and generated 1000 time-series $\{x_{s,0},\ldots x_{s,100}\}$ for $s=1,\ldots,1000$ with different initial values distributed uniformly on $[-1,1]\times[-1,1]$.
We added a random noise, which was drawn from the normal distribution of mean 0 and standard deviation 0.01, to each $x_{s,j}$ and set it as $\tilde{x}_{s,j}$.
For training, we used the pairs $\{\tilde{x}_{s,0},\tilde{x}_{s,100}\}$ for $s=1,\ldots,1000$.
Then, we trained deep Koopman-layered models on $\mathbb{T}^{3}$ by minimizing the loss $\sum_{s=1}^{1000}\|Q_{N}\mathrm{e}^{\mathbf{L}_{1}}\cdots\mathrm{e}^{\mathbf{L
}_{J}}Q_{N}^{*}v(\tilde{x}_{s,0})-\tilde{x}_{s,100}\|^{2}$ using the Adam optimizer (Kingma & Ba, ) with the learning rate $0.001$.
We created data for testing in the same manner as the training dataset.
We set $v(x,y)=\sin(y)x_{1}+\cos(y)x_{2}$ for $x=[x_{1},x_{2}]\in\mathbb{T}^{2}$ and $y\in\mathbb{T}$.
Note that based on Remark 1, we constructed Kooman-layers on $\mathbb{T}^{d+1}$ for the input dimension $d$, and we designed the function $v$ so that it recovers $x_{1}$ by $v(x,\pi/2)$ and $x_{2}$ by $v(x,0)$.
We used the sine and cosine functions for designing $v$ since the representation space is constructed with the Fourier functions.
We set $N=\{n=[n_{1},n_{2},n_{3}]\in\mathbb{Z}^{3}\,\mid\,-5\leq n_{1},n_{2},n_{3}\leq
5
\}\setminus\{0\}$, $R=1$, and $M_{1}=\{n=[n_{1},n_{2},n_{3}]\in\mathbb{Z}^{3}\,\mid\,-2\leq n_{1},n_{2}\leq 2
,-1\leq n_{3}\leq 1\}\setminus\{0\}$ for all the layers.
We applied the Arnoldi method (Gallopoulos & Saad, ) to compute the exponential of $\mathbf{L}_{j}$.

Figure 1 (a) shows the test error for $J=1$ and $J=2$.
We can see that the performance becomes higher when $J=2$ than $J=1$.
Note that Theorem 5 is a fundamental result for autonomous systems, and according to Theorem 5, we may need more than one layer even for the autonomous systems.
The result reflects this theoretical result.
In addition, based on Remark 15, we added the regularization term $10^{-5}(\|\mathrm{e}^{\mathbf{L}_{1}}\|+\cdots+\|\mathrm{e}^{\mathbf{L}_{J}}\|)$ and observed the behavior.
We consider the case where the training data is noisy, and its sample size is small.
We generated training data as above, but the sample size was 30, and the standard deviation of the noise was 0.03.
We used the test data without the noise.
The sample size of the test data was 1000.
We set $J=3$ to consider the case where the number of parameters is large.
The result is illustrated in Figure 1 (b).
We can see that with the regularization, we can achieve smaller test errors than without the regularization, which implies that with the regularization, the model generalizes well.

Figure: (a) Without the regularization
Refer to caption: extracted/5897765/layers.png

### 6.2 Eigenvalues of the Koopman-layers for nonautonomous systems

To confirm that we can extract information about the underlined nonautonomous dynamical systems of time-series data using the deep Koopman-layered model, we observed the eigenvalues of the Koopman-layers.

#### 6.2.1 Measure-preserving dynamical system

Consider the nonautonomous dynamical system on $\mathbb{T}^{2}$

$$ $\bigg{(}\frac{\mathrm{d}x_{1}(t)}{\mathrm{d}t},\frac{\mathrm{d}x_{2}(t)}{ \mathrm{d}t}\bigg{)}=\bigg{(}-\frac{\partial\zeta}{\partial x_{2}}(t,x(t)), \frac{\partial\zeta}{\partial x_{1}}(t,x(t))\bigg{)}=:f(t,x),$ (6) $$

where $\zeta(t,[x_{1},x_{2}])=\mathrm{e}^{\kappa(\cos(x_{1}-t)+\cos x_{2})}$.
Since the dynamical system $f(t,\cdot)$ is measure-preserving for any $t\in\mathbb{R}$, the corresponding Koopman operator $K^{t}$ is unitary for any $t\in\mathbb{R}$.
Thus, the spectrum of $K^{t}$ is on the unit disk in the complex plane.
We discretized Eq. (6) with the time-interval $\Delta t=0.01$, and generated 1000 time-series $\{x_{s,0},\ldots x_{s,119}\}$ for $s=1,\ldots,1000$ for training with different initial values distributed uniformly on $[-1,1]\times[-1,1]$.
We split the data into 6 subsets $S_{t}=\{x_{s,j}\,\mid\,s\in\{1,\ldots,1000\},\ j\in\{20t,\ldots,20(t+1)-1\}\}$ for $t=0,\ldots,5$.
Then, we trained the model with 5 Kooman-layers on $\mathbb{T}^{3}$ by minimizing the loss $\sum_{j=1}^{5}\sum_{s=1}^{1000}\sum_{l=0}^{19}\|Q_{N}\mathrm{e}^{\mathbf{L}_{j
}}\cdots\mathrm{e}^{\mathbf{L}_{5}}Q_{N}^{*}v(x_{s,20(j-1)+l})-x_{s,100+l}\|^{2}$ using the Adam optimizer with the learning rate $0.001$.
In the same manner as Subsection 6.1, we set $v(x,y)=\sin(y)x_{1}+\cos(y)x_{2}$ for $x=[x_{1},x_{2}]\in\mathbb{T}^{2}$ and $y\in\mathbb{T}$.
Note that we trained the model so that $Q_{N}\mathrm{e}^{\mathbf{L}_{j}}\cdots\mathrm{e}^{\mathbf{L}_{5}}Q_{N}^{*}v$ maps samples in $S_{j-1}$ to $S_{5}$.
We set $N=\{n=[n_{1},n_{2},n_{3}]\in\mathbb{Z}^{3}\,\mid\,-5\leq n_{1},n_{2}\leq 5,-2
\leq n_{3}\leq 2\}$, $R=1$, and $M_{1}=\{n=[n_{1},n_{2},n_{3}]\in\mathbb{Z}^{3}\,\mid\,-2\leq n_{1},n_{2}\leq 2
,-1\leq n_{3}\leq 1\}$ for all the layers.
We applied the Arnoldi method to compute the exponential of $\mathbf{L}_{j}$.
In addition, we assumed the continuity of the flow of the nonautonomous dynamical system and added a regularization term $0.01\sum_{j=2}^{5}\|\mathrm{e}^{\mathbf{L}_{j}}-\mathrm{e}^{\mathbf{L}_{j-1}}\|$ to make the Koopman layers
next to each other become close.
After training the model sufficiently (after 3000 epochs), we computed the eigenvalues of the approximation $\mathrm{e}^{\mathbf{L}_{j}}$ of the Koopman operator for each layer $j=1,\ldots,5$.
For comparison, we estimated the Koopman operator $K^{t_{j}}_{j}$ using EDMD and KDMD (Kawahara, ) with the dataset $S_{j-1}$ and $S_{j}$ separately for $j=1,\ldots,5$.
For EDMD, we used the same Fourier functions $\{q_{j}\,\mid\,j\in N\}$
as the deep Koopman-layered model for the dictionary functions.
For KDMD, we transformed $[x_{1},x_{2}]\in\mathbb{T}^{2}$ into $\tilde{x}=[\mathrm{e}^{\mathrm{i}x_{1}},\mathrm{e}^{\mathrm{i}x_{2}}]\in
\mathbb{C}^{2}$ and applied the Gaussian kernel $k(x,y)=\mathrm{e}^{-0.1\|\tilde{x}-\tilde{y}\|^{2}}$.
For estimating $K^{t_{j}}_{j}$, we applied the principal component analysis to the space spanned by $\{k(\cdot,x)\,\mid\,x\in S_{j-1}\}$ to obtain $|N|$ principal vectors $p_{1},\ldots,p_{|N|}$.
We estimated $K^{t_{j}}_{j}$ by constructing the projection onto the space spanned by $p_{1},\ldots,p_{|N|}$.
Figure 2 illustrates the results.
We can see that the eigenvalues of the estimated Koopman operators by the deep Koopman-layered model are distributed on the unit circle for $j=1,\ldots,5$, which enables us to observe that the dynamical system is measure-preserving for any time.
On the other hand, the eigenvalues of the estimated Koopman operators with EDMD and KDMD are not on the unit circle, which implies that the separately applying EDMD and KDMD failed to capture the property of the dynamical system since the system is nonautonomous.

Figure: (a) Deep Koopman-layerd model
Refer to caption: extracted/5897765/vol_1.png

#### 6.2.2 Damping oscillator with external force

Consider the nonautonomous dynamical system regarding a damping oscillator on a compact subspace of $\mathbb{R}$

$$ $\frac{\mathrm{d}^{2}x(t)}{\mathrm{d}t^{2}}=-\alpha\frac{\mathrm{d}x(t)}{ \mathrm{d}t}-x(t)-a\sin(bt),$ (7) $$

where $\alpha=0.1$, $a=b=1$.
By setting $\mathrm{d}x/\mathrm{d}t$ as a new variable, we regard Eq. (7) as a first-ordered system on the two-dimensional space.
We generated data, constructed the deep Koopman-layered model, and applied EDMD and KDMD for comparison in the same manner as Subsection 6.2.1.
Figure 3 illustrates the results.
In this case, since the dynamical system is not measure preserving, it is reasonable that the estimated Koopman operators have eigenvalues inside the unit circle.
We can see that many eigenvalues for the deep Koopman-layered model are distributed inside the unit circle, and the distribution changes along the layers.
Since the external force becomes large as $t$ becomes large, the damping effect becomes small as $t$ becomes large (corresponding to $j$ becoming large).
Thus, the number of eigenvalues distributed inside the unit circle becomes small as $j$ becomes large.
On the other hand, we cannot obtain this type of observation from the separate estimation of the Koopman operators by EDMD and KDMD.

Figure: (a) Deep Koopman-layerd model
Refer to caption: extracted/5897765/force_1.png

## 7 Connection with other methods

### 7.1 Deep Koopman-layered model as a neural ODE-based model

The model (1) can also be regarded as a model with multiple neural ODEs (Teshima et al., ; Li et al., , Section 3.3).
From this perspective, we can also apply the model to standard tasks with ResNet.
For existing Neural ODE-based models, we use numerical methods, such as Runge-Kutta methods, to solve the ODE (Chen et al., ).
In our framework, solving the ODE corresponds to computing $\mathrm{e}^{\mathbf{L}_{j}}u$ for a matrix $\mathbf{L}_{j}$ and a vector $u$.
As we stated in Subsection 3.2, we use Krylov subspace methods to compute $\mathrm{e}^{\mathbf{L}_{j}}u$.
In this sense, our framework provides numerical linear algebraic way to train Neural ODE-based models by virtue of introducing Koopman generators and operators.

### 7.2 Connection with neural network-based Koopman approaches

In the framework of neural network-based Koopman approaches, we train an encoder $\phi$ and a decoder $\psi$ that minimizes $\|x_{t+1}-\psi(K\phi(x_{t}))\|$ for the given time-series
$x_{0},x_{1},\ldots$ (Lusch et al., ; Azencot et al., ; Shi & Meng, ).
Here, $K$ is a linear operator, and we can construct $K$ using EDMD or can train $K$ simultaneously with $\phi$ and $\psi$.
For neural network-based Koopman approaches, since the encoder $\phi$ changes along the learning process, the representation space of the operator $K$ also changes.
Thus, the theoretical analysis of these approaches is challenging.
On the other hand, our deep Koopman-layered approach fixes the representation space using the Fourier functions and learns only the linear operators corresponding to Koopman generators by restricting the linear operator to a form based on the Koopman operator.

## 8 Conclusion and discussion

In this paper, we proposed deep Koopman-layered models based on the Koopman operator theory combined with Fourier functions and Toeplitz matrices.
We showed that the Fourier basis forms a proper representation space of the Koopman operators in the sense of the universal and generalization property of the model.
In addition to the theoretical solidness, the flexibility of the proposed model allows us to train the model to fit time-series data coming from nonautonomous dynamical systems.

According to Lemma 10 and Theorem 5, to represent any function, we need more than one Koopman layer.
Investigating how many layers we need and how the representation power grows as the number of layers increases theoretically remains for future work.
In addition, we applied Krylov subspace methods to approximate the actions of the Koopman operators to vectors.
Since the Krylov subspace methods are iterative methods, we can control the accuracy of the approximation by controlling the iteration number.
How to decide and change the iteration number throughout the learning process for more efficient computations is also future work.

### Acknowledgements

We would like to thank Dr. Isao Ishikawa for a constructive discussion.

## Appendix

## Appendix A Proofs

We provide the proofs of statements in the main text.

Proposition 2
The $(n,l)$-entry of the representation matrix $Q_{N}^{*}L_{j}Q_{N}$ of the approximated operator is

$$ $\displaystyle\sum_{k=1}^{d}\sum_{n_{R_{j}}-l\in M^{j}_{R_{j}}}\sum_{n_{R_{j}-1 }-n_{R_{j}}\in M^{j}_{R_{j}-1}}\cdots\sum_{n_{2}-n_{3}\in M^{j}_{2}}\sum_{n-n_ {2}\in M^{j}_{1}}$ $\displaystyle\qquad\qquad\qquad\qquad a^{j,k}_{n_{R_{j}}-l,R_{j}}a^{j,k}_{n_{R _{j}-1}-n_{R_{j}},R_{j-1}}\cdots a^{j,k}_{n_{2}-n_{3},2}a^{j,k}_{n-n_{2},1} \mathrm{i}l_{k},$ $$

where $l_{k}$ is the $k$th element of the index $l\in\mathbb{Z}^{d}$.
Moreover, we set $n_{r}=m_{R_{j}}+\cdots+m_{r}+l$, thus $n_{1}=n$, $m_{r}=n_{r}-n_{r+1}$ for $r=1,\ldots,R_{j}-1$, and $m_{R_{j}}=n_{R_{j}}-l$.

Proof  We have

$$ $\displaystyle\left\langle q_{n},L_{j}q_{l}\right\rangle$ $\displaystyle=\bigg{\langle}q_{n},\sum_{k=1}^{d}\sum_{m_{R_{j}}\in M^{j}_{R_{j }}}a^{j,k}_{m_{R_{j}},R_{j}}q_{m_{R_{j}}}\cdots\sum_{m_{1}\in M^{j}_{1}}a^{j,k }_{m_{1},1}q_{m_{1}}\mathrm{i}l_{k}q_{l}\bigg{\rangle}$ $\displaystyle=\bigg{\langle}q_{n},\sum_{k=1}^{d}\sum_{m_{R_{j}}\in M^{j}_{R_{j }}}\cdots\sum_{m_{1}\in M^{j}_{1}}a^{j,k}_{m_{R_{j}},R_{j}}\cdots a^{j,k}_{m_{ 1},1}q_{m_{R_{j}}+\cdots+m_{1}+l}\mathrm{i}l_{k}\bigg{\rangle}$ $\displaystyle=\sum_{k=1}^{d}\sum_{\begin{subarray}{c}m_{R_{j}}+\cdots+m_{1}+l= n\\ m_{R_{j}}\in M^{j}_{R_{j}}\cdots m_{1}\in M^{j}_{1}\end{subarray}}a^{j,k}_{m_{ R_{j}},R_{j}}\cdots a^{j,k}_{m_{1},1}\mathrm{i}l_{k}$ $\displaystyle=\sum_{k=1}^{d}\sum_{n_{R_{j}}-l\in M^{j}_{R_{j}}}\sum_{n_{R_{j}- 1}-n_{R_{j}}\in M^{j}_{R_{j}-1}}\cdots\sum_{n_{2}-n_{3}\in M^{j}_{2}}\sum_{n-n _{2}\in M^{j}_{1}}$ $\displaystyle\qquad\qquad\qquad\qquad a^{j,k}_{n_{R_{j}}-l,R_{j}}a^{j,k}_{n_{R _{j}-1}-n_{R_{j}},R_{j-1}}\cdots a^{j,k}_{n_{2}-n_{3},2}a^{j,k}_{n-n_{2},1} \mathrm{i}l_{k}.$ $$

$\Box$

Corollary 6
Assume $v\in L_{0}^{2}(\mathbb{T}^{d})$ and $v\neq 0$.
For any sequence $g_{1}(t_{1},\cdot),\ldots,g_{\tilde{J}}(t_{\tilde{J}},\cdot)$ of flows that satisfies $v\circ g_{\tilde{J}}(t_{\tilde{J}},\cdot)\circ\cdots\circ g_{j}(t_{j},\cdot)
\in L^{2}_{0}(\mathbb{T}^{d})$ and $v\circ g_{\tilde{J}}(t_{\tilde{J}},\cdot)\circ\cdots\circ g_{j}(t_{j},\cdot)\neq
0$ for $j=1\ldots\tilde{J}$, and for any $\epsilon>0$, there exist a finite set $N\subset\mathbb{Z}\setminus\{0\}$, integers $0<J_{1}<\cdots<J_{\tilde{J}}$, and matrices $\mathbf{L}_{1},\ldots,\mathbf{L}_{J_{\tilde{J}}}\in T(N,\mathbb{C})$ such that $\|v\circ g_{\tilde{J}}(t_{\tilde{J}},\cdot)\circ\cdots\circ g_{j}(t_{j},\cdot)
-\mathbf{G}_{j}\|\leq\epsilon$ and $\mathbf{G}_{j}=Q_{N}\mathrm{e}^{\mathbf{L}_{J_{j-1}+1}}\cdots\mathrm{e}^{
\mathbf{L}_{J_{\tilde{J}}}}Q_{N}^{*}v$ for $j=1,\ldots,\tilde{J}$, where $J_{0}=1$.

Proof  Since $v\circ g_{\tilde{J}}(t_{\tilde{J}},\cdot)\circ\cdots\circ g_{j}(t_{j},\cdot)
\in L^{2}_{0}(\mathbb{T}^{d})$ and $v\circ g_{\tilde{J}}(t_{\tilde{J}},\cdot)\circ\cdots\circ g_{j}(t_{j},\cdot)\neq
0$, there exist finite $N_{j}\subset\mathbb{Z}^{d}\setminus\{0\}$ and $\mathbf{G}_{j}\in V_{N_{j}}$, $\mathbf{G}_{j}\neq 0$ such that $\|v\circ g_{\tilde{J}}(t_{\tilde{J}},\cdot)\circ\cdots\circ g_{j}(t_{j},\cdot)
-\mathbf{G}_{j}\|\leq\epsilon$ for $j=1,\ldots,\tilde{J}$.
Since $v\in L_{0}^{2}(\mathbb{T}^{d})$ and $v\neq 0$, there exist finite $N_{\tilde{J}+1}\subset\mathbb{Z}^{d}\setminus\{0\}$ such that $Q_{N_{\tilde{J}+1}}^{*}v\neq 0$.
Let $N=\bigcup_{j=1}^{\tilde{J}+1}N_{j}$.
By Lemma 11, since $Q_{N}^{*}v\neq 0$, there exist $J_{\tilde{J}-1},J_{\tilde{J}}\in\mathbb{N}$ and $\mathbf{L}_{J_{\tilde{J}-1}+1},\ldots,\mathbf{L}_{J_{\tilde{J}}}\in T(N,
\mathbb{C})$ such that $\mathbf{G}_{\tilde{J}}=Q_{N}\mathrm{e}^{\mathbf{L}_{J_{\tilde{J}-1}+1}}\cdots
\mathrm{e}^{\mathbf{L}_{J_{\tilde{J}}}}Q_{N}^{*}v$.
Since $\mathbf{G}_{\tilde{J}}\neq 0$, again by Lemma 11,
there exist
$J_{\tilde{J}-2}\in\mathbb{N}$ and $\mathbf{L}_{J_{\tilde{J}-2}+1},\ldots,\mathbf{L}_{J_{\tilde{J}-1}}\in T(N,
\mathbb{C})$ such that $\mathbf{G}_{\tilde{J}-1}=Q_{N}\mathrm{e}^{\mathbf{L}_{J_{\tilde{J}-2}+1}}
\cdots\mathrm{e}^{\mathbf{L}_{J_{\tilde{J}-1}}}\mathrm{e}^{\mathbf{L}_{J_{
\tilde{J}-1}+1}}\cdots\mathrm{e}^{\mathbf{L}_{J_{1}}}Q_{N}^{*}v=Q_{N}\mathrm{e
}^{\mathbf{L}_{J_{\tilde{J}-2}+1}}\cdots\mathrm{e}^{\mathbf{L}_{J_{\tilde{J}-1
}}}Q_{N}^{*}\mathbf{G}_{\tilde{J}}$.
We continue to apply Lemma 11 to obtain the result.
$\Box$

Lemma 9
Assume $N\subset\mathbb{Z}^{d}\setminus\{0\}$.
Then, we have $\mathbb{C}^{N\times N}=T(N,\mathbb{C})$.

Proof  We show $\mathbb{C}^{N\times N}\subseteq T(N,\mathbb{C})$.
The inclusion $\mathbb{C}^{N\times N}\supseteq T(N,\mathbb{C})$ is trivial.
Since $N\subset\mathbb{Z}^{d}\setminus\{0\}$, for any $n=[n_{1},\ldots,n_{d}]\in N$, there exists $k\in\{1,\ldots,d\}$ such that $\mathrm{i}n_{k}=(D_{k})_{n,n}\neq 0$.
We denote by $k_{\operatorname{min}}(n)$ the minimal index $k\in\{1,\ldots,d\}$ that satisfies $(D_{k})_{n,n}\neq 0$.
Let $B\in\mathbb{C}^{N\times N}$.
We decompose $B$ as $B=B_{1}+\ldots+B_{d}$, where $(B_{k})_{:,n}=B_{:,n}$ if $k=k_{\operatorname{min}}(n)$ and $(B_{k})_{:,n}=\mathbf{0}$ otherwise.
Here, $(B_{k})_{:,n}$ is the $n$th column of $B_{k}$.
Then, we have $(B_{k})_{:,n}=\mathbf{0}$ if $(D_{k})_{n,n}=0$.
Let $D_{k}^{+}$ be the diagonal matrix defined as $(D_{k}^{+})_{n,n}=1/(D_{k})_{n,n}$ if $(D_{k})_{n,n}\neq 0$ and $(D_{k}^{+})_{n,n}=0$ if $(D_{k})_{n,n}=0$.
In addition, let $C_{k}=B_{k}D_{k}^{+}$.
Then, we have $B=\sum_{k=1}^{d}C_{k}D_{k}$.
Applying Propostion 8, we have $B\in T(N,\mathbb{C})$, and obtain $\mathbb{C}^{N\times N}\subseteq T(N,\mathbb{C})$.
$\Box$

Lemma 11
For any $\mathbf{u},\mathbf{v}\in\mathbb{C}^{N}\setminus\{0\}$, there exists $A\in GL(N,\mathbb{C})$ such that $\mathbf{u}=A\mathbf{v}$.

Proof  Let $n_{0}\in N$ and let
$B\in\mathbb{N\times N}$ be defined as $B_{n,:}=1/\|\mathbf{v}\|^{2}\mathbf{v}^{*}$ for $n=n_{0}$ and so that $B_{n,:}$ and $B_{m,:}$ becoming orthogonal if $n\neq m$.
Then, the $n$th element of $B\mathbf{v}$ is $1$ for $n=n_{0}$ and is 0 0 for $n\neq n_{0}$.
Let $C\in\mathbb{N\times N}$ be defined as $C_{n,:}=\mathbf{u}$ for $n=n_{0}$ and so that $C_{n,:}$ and $C_{m,:}$ becoming orthogonal if $n\neq m$.
Then, $B,C\in GL(N,\mathbb{C})$ and $CB\mathbf{v}=\mathbf{u}$.
$\Box$

Lemma 13
We have

$$ $\displaystyle\hat{R}_{S}(\mathcal{G}_{N})\leq\frac{\alpha}{\sqrt{S}}\max_{j\in N }\mathrm{e}^{\tau\|j\|_{1}}\sup_{\mathbf{L}_{1},\ldots,\mathbf{L}_{J}\in T(N, \mathbb{C})}\|\mathrm{e}^{\mathbf{L}_{1}}\|\cdots\|\mathrm{e}^{\mathbf{L}_{J}} \|\,\|v\|,$ $$

where $\alpha=\sum_{j\in\mathbb{Z}^{d}}\mathrm{e}^{-2\tau\|j\|_{1}}$.

Proof

$$ $\displaystyle\hat{R}_{S}(\mathcal{G}_{N})$ $\displaystyle=\frac{1}{S}\mathrm{E}\bigg{[}\sup_{\mathbf{G}\in\mathcal{G}_{N}} \sum_{s=1}^{S}\mathbf{G}(x_{s})\sigma_{s}\bigg{]}=\frac{1}{S}\mathrm{E}\bigg{[ }\sup_{\mathbf{G}\in\mathcal{G}_{N}}\sum_{s=1}^{S}\iota_{N}\mathbf{G}(x_{s}) \sigma_{s}\bigg{]}$ $\displaystyle=\frac{1}{S}\mathrm{E}\bigg{[}\sup_{\mathbf{G}\in\mathcal{G}_{N}} \bigg{\langle}\sum_{s=1}^{S}\sigma_{s}\phi(x_{s}),\iota_{N}\mathbf{G}\bigg{ \rangle}\bigg{]}\leq\frac{1}{S}\sup_{\mathbf{G}\in\mathcal{G}_{N}}\|\iota_{N} \mathbf{G}\|_{\mathcal{H}_{K}}\bigg{(}\sum_{s=1}^{S}K(x_{s},x_{s})\bigg{)}^{1/2}$ $\displaystyle\leq\frac{\alpha}{\sqrt{S}}\sup_{\mathbf{G}\in\mathcal{G}_{N}}\| \iota_{N}\|\|\mathbf{G}\|_{L^{2}(\mathbb{T}^{d})}\leq\frac{\alpha}{\sqrt{S}} \max_{j\in N}\mathrm{e}^{\tau\|j\|_{1}}\sup_{\mathbf{L}_{1},\ldots,\mathbf{L}_ {J}\in T(N,\mathbb{C})}\|Q_{N}\mathrm{e}^{\mathbf{L}_{1}}\cdots\mathrm{e}^{ \mathbf{L}_{J}}Q_{N}^{*}v\|$ $\displaystyle\leq\frac{\alpha}{\sqrt{S}}\max_{j\in N}\mathrm{e}^{\tau\|j\|_{1} }\sup_{\mathbf{L}_{1},\ldots,\mathbf{L}_{J}\in T(N,\mathbb{C})}\|\mathrm{e}^{ \mathbf{L}_{1}}\|\cdots\|\mathrm{e}^{\mathbf{L}_{J}}\|\,\|v\|,$ $$

where $\alpha=\sum_{j\in\mathbb{Z}^{d}}\mathrm{e}^{-2\tau\|j\|_{1}}$.
$\Box$

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 5
presentation: 2
questions: I) What are some potential real-world applications of the deep Koopman-layered model? 

II)  The paper discusses the concept of multiple Koopman layers. What criteria should be used to determine the **optimal number of layers** in the model? Could you provide guidelines or heuristics for selecting the number of layers for different types of problems or datasets?

Please also see "Weaknesses".
soundness: 3
strengths: 1. **Originality**
The paper introduces an innovative integration of deep learning techniques with Koopman operator theory. By proposing deep Koopman-layered models that utilize learnable Toeplitz matrices and Fourier functions, the authors offer a fresh perspective for analyzing nonautonomous systems, which are often challenging for traditional methods. The approach of simultaneously estimating multiple Koopman operators is a notable contribution that broadens the model's applicability.

2. **Quality**
The authors provide a solid mathematical framework that substantiates their claims regarding the universality and generalization properties of the proposed model.

3. **Clarity**
The paper is well-structured and almost clearly written. The authors effectively convey the motivations behind their approach and highlight the significance of Koopman operator theory.

4. **Significance**
The significance of the paper lies in its potential impact on the fields of dynamical systems and machine learning. By bridging these two areas, the proposed deep Koopman-layered models could advance the analysis of time-series data from nonautonomous dynamical systems.
summary: The paper presents a novel deep Koopman-layered framework for modeling dynamical systems, particularly suited for nonautonomous time-series data. This approach integrates Koopman operator theory with Fourier functions and learnable Toeplitz matrices, enabling the simultaneous estimation of multiple Koopman operators to capture temporal dynamics. The model is theoretically robust and flexible, leveraging the universal and reproducing properties of Toeplitz matrices, which enhance its generalization capabilities. Efficient training is achieved through Krylov subspace methods.
weaknesses: 1. To my understanding, the choice of the Fourier basis in the proposed deep Koopman-layered model is motivated by its universality in function representation, desirable theoretical properties for analyzing Koopman operators, flexibility in learning multiple operators simultaneously, and compatibility with efficient Krylov subspace methods for low computational cost. However, it is still unclear:

- What specific properties of the Fourier basis make it particularly suitable for capturing the dynamics of nonautonomous systems compared to wavelet bases or polynomials?

- How does the use of the Fourier basis influence the convergence rates of the learning algorithms employed in the deep Koopman-layered model?

- Can the model's performance be enhanced by incorporating additional basis functions alongside the Fourier basis, and if so, how would this be implemented?

- How does the choice of the Fourier basis compare to other potential basis functions in terms of computational efficiency and accuracy when modeling complex nonautonomous dynamical systems?

2. The theoretical and numerical results in the paper are restricted to the torus. Although this has several advantages theoretically, it may limit the applicability of the theoretical results for many cases, especially for real-world scenarios. It would be great to know how the insights gained from the toroidal analysis may provide a foundational understanding to explore and develop models for more general cases, ultimately enhancing the applicability of the theoretical results to real-world scenarios. 

3. There are no experiments for chaotic dynamics. So, what happens if the training dynamics are more complex, such as in the case of (autonomous) chaotic dynamics with a **mixed Koopman spectrum**?

4. The sentence "By computing the eigenvalues of Koopman operators, we can understand the long-term behavior of the undelined dynamical systems" in the Introduction, line 29, may be better if mentioned specifically for systems with a discrete Koopman spectrum, but not in general (for example, not for systems that exhibit a continuous spectrum rather than a purely discrete set of eigenvalues).

---

confidence: 2
contribution: 2
initial_rating: 3
presentation: 1
questions: - In my understanding of your work, you learn an approximation to a time-series dataset over the time-window $t \in [t_1, t_J]$. In most time-series forecasting problems, one goal is to generate forecasts beyond the time-window in which the observations were collected. With your approach how do you generate forecasts for $t \geq t_J$?
- It wasn't clear to me under which conditions the Toeplitz matrix will be sparse in practice. Can you discuss in which cases this will be true? How will your approach scale in cases where you cannot make this assumption?
- My understanding is that you haven't placed requirements on your Fourier functions that they be invertible. In this case, is it possible to reconstruct the original space once you've made a forecast in the observable space?
soundness: 2
strengths: - Transforming the original space into the space of observables using fixed basis functions (as opposed to learning the space of observables like in Lusch 2017) allows the authors to develop a strong theoretical backing for their approach including showing universality and a generalization bound. I think the theoretical foundation of the work is its strongest feature.
- A nice feature of the approach is that one can analyze the eigenvalues of the trained model to interpret the underlying system (for example whether it is measure preserving or not).
- The idea of using Krylov subspace methods to approximate the matrix exponential matrix product in this context is novel as far as I'm aware.
summary: The authors present an approach for learning Koopman representations of time-series datasets. The main idea is to transform the data into a new space using fourier functions where the dynamics can be approximated by a switching linear dynamical system. They also propose to use Krylov subspace methods to approximate the matrix exponentials which define the time evolution of the linear system.
weaknesses: - Overall, I found it quite challenging to understand the details of your proposed learning algorithm and I am not confident that I would be able to reproduce your approach using your paper alone (that said, I appreciate the authors providing their code in the supplementary materials). It would have been helpful for me if you had included a section which provides a step-by-step summary or an algorithm block showing how your approach works for a given time-series dataset.
- I think the claims you make in S7 (that your approach provides an alternative to neural ODE-based models or other deep Koopman based approaches) is not well-supported by your numerical studies. To put these claims on solid footing you need to include numerical studies showing your approach can provide equivalent performance on some common time-series benchmarks; for example see the numerical studies [1].
- It wasn't clear to me how your approach can be used in time-series forecasting (i.e. making predictions for $t>t_J$.

[1] Wang, R., Dong, Y., Arik, S. Ö., & Yu, R. (2022). Koopman neural forecaster for time series with temporal distribution shifts. arXiv preprint arXiv:2210.03675.

Overall because the presentation of your approach was challenging to parse and because your numerical studies did not do a good job supporting some of your claims, I feel that I cannot recommend acceptance of your work at this time. 

*Minor comments which will not impact my decision*
- I noticed a couple of spelling mistakes -- please run your draft through a spell checker for the next submission.

---

confidence: 5
contribution: 2
initial_rating: 3
presentation: 2
questions: * I do not understand the sentence in the introduction: "In addition, since each Koopman operator for a time window is estimated individually in these frameworks, we cannot take the information of other Koopman operators into account.". A single system only has a single Koopman operator (family, if one considers one operator per time t). What do the authors mean by "other Koopman operators"? The family?

 * For proposition 5.1: is it not possible to use simple Monte-Carlo arguments to show that the variance of $h(x,y)$ goes to zero with $1/\sqrt{S}$ (i.e., $1/S \sum_i(h(x_i,y_i)$) converges to $E[h(x,y)]$ with the Monte-Carlo rate of $1/\sqrt{S})$? The current bound involves operator norms which may be unbounded, correct?

 * l478 "our framework provides numerical linear algebraic way to train Neural ODE-based models" which part of the "deep Koopman" model represents the neural network from the Neural ODE? The neural ODE framework was introduced by Duvenaud and coauthors (Chen et al. (2018)) to mitigate the memory consumption of back-propagation through classical solvers, they propose instead to use the adjoint method. Where does the memory footprint come in with the new method? "Finding the vector field of an ODE" is not the novelty of Chen et al. (2018), it is finding it *efficiently* (without a lot of memory cost).

 * How many internal layers would one need to approxiamte the simple non-autonomous system "\dot{x} = -x+t"?
soundness: 2
strengths: Constructing models for non-autonomous, non-linear dynamical systems is a very important and very challenging problem. On the n-torus, Fourier series provide a very efficient and useful dictioanry to approximate functions in L2, so the setting the authors chose to present their Koopman approximation is reasonable. The theoretical contributions cover most of the important, first aspects that should be treated when introducing a new model, and seem to be proven appropriately. The numerical examples are chosen well, and illustrate the efficacy of the new method. The connection to neural ODEs (or rather, approximation of ODE vector fields using neural networks) is interesting and should be explored further.
summary: The authors introduce a concatenation of Koopman operator approximations ("Deep Koopman model") defined on a Fourier approximation of the L^2 space on an n-torus, for the approximation of nonlinear, nonautonomous dynamical systems.
They propose to learn the generator in each of the Koopman operators (i.e., in each "layer" of the model) as Toeplitz matrices, and then map the learned generators to the corresponding operator using matrix exponentials.
The model is analyzed regarding universal approximation and generalization, and compared to KDMD and EDMD on several numerical examples on the torus.
weaknesses: The authors chose an extremely specific setting (systems on n-tori) and can therefore avoid the most challenging question in Koopman operator approximation - how to choose the proper dictionary for a given problem? The challenging aspect in their setting thus only comes from the fact that the chosen systems are non-autonomous, not that they are non-linear. Unfortunately, the authors do not discuss any existing methods for non-autonomous systems or compare their method in the numerical setting against them. In fact, the "deep Koopman model" proposed by the authors may not even be appropriate for general non-autonomous systems, because the number of internal "layers" may need to go to infinity for such systems, e.g. if the time-dependent part is not periodic (but, e.g., linear: $\dot{x}= -x+t$). I list those concerns in more detail below.

 * One of the most challenging aspects of numerical Koopman operator approximation is the choice of function space and corresponding truncation. The authors only work on the n-torus, for which the Fourier basis is a good choice - but this also means that they do not discuss (at all) the issues arising when the base space is not an n-torus (which, I would argue, is a very common case in practice).

 * The authors should incorporate and discuss existing literature on learning dictionary functions and switched systems:
   - Li, Qianxiao, Felix Dietrich, Erik M. Bollt, and Ioannis G. Kevrekidis. “Extended Dynamic Mode Decomposition with Dictionary Learning: A Data-Driven Adaptive Spectral Decomposition of the Koopman Operator.” Chaos: An Interdisciplinary Journal of Nonlinear Science 27, no. 10 (October 2017): 103111. https://doi.org/10.1063/1.4993854.
   - Peitz, Sebastian, and Stefan Klus. “Koopman Operator-Based Model Reduction for Switched-System Control of PDEs.” Automatica 106 (August 2019): 184–91. https://doi.org/10.1016/j.automatica.2019.05.016.

 * Experiment 6.2.1 is not using proper dictionary functions for KDMD and EDMD, and hence the numerical results are unfavourable for these methods compared to the presented method. The authors should either learn Koopman operators separately for each subset, or use time-dependent dictionary functions for EDMD - the true model is periodic with period $2\pi$, so it should be easily possible to properly encode the dynamics even with classical EDMD (and time-dependent dictionary).

---

confidence: 3
contribution: 3
initial_rating: 8
presentation: 2
questions: - Page 7, line 358 - 359: "according to Theorem 4.1, we may need more than one layer even for the autonomous systems." -- Indeed, the theorem asserts that "there exists some positive integer $J$", but it doesn't align w/ our understanding for autonomous systems - for which the Koopman operator (or generator) should be time-invariant as defined in section 2.3? Is the statement true? How do you explain the discrepancy?

- Even for non-autonomous systems, do you have suggestions on how to choose $J$ (page 8, line 397-398) and the function $v$ beforehand? Any general guidelines? I can see because we construct the representation space  with the Fourier functions, so we could use sinusoidal functions for representing $v$ for simplicity, etc. But are there deeper reasons for why you choose $v(x,y)$ that way (page 7, line 350)?

- Page 8, 416-417: what's the fundamental reason of seeing eigenvalues distributed on the unit circle in the deep Koopman-layered model? Why EDMD/KDMD failed? Since you are benchmarking against DMD methods on i) non-autonomous system (transient dynamics) and the goal is to ii) find stable eigenvalues, I strongly recommend (NOT required) you to compare against more relevant DMD variants below:
  - "Multi-Resolution Dynamic Mode Decomposition" by J Kutz et al.
  - "Forward-Backward Extended DMD with an Asymptotic Stability Constraint" by Forbes et al.
soundness: 3
strengths: The idea of using sequence of Toeplitz matrices to parameterize Koopman generator and leverage its universal property to prove the universality and generalization of the approximation scheme is novel and original. The work is solid in the sense that the framework allows explicit error control which is often missed in other Koopman operator learning literatures. Experimental results also confirms the importance of regularization as derived in the generalization bound and shows superior performance against some other approaches, in terms of the distribution of the identified eigenvalues (around the unit circle). Although the idea of “using a moving stencils method to compute the time-dependent Koopman operator (and in each local stencil, the dynamical system is assumed time invariant)” is not new (see eg. “Data-driven reduced-order modeling for nonautonomous dynamical systems in multiscale media” by Mengnan Li, Lijian Jiang), this work makes decent contribution to the workstream of Koopman operator learning for non-autonomous systems with theoretical basis.
summary: The paper proposes a deep Koopman-layered model designed to analyze the dynamics of time-series data, particularly focusing on nonautonomous dynamical systems, where the system's behavior changes over time. The model leverages Koopman operator theory, which provides a linear framework for analyzing nonlinear dynamical systems, and it is structured as a series of "Koopman layers", each approximating a Koopman operator at a specific time point. These layers allow the model to capture the evolving dynamics of the system. Toeplitz matrices are utilized to construct the Koopman layers. Krylov subspace methods are utilized for efficient computation of the Koopman operators. The model further uses the Fourier basis functions as the representation space for the Koopman operators. This choice is justified by the theoretical properties of the Fourier basis, including its universality and the ability to derive generalization bounds for the model. The paper establishes the universality of the proposed model, meaning that it can approximate any function within the chosen representation space with arbitrary accuracy, given enough layers. Additionally, the authors derive a generalization bound, indicating the model's ability to perform well on unseen data.
weaknesses: - I don’t agree with the statement “For existing Neural ODE-based models, we use numerical methods, such as Runge-Kutta methods, to solve the ODE  (Chen et al., 2018)” in section 7.1. In fact, one of the innovations in that work was the authors solved the adjoint equation backwards in time to obtain the gradients automatically for parameter updates, so that one can circumvent the memory issues caused by using other numerical methods, such as Runge-Kutta, although I do agree that using Krylov subspace methods to compute $e^{L_j u}$ in the work has similar potential too. 

- The notations are sometimes confusing or overloaded to digest. For example, page 3 line 161, it seems the notation of $x_k$ was not properly defined before its appearance? The superscripts and subscripts in equation (2), (3) are heavy, the authors may consider simplifying it with a bit abuse of notation OR show how it works 1D case and with 1 system so that j=1 and k=1 OR provide a more concrete examples to guide the readers through in appendix. 

- In addition, consider citing two papers from Mitsubishi lab
  - "Physics-Informed Koopman Network"
  - "Physics-informed neural ODE (PINODE): embedding physics into models using collocation points"

The spirit of this work aligns quite well with above two in that both proposes to learn the Koopman generator $L$. The difference lies in that this work is “data-driven” whereas the other two are “physics/equations-driven”.

Meta-review: This paper proposes to compose a series of learnable Koopman operator approximations to construct "deep Koopman-layered models", for the purpose of approximating nonlinear, nonautonomous differential equations. Each Koopman operator (i.e. each "layer") has its generator represented by Toeplitz matrices so that it can represent a Fourier approximation on an n-torus. Krylov subspace methods were also employed to improve the efficiency of training. While reviewers and I feel this is an interesting approach, there were major concerns unresolved, such as applicability beyond the quasiperiodic cases, fairness in the comparison with existing approaches, and the lack of empirical demonstration on more challenging problems. Recognizing the potential of the idea, I recommend the authors take the discussions into consideration and re-submit a substantially revised version.

Decision: Reject