Koopman operator framework for spectral analysis and identification of infinite-dimensional systems

We consider Koopman operator theory in the context of nonlinear infinite-dimensional systems, where the operator is defined over a space of bounded continuous functionals. The properties of the Koopman semigroup are described and a finite-dimensional projection of the semigroup is proposed, which provides a linear finite-dimensional approximation of the underlying infinite-dimensional dynamics. This approximation is used to obtain spectral properties from data, a method which can be seen as a generalization of Extended Dynamic Mode Decomposition for infinite-dimensional systems. Finally, we exploit the proposed framework to identify (a finite-dimensional approximation of) the Lie generator associated with the Koopman semigroup. This approach yields a linear method for nonlinear PDE identification, which is complemented with theoretical convergence results.


INTRODUCTION
Koopman operator theory is a powerful framework which provides an alternative approach to dynamical systems. Through the so-called Koopman (or composition) operator (Koopman (1931)), nonlinear dynamical systems are approximated by linear, higher-dimensional systems that are amenable to systematic analysis. Under the impulse of the seminal work Mezić (2005), the Koopman operator framework has grown in popularity over the last decade in dynamical systems theory, and more recently has attracted attention in control theory (see  and references therein). However, the research effort has focused on finite-dimensional systems, and little work has been devoted to infinite-dimensional dynamical systems, such as nonlinear partial differential equations (PDEs). In this context, one can mention the early work by Banks (1985), where (an equivalent of) the Koopman operator was defined on a separable L 2 space of functionals (themselves defined on a Hilbert space). Later on, Dorroh and Neuberger (1996) followed a different path, defining the composition operator in the space of bounded continuous functionals and investigating the properties of the associated Lie generator. In the recent work by Farkas and Kreidler (2020), this approach was pushed further and leveraged in the framework of strongly continuous semigroups by considering a space of bounded continuous functionals equipped with a mixed topology. Finally, it is also recently that the Koopman operator framework has been considered to study the spectral properties of nonlinear PDEs (Mezić (2020); Nakao and Mezić (2020)).
In this paper, we further exploit and investigate the Koopman operator framework for systems described by infinite-dimensional differential equations. We do not rely on strongly continuous semigroup theory, but rather adopt the approach proposed by Dorroh and Neuberger (1996). In this context, the Lie generator is shown to be related to a Gâteaux derivative and a finite-dimensional approximation of the Koopman operator is presented, which allows to approximate a nonlinear, infinite-dimensional system by a linear, finite-dimensional one. Based on this framework, our main contributions are twofold. First, we complement the spectral analysis in Nakao and Mezić (2020) by developing a data-driven method to compute the spectral properties of the Koopman operator. This method is a generalization of the Extended Dynamic Mode Decomposition (EDMD) method (Williams et al. (2015)) to infinitedimensional systems and is more general than the mere application of standard EDMD to spatially discretized PDEs. In fact, while the method is developed for general basis functionals, it narrows down to the EDMD method only for a specific choice of basis functionals. Second, combining our previous work (Mauroy and Goncalves (2020)) with the above framework, we propose a novel method for nonlinear identification of infinite-dimensional systems. This method relies on a linear estimation of the Koopman generator, similarly to Kaiser (Kutz);  in the context of finite-dimensional (possibly stochastic) systems, but in contrast our proposed method does not require to evaluate time derivatives. In this sense, it is an indirect alternative method to recent direct methods for data-driven discovery of nonlinear PDEs (Rudy et al. (2017); Long et al. (2018); Li et al. (2019); Gurevich et al. (2019)).
The rest of the paper is organized as follows. In Section 2, we introduce the Koopman operator framework for infinite-dimensional systems, with a focus on the Lie generator and finite-dimensional approximation. Section 3 is devoted to spectral analysis and presents the generalized Extended Dynamic Mode Decomposition (EDMD) method. The identification method for infinite-dimensional systems is presented in Section 4 along with theoretical convergence results and numerical examples. Finally, concluding remarks are given in Section 5.

