There is an increasing interest in recent years in transferring the operator theoretic techniques such as the Koopman operator [1], [2] to the analysis of dynamical systems. Such operator based methods differ from the classical approaches, in that they define the evolution of observable functions in a function space rather than using state vectors in the state space. The power of these operator theoretic methods is that it provides linear representations of nonlinear timeinvariant systems, albeit in higher dimensional spaces that are countable or uncountable. Various numerical approaches, such as dynamic mode decomposition(DMD), Hankel-DMD, extended dynamic mode decomposition (E-DMD) and structured dynamic mode decomposition (S-DMD), have been proposed for discovering the Koopman operator of a nonlinear system, using a series of dictionary functions with spanning or universal function approximation properties [2]– [6]. Researchers have recently shown that it is possible to integrate machine-learned representations with dynamic mode decomposition algorithms, using variational autoencoders to achieve phase-dependent representations of spectra [7] or delay embeddings [8], shallow neural networks [3], linearly recurrent neural networks for balancing expressiveness and overfitting [9], and deep RELU feedforward networks for predictive modeling in biological and transmission systems [10]. E-DMD [3] and Deep-DMD [10] have been utilized in various domains, including nonlinear system identification [10]–[13], image processing [4], [14] and robotic control [15], [16].

Z. Liu, G. Ding and L. Chen are with the Department of Computer Science, University of Colorado, Boulder, CO 80309, USA (emails: {zhiyuan.liu, duohui.ding, lijun.chen}@colorado.edu).

E. Yeung is with the Department of Mechanical Engineering, the Center for Control, Dynamical Systems, and the Center for Biological Engineering, University of California, Santa Barbara, CA 93106, USA (email: eyeung@ucsb.edu).

Generally speaking, the learning especially the training phase of the Koopman operator tries to minimize the empirical loss based on the training set, e.g., the data sampled from the real trajectory of dynamic system. Compared to the traditional machine learning problem which learns the unknown mapping from input to output, the Koopman learning has two tasks: 1) Learning the function space that lifts state space to a high even infinite dimensional space, and 2) learning a linear mapping within that function space. These two tasks are highly intertwined, e.g., inappropriate function space learned will lead to poor learning performance even the linear mapping is perfect. However, to the best of our knowledge, the method of Koopman training has not gotten enough attention up to now. Another challenge is that, when parameterized function approximation such as neural network is used, the learning problem is nonconvex. For instance, even for a single layer neural network, it is NPcomplete to find the global optimal [17]. However, recent work [18]–[20] shows that for over-parameterized (wide) shallow neural networks, local optima provide satisfactory performance. Specifically, they show that every local optimum is global optimum if the hidden layer is non-singular; every local minimum of the simplified objective is close to the global minimum. In this paper, we contribute a proof of convergence for Koopman learning algorithms utilizing shallow neural networks, and derive conditions for first-order optimality, the properties of the so-called dictionary functions used in deep and E-DMD that guarantee convergence. We propose alternating optimization algorithm with an optimal tracker for training the Koopman operator. By proving the objective function’s smoothness property, we show that our algorithm admits O(1/T) convergence rate for chosen constant learning rate and O(1/ log T) for diminishing learning rate. We illustrate convergence of the alternating optimization algorithm for single-node training (non-distributed) on two nonlinear systems with oscillatory dynamics.

A second major contribution of this paper is the development of a distributed Koopman operator learning algorithm. Most Koopman operator learning algorithms operate under the assumption of full-state measurements. However, in engineered and natural systems represented by data, full-state measurements are often not available, or are too expensive to collect. For example, power distribution networks consisting of hundreds of thousands of nodes exhibit real-time dynamics on systems that are poorly modeled, calibrated, or dated. Biological networks operate on thousands of genes to generate transcriptomic reponse profiles as a function of time; full-state measurement via deep sequencing is prohibitively expensive. In many instances, it is much more feasible to collect measurements from select locations, via strategic placement of observers [21], which gives rise to a different form of data — time-series data that is spatially distributed or fragmented across the whole network. We address the challenge of training distributed representations of Koopman operators and develop a distributed Koopman learning algorithm, proving its asymptotic convergence, and illustrating predictive accuracy and convergence on several simulated examples.

The rest of the paper is organized as follows. Section II introduces the Koopman operator learning problem. Section III describes our alternating optimization algorithm for the Koopman learning and proves the convergence. Section IV extends the algorithm to a distributed setting and shows its convergence. Section V evaluates the performance of the two algorithms on two nonlinear systems. Section VI concludes the paper.

We consider a discrete time open-loop nonlinear dynamic system of the following form:

where and are continuously differentiable. The function f is the state-space model and the function h maps current state to a vector of observables or output . The Koopman operator K of system (1), if it exists, is a linear operator that acts on observable functions and forward propagates them in time. To be more specific, the Koopman operator for this system must satisfy the equations:

where is a basis function that defines the lifted space of observables and is a constant matrix. Based on the Koopman operator theory, is the basis function of observables under which is K-invariant for all n. This implies that the Koopman operator comprehensively captures the flow of the observable trajectory .

Based on the data-driven method [10] [22], a general model for approximating Koopman operator given the data trajectory can be formulated as follows:

The above model aims to minimize the empirical loss from the learning perspective. One can slightly change the objective function by adding certain regularized term, e.g., for sparse operator or for avoiding large training lost, to make the tradeoff between the training and generalization errors.

While there has been a surge of interest in using neural networks to perform Koopman learning, little is known regarding the convergence and numerical stability of the training process. This motivates us to investigate the property of Koopman learning during its training phase. There are two challenges in solving the optimization problem (5) in practice. First, the basis function is unknown. This makes it difficult to ascertain what functions and how many functions to include, let alone the minimal number of functions, to ensure K-invariant. Recently, EDMD [22] uses an expansive set of orthonomal polynomial basis functions, but this approach does not scale well and suffers from overfitting with an increasing number of basis functions. Deep-DMD [10] adopts the neural networks to approximate the basis function based on universal approximation theorem, but it lacks the theoretical guarantee, e.g., the stability and convergence. Second, the objective function is nonconvex. Therefore it is unrealistic to expect an algorithm to converge to global minima.

Here we focus on the basis function based on parametric method. Specifically, we redefine . The term means matrix product. A typical example is a fully connected one-layer neural network since for wide shallow neural network, local optima provide satisfactory under some mild conditions [18]–[20], where W is the layer parameter and is activation function. With the parametric basis method, problem (5) becomes

Although this problem is nonconvex, there are some interesting structures. For example, if we fix the parameter W of the basis function, optimizing K is a quadratic problem that finds the linear mapping from to . On the other hand, with fixed K, optimizing W is to adjust the parameter W to find the function space that satisfies the linear transformation mapping (this is still highly nonconvex but will reduce the complexity a lot). We thus consider the algorithm that alternatively optimize over W and K.

In this section, we first state our alternating algorithm and then investigate its convergence properties. Let and denote by the Frobenius norm. The detail of the algorithm is shown in Algorithm 1. Here E measures how far the gradient is from that at the critical point and track the best parameters so far. We make the following assumptions.

Assumption 1: The function is bounded and has a bounded gradient and Hessian.

Assumption 2: The parameters K and W are bounded, i.e., there exist two constant and such that and .

Assumption 1 looks strong. However, it holds for several popular activation functions such as logistic function

( ), hyperbolic tangent (tanh(x)), and inverse hyper- bolic tangent (arctan(x)). By Assumptions 1 and 2, one can verify that the objective function F is bounded, i.e., there exists a constant R such that . We can show that F has Lipschitz-continuous gradient with respect to the parameter W of basis functions.

Lemma 1: Under Assumptions 1 and 2 and given the data trajectory , we have

with

where and is the Lipschitz constant for the function .

Proof: First, denote by K[:, i] the i-th column of matrix the j-th row of matrix W, and the k-th dimension of . We can compute the element [j, k] of as:

and as:

where denotes the element-wise production. We can then write the gradient difference with respect to and as

So if we can show that is Lipschitz-continuous, the proof is done. We have

Consider function in R and its gradient By

Assumption 1, is bounded by some constant, denoted by . Let , we can bound as follows:

∥≤ − W∥W− W

where the last inequality is due to Cauchy-Schwarz inequality. Similarly, we can bound as follows:

∥≤

where the second inequality is by Assumption 2. Combining the above results, we have

Lemma 2: Under Assumption 1 and assume that the basis function is bounded by h, we have

and

With Lemmas 1 and 2, we now show that Algorithm 1 will converge to a critical point with convergence rate or . Theorem 1: Under Assumptions 1 and 2, Algorithm 1 for the Koopman operator learning will converge to a critical point. With constant learning rate , its convergence rate is ; and with diminishing learning rate , its convergence rate is . Proof: Since the objective function is Lipschitz gradient continuous with respect to K, the descent lemma [23] can be applied and we have

where tr denotes the trace of the matrix. The first equality is due to the gradient update of and the second equality is by the fact that tr.

As for the basis function’s parameter W, we can have the similar result since the objective function is Lipschitz gradient continuous with respect to W:

So by equations (7) and (8), we have the following for each complete update from :

