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

Thank you very much for your constructive comments. We addressed you comments and questions and revised the paper. We summarize the answers below.

[W1] About 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)”
We are sorry for the confusion and thank you for pointing out. We replaced the sentence with "For existing Neural ODE-based models, we solve ODEs for the forward computation and solve adjoint equations for backward computation (Chen et al., 2018)". In addition, we replaced the last sentence with "our framework provides numerical linear algebraic way to solve Neural ODE-based models".

[W2] Notations and readability
Thank you for pointing out. We added the definition of $x_k$ in line 161. We also omit the sub- and super-script of $j$ in Section 3.2 and added the explanation "For the remaining part of this section, we omit the subscript $j$ for simplicity. However, in practice, the approximation is computed for the generator $L_j$ for each layer $j=1,\ldots,J$."

[W3] Additional references
We added the reference to section 7. They are deeply related to this paper. Thank you for the references.

[Q1] Necessity of more than one layer even for the autonomous systems
This is an effect of the approximation of the generator. If we can use the true Koopman generator, then we only need one layer for autonomous systems. However, since we approximated the generator using matrices, we may need more than one layer. We added the above explanation right after the sentence.

[Q2] Choice of $J$ (number of layers) and $v$ (fixed nonlinear function in the proposed model)
Regarding the choice of $J$, heuristically, we can use validation data to determine an optimal number of layers. For example, we begin by one layer and compute the validation loss. Then, we set two layers and compute the validation loss, and continue with more layers. We can set the number of layers as the number that achieves the minimal validation loss. Another way is to set a sufficiently large number of layers and train the model with the validation data. As we discussed in Section 6.2, we can add a regularization term to the loss function so that the Koopman layers next to each other become close. After the training, if there are Koopman layers next to each other and sufficiently close, then we can regard them as one Koopman layer and determine an optimal number of layers. We added Appendix G for the above explanation.

Regarding the choice of $v$, as you pointed out, we chose $v$ so that it is suitable for the Fourier functions. Heuristically, we can parameterize $v$ and choose it using validation data or learn it simultaneously with the learnable parameters in the proposed model with training data. More deeper investigation is left to future work.

[Q3] Fundamental reason of seeing eigenvalues distributed on the unit circle in the deep Koopman-layered model and comparison against more relevant DMD variants
We guess that since EDMD and KDMD are for autonomous systems and we applied them individually for each time-window, they tried to interpret the underlined dynamical system as an autonomous system, and fail to capture nonautonomous behavior of the system. Based on other reviewers' comments, we added numerical results regarding the damping oscillator with external force (Section 6.3.2). We conducted experiments with the Koopman-based approach with learned dictionary functions. We considered the following two settings.

1. Learn a set of dictionary functions to construct the representation space of five Koopman generators. (learning a common set of dictionary functions is also considered in [1,2])
2. Learn five sets of dictionary functions each of which is for each Koopman generator.

We used a 3-layered ReLU neural network to learn the dictionary functions. The result is illustrated in Appendix E. We cannot capture the transition of the distribution of the eigenvalues through $j=1,\ldots,5$ even though we learned the dictionary functions. We can also see that there are some eigenvalues equally spaced on the unit circle. This behavior is typical for autonomous systems with a constant frequency. Since the dynamical system is nonautonomous and the frequency of the system changes over time, the above behavior is not suitable for this example. We conclude that the proposed method can capture the transition of the dynamics more properly than DMD-based method regardless of whether the dictionary functions are fixed or not.

We also implemented the forward-backward extended DMD (the second reference) with learned dictionary functions, based on your comment (The first one also seems to be deeply related to the proposed method, but since it is mainly for the standard DMD without Koopman operators, we focused more on the second reference). The result is also shown in Appendix E, and it is similar to the above two settings.

Thank you very much for your constructive comments. We addressed you comments and questions and revised the paper. We summarize the answers below.