Koopman semigroup
We consider (infinite-dimensional) dynamical systems of the formu = W (u) , u ∈ U (1) where U is a separable Hilbert space and W : D(W ) → U is a nonlinear operator, with D(W ) the domain of W . If the system is described by a PDE, then W is typically a differential operator. Moreover, we assume that W generates a (possibly nonlinear) semiflow (ϕ t ) t≥0 : U → U, i.e. u(t) = ϕ t (u) is a classical solution to the abstract differential equation (1) associated with the initial condition u 0 ∈ D(W ). We also make the standing assumption that each map ϕ t : U → U, t ≥ 0, is continuous and that the mapping t → ϕ t (u) is continuous from R + into U (strong continuity).
The semigroup of Koopman operators (or Koopman semigroup in short) associated with (1) is defined on a space of "observable-functionals." Definition 1. (Koopman semigroup). Consider the space E of complex-valued functionals ζ : U → C, where the definition domain U ⊂ D(W ) is invariant under ϕ t . The semigroup of Koopman operators (K t ) t≥0 associated with the semiflow (ϕ t ) t≥0 is defined by K t ζ = ζ • ϕ t , ζ ∈ E.
In the rest of the paper, E is the space of bounded continuous functionals, endowed with the supremum norm ζ = sup u∈U |ζ(u)|, i.e. E = C(U).

Lie generator
Following the work by Dorroh and Neuberger (1996), we define the Lie generator of the Koopman semigroup. Definition 2. (Lie generator). The Lie generator of the semigroup (K t ) t≥0 is the linear operator L : D(L) → E that satisfies for all u ∈ U. Remark 1. (Infinitesimal generator). The Lie generator is reminiscent of the infinitesimal generator of strongly continuous semigroups (Engel and Nagel (1999)). However, the limit in (2) is defined pointwise, while it is defined in the strong sense in the case of the infinitesimal generator.
In fact, the infinitesimal generator of the Koopman semigroup cannot be defined unless the Koopman semigroup is strongly continuous (i.e. lim t↓0 K t ζ − ζ = 0 ∀ζ ∈ E). This property does not hold in our setting, but can be satisfied with the mixed topology on C(U), as shown in Farkas and Kreidler (2020).
The Lie generator enjoys a few properties (e.g. dense domain, bounded resolvent) and we refer to Dorroh and Neuberger (1996) for more details. In the case of a Koopman semigroup associated with a semiflow generated by the abstract differential equation (1), it follows from the chain rule property that the Lie generator is given by where D W (u) ζ(u) denotes the Gâteaux derivative of ζ at u in the direction W (u): Note that this can be interpreted as the Lie derivative associated with the infinite-dimensional vector field W (·).

Finite-dimensional representation
It is convenient to approximate the Koopman semigroup K t or the Lie generator L in a finite-dimensional subspace of E. Toward that end, we can consider the compressions K t n = P n K t | En and L n = P n L| En , where E n ⊂ D(L) is a ndimensional subspace of E and P n : E → E n is a projection operator. Suppose that E n is spanned by the basis of functionals {ζ i } n i=1 . In this basis, the finite-dimensional operators K t n and L n can be represented by the matrices K ∈ R n×n and L ∈ R n×n , respectively, which are defined so that The choice of the basis functions is crucial, as it affects the accuracy of the approximation and the performance of the methods based on this approximation (see below). However, finding the optimal set of basis functions is not trivial, and may require a priori knowledge on the system.

