Learning Koopman Invariant Subspaces for Dynamic Mode Decomposition

Naoya Takeishi, Yoshinobu Kawahara, Takehisa Yairi

Introduction

A variety of time-series data are generated from nonlinear dynamical systems, in which a state evolves according to a nonlinear map or differential equation. In summarization, regression, or classification of such time-series data, precise analysis of the underlying dynamical systems provides valuable information to generate appropriate features and to select an appropriate computation method. In applied mathematics and physics, the analysis of nonlinear dynamical systems has received significant interest because a wide range of complex phenomena, such as fluid flows and neural signals, can be described in terms of nonlinear dynamics. A classical but popular view of dynamical systems is based on state space models, wherein the behavior of the trajectories of a vector in state space is discussed (see, e.g., ). Time-series modeling based on a state space is also common in machine learning. However, when the dynamics are highly nonlinear, analysis based on state space models becomes challenging compared to the case of linear dynamics.

Recently, there is growing interest in operator-theoretic approaches for the analysis of dynamical systems. Operator-theoretic approaches are based on the Perron–Frobenius operator or its adjoint, i.e., the Koopman operator (composition operator) . The Koopman operator defines the evolution of observation functions (observables) in a function space rather than state vectors in a state space. Based on the Koopman operator, the analysis of nonlinear dynamical systems can be lifted to a linear (but infinite-dimensional) regime. Consequently, we can consider modal decomposition, with which the global characteristics of nonlinear dynamics can be inspected . Such modal decomposition has been intensively used for scientific purposes to understand complex phenomena (e.g., ) and also for engineering tasks, such as signal processing and machine learning. In fact, modal decomposition based on the Koopman operator has been utilized in various engineering tasks, including robotic control , image processing , and nonlinear system identification .

One of the most popular algorithms for modal decomposition based on the Koopman operator is dynamic mode decomposition (DMD) . An important premise of DMD is that the target dataset is generated from a set of observables that spans a function space invariant to the Koopman operator (referred to as Koopman invariant subspace). However, when only the original state vectors are available as the dataset, we must prepare appropriate observables manually according to the underlying nonlinear dynamics. Several methods have been proposed to utilize such observables, including the use of basis functions and reproducing kernels . Note that these methods work well only if appropriate basis functions or kernels are prepared; however, it is not always possible to prepare such functions if we have no a priori knowledge about the underlying dynamics.

In this paper, we propose a fully data-driven method for modal decomposition via the Koopman operator based on the principle of learning Koopman invariant subspaces (LKIS) from scratch using observed data. To this end, we estimate a set of parametric functions by minimizing the residual sum of squares (RSS) of linear least-squares regression, so that the estimated set of functions transforms the original data into a form in which the linear regression fits well. In addition to the principle of LKIS, an implementation using neural networks is described. Moreover, we introduce empirical performance of DMD based on the LKIS framework with several nonlinear dynamical systems and applications, which proves the feasibility of LKIS-based DMD as a fully data-driven method for modal decomposition via the Koopman operator.

Background

We focus on a (possibly nonlinear) discrete-time autonomous dynamical system

Here, the value of gg is decomposed into a sum of Koopman modes wi=φi(x0)ciw_{i}=\varphi_{i}(\bm{x}_{0})c_{i}, each of which evolves over time with its frequency and decay rate respectively given by ∠λi\angle\lambda_{i} and ∣λi∣|\lambda_{i}|, since λi\lambda_{i} is a complex value. The Koopman modes and their eigenvalues can be investigated to understand the dominant characteristics of complex phenomena that follow nonlinear dynamics. The above discussion can also be applied straightforwardly to continuous-time dynamical systems .

Modal decomposition based on K\mathcal{K}, often referred to as Koopman spectral analysis, has been receiving attention in nonlinear physics and applied mathematics. In addition, it is a useful tool for engineering tasks including machine learning and pattern recognition; the spectra (eigenvalues) of K\mathcal{K} can be used as features of dynamical systems, the eigenfunctions are a useful representation of time-series for various tasks, such as regression and visualization, and K\mathcal{K} itself can be used for prediction and optimal control. Several methods have been proposed to compute modal decomposition based on K\mathcal{K}, such as generalized Laplace analysis , the Ulam–Galerkin method , and DMD . DMD, which is reviewed in more detail in the next subsection, has received significant attention and been utilized in various data analysis scenarios (e.g., ).