[W1] Issues arising when the base space is not an n-torus
We do not think the analysis in the n-torus is a specific setting. In many practical cases, we are interested in dynamics in a bounded domain $\Omega$ in $\mathbb{R}^d$. For example, dynamics in a space around a certain object (e.g., heat source). In more physical settings, we are often interested in the dynamics around a limit cycle or an equilibrium point. Let $S^{d-1}$ be the $(d-1)$-dimensional sphere. If $\Omega$ is diffeomorphic to $S^{d-1}$, then we can construct a dynamical system $\check{f}_j$ on the torus $\mathbb{T}^d$ that satisfies $\check{f}_j(x)=\tilde{f}_j(x)$ for $x\in S^{d-1}$, where $\tilde{f}_j$ is the equivalent dynamical system on $S^{d-1}$ with $f_j$. For more details, please see Appendix B. We added Remark 3.2 for the above argument.

We agree that in some cases, it would be reasonable to consider $S^d$ directly or other types of domains (e.g. $\mathbb{R}^d$). In these cases, we can use the spherical harmonics or wavelet basis functions. However, the theoretical analysis in these cases is challenging and is left to future work.

[W2] Discussion of existing literature on learning dictionary functions and switched systems
Thank you for the references. We added the references. Regarding the first reference, we added it to Section 7.2 since in our understandings, it is mainly for autonomous systems. Regarding the second reference, in our understandings, each Koopman operator for each time window is estimated individually. We added it to Section 1.

We emphasize that an advantage of the proposed method over these approaches is that it has both theoretical solidness and flexibility. By fixing the dictionary functions, we can derive a strong theoretical guarantee; we can represent any function using the proposed model. At the same time, we can learn multiple Koopman operators simultaneously, which enables us to extract the transition of properties of dynamical systems along time.

[W3] Dictionary functions for EDMD and KDMD in the experiment
We used the same Fourier dictionary functions as the proposed deep Koopman-layered model for EDMD. Also, we applied the principal component analysis to select better dictionary functions for KDMD. Thus, we believe they are proper dictionary functions for the comparison with the proposed model. However, we agree that comparing the proposed method only to the method with fixed dictionary functions was not enough. Thus, we conducted additional experiments with the Koopman-based approach with learned dictionary functions. We considered the following two settings.

1. Learn a set of dictionary functions to construct the representation space of five Koopman generators. (learning a common set of dictionary functions is also considered in [1,2])
2. Learn five sets of dictionary functions each of which is for each Koopman generator.

We used a 3-layered ReLU neural network to learn the dictionary functions. The result is illustrated in Appendix E. We cannot capture the transition of the distribution of the eigenvalues through $j=1,\ldots,5$ even though we learned the dictionary functions. We can also see that there are some eigenvalues equally spaced on the unit circle. This behavior is typical for autonomous systems with a constant frequency. Since the dynamical system is nonautonomous and the frequency of the system changes over time, the above behavior is not suitable for this example. We conclude that the proposed method can capture the transition of the dynamics more properly than DMD-based method regardless of whether the dictionary functions are fixed or not.

[Q1] About the sentence ""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."
We are sorry for the confusion. This is about the existing works based on EDMD for nonautonomous systems with fixed dictionary functions. In these methods, each Koopman operator for each time-window is estimated separately. We revised the sentence.

[1] Rui Wang, Yihe Dong, Sercan O. Arik, and Rose Yu, "Koopman Neural Operator Forecaster for Time-series with Temporal Distributional Shifts", ICLR2023.
[2] Yong Liu, Chenyu Li, Jianmin Wang, and Mingsheng Long, "Koopa: Learning non-stationary time series dynamics with Koopman predictors", NeurIPS 2023.

Thank you very much for your response and additional comments, and thank you for updating the score.