We sum both sides of inequality (9) from and obtain

One can see that, each term on the right is non-negative and their summation is bounded by some constant. We can conclude that the alternating optimization algorithm will converge asymptotically to one critical point even without optimal tracker and when

Based on inequality (11), one can bound the minimum gradients up to T for Algorithm 1 as follows:

∥∇∥∇≤− F≤ 2RST .

(2) Diminishing learning rate If we choose the diminishing learning rate, e.g., , the result becomes

∥∇∥∇≤− F−

We know

So for diminishing stepsize, we can obtain

We see that the constant stepsize has the better convergence rate than the diminishing stepsize with the help of optimal tracker. Both cases show that

We now develop an algorithm to handle the learning problem for the Koopman operator of high dimensional nonlinear dynamical systems. Even if there are only a thousand states in the underlying nonlinear system, the dimension of the dictionary functions explodes exponentially with the number of states. Memory constraints thus make it infeasible to train a Koopman operator using a centralized or stand-alone computing node. This motivates the derivation of a scalable, distributed approximation algorithm to relieve this problem.

Assumption 3: The basis function can be decomposed or approximated by where is the new basis function for and is a subset of x with .

Based on Assumption 3, we can reformulate the centralized Koopman objective function as

lows. Denote by the set of computation nodes which can communicate with each other. For each computation node , it only store part of the data set , its corresponding row and column of Koopman operator and its basis function .

part can be calculated based on its own knowledge. Another part needs the information from other nodes. We first define by the error term for node i with data point , where

and define by the Jacobi matrix of function , we then have the following distributed Koopman learning algorithm shown in Algorithm 2.

1 Initialization:

2 node i randomly initilizes its . 3 while Not Converge do

4 = 0= 0= 05 for j = 1; + 1 do 6 node i calculates

13 for i ∈ Q do 14 W← W. 15 · · · · · · 16 end

17 end

For our distributed Koopman operator learning Algorithm 2, line 6-8 is the communication stage, each computation node i calculates its result in the lifted dimensional space and then broadcast it over a communication network. After the communication, the information is enough to compute local error term , and node i send to node v (line 8). Here the communication stage ends and computation stage (line 9-11) begins. will sum up all the information for each data point. The last is update stage with gradient descent method (line 14-15). Based on this distributed algorithm and Assumption 3, we can prove this is equivalent to the centralized gradient descent algorithm.

Lemma 3: Under Assumption 3, the distributed Koopman learning in Algorithm 2 is equivalent to the following update:

verify the following after updating all the data points:

Compared to gradients of F(W, K), we can find that

So the update stage (line 14-15) is the same with equation (13), which finishes the proof.

Remark 1: Our alternating Koopman operator learning (Algorithm 1) can be regarded as nonlinear Gauss-Seidel iterations [24], while our distributed Koopman operator learning lies in the model of nonlinear Jacobi iteration [24]. Here we choose nonlinear Jacobi iteration for distributed Koopman operator learning due to that nonlinear Jacobi iteration is (1) suitable for parallel computation and (2) with less communication overhead.

By lemma 3, we have the following convergence result for our distributed Koopman operator learning.

Theorem 2: Under Assumption 1, 2 and 3, the distributed Koopman operator learning based on Algorithm 2 will converge to a critical point asymptotically. Proof: The equation (13) is the case without optimal tracker of Algorithm 1. The proof of theorem 1 can be applied directly here (equation (11)).

The advantages of the distributed Koopman learning over the centralized one are not only the scalability, e.g., the ability to handle the high dimensional nonlinear dynamics, but also the feasibility to adjust to different complexity of the partial observations. For example, if one partial observation is with complexity dynamic, we can increase the number of basis function.

On the other hand, Algorithm 2 for distributed Koopman learning is under the ideal synchronous model. Although the computation of each node i is parallel, the computation will not start until the broadcast process finishes. This can lead to significant inefficiency in the distributed Koopman learning when, e.g., one node has very limited communication capability so that all other nodes need wait for this node. Also, one packet loss will lead to all nodes waiting for the resending. However, it is easily to extend Algorithm 2 to handle asynchronous mode as shown in Algorithm 3. Each node will store the received information in its memory, updating it when new information arrives. Once the computation node comes to computation stage, it will directly use the information stored in the memory instead of waiting for the newest information.

For the asynchronous version of gradient descent algorithm, existing work such as [25] (Theorem 2), [26] (Proposition 2.1), and [27] (Theorem 7) shows that the synchronous and asynchronous algorithms will converge to the same point once as long as communication delay is bounded. Their proof applies to our problem with slight modification.