Note that the Koopman operator and modal decomposition based on it can be extended to random dynamical systems actuated by process noise . In addition, Proctor et al. discussed Koopman analysis of systems with control signals. In this paper, we primarily target autonomous deterministic dynamics (e.g., Eq. (1)) for the sake of presentation clarity.

2 Dynamic mode decomposition and Koopman invariant subspace

where m+1m+1 is the number of snapshots in the dataset. The core functionality of DMD algorithms is computing the eigendecomposition of matrix A=Y1Y0†\bm{A}=\bm{Y}_{1}\bm{Y}_{0}^{\dagger} , where Y0†\bm{Y}_{0}^{\dagger} is the Moore–Penrose pseudoinverse of Y0\bm{Y}_{0}. The eigenvectors of A\bm{A} are referred to as dynamic modes, and they coincide with the Koopman modes if the corresponding eigenfunctions of K\mathcal{K} are in span⁡{g1,…,gn}\operatorname{span}\{g_{1},\dots,g_{n}\} . Alternatively (but nearly equivalently), the condition under which DMD works as a numerical realization of Koopman spectral analysis can be described as follows.

Here, an important problem in the practice of DMD arises, i.e., we often have no access to g\bm{g} that spans a Koopman invariant subspace G\mathcal{G}. In this case, for nonlinear dynamics, we must manually prepare adequate observables. Several researchers have addressed this issue; Williams et al. leveraged a dictionary of predefined basis functions to transform original data, and Kawahara defined Koopman spectral analysis in a reproducing kernel Hilbert space. Brunton et al. proposed the use of observables selected in a data-driven manner from a function dictionary. Note that, for these methods, we must select an appropriate function dictionary or kernel function according to the target dynamics. However, if we have no a priori knowledge about them, which is often the case, such existing methods do not have to be applied successfully to nonlinear dynamics.

Learning Koopman invariant subspaces

In this paper, we propose a method to learn a set of observables {g1,…,gn}\{g_{1},\dots,g_{n}\} that spans a Koopman invariant subspace G\mathcal{G}, given a sequence of measurements as the dataset. In the following, we summarize desirable properties for such observables, upon which the proposed method is constructed.

If ∀x∈M, g(f(x))=Gg(x)\forall\bm{x}\in\mathcal{M},~{}\bm{g}(\bm{f}(\bm{x}))=G\bm{g}(\bm{x}), then for any g^=∑i=1naigi∈span⁡{g1,…,gn}\hat{g}=\sum_{i=1}^{n}a_{i}g_{i}\in\operatorname{span}\{g_{1},\dots,g_{n}\},

According to Theorem 1, we should obtain g\bm{g} that makes g∘f−Gg\bm{g}\circ\bm{f}-G\bm{g} zero. However, such problems cannot be solved with finite data because g\bm{g} is a function. Thus, we give the corresponding empirical risk minimization problem based on the assumption of ergodicity of f\bm{f} and the convergence property of the empirical matrix as follows.

Define Y0\bm{Y}_{0} and Y1\bm{Y}_{1} by Eq. (4) and suppose that Assumption 1 holds. If all modes are sufficiently excited in the data (i.e., rank⁡(Y0)=n\operatorname{rank}(\bm{Y}_{0})=n), then matrix A=Y1Y0†\bm{A}=\bm{Y}_{1}\bm{Y}_{0}^{\dagger} almost surely converges to the matrix form of linear operator GG in m→∞m\to\infty.

Since A=Y1Y0†\bm{A}=\bm{Y}_{1}\bm{Y}_{0}^{\dagger} is the minimum-norm solution of the linear least-squares regression from the columns of Y0\bm{Y}_{0} to those of Y1\bm{Y}_{1}, we constitute the learning problem to estimate a set of function that transforms the original data into a form in which the linear least-squares regression fits well. In particular, we minimize RSS, which measures the discrepancy between the data and the estimated regression model (i.e., linear least-squares in this case). We define the RSS loss as follows:

which becomes zero when g\bm{g} spans a Koopman invariant subspace. If we implement a smooth parametric model on g\bm{g}, the local minima of LRSS\mathcal{L}_{\text{RSS}} can be found using gradient descent. We adopt g\bm{g} that achieves a local minimum of LRSS\mathcal{L}_{\text{RSS}} as a set of observables that spans (approximately) a Koopman invariant subspace.

2 Linear delay embedder for state space reconstruction

In this paper, we propose to surrogate the parameter selection of the delay-coordinate embedding by learning a linear delay embedder from data. Formally, we learn embedder ϕ\bm{\phi} such that

3 Reconstruction of original measurements

and, if h\bm{h} is a smooth parametric model, this term can also be reduced using gradient descent. Finally, the objective function to be minimized becomes

where α\alpha is a parameter that controls the balance between LRSS\mathcal{L}_{\text{RSS}} and Lrec\mathcal{L}_{\text{rec}}.

4 Implementation using neural networks

In Sections 3.1–3.3, we introduced the main concepts for the LKIS framework, i.e., RSS loss minimization, learning the linear delay embedder, and reconstruction of the original measurements. Here, we demonstrate an implementation of the LKIS framework using neural networks.

After estimating the parameters of ϕ\bm{\phi}, g\bm{g}, and h\bm{h}, DMD can be performed normally by using the values of the learned g\bm{g}, defining the data matrices in Eq. (4), and computing the eigendecomposition of A=Y1Y0†\bm{A}=\bm{Y}_{1}\bm{Y}_{0}^{\dagger}; the dynamic modes are obtained by w\bm{w}, and the values of the eigenfunctions are obtained by φ=zHg\varphi=\bm{z}^{\mathsf{H}}\bm{g}, where w\bm{w} and z\bm{z} are the right- and left-eigenvectors of A\bm{A}. See Section 2.2 for details.

In the numerical experiments described in Sections 5 and 6, we performed optimization using first-order gradient descent. To stabilize optimization, batch normalization was imposed on the inputs of hidden layers. Note that, since RSS loss function (5) is not decomposable with regard to data points, convergence of stochastic gradient descent (SGD) cannot be shown straightforwardly. However, we empirically found that the non-decomposable RSS loss was often reduced successfully, even with mini-batch SGD. Let us show an example; the full-batch RSS loss (denoted LRSS⋆\mathcal{L}^{\star}_{\text{RSS}}) under the updates of the mini-batch SGD are plotted in the rightmost panel of Figure 4. Here, LRSS⋆\mathcal{L}^{\star}_{\text{RSS}} decreases rapidly and remains small. For SGD on non-decomposable losses, Kar et al. provided guarantees for some cases; however, examining the behavior of more general non-decomposable losses under mini-batch updates remains an open problem.

Related work

The proposed framework is motivated by the operator-theoretic view of nonlinear dynamical systems. In contrast, learning a generative (state-space) model for nonlinear dynamical systems directly has been actively studied in machine learning and optimal control communities, on which we mention a few examples. A classical but popular method for learning nonlinear dynamical systems is using an expectation-maximization algorithm with Bayesian filtering/smoothing (see, e.g., ). Recently, using approximate Bayesian inference with the variational autoencoder (VAE) technique to learn generative dynamical models has been actively researched. Chung et al. proposed a recurrent neural network with random latent variables, Gao et al. utilized VAE-based inference for neural population models, and Johnson et al. and Krishnan et al. developed inference methods for structured models based on inference with a VAE. In addition, Karl et al. proposed a method to obtain a more consistent estimation of nonlinear state space models. Moreover, Watter et al. proposed a similar approach in the context of optimal control. Since generative models are intrinsically aware of process and observation noises, incorporating methodologies developed in such studies to the operator-theoretic perspective is an important open challenge to explicitly deal with uncertainty.

We would like to mention some studies closely related to our method. After the first submission of this manuscript (in May 2017), several similar approaches to learning data transform for Koopman analysis have been proposed . The relationships and relative advantages of these methods should be elaborated in the future.

Numerical examples