Although we focus on the fundamental case of $\mathbb{T}^d$, the analysis on $\mathbb{T}^d$ opens up methods for more general cases. Indeed, $\mathbb{T}^d$ is the simplest example of locally compact groups, and the Fourier functions are generalized to the irreducible representations [1]. For a group $G$, a representation $\rho$ is a map $\rho:G\to GL(V)$ for some vector space $V$ that satisfies $\rho(x)\rho(y)=\rho(x\cdot y)$ for $x,y\in G$. Here, $GL(V)$ is the space of all bijective linear transformations from $V$ to $V$. We note that if $G$ is abelian, then $V$ is always one-dimensional. For example, for $\mathbb{R}^d$, the irreducible representations are $\rho\_\{\xi\}(x):=\mathrm{e}^{\mathrm{i}\xi \cdot x}$ for $x\in\mathbb{R}^d$ and $\xi\in\mathbb{R}^d$. As in Section 3.2, a crucial property of the Fourier function $q_n$ is the product of two Fourier functions are also a Fourier function, i.e., $q_n\cdot q_m=q_{m+n}$. This property is valid also for general irreducible representations since we have $\rho(x)\rho(y)=\rho(x\cdot y)$. Thus, we can consider a generalized version of Toeplitz matrices by using $\rho$ [2]. Although showing the universal property of the generalized Toeplitz matrices is challenging and future work, this type of argument gives us a promising way of generalizing our framework. We added this argument in Remark 3.2 and Appendix B.2.

We would like to emphasize that this work is the first step for constructing a method both with the theoretical guarantee and the flexibility. The torus is the most fundamental example of locally compact groups, and although there are some mathematically challenging problems, we have a clear path for the generalization of the proposed method.

We also added some explanations in the main text to improve the readability and corrected the typo. The revised parts are colored in blue. We would like to note that based on other reviewers' comments, we improved Subsection 6.1 and added Appendix E to help readers to understand the algorithm well. Also, we omitted the sub- and super-script $j$ in Section 3.2 to simplify the notation.

[1] William Fulton and Joe Harris. Representation Theory –A First Course–, Springer, 2004.
[2] Marc A. Rieffel. Matrix algebras converge to the sphere for quantum Gromov–Hausdorff distance. Memoirs of the American Mathematical Society, 168:67-91, 2004.

Thank you very much for your constructive comments. We addressed you comments and questions and revised the paper. We summarize the answers below.

[W1] Readability to understand the details of the proposed learning algorithm
We moved the paragraph that explains the details of the training the deep Koopman-layered model in Section 4 to the beginning of Section 6 as a new section 6.1. We also added more detailed explanation of how to compute the output of the model so that readers can understand the details of the proposed learning algorithm. Moreover, we added the pseudocode and more details in Appendix D. We hope the revised version is helpful for readers to understand the details of the algorithm and reproduce it.

[W2] Numerical studies to show the proposed approach can provide equivalent performance on some common time-series benchmarks
An advantage of the Koopman-based approaches over neural ODE-based approaches is that we can extract the information of the transition of the dynamics by computing the eigenvalues of the approximated Koopman operators. Therefore, we compared the proposed method to other Koopman-based approaches in the sense that how we can extract the information of the transition of the dynamics. We admit that comparing the proposed method only to the method with fixed dictionary functions was not enough. Thus, we conducted additional experiments with the Koopman-based approach with learned dictionary functions. We considered the following two settings.

1. Learn a set of dictionary functions to construct the representation space of five Koopman generators. (learning a common set of dictionary functions is also considered in [1,2])
2. Learn five sets of dictionary functions each of which is for each Koopman generator.

We used a 3-layered ReLU neural network to learn the dictionary functions. The result is illustrated in Appendix E. We cannot capture the transition of the distribution of the eigenvalues through $j=1,\ldots,5$ even though we learned the dictionary functions. We can also see that there are some eigenvalues equally spaced on the unit circle. This behavior is typical for autonomous systems with a constant frequency. Since the dynamical system is nonautonomous and the frequency of the system changes over time, the above behavior is not suitable for this example. We conclude that the proposed method can capture the transition of the dynamics more properly than DMD-based method regardless of whether the dictionary functions are fixed or not.

[W3, Q1] Application to time-series forecasting for $t>t_J$
As we also insist above, an advantage of the proposed method is that we can extract the information of the transition of the dynamics by computing the eigenvalues of the approximated Koopman operators. We can also forecast the time-series starting with a different initial value for $t\le t_J$. Our main goal in this paper is not the time-series forecasting for $t>t_J$. However, as you pointed out, applying the proposed method to time-series forecasting is also interesting. Applying the idea of [1,2], we can decompose the Koopman operators into time-invariant and time-variant parts. By extracting time-invariant features of the dynamics using the approximated Koopman operators (e.g., time-invariant eigenvectors or singular vectors), we can combine it with time-variant Koopman operators constructed by local time-series to construct the forecast.