SPECTRAL ANALYSIS AND EXTENDED DYNAMIC MODE DECOMPOSITION
The spectral properties of the Koopman operator reveal important geometric properties of the underlying dynamics (see e.g. Mezić (2005), and Mauroy and Mezić (2012); Nakao and Mezić (2018) in the context of phase reduction).
In this section, we exploit the proposed framework for infinite-dimensional systems and compute the spectrum of the Koopman operator from data. This yields a generalization of the Extended Dynamical Mode Decomposition method for infinite-dimensional systems.
3.1.0.1. Case of linear systems. It is well-known that the spectrum of linear finite-dimensional systems is contained in the spectrum of the related Koopman operator. As shown in Mezić (2020) and Nakao and Mezić (2020), this result also holds for infinite-dimensional systems. Consider a linear systemu = Au, u ∈ U, and suppose that λ is an eigenvalue of A, so that there exists an eigenfunction w λ ∈ D(A * ) with A * w λ =λw λ , where A * denotes the adjoint operator of A andλ is the complex conjugate of λ.
Then the functional ζ λ (·) = ·, w λ satisfies Lξ λ (u) = D Au u, w λ = Au, w λ = u, A * w λ = u,λw λ = λ ξ λ (u) , so that λ is a Koopman eigenvalue. Moreover, it is easy to verify that (ζ λ ) α is a Koopman eigenfunctional associated with the Koopman eigenvalue αλ, provided that it belongs to the space E. If A generates a strongly continuous (linear) semigroup (ϕ t ) t≥0 = (T t ) t≥0 , then the spectral mapping theorem implies that e λt is an eigenvalue of the operator T t (see e.g. Engel and Nagel (1999), Chapter IV). It follows that and e λt is also an eigenvalue of K t .

Extended Dynamic Mode Decomposition for infinitedimensional systems
Extended Dynamic Mode Decomposition (EDMD) is a data-driven method that builds a finite-dimensional approximation of the Koopman semigroup and computes the approximate spectral properties of the operator (Williams et al. (2015)). It can be easily extended to the infinitedimensional framework that we consider here.
Suppose we have access to a set of m pairs (u k , ϕ ts (u k )) ∈ U 2 , where t s is a sampling time, in such a way that we can The generalized EDMD method proceeds with the following steps.
1. Compute the data matrices and 2. Provided that m ≥ n, a matrix approximation of K ts n = P n K ts | En is given by the least squares solution K = Θ † 1 Θ 2 , where Θ † denotes the Moore-Penrose pseudoinverse of Θ. Note that, in this case, P n is the discrete orthogonal projection 3. The eigenvalues λ K of K ts are approximated by the eigenvalues of K and estimates of the eigenvalues λ L of the Lie generator are given by λ L = log(λ K )/t s . Moreover, Koopman eigenfunctionals are approximated in the basis of functionals by the components of the corresponding (right) eigenvectors of K.
The method only requires to know the samples u k in a weak sense, i.e. through the values ζ j (u k ) of a finite number of functionals. For instance, these values could be the weighted averages of u k : X → R over the definition domain X. For the specific choice of evaluation functionals ζ j (u) = u(x j ) with x j ∈ X, we recover the classical DMD method (Tu et al. (2014)) applied to a discretized version of the infinite-dimensional system (1). Remark 2. The method is more general than the standard EDMD method in that it only requires to know the samples u k in a weak sense, i.e. through the values ζ j (u k ) of a finite number of functionals. For instance, these values could be the weighted averages of u k : X → R over the definition domain X. For the specific choice of evaluation functionals ζ j (u) = u(x j ) with the sample points x j ∈ X, we recover the classical DMD method (Tu et al. (2014)) applied to a spatially discretized version of the infinitedimensional system (1). Similarly, evaluation functionals of the form ζ j (u) = ψ(u(x j )), where ψ is a basis function, yield the EDMD method Williams et al. (2015) applied to the discretized infinite-dimensional system. Remark 3. (Convergence properties). The EDMD method must be used with some care since it might yield spurious eigenvalues. Similarly to Korda and Mezić (2018), convergence properties should be characterized as m, n → ∞ and the general validity of the spectral mapping theorem should be investigated. We leave these questions for future research.