In this section, we provide numerical examples of DMD based on the LKIS framework (LKIS-DMD) implemented using neural networks. We conducted experiments on three typical nonlinear dynamical systems: a fixed-point attractor, a limit-cycle attractor, and a system with multiple basins of attraction. We show the results of comparisons with other recent DMD algorithms, i.e., Hankel DMD , extended DMD , and DMD with reproducing kernels . The detailed setups of the experiments discussed in this section and the next section are described in the appendix.

Consider a two-dimensional nonlinear map on xt=[x1,tx2,t]T\bm{x}_{t}=\begin{bmatrix}x_{1,t}&x_{2,t}\end{bmatrix}^{\mathsf{T}}:

which has a stable equilibrium at the origin if λ,μ<1\lambda,\mu<1. The Koopman eigenvalues of system (9) include λ\lambda and μ\mu, and the corresponding eigenfunctions are φλ(x)=x1\varphi_{\lambda}(\bm{x})=x_{1} and φμ(x)=x2−x12\varphi_{\mu}(\bm{x})=x_{2}-x_{1}^{2}, respectively. λiμj\lambda^{i}\mu^{j} is also an eigenvalue with corresponding eigenfunction φλiφμj\varphi_{\lambda}^{i}\varphi_{\mu}^{j}. A minimal Koopman invariant subspace of system (9) is span⁡{x1,x2,x12}\operatorname{span}\{x_{1},x_{2},x_{1}^{2}\}, and the eigenvalues of the Koopman operator restricted to such subspace include λ\lambda, μ\mu and λ2\lambda^{2}. We generated a dataset using system (9) with λ=0.9\lambda=0.9 and μ=0.5\mu=0.5 and applied LKIS-DMD (n=4n=4), linear Hankel DMD (delay 2), and DMD with basis expansion by {x1,x2,x12}\{x_{1},x_{2},x_{1}^{2}\}, which corresponds to extended DMD with a right and minimal observable dictionary. The estimated Koopman eigenvalues are shown in Figure 3, wherein LKIS-DMD successfully identifies the eigenvalues of the target invariant subspace. In Figure 3, we show eigenvalues estimated using data contaminated with white Gaussian observation noise (σ=0.1\sigma=0.1). The eigenvalues estimated by LKIS-DMD coincide with the true values even with the observation noise, whereas the results of DMD with basis expansion (i.e., extended DMD) are directly affected by the observation noise.

We generated data from the limit cycle of the FitzHugh–Nagumo equation

where a=0.7a=0.7, b=0.8b=0.8, c=0.08c=0.08, and I=0.8I=0.8. Since trajectories in a limit-cycle are periodic, the (discrete-time) Koopman eigenvalues should lie near the unit circle. Figure 4 shows the eigenvalues estimated by LKIS-DMD (n=16n=16), linear Hankel DMD (delay 8), and DMDs with reproducing kernels (polynomial kernel of degree 4 and RBF kernel of width 1). The eigenvalues produced by LKIS-DMD agree well with those produced by kernel DMDs, whereas linear Hankel DMD produces eigenvalues that would correspond to rapidly decaying modes.

where α=1\alpha=1, β=−1\beta=-1, and δ=0.5\delta=0.5. States x\bm{x} following (11) evolve toward [10]T\begin{bmatrix}1&0\end{bmatrix}^{\mathsf{T}} or [−10]T\begin{bmatrix}-1&0\end{bmatrix}^{\mathsf{T}} depending on which basin of attraction the initial value belongs to unless the initial state is on the stable manifold of the saddle. Generally, a Koopman eigenfunction whose continuous-time eigenvalue is zero takes a constant value in each basin of attraction ; thus, the contour plot of such an eigenfunction shows the boundary of the basins of attraction. We generated 1,000 episodes of time-series starting at different initial values uniformly sampled from 2^{2}. The left plot in Figure 5 shows the continuous-time Koopman eigenvalues estimated by LKIS-DMD (n=100n=100), all of which correspond to decaying modes (i.e., negative real parts) and agree with the property of the data. The center plot in Figure 5 shows the true basins of attraction of (11), and the right plot shows the estimated values of the eigenfunction corresponding to the eigenvalue of the smallest magnitude. The surface of the estimated eigenfunction agrees qualitatively with the true boundary of the basins of attractions, which indicates that LKIS-DMD successfully identifies the Koopman eigenfunction.