More precisely, we can decompose the Koopman operator $K^t$ for time $t$ as $K^t=K_{inv}+K^t_{var}$, where $K_{inv}=\sum_{i=1}^n\sigma_iv_iu_i^*$ and $K_{var}=\sum_{i=1}^m\tilde{\sigma}_i\tilde{v}_i\tilde{u}_i^*$, $\sigma_i$, $v_i$, $u_i$ are time-invariant singular values and the corresponding singular vectors of the approximated Koopman operators for $j=1,\ldots, J$, $\tilde{v}_i$ are the singular vectors of the local Koopman operator that is orthogonal to $v_i$, and $\tilde{\sigma}_i$ and $\tilde{u}_i$ are singular values and singular vectors corresponding to $\tilde{v}_i$. Since we can use the time-invariant property of $t\le t_J$, we can forecast time-series well even for $t>t_J$. We added Appendix F for explaining the above point.

[1] Rui Wang, Yihe Dong, Sercan O. Arik, and Rose Yu, "Koopman Neural Operator Forecaster for Time-series with Temporal Distributional Shifts", ICLR2023.
[2] Yong Liu, Chenyu Li, Jianmin Wang, and Mingsheng Long, "Koopa: Learning non-stationary time series dynamics with Koopman predictors", NeurIPS 2023.

[W2] Restriction of the theoretical and numerical results to the torus
We do not think the analysis in the torus limits the applicability of the proposed model for many cases. In many practical cases, we are interested in dynamics in a bounded domain $\Omega$ in $\mathbb{R}^d$. For example, dynamics in a space around a certain object (e.g., heat source). In more physical settings, we are often interested in the dynamics around a limit cycle or an equilibrium point. Let $B_d$ be the unit ball in $\mathbb{T}^d$. If $\Omega$ is diffeomorphic to $B_d$, then we can construct a dynamical system $\check{f}_j$ on the torus $\mathbb{T}^d$ that satisfies $\check{f}_j(x)=\tilde{f}_j(x)$ for $x\in B_d$, where $\tilde{f}_j$ is the equivalent dynamical system on $B_d$ with $f_j$. For more details, please see Appendix B. We added Remark 3.2 for the above argument.

We agree that in some cases, it would be reasonable to consider other types of domains (e.g. $\mathbb{R}^d$ and $S^{d}$). In these cases, we can use the spherical harmonics or wavelet basis functions, as you pointed out. However, the theoretical analysis in these cases is challenging and is left to future work.
(We corrected typos in this comments besed on the pointing out by Reviewer naMD.)

[W2] Experiments for chaotic dynamics
For chaotic dynamical systems, the Koopman operator has continuous and residual spectra. Thus, it is difficult to capture the behavior of the system only with the eigenvalues. The situation is the same for the standard DMD-based methods, and not specific for the proposed method. Dealing with more complex dynamics is future work.

[W4] About the sentence "By computing the eigenvalues of Koopman operators, we can understand the long-term behavior of the undelined dynamical systems"
Thank you for the suggestion. We are sorry for the confusion, and we agree that the explanation is only for the Koopman operators with discrete spectra. We revised the sentence.

[Q1] Potential real-world applications
As we discussed above, for example, we can analyze dynamics around a certain object (e.g. heat source) in $\mathbb{R}^d$ using the proposed model. For the example of a heat source, the amount of heat created by the source may switch or change over time. In addition, the environment around the source may also switch or change over time (e.g. by wind, by humans). In this case, we can use the proposed model to analyze the transition of the temperature around the heat source.

[Q2] Optimal number of layers
Heuristically, we can use validation data to determine an optimal number of layers. For example, we begin by one layer and compute the validation loss. Then, we set two layers and compute the validation loss, and continue with more layers. We can set the number of layers as the number that achieves the minimal validation loss. Another way is to set a sufficiently large number of layers and train the model with the validation data. As we discussed in Section 6.2, we can add a regularization term to the loss function so that the Koopman layers next to each other become close. After the training, if there are Koopman layers next to each other and sufficiently close, then we can regard them as one Koopman layer and determine an optimal number of layers. We added Appendix G for the above explanation.