Numerical example
We illustrate the generalized EDMD method with the Burgers equatioṅ associated with homogeneous Dirichlet boundary conditions u(−1) = u(1) = 0. The dynamics (7) is conjugated to the linear diffusion dynamicsv = ∂ 2 v/∂t 2 through the so-called Cole-Hopf transformation (Hopf (1950)). As explained in Nakao and Mezić (2020); Page and Kerswell (2018), this implies that the Koopman spectrum associated with the Burgers dynamics coincides with the Koopman spectrum associated with the linear diffusion dynamics. It follows that Koopman eigenfunctionals are related to eigenfunctions of the diffusion operator (see Section 3.1) and, in particular, they are of the form ζ(v) = sin(kπx/2), v(x) α (in the new variable v). The associated Koopman eigenvalues are given by λ = α(kπ/2) 2 .
We use the dynamics (7) to generate m = 50 datapairs taken from 10 trajectories (with random, arbitrarily chosen initial conditions of the form (x 2 −1) cos(aπx+bπ), a, b ∈ [0, 1]). The sampling time is t s = 0.2. Estimates of the Koopman eigenvalues are computed through our generalized EDMD method, with n = 27 basis functionals ζ j,k,l = cos(a j (πx/2) + b j π/2), (u(x)) k l , with (j, k, l) ∈ {1, 2, 3} 3 and where a j , b j are randomly chosen over the interval [0, 1]. As shown in Fig. 1, dominant eigenvalues −α(π/2) 2 , α ∈ N, are captured. This is consistent with the results presented in Page and Kerswell (2018) and shows the importance of selecting an augmented basis of nonlinear functionals to capture other eigenvalues than the principal ones of the form (kπ/2) 2 , k ∈ N. Note that other (complex) eigenvalues may appear over different numerical tests. They should be considered in light of future theoretical analysis (see Remark 3), and more advanced methods should be proposed to distinguish true eigenvalues from spurious ones.  (7). Dominant eigenvalues −α(π/2) 2 , α ∈ N (blue circles) are correctly estimated (red crosses).

IDENTIFICATION OF INFINITE-DIMENSIONAL SYSTEMS
In this section, we use the Koopman operator framework for infinite-dimensional systems in the context of identification. Our goal is to identify the coefficients c i ∈ R of the infinite-dimensional dynamicṡ with u ∈ U = L 2 (X) (where X ⊂ R p is a compact set), using m data pairs (u k , ϕ ts (u k )) generated by the dynamics (8). We assume that the operators W i : D(W i ) → U are known a priori, so that this can be seen as a parameter estimation problem. Note that (8) may be described by a partial differential equation, but the proposed method is not limited to that case (see Section 4.4).

Lifting identification method
The lifting identification method proposed in our previous work (Mauroy and Goncalves (2020)) is generalized to the case of infinite-dimensional systems. This method consists in three steps.
1. Lifting of the data. We compute the data matrices (4) and (5) with the basis of functionals where ·, · denotes the inner product in L 2 (X) and w ∈ L 2 (X) is a weighting function. We suppose without loss of generality that W 1 (u) = u (possibly with the coefficient c 1 = 0 in (8)), so that the linear functional ζ 1 (u) = u, w belongs to E n . 2. Identification of the Lie generator. We compute the matrix representation K = Θ † 1 Θ 2 of the compression K ts n in the subspace E n = span(ζ 1 , . . . , ζ n ) (step 2 in Section 3.2). Then we obtain a finite-dimensional approximation L (ts) of the Lie generator by taking the matrix logarithm Note that this approximation is not equal to the matrix representation L of the L n (see (3)). 3. Identification of the coefficients. Estimatesĉ i of the coefficients c i are given by the entries of the first column of L (ts) , i.e.ĉ i = L (ts) i1 . Remark 4. For the specific choice of basis functionals of the form ζ j (u) = ψ(u(x j )), with x j ∈ X and where ψ is a basis function, one recovers the original lifting identification method Mauroy and Goncalves (2020) applied to a spatially discretized version of the infinite-dimensional system. The proposed method is more general since it allows any basis functionals.
The lifting identification method does not require to compute time derivatives and is therefore an alternative to direct methods for PDE identification such as those proposed in Rudy et al. (2017) 2019), which makes a clever use of a weak formulation of the data in space and time. We note that this method requires sufficiently long time-series to allow accurate time integration, while our method can deal with data pairs belonging to different trajectories.