Applications

The numerical experiments in the previous section demonstrated the feasibility of the proposed method as a fully data-driven method for Koopman spectral analysis. Here, we introduce practical applications of LKIS-DMD.

Prediction of a chaotic time-series has received significant interest in nonlinear physics. We would like to perform the prediction of a chaotic time-series using DMD, since DMD can be naturally utilized for prediction as follows. Since g(xt)\bm{g}(\bm{x}_{t}) is decomposed as ∑i=1nφi(xt)ci\sum_{i=1}^{n}\varphi_{i}(\bm{x}_{t})\bm{c}_{i} and φ\varphi is obtained by φi(xt)=ziHg(xt)\varphi_{i}(\bm{x}_{t})=\bm{z}_{i}^{\mathsf{H}}\bm{g}(\bm{x}_{t}) where zi\bm{z}_{i} is a left-eigenvalue of K\bm{K}, the next step of g\bm{g} can be described in terms of the current step, i.e., g(xt+1)=∑i=1nλi(ziHg(xt))ci\bm{g}(\bm{x}_{t+1})=\sum_{i=1}^{n}\lambda_{i}(\bm{z}_{i}^{\mathsf{H}}\bm{g}(\bm{x}_{t}))\bm{c}_{i}. In addition, in the case of LKIS-DMD, the values of g\bm{g} must be back-projected to y\bm{y} using the learned h\bm{h}. We generated two types of univariate time-series by extracting the {x}\{x\} series of the Lorenz attractor and the Rossler attractor . We simulated 25,000 steps for each attractor and used the first 10,000 steps for training, the next 5,000 steps for validation, and the last 10,000 steps for testing prediction accuracy. We examined the prediction accuracy of LKIS-DMD, a simple LSTM network, and linear Hankel DMD , all of whose hyperparameters were tuned using the validation set. The prediction accuracy of every method and an example of the predicted series on the test set by LKIS-DMD are shown in Figure 7. As can be seen, the proposed LKIS-DMD achieves the smallest root-mean-square (RMS) errors in the 30-step prediction.

One of the most popular applications of DMD is the investigation of the global characteristics of dynamics by inspecting the spatial distribution of the dynamic modes. In addition to the spatial distribution, we can investigate the temporal profiles of mode activations by examining the values of corresponding eigenfunctions. For example, assume there is an eigenfunction φλ≪1\varphi_{\lambda\ll 1} that corresponds to a discrete-time eigenvalue λ\lambda whose magnitude is considerably smaller than one. Such a small eigenvalue indicates a rapidly decaying (i.e., unstable) mode; thus, we can detect occurrences of unstable phenomena by observing the values of φλ≪1\varphi_{\lambda\ll 1}. We applied LKIS-DMD (n=10n=10) to a time-series generated by a far-infrared laser, which was obtained from the Santa Fe Time Series Competition Data . We investigated the values of eigenfunction φλ≪1\varphi_{\lambda\ll 1} corresponding to the eigenvalue of the smallest magnitude. The original time-series and values of φλ≪1\varphi_{\lambda\ll 1} obtained by LKIS-DMD are shown in Figure 7. As can be seen, the activations of φλ≪1\varphi_{\lambda\ll 1} coincide with sudden decays of the pulsation amplitudes. For comparison, we applied the novelty/change-point detection technique using one-class support vector machine (OC-SVM) and direct density-ratio estimation by relative unconstrained least-squares importance fitting (RuLSIF) . We computed AUC, defining the sudden decays of the amplitudes as the points to be detected, which were 0.924, 0.799, and 0.803 for LKIS, OC-SVM, and RuLSIF, respectively.

Conclusion

In this paper, we have proposed a framework for learning Koopman invariant subspaces, which is a fully data-driven numerical algorithm for Koopman spectral analysis. In contrast to existing approaches, the proposed method learns (approximately) a Koopman invariant subspace entirely from the available data based on the minimization of RSS loss. We have shown empirical results for several typical nonlinear dynamics and application examples.