Thank you very much for your constructive comments. We addressed you comments and questions and revised the paper. We summarize the answers below.

[W1-1] Specific properties of the Fourier basis making it particularly suitable for capturing the dynamics of nonautonomous systems compared to wavelet bases or polynomials
The following two properties of Fourier functions $q_n(x)=\mathrm{e}^{\mathrm{i}n \cdot x}$ are suitable for the deep structure of the Koopman-layered model;

1. They form an orthonormal basis of $L^2(\mathbb{T}^d)$ and the RKHS on $\mathbb{T}^d$ considered in Section 5.
2. $\frac{d q_n}{dx}=\mathrm{i}nq_n$ and $q_n\cdot q_m=q_{n+m}$ (for simplicity, we consider the case of $d=1$, but same properties are valid for $d\ge 2$).

The first property enables us to construct the model in a dense subspace of $L^2(\mathbb{T}^d)$ and the RKHS and to derive generalization bound (since they are also an orthonormal basis of an RKHS). By the second property, if a Koopman generator acts on a Fourier function, the resulting function is also represented by the weighted sum of the Fourier functions. This feature is useful for the deep structure of the Koopman-layered model. In addition, the second property enables us to approximate Koopman operators by Toeplitz and diagonal matrices, which have nice properties such as the universality and the computational efficiency (please see below for further details). By virtue of the above two properties of the Fourier functions, we have the theoretical advantages (universality and generalization bound) and the computational advantage.

[W1-2] Influence of the use of the Fourier basis on the convergence rates of the learning algorithms
In the same manner as Theorem 4.1, we can show that we can represent any function in $V_N=\operatorname{Span}\\{q_n\,\mid\,n\in N\\}$ exactly using the deep Koopman-layered model, where $N$ is the index set characterizes the approximation space of the Koopman operators. Thus, if the decay rate of the Fourier transform of the target function $h$ is $\alpha$, then the convergence rate with respect to $N$ is $O(({1-\alpha^2})^{-d/2})$. We added Remark 4.4 and Appendix C regarding the above argument. Please see Appendix C for more details.

[W1-3] Performance enhancement by incorporating additional basis functions alongside the Fourier basis
Thank you for the suggestion. Formally, we can replace some of the Fourier basis functions with other functions. However, in that case, the approximated Koopman operators have more complicated structures and it is challenging to analyze it theoretically. Moreover, the computational cost may become expensive since the operators are not represented by Toeplitz matrices any more. On the other hand, by adding other functions than the Fourier functions, the convergence rate may become better. This paper is the first step to propose a theoretically guaranteed method for analyzing nonautonomous dynamical systems using Koopman operators. Investigating more sophisticated way of choosing basis functions is future work.

[W1-4] Comparison of Fourier basis to other potential basis functions in terms of computational efficiency and accuracy
In terms of the computational efficiency, the product of an $n$ by $n$ Toeplitz matrix and a vector is efficiently computed with the cost of $O(n\log n)$ by using the fast Fourier transform (FFT). As we discussed before Remark 3.4, if we set the index set $M_r$ much smaller than $N$, then the Toeplitz matrix is sparse. However, even if the Toeplitz matrices are dense, the computational cost of one iteration of the Krylov subspace method for computing $\mathrm{e}^{\mathbf{L}}v$ for a vector $v$ is $O(\vert N\vert \log\vert N\vert)$. If we use other basis functions such as wavelet and polynomial functions, then the approximation of a Koopman generator do not have a structure, and we have to deal with it as a standard dense matrix. Then, the computational cost is $O(\vert N\vert^2)$. Thank you for pointing out, and we added the explanation before Remark 3.4. In terms of the accuracy, as we insisted in the answer of the second question (about the convergence rate), we can represent any function in $\operatorname{Span}\\{q_n\,\mid\,n\in N\\}$ exactly using the deep Koopman-layered model.