Convergence results
Now we show that the estimated coefficients converge to the true coefficients as t s → 0. This is summarized in the following proposition. Proposition 1. Let ϕ t be the continuous flow generated by (8) and let E n ⊂ D(L) be the subspace spanned by the basis functionals ζ i (·) = W i (·), w ∈ C(U), with ζ 1 (·) = ·, w . Assume that the pairs (u k , ϕ ts (u k )) ∈ U 2 are such that the n ≤ m vectors (ζ i (u 1 ), . . . , ζ i (u m )), i = 1, . . . , n, are linearly independent. Then, lim ts→0 L (ts) i1 = c i where L (ts) is given by (10).
Proof. The discrete orthogonal projection P n (6) is welldefined since the vectors (ζ i (u 1 ), . . . , ζ i (u m )), i = 1, . . . , n, are linearly independent. It is clear that and we also have Since P n ζ i = ζ i for all i = 1, . . . , n, it follows that For t s small enough, one has log A ts n = P n LP n , where A ts n = E n → E n is the finite-dimensional operator A ts n := e tsPnLPn , so that Since the basis functionals ζ i are linearly independent, lim ts→0 1 t s |(log A ts n − log K ts n )ζ 1 (u)| = 0 ∀u (13) implies that lim ts→0 L (ts) i1 − c i = 0 . Thus it remains to show that (13) holds.
Since lim ts→0 A ts n − I = 0, it is easy to check that lim [note that lim t→0 (log f (t)−(f (t)−1))/t = 0 for all f ∈ C 1 such that lim t→0 f (t) = 1] and it follows similarly that lim ts→0 1 t s log K ts n − (K ts n − I) = 0 .
Since the mapping t → (A t n − K t n )ζ(u) is differentiable, the mean value theorem implies that 1 for all ζ ∈ E n and u ∈ U, and for some τ ∈ [0, t s ]. Then we can write 1 t s (A ts n − K ts n )ζ(u) ))| where we used the fact that P n ζ = ζ and that K t and L commute. Since the functionals A τ n P n LP n ζ and A τ n P n Lζ are continuous and the flow ϕ t is continuous in t, we have lim ts→0 |A τ n P n LP n (ζ(ϕ ts−τ (u)) − ζ(u))| = 0 lim ts→0 |A τ n P n L(ζ(u) − ζ(ϕ ts−τ (u)))| = 0 so that lim Finally (14), (15) and (16) yield This concludes the proof.
Although the result requires that the sampling time t s tend to zero, accurate results can be obtained even if t s is not so small (see Section 4.4). However, large sampling times might affect the method performance. In particular, if the sampling time is too large, the principal branch of the matrix logarithm 1/t s log A ts n = 1/t s log e tsPnLPn is not equal to P n LP n . This is related to the so-called system aliasing issue Yue et al. (2016) (see also Mauroy and Goncalves (2020) for more details).
Corollary 1. Under the assumptions of Proposition 1, but with ζ i (·) = V i (·), w , V i : D(V i ) → U (and V 1 = I), the estimated nonlinear operatorŴ (ts) Proof. It follows from (11), (12) and (13) that The result follows from the definition of the discrete orthogonal projection (6). Remark 5. The proof of Proposition 1 is adapted from a proof in Mauroy and Goncalves (2020), but does not rely on semigroup theory. In particular, only weak (i.e. pointwise) convergence properties are used in the proof of Proposition 1, while strong convergence is considered in the proof in Mauroy and Goncalves (2020). For this reason, strong convergence results can be obtained in Mauroy and Goncalves (2020) in the case of arbitrary bases, but a weaker result is obtained here (see Corollary 1). In this context, considering a space of bounded continuous functionals equipped with a mixed topology would allow to exploit the strong continuity property of the semigroup, as shown in Farkas and Kreidler (2020), and possibly recover stronger convergence results.

Case of linearly dependent basis functionals
The basis functionals (9) might not be linearly independent, even if the operators W i are linearly independent (see Section 4.4.2). In this case, the result of Proposition 1 does not hold and in particular the estimated coefficientŝ c i L (ts) i1 do not approximate the exact coefficients c i . However, it follows from (12)-(13) that these coefficients satisfy the equality as t s goes to zero. Considering several weighting functions w (j) , with j = 1, . . . , J, we can use the proposed identification method with J different sets of basis functionals ζ Considering the above set of equations at u = u k , for k = 1, . . . , m, we obtain the matrix equality with is the data matrix (4) obtained with basis functionals ζ (j) i . Provided that J is large enough and the choice of weighting functions w (j) is appropriate, the matrix Θ f ull can be full rank so that (17) admits a unique solution Θ † f ull b. In this case, the coefficients c i are recovered despite the fact that every set of basis functionals is linearly dependent. Remark 6. Several weighting functions can also be used when the basis functionals are linearly independent (provided that the observations are available in practice). Then, the sets of estimated coefficientsĉ (j) i can be averaged to improve the accuracy of the results. Inconsistent results can also be discarded by comparing the different sets of coefficients.