We have also introduced an implementation using multi-layer perceptrons; however, one possible drawback of such an implementation is the local optima of the objective function, which makes it difficult to assess the adequacy of the obtained results. Rather than using neural networks, the observables to be learned could be modeled by a sparse combination of basis functions as in but still utilizing optimization based on RSS loss. Another possible future research direction could be incorporating approximate Bayesian inference methods, such as VAE . The proposed framework is based on a discriminative viewpoint, but inference methodologies for generative models could be used to modify the proposed framework to explicitly consider uncertainty in data.

This work was supported by JSPS KAKENHI Grant No. JP15J09172, JP26280086, JP16H01548, and JP26289320.

References

Appendix A Algorithm of dynamic mode decomposition

Dynamic mode decomposition (DMD) was originally invented as a tool for inspecting fluid flows , and it has been utilized in various fields other than fluid dynamics. An output of DMD coincides with Koopman spectral analysis if we have g\bm{g} that spans a Koopman invariant subspace. The popular algorithm of DMD , which is based on the singular value decomposition (SVD) of a data matrix, is defined as follows.

Given a sequence of g(x)\bm{g}(\bm{x}), define data matrices Y0=[g(x0)⋯g(xm−1)]\bm{Y}_{0}=\begin{bmatrix}\bm{g}(\bm{x}_{0})&\cdots&\bm{g}(\bm{x}_{m-1})\end{bmatrix} and Y1=[g(f(x0))⋯g(f(xm−1))]\bm{Y}_{1}=\begin{bmatrix}\bm{g}(\bm{f}(\bm{x}_{0}))&\cdots&\bm{g}(\bm{f}(\bm{x}_{m-1}))\end{bmatrix}.

Calculate the compact SVD of Y0\bm{Y}_{0} as Y0=UrSrVrH\bm{Y}_{0}=\bm{U}_{r}\bm{S}_{r}\bm{V}_{r}^{\mathsf{H}}, where rr is the rank of Y0\bm{Y}_{0}.

Normalize w\bm{w} and z\bm{z} such that wiHzj=δi,j\bm{w}_{i}^{\mathsf{H}}\bm{z}_{j}=\delta_{i,j}, where δi,j=1\delta_{i,j}=1 if i=ji=j and δi,j=0\delta_{i,j}=0 otherwise.

Return dynamic modes w\bm{w} and corresponding eigenvalues λ\lambda. In addition, return the values of corresponding eigenfunctions by φ=zHg\varphi=\bm{z}^{\mathsf{H}}\bm{g}.

The eigenvalues computed using the algorithm above are “discrete-time” ones, in the sense that they represent frequencies and decay rates in terms of discrete-time dynamical systems. The “continuous-time” counterparts can be computed easily by λc=log⁡(λ)/Δt\lambda_{c}=\log(\lambda)/\Delta t, where Δt\Delta t is the time interval in the discrete-time setting.

Note that the definition above presumes the access to an appropriate observable g\bm{g}. In contrast, in the proposed method, the above-mentioned algorithm is run after applying the proposed framework to learn g\bm{g} from observed data.

Appendix B Detailed experimental setup

In this appendix section, the configurations of the numerical examples and applications, which were omitted in the main text, are described.

One must not subtract the mean from the original data because subtracting something from the data may change the spectra of the underlying dynamical systems (see, e.g., ). If the absolute values of the data were too large, we simply divided the data by the maximum absolute value for each series.

In optimization, we found that the adaptive learning rate by SMORMS3 achieved fast convergence compared to a fixed learning rate and other adaptation techniques. The maximum learning rate of SMORMS3 was selected from 10−310^{-3} to 10−210^{-2} in each experiment according to the amount of data. In some cases, optimization was performed in two stages: the parameters of ϕ\bm{\phi}, g\bm{g}, and h\bm{h} were updated in the first stage, and, in the second stage, the parameters of ϕ\bm{\phi} and g\bm{g} were fixed and only h\bm{h} was updated. This two-stage optimization was particularly useful for the application of prediction, where a precise reconstruction of the original measurements was necessary. Moreover, when the original states x\bm{x} of the dynamical system were available and used without delay (i.e., k=1k=1 and p=rp=r), parameter Wϕ\bm{W}_{\phi} of the linear embedder was fixed to be an identity matrix (i.e., no embedder was used). Also, we set the mini-batch size from 100 to 500 because smaller mini-batches often led to an unstable computation of pseudo-inverse.