Lemma 4: ( [25], [26], [27]) If the communication delay is bounded by some constant, with small enough stepsize, the asynchronous algorithm 3 will asymptotically converge to the same point as the synchronous one.

We now evaluate the performance of the proposed alternating optimization and distributed algorithms. For each experiment, we sample some points from the real trajectory to prepare the training set and prediction set. Note that our prediction phase is multi-step prediction, i.e., given one initial state, our algorithm will predict the following trajectory using the K-invariant property:

To evaluate the performance of alternating optimization, we consider Van der Pol oscillator shown in Example 1.

Example 1: Van der Pol oscillator

(15) In this example, we choose . The number of date points we sampled is 600 with 400 for training and 200 for prediction. We construct a very simple network with one layer and 3 dimensions to learn the pattern of Van der Pol oscillator. The total training time is around 1.08s (i7-8700 CPU @ 3.20 GHz, 8 GB RAM) with 500 iterations and constant stepsize is 0.23. Fig. 1 shows the multi-step prediction result with alternating optimization method. One step prediction error is around 0.16% and 200 step prediction error is around 1.89%.

Fig. 1: Alternating optimization for centralized Koopman operator learning with Van der Pol oscillator. In this experiments, only the points at time 0 are given. All the data points [1-200] are our predictions with Koopman learning.

Example 2: Glycolytic pathway

mented on a larger nonlinear dynamical system shown in Example 2, namely the glycolysis network from cellular biology [28]. We adopt the parameter setting: J = 2.5, a =

from [28]. 1000 data points are sampled from the real trajectory with 600 points for training and 400 for prediction. We create 7 threads to simulate the distributed learning and each thread only learn the dynamic pattern of one state by a simple 3-layer neural network with 15 dimensions. The total training time is 400.4s with 10000 iterations. Results for each state is shown in Fig. 2. One step error is around 0.02% and 400 step prediction error is around 2.7%. We see that our alternating optimization and distributed algorithms both achieve good performance with multi-step prediction. Even though partial state measurements are provided for training, the trained distributed Koopman operator is able to predict the behavior of the glycolysis network over 400 steps. Further, these results provide a glimmer of hope for whole cell network modeling, using strategically placed reporter libraries that provide partial measurements of an entire transcriptome [21], [29].

We have proposed an alternating optimization algorithm to the nonconvex Koopman operator learning problem for nonlinear dynamic systems. We prove that the proposed algorithm converges to a critical point with rate O(1/T) and for the constant and diminishing learning rates, respectively, under some mild conditions. To cope with the high dimensional nonlinear dynamical systems, we have further proposed a distributed Koopman operator learning algorithm with an appropriate communication mechanism. We show that the distributed Koopman operator learning is of the same convergence property with the centralized one if the basis functions are decomposable. Numerical experiments are provided to complement our theoretical results.

We thank Igor Mezic, Nathan Kutz, Robert Egbert, Bassam Bamieh, Sai Nandanoori Pushpak, Sean Warnick, Jongmin Kim, Umesh Vaidya, and Erik Bollt for stimulating conversations. Any opinions, findings and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the Defense Advanced Research Projects Agency (DARPA), the Department of Defense, or the United States Government. This work was supported partially by a Defense Advanced Research Projects Agency (DARPA) Grant No. DEAC0576RL01830 and an Institute of Collaborative Biotechnologies Grant.

[1] Bernard O Koopman. Hamiltonian systems and transformation in hilbert space. Proceedings of the National Academy of Sciences, 17(5):315–318, 1931.

[2] Igor Mezi´c. Spectral properties of dynamical systems, model reduction and decompositions. Nonlinear Dynamics, 41(1-3):309–325, 2005.

[3] Qianxiao Li, 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(10):103111, 2017.

[4] J Nathan Kutz, Xing Fu, and Steven L Brunton. Multiresolution dynamic mode decomposition. SIAM Journal on Applied Dynamical Systems, 15(2):713–735, 2016.

[5] Clarence W Rowley, Igor Mezi´c, Shervin Bagheri, Philipp Schlatter, and Dan S Henningson. Spectral analysis of nonlinear flows. Journal of fluid mechanics, 641:115–127, 2009.

[6] Subhrajit Sinha, Umesh Vaidya, and Enoch Yeung. On computation of koopman operator from sparse data. In 2019 American Control Conference (ACC), pages 5519–5524. IEEE, 2019.

[7] Bethany Lusch, J Nathan Kutz, and Steven L Brunton. Deep learning for universal linear embeddings of nonlinear dynamics. Nature communications, 9(1):4950, 2018.