Numerical examples
We can now use the lifting identification method with two illustrating examples: a nonlinear partial differential equation and a nonlinear diffusive dynamics on a graphon.
In particular, we can compute the basis functionals ζ i (·) = W i (·), w , with i = 5, . . . , 8 and we obtain with C 1 = 1 0 w(x)dx, C 2 = −1, C 3 = 1 0 x w(x)dx, and C 4 = −1/2. Then it is easy to see that the functionals ζ i , with i = 5, . . . , 8, are linearly dependent for any weight function w. We therefore follow the procedure described in Section 4.3, using the weighting functions w (j) (x) = x j , with j = 1, . . . , 4. As shown in Fig. 3, the coefficients are correctly estimated and, in particular, the graphon is identified.

Numerical performance
In this section, we compare the lifting identification method with a simple direct identification method. For this latter method, the time derivativeu k at u k is estimated through (forward) finite differencesu The estimated coefficientsĉ i are obtained through least squares regression of u, w over basis functionals of the form (9), i.e.   ĉ where Θ 1 is the data matrix (4).
We apply the two methods on data generated by the PDE (18), with the same setting and set of parameters as in Section 4.4.1. However several values of the sampling time t s are considered and two weighting functions w (1) (x) = exp(−1/(1 − (x/L) 2 )) and w (2) (x) = exp(−0.5/(1 − (x/L) 2 )) are also used (see Remark 6). For each method, we compute two sets of coefficientsĉ (1) i and c (2) i with the two weighting functions and take the mean i )/2. We finally compute the root mean square error (RMSE) which is shown in Figure 4 for both methods as a function of the sampling time. We can see that the lifting method outperforms the direct method and, in particular, is characterized by a small RMSE even for large sampling times. However, it is characterized by a higher variability.
We have also considered the effect of measurement noise.
To do so, we have considered the same setting as above (with t s = 0.5) and we have added to the data a Gaussian noise with zero mean and standard deviation σ · std(data), where std(data) stands for the standard deviation of the data. As shown in Figure 5, the lifting method still outperforms the direct method for low noise level (i.e. 0.001% and 0.01%), but is not robust to larger noise level in which case the direct method yields better results. The variability of the results is also higher with the lifting method. This can be explained by the fact that our proposed parameter estimation method is biased and not consistent, due to the lifting of (noisy) data. In future work, noise robustness of the method should be improved. For instance, the weak formulation of the method could be further exploited to enhance noise robustness (see e.g. Gurevich et al. (2019) i | > 1 for some i are discarded.

CONCLUSION
In this paper, the Koopman operator framework has been leveraged in the case of infinite-dimensional nonlinear systems. Building on previous theoretical works, we have considered a semigroup of composition operators and the associated Lie generator on a space of continuous functionals. A finite-dimensional representation of these operators has been proposed and used in the context of spectral analysis. This approach yields a generalization of the data-driven EDMD method for infinite-dimensional systems. We have also developed a novel identification method, which allows to identify nonlinear PDEs although it relies solely on linear techniques. This method has been complemented with convergence results.
The Koopman operator framework for infinite-dimensional systems is still in its infancy. Convergence properties of the finite-dimensional approximation of the Koopman operator should be thoroughly studied, in particular in light of the recent results by Farkas and Kreidler (2020) in semigroup theory. This could provide some insight into the results obtained with the generalized EDMD method. The possibility to consider several weights functions should be further investigated and exploited. In particular, some guidelines to carefully select the weight functions could be provided. Finally, robustness to noise should be improved.