B.2 Fixed-point attractor experiment

In the experiment using the fixed-point attractor, the data were generated with four initial values: [55]T\begin{bmatrix}5&5\end{bmatrix}^{\mathsf{T}}, [−55]T\begin{bmatrix}-5&5\end{bmatrix}^{\mathsf{T}}, [5−5]T\begin{bmatrix}5&-5\end{bmatrix}^{\mathsf{T}}, and [−55]T\begin{bmatrix}-5&5\end{bmatrix}^{\mathsf{T}}, with the length of each episode being 3030. In the case of noisy dataset, the standard deviation of the observation noise was set to 0.10.1. In both experiments (with and without observation noise), we set k=2k=2 and n=4n=4 to cover the minimal three-dimensional Koopman invariant subspace.

B.3 Limit-cycle attractor experiment

The data were generated using MATLAB’s ode45 function , which was run with time-step Δt=0.1\Delta t=0.1 and initial value x0=[11.6]T\bm{x}_{0}=\begin{bmatrix}1&1.6\end{bmatrix}^{\mathsf{T}} for 2,000 steps. The hyperparameters of LKIS-DMD, linear Hankel DMD, and kernel DMDs were set such that they produced 16 eigenvalues, i.e., k=8k=8 and n=16n=16 for LKIS-DMD, and POD modes whose singular value was less than ε\varepsilon were disposed in kernel DMDs (ε=0.0001\varepsilon=0.0001 for the polynomial kernel and ε=0.05\varepsilon=0.05 for the RBF kernel).

B.4 Multiple basins of attraction experiment

The data were generated using the settings provided in the literature ; 1,000 initial values were drawn from the uniform distribution on ×\times and each initial value was proceeded in time for 11 steps with Δt=0.25\Delta t=0.25. We used MATLAB’s ode45 function for numerical integration. For LKIS-DMD, we set k=1k=1 and n=100n=100. Note that the values of the estimated eigenfunction were evaluated and plotted in consideration of each data point.

B.5 Chaotic time-series prediction experiment

The data were generated from the Lorenz attractor (parameters β=\nicefrac83\beta=\nicefrac{{8}}{{3}}, σ=10\sigma=10, and ρ=28\rho=28) and the Rossler attractor (parameters a=0.2a=0.2, b=0.2b=0.2, and c=5.7c=5.7). We generated 25,000 steps for each attractor and divided them into training, validation, and test sets. For all methods, the delay dimension was fixed at 77, i.e., k=7k=7 for LKIS-DMD and linear Hankel DMD, and backpropagation was truncated to length 77 to learn the LSTM network. We tuned nn of LKIS-DMD and the dimensionality of LSTM’s hidden state (denoted nhn_{h}) according to the 30-step prediction accuracies obtained using the validation set. Here, we obtained n=5n=5 and nh=5n_{h}=5 for the Lorenz data and n=6n=6 and nh=3n_{h}=3 for the Rossler data.

In this experiment, LSTM was applied because it had been utilized for various nonlinear time-series, and Hankel DMD was used because it had been successfully utilized for analysis of chaotic systems .

B.6 Unstable phenomena detection experiment

The dataset was obtained from the Santa Fe Time Series Competition Data . Note that the author’s original web page was not available on the date of submission of this manuscript (May 2017); however, the dataset itself was still available online. The length of delay (or sliding window) was fixed to 1010 for all methods applied in this experiment. In addition, no intensive tuning of the other hyperparameters was conduct because the purpose was qualitative. The default settings of libsvm were used for the one-class SVM (except for ν=0.05\nu=0.05). For the density-ratio estimation by RuLSIF, the default values of the implementation by the authors of were used.

In this experiment, OC-SVM was applied because it was a kind of de facto standard for novelty/change-point detection, and RuLSIF ws used because it had achieved the best performance among methods based on density-ratio estimation .