[8] Naoya Takeishi, Yoshinobu Kawahara, and Takehisa Yairi. Learning koopman invariant subspaces for dynamic mode decomposition. In Advances in Neural Information Processing Systems, pages 1130– 1140, 2017.

[9] Samuel E Otto and Clarence W Rowley. Linearly recurrent autoencoder networks for learning dynamics. SIAM Journal on Applied Dynamical Systems, 18(1):558–593, 2019.

[10] Enoch Yeung, Soumya Kundu, and Nathan Hodas. Learning deep neural network representations for koopman operators of nonlinear dynamical systems. In 2019 American Control Conference (ACC), pages 4832–4839. IEEE, 2019.

[11] Charles A Johnson and Enoch Yeung. A class of logistic functions for approximating state-inclusive koopman operators. In 2018 Annual American Control Conference (ACC), pages 4803–4810. IEEE, 2018.

[12] Zhiyuan Liu, Soumya Kundu, Lijun Chen, and Enoch Yeung. Decomposition of nonlinear dynamical systems using koopman gramians. In 2018 Annual American Control Conference (ACC), pages 4811–4818. IEEE, 2018.

[13] Prashant G Mehta and Umesh Vaidya. On stochastic analysis approaches for comparing complex systems. In Proceedings of the 44th IEEE Conference on Decision and Control, pages 8082–8087. IEEE, 2005.

[14] J Nathan Kutz, Xing Fu, Steve L Brunton, and N Benjamin Erichson. Multi-resolution dynamic mode decomposition for foreground/background separation and object tracking. In 2015 IEEE International Conference on Computer Vision Workshop (ICCVW), pages 921–929. IEEE, 2015.

[15] Ian Abraham, Gerardo De La Torre, and Todd D Murphey. Model-based control using koopman operators. arXiv preprint arXiv:1709.01568, 2017.

[16] Erik Berger, Mark Sastuba, David Vogt, Bernhard Jung, and Heni Ben Amor. Estimation of perturbations in robotic behavior using dynamic mode decomposition. Advanced Robotics, 29(5):331–343, 2015.

[17] Avrim Blum and Ronald L Rivest. Training a 3-node neural network is np-complete. In Advances in neural information processing systems, pages 494–501, 1989.

Fig. 2: Distributed Koopman learning for Glycolytic pathway.

[18] Mahdi Soltanolkotabi, Adel Javanmard, and Jason D Lee. Theoretical insights into the optimization landscape of over-parameterized shallow neural networks. IEEE Transactions on Information Theory, 65(2):742–769, 2018.

[19] Digvijay Boob and Guanghui Lan. Theoretical properties of the global optimizer of two layer neural network. arXiv preprint arXiv:1710.11241, 2017.

[20] Anna Choromanska, Mikael Henaff, Michael Mathieu, G´erard Ben Arous, and Yann LeCun. The loss surfaces of multilayer networks. In Artificial Intelligence and Statistics, pages 192–204, 2015.

[21] Aqib Hasnain, Nibodh Boddupalli, and Enoch Yeung. Optimal reporter placement in sparsely measured genetic networks using the koopman operator. arXiv preprint arXiv:1906.00944, 2019.

[22] Matthew O Williams, Ioannis G Kevrekidis, and Clarence W Rowley. A data–driven approximation of the koopman operator: Extending dynamic mode decomposition. Journal of Nonlinear Science, 25(6):1307–1346, 2015.

[23] Yurii Nesterov. Introductory lectures on convex optimization: A basic course, volume 87. Springer Science & Business Media, 2013.

[24] MN Vrahatis, GD Magoulas, and VP Plagianakos. From linear to nonlinear iterative methods. Applied Numerical Mathematics, 45(1):59–77, 2003.

[25] Steven H Low and David E Lapsley. Optimization flow controli: basic algorithm and convergence. IEEE/ACM Transactions on Networking (TON), 7(6):861–874, 1999.

[26] Dimitri P Bertsekas and John N Tsitsiklis. Parallel and distributed computation: numerical methods, volume 23. Prentice hall Englewood Cliffs, NJ, 1989.

[27] Zhiyuan Liu and Lijun Chen. Proportional control applied to dynamic network resource allocation. In American Control Conference (ACC), 2017, pages 1948–1953. IEEE, 2017.

[28] Bryan C Daniels and Ilya Nemenman. Efficient inference of parsimonious phenomenological models of cellular dynamics using s-systems and alternating regression. PloS one, 10(3):e0119821, 2015.

[29] C Ward, E Yeung, T Brown, B Durtschi, S Weyerman, R Howes, and Jorge Goncalves. A comparison of network reconstruction methods for chemical reaction networks.