Koopman Operator Framework for Spectral Analysis and Identification of Infinite-Dimensional Systems

Abstract

We consider the 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 the data, a method which can be seen as a generalization of the 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.

1. Introduction

The Koopman operator theory is a powerful framework which provides an alternative approach to dynamical systems. Through the so-called Koopman (or composition) operator [1], nonlinear dynamical systems are approximated by linear, higher-dimensional systems, that are amenable to systematic analysis. Under the inspiration of the seminal work by [2], 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 [3], 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 [4], 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, ref. [5] 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 recent work by [6], 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 for studying the spectral properties of nonlinear PDEs [7,8].

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 the strongly continuous semigroup theory, but rather adopt the approach proposed by [5]. 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 us 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 [8] 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 [9] to infinite-dimensional systems and is more general than the mere application of standard EDMD to spatially discretized PDEs. In fact, while the method was 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 [10] with the above framework, we propose a novel method for the nonlinear identification of infinite-dimensional systems. This method relies on a linear estimation of the Koopman generator, similar to [11,12,13] in the context of finite-dimensional (possibly stochastic) systems; in contrast, however, our proposed method does not require the evaluation of time derivatives. In this sense, it is an indirect alternative method for the data-driven discovery of nonlinear PDEs [14,15,16,17].

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 approximations. 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.

2. Koopman Operator Theory for Infinite-Dimensional Systems

2.1. Koopman Semigroup

We consider (infinite-dimensional) dynamical systems of the form:

$$\dot{u}=W\left(u\right),u\in \mathcal{U}$$

(1)

where $\mathcal{U}$ is a separable Hilbert space and $W:\mathcal{D}\left(W\right)\to \mathcal{U}$ is a nonlinear operator, with $\mathcal{D}\left(W\right)$ being 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 ${\left({\phi }^t\right)}_{t\ge 0}:\mathcal{U}\to \mathcal{U}$, i.e., $u\left(t\right)={\phi }^t\left(u\right)$, that is a classical solution to the abstract differential Equation (1) associated with the initial condition $u_0\in \mathcal{D}\left(W\right)$. We also make the standing assumption that each map ${\phi }^t:\mathcal{U}\to \mathcal{U}$, $t\ge 0$, is continuous and that the mapping $t\mapsto {\phi }^t\left(u\right)$ is continuous from ${\mathbb{R}}^{+}$ into $\mathcal{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 $\mathcal{E}$ of complex-valued functionals $\zeta :\overline{\mathcal{U}}\to \mathbb{C}$, where the definition domain $\overline{\mathcal{U}}\subset \mathcal{D}\left(W\right)$ is invariant under ${\phi }^t$. The semigroup of Koopman operators ${\left(K^t\right)}_{t\ge 0}$ associated with the semiflow ${\left({\phi }^t\right)}_{t\ge 0}$ is defined by $K^t\zeta =\zeta \circ {\phi }^t$, $\zeta \in \mathcal{E}$.

In the rest of the paper, $\mathcal{E}$ is the space of bounded continuous functionals, endowed with the supremum norm $\parallel \zeta \parallel ={sup}_{u\in \overline{\mathcal{U}}}\left|\zeta \left(u\right)\right|$, i.e., $\mathcal{E}=C\left(\overline{\mathcal{U}}\right)$.

2.2. Lie Generator

Following the work by [5], we define the Lie generator of the Koopman semigroup.

Definition 2

(Lie generator). The Lie generator of the semigroup ${\left(K^t\right)}_{t\ge 0}$ is the linear operator $L:\mathcal{D}\left(L\right)\to \mathcal{E}$ that satisfies:

$$L\zeta \left(u\right)=\underset{t\downarrow 0}{lim}\frac{K^t\zeta \left(u\right)-\zeta \left(u\right)}{t}\zeta \in \mathcal{D}\left(L\right)$$

(2)

for all $u\in \overline{\mathcal{U}}$.

Remark 1

(Infinitesimal generator). The Lie generator is reminiscent of the infinitesimal generator of strongly continuous semigroups [18]. 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\downarrow 0}\parallel K^t\zeta -\zeta \parallel =0$$\forall \zeta \in \mathcal{E}$). This property does not hold in our setting, but can be satisfied with the mixed topology on $C\left(\overline{\mathcal{U}}\right)$, as shown in [6].

The Lie generator enjoys a few properties (e.g., a dense domain and bounded resolvent) and we refer to [5] 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:

$$L\zeta \left(u\right)=\underset{t\downarrow 0}{lim}\frac{\zeta \left({\phi }^t\left(u\right)\right)-\zeta \left(u\right)}{t}=D_{D_t{\phi }^t\left(u\right)}\zeta \left(u\right)=D_{W\left(u\right)}\zeta \left(u\right)$$

where $D_{W\left(u\right)}\zeta \left(u\right)$ denotes the Gâteaux derivative of $\zeta$ at u in the direction $W\left(u\right)$:

$$D_{W\left(u\right)}\zeta \left(u\right)=\underset{\lambda \to 0}{lim}\frac{\zeta (u+\lambda W(u\left)\right)-\zeta \left(u\right)}{\lambda }.$$

Note that this can be interpreted as the Lie derivative associated with the infinite-dimensional vector field $W(·)$.

2.3. Finite-Dimensional Representation

It is convenient to approximate the Koopman semigroup $K^t$ or the Lie generator L in a finite-dimensional subspace of $\mathcal{E}$. Toward that end, we can consider the compressions $K_n^t=P_n{K}^t{|}_{{\mathcal{E}}_n}$ and $L_n=P_n{L|}_{{\mathcal{E}}_n}$, where ${\mathcal{E}}_n\subset \mathcal{D}\left(L\right)$ is a n-dimensional subspace of $\mathcal{E}$ and $P_n:\mathcal{E}\to {\mathcal{E}}_n$ is a projection operator. Suppose that ${\mathcal{E}}_n$ is spanned by the basis of functionals ${\left\{{\zeta }_i\right\}}_{i=1}^n$. In this basis, the finite-dimensional operators $K_n^t$ and $L_n$ can be represented by the matrices $\mathbf{K}\in {\mathbb{R}}^{n\times n}$ and $\mathbf{L}\in {\mathbb{R}}^{n\times n}$, respectively, which are defined so that:

$$K_n^t{\zeta }_j=\sum _{i=1}^n{\mathbf{K}}_{ij}{\zeta }_i\mathrm{and}{L}_n{\zeta }_j=\sum _{i=1}^n{\mathbf{L}}_{ij}{\zeta }_i.$$

(3)

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 of the system.

3. 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., [2], and [19,20], 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 the data. This yields a generalization of the Extended Dynamical Mode Decomposition method for infinite-dimensional systems.

3.1. Spectrum of the Koopman Operator

We consider the spectrum of the Lie generator (2), i.e., the set of (Koopman) eigenvalues $\lambda$ such that $L{\zeta }_{\lambda }=\lambda {\zeta }_{\lambda }$ for some (Koopman) eigenfunctional ${\zeta }_{\lambda }\in \mathcal{E}$.

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 [7,8], this result also holds for infinite-dimensional systems. Consider a linear system $\dot{u}=Au$, $u\in \mathcal{U}$, and suppose that $\lambda$ is an eigenvalue of A, so that there exists an eigenfunction $w_{\lambda }\in \mathcal{D}\left(A^{*}\right)$ with $A^{*}{w}_{\lambda }=\overline{\lambda }{w}_{\lambda }$, where $A^{*}$ denotes the adjoint operator of A and $\overline{\lambda }$ is the complex conjugate of $\lambda$. Then the functional ${\zeta }_{\lambda }(·)=\langle ·,w_{\lambda }\rangle$ satisfies:

$$L{\zeta }_{\lambda }\left(u\right)=D_{Au}\langle u,w_{\lambda }\rangle =\langle Au,w_{\lambda }\rangle =\langle u,A^{*}{w}_{\lambda }\rangle =\langle u,\overline{\lambda }{w}_{\lambda }\rangle =\lambda {\zeta }_{\lambda }\left(u\right),$$

so that $\lambda$ is a Koopman eigenvalue. Moreover, it is easy to verify that ${\left({\zeta }_{\lambda }\right)}^{\alpha }$ is a Koopman eigenfunction associated with the Koopman eigenvalue $\alpha \lambda$, provided that it belongs to the space $\mathcal{E}$. If A generates a strongly continuous (linear) semigroup ${\left({\phi }^t\right)}_{t\ge 0}={\left(T^t\right)}_{t\ge 0}$, then the spectral mapping theorem implies that $e^{\lambda t}$ is an eigenvalue of the operator $T^t$ (see, e.g., [18], Chapter IV). It follows that:

$$K^t{\zeta }_{\lambda }=\langle T^{t}u,w_{\lambda }\rangle =\langle u,{\left(T^t\right)}^{*}{w}_{\lambda }\rangle =\langle u,e^{\overline{\lambda }t}{w}_{\lambda }\rangle =e^{\lambda t}{\zeta }_{\lambda }$$

and $e^{\lambda t}$ is also an eigenvalue of $K^t$.

3.2. Extended Dynamic Mode Decomposition for Infinite-Dimensional 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 [9]. It can be easily extended to the infinite-dimensional framework that we consider here.

Suppose we have access to a set of m pairs $(u_k,{\phi }^{t_s}\left(u_k\right))\in {\overline{\mathcal{U}}}^2$, where $t_s$ is a sampling time, in such a way that we can measure the values of n functionals $({\zeta }_i\left(u_k\right),{\zeta }_i\left({\phi }^{t_s}\left(u_k\right)\right))\in {\mathbb{R}}^2$ for all $i\in \{1,\cdots ,n\}$. The generalized EDMD method proceeds with the following steps.

1.

Compute the data matrices:

$${\Theta }_1=\left(\begin{array}{ccc}{\zeta }_1\left(u_1\right)& \cdots & {\zeta }_n\left(u_1\right)\\ ⋮& & ⋮\\ {\zeta }_1\left(u_m\right)& \cdots & {\zeta }_n\left(u_m\right)\end{array}\right)$$

(4)

and:

$${\Theta }_2=\left(\begin{array}{ccc}{\zeta }_1\left({\phi }^{t_s}\left(u_1\right)\right)& \cdots & {\zeta }_n\left({\phi }^{t_s}\left(u_1\right)\right)\\ ⋮& & ⋮\\ {\zeta }_1\left({\phi }^{t_s}\left(u_m\right)\right)& \cdots & {\zeta }_n\left({\phi }^{t_s}\left(u_m\right)\right)\end{array}\right).$$

(5)

2.

Provided that $m\ge n$, a matrix approximation of $K_n^{t_s}=P_n{K}^{t_s}{|}_{{\mathcal{E}}_n}$ is given by the least squares solution $\mathbf{K}={\Theta }_1^{†}{\Theta }_2$, where ${\Theta }^{†}$ denotes the Moore-Penrose pseudoinverse of $\Theta$. Note that, in this case, $P_n$ is the discrete orthogonal projection:

$$P_n\zeta =\underset{\tilde{\zeta }\in {\mathcal{E}}_n}{\mathrm{argmin}}\sum _{k=1}^m{|\zeta \left(u_k\right)-\tilde{\zeta }\left(u_k\right)|}^2.$$

(6)

3.

The eigenvalues ${\lambda }_K$ of $K^{t_s}$ are approximated by the eigenvalues of $\mathbf{K}$ and estimates of the eigenvalues ${\lambda }_L$ of the Lie generator are given by ${\lambda }_L=log\left({\lambda }_K\right)/t_s$. Moreover, Koopman eigenfunctionals are approximated in the basis of functionals by the components of the corresponding (right) eigenvectors of $\mathbf{K}$.

Remark 2.

The method is more general than the standard EDMD method in that it only requires the samples $u_k$ be known in a weak sense, i.e., through the values ${\zeta }_j\left(u_k\right)$ of a finite number of functionals. For instance, these values could be the weighted averages of $u_k:X\to \mathbb{R}$ over the definition domain X. For the specific choice of evaluation functionals ${\zeta }_j\left(u\right)=u\left(x_j\right)$ with the sample points $x_j\in X$, we recover the classical DMD method [21] applied to a spatially discretized version of the infinite-dimensional system (1). Similarly, evaluation functionals of the form ${\zeta }_j\left(u\right)=\psi \left(u\left(x_j\right)\right)$, where ψ is a basis function, yield the EDMD method [9] 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. Similar to [22], convergence properties should be characterized as $m,n\to \infty$ and the general validity of the spectral mapping theorem should be investigated. We leave these questions for future research.

3.3. Numerical Example

We illustrate the generalized EDMD method with the Burgers equation:

$$\dot{u}=-u\frac{\partial u}{\partial x}+\frac{{\partial }^{2}u}{\partial x^2}u\in L^2\left([-1,1]\right)$$

(7)

associated with homogeneous Dirichlet boundary conditions $u(-1)=u\left(1\right)=0$. The dynamics (7) are conjugated to the linear diffusion dynamics $\dot{v}={\partial }^{2}v/\partial t^2$ through the so-called Cole-Hopf transformation [23]. As explained in [8,24], this implies that the Koopman spectrum, associated with Burgers dynamics, coincides with the Koopman spectrum associated with linear diffusion dynamics. Thus, it follows that Koopman eigenfunctions are related to the eigenfunctions of the diffusion operator (see Section 3.1) and, in particular, they are of the form $\zeta \left(v\right)={\langle sin(k\pi x/2),v\left(x\right)\rangle }^{\alpha }$ (in the new variable v). The associated Koopman eigenvalues are given by $\lambda =\alpha {(k\pi /2)}^2$.

We use the dynamics (7) to generate $m=50$ data-pairs taken from 10 trajectories (with random, arbitrarily chosen initial conditions of the form $(x^2-1)cos(a\pi x+b\pi )$, $a,b\in [0,1]$, which satisfy the boundary conditions). 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 ${\zeta }_{j,k,l}={\langle cos(a_j(\pi x/2)+b_j\pi /2),{\left(u\left(x\right)\right)}^k\rangle }^l$, with $(j,k,l)\in {\{1,2,3\}}^3$ and where $a_j,b_j$ are randomly chosen over the interval $[0,1]$. As shown in Figure 1, dominant eigenvalues $-\alpha {(\pi /2)}^2$, $\alpha \in \mathbb{N}$, are captured. This is consistent with the results presented in [24] and shows the importance of selecting an augmented basis of nonlinear functionals to capture eigenvalues other than the principal ones of the form ${(k\pi /2)}^2$, $k\in \mathbb{N}$. Note that other (complex) eigenvalues may appear over different numerical tests. They should be considered in light of future theoretical analyses (see Remark 3), and more advanced methods should be proposed to distinguish true eigenvalues from spurious ones.

4. 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\in \mathbb{R}$ of the infinite-dimensional dynamics,

$$\dot{u}=W\left(u\right)=\sum _i{c}_i{W}_i\left(u\right),$$

(8)

with $u\in \mathcal{U}=L^2\left(X\right)$ (where $X\subset {\mathbb{R}}^p$ is a compact set), using m data pairs $(u_k,{\phi }^{t_s}\left(u_k\right))$ generated by the dynamics (8). We assume that the operators $W_i:\mathcal{D}\left(W_i\right)\to \mathcal{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).

4.1. Lifting Identification Method

The lifting identification method proposed in our previous work [10] is generalized to the case of infinite-dimensional systems. This method consists of three steps.

1.

Lifting of the data. We compute the data matrices (4) and (5) with the basis of functionals:

$${\zeta }_i\left(u\right)=\langle W_i\left(u\right),w\rangle$$

(9)

where $\langle ·,·\rangle$ denotes the inner product in $L^2\left(X\right)$ and $w\in L^2\left(X\right)$ is a weighting function. We suppose without loss of generality that $W_1\left(u\right)=u$ (possibly with the coefficient $c_1=0$ in (8)), so that the linear functional ${\zeta }_1\left(u\right)=\langle u,w\rangle$ belongs to ${\mathcal{E}}_n$.

2.

Identification of the Lie generator. We compute the matrix representation $\mathbf{K}={\Theta }_1^{†}{\Theta }_2$ of the compression $K_n^{t_s}$ in the subspace ${\mathcal{E}}_n=\mathrm{span}({\zeta }_1,\cdots ,{\zeta }_n)$ (step 2 in Section 3.2). Then we obtain a finite-dimensional approximation ${\tilde{\mathbf{L}}}^{\left(t_s\right)}$ of the Lie generator by taking the matrix logarithm:

$${\tilde{\mathbf{L}}}^{\left(t_s\right)}=\frac{1}{t_s}log\left({\Theta }_1^{†}{\Theta }_2\right).$$

(10)

Note that this approximation is not equal to the matrix representation $\mathbf{L}$ of the $L_n$ (see (3)).

3.

Identification of the coefficients. Estimates ${\widehat{c}}_i$ of the coefficients $c_i$ are given by the entries of the first column of ${\tilde{\mathbf{L}}}^{\left(t_s\right)}$, i.e., ${\widehat{c}}_i={\tilde{\mathbf{L}}}_{i1}^{\left(t_s\right)}$.

Remark 4.

For the specific choice of basis functionals of the form ${\zeta }_j\left(u\right)=\psi \left(u\left(x_j\right)\right)$, with $x_j\in X$ and where ψ is a basis function, one recovers the original lifting identification method [10] 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 [14,15,16,17]. It is also noticeable that the use of basis functionals of the form (9) bears similarity to the method developed in [14], 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.

4.2. Convergence Results

Now we show that the estimated coefficients converge to the true coefficients as $t_s\to 0$. This is summarized in the following proposition.

Proposition 1.

Let ${\phi }^t$ be the continuous flow generated by (8) and let ${\mathcal{E}}_n\subset \mathcal{D}\left(L\right)$ be the subspace spanned by the basis functionals ${\zeta }_i(·)=\langle W_i(·),w\rangle \in C\left(\overline{\mathcal{U}}\right)$, with ${\zeta }_1(·)=\langle ·,w\rangle$. Assume that the pairs $(u_k,{\phi }^{t_s}\left(u_k\right))\in {\overline{\mathcal{U}}}^2$ are such that the $n\le m$ vectors $({\zeta }_i\left(u_1\right),\cdots ,{\zeta }_i\left(u_m\right))$, $i=1,\cdots ,n$, are linearly independent. Then,

$$\underset{t_s\to 0}{lim}{\tilde{\mathbf{L}}}_{i1}^{\left(t_s\right)}=c_i$$

where ${\tilde{\mathbf{L}}}^{\left(t_s\right)}$ is given by (10).

Proof. 

The discrete orthogonal projection $P_n$ (6) is well-defined since the vectors

$({\zeta }_i\left(u_1\right),\cdots ,{\zeta }_i\left(u_m\right))$, $i=1,\cdots ,n$, are linearly independent. It is clear that:

$$\frac{1}{t_s}log{K}_n^{t_s}{\zeta }_1\left(u\right)=\sum _{i=1}^n{\tilde{\mathbf{L}}}_{i1}^{\left(t_s\right)}{\zeta }_i\left(u\right)$$

and we also have:

$$L{\zeta }_1\left(u\right)=D_{W\left(u\right)}\langle u,w\rangle =\langle W\left(u\right),w\rangle =\sum _{i=1}^n{c}_i{\zeta }_i\left(u\right).$$

(11)

Since $P_n{\zeta }_i={\zeta }_i$ for all $i=1,\cdots ,n$, it follows that:

$$\left|\sum _{i=1}^n\left({\tilde{\mathbf{L}}}_{i1}^{\left(t_s\right)}-c_i\right){\zeta }_i\left(u\right)\right|=\left|\left(\frac{1}{t_s}log{K}_n^{t_s}-P_{n}L{P}_n\right){\zeta }_1\left(u\right)\right|.$$

For $t_s$ small enough, one has $log{A}_n^{t_s}=P_{n}L{P}_n$, where $A_n^{t_s}={\mathcal{E}}_n\to {\mathcal{E}}_n$ is the finite-dimensional operator $A_n^{t_s}:=e^{t_s{P}_{n}L{P}_n}$, so that:

$$\left|\sum _{i=1}^n\left({\tilde{\mathbf{L}}}_{i1}^{\left(t_s\right)}-c_i\right){\zeta }_i\left(u\right)\right|=\frac{1}{t_s}\left|\left(log{K}_n^{t_s}-log{A}_n^{t_s}\right){\zeta }_1\left(u\right)\right|.$$

(12)

Since the basis functionals ${\zeta }_i$ are linearly independent,

$$\underset{t_s\to 0}{lim}\frac{1}{t_s}\left|(log{A}_n^{t_s}-log{K}_n^{t_s}){\zeta }_1\left(u\right)\right|=0\forall u$$

(13)

implies that:

$$\underset{t_s\to 0}{lim}\left|{\tilde{\mathbf{L}}}_{i1}^{\left(t_s\right)}-c_i\right|=0.$$

Thus it remains to show that (13) holds. Since ${lim}_{t_s\to 0}\parallel A_n^{t_s}-I\parallel =0$, it is easy to check that:

$$\underset{t_s\to 0}{lim}\frac{1}{t_s}∥log{A}_n^{t_s}-(A_n^{t_s}-I)∥=0$$

(14)

[note that ${lim}_{t\to 0}(logf\left(t\right)-(f\left(t\right)-1))/t=0$ for all $f\in C^1$ such that ${lim}_{t\to 0}f\left(t\right)=1$] and it follows similarly that:

$$\underset{t_s\to 0}{lim}\frac{1}{t_s}∥log{K}_n^{t_s}-(K_n^{t_s}-I)∥=0.$$

(15)

Since the mapping $t\mapsto (A_n^t-K_n^t)\zeta \left(u\right)$ is differentiable, the mean value theorem implies that:

$$\frac{1}{t_s}\left|(A_n^{t_s}-K_n^{t_s})\zeta \left(u\right)\right|=\left|\frac{d}{dt}\left(A_n^t{K}^{t_s-t}\zeta \left(u\right)\right){|}_{t=\tau }\right|$$

for all $\zeta \in {\mathcal{E}}_n$ and $u\in \overline{\mathcal{U}}$, and for some $\tau \in [0,t_s]$. Then we can write:

$$\begin{array}{cc}\hfill \frac{1}{t_s}\left|(A_n^{t_s}-K_n^{t_s})\zeta \left(u\right)\right|& =\left|A_n^{\tau }(P_{n}L{P}_n-P_{n}L)K^{t_s-\tau }\zeta \left(u\right)\right|\hfill \\ \hfill & \le |A_n^{\tau }{P}_{n}L{P}_n(K^{t_s-\tau }\zeta \left(u\right)-\zeta \left(u\right))|+|A_n^{\tau }{P}_{n}L(P_n\zeta \left(u\right)-\zeta \left(u\right))|\hfill \\ \hfill & +|A_n^{\tau }{P}_{n}L(\zeta \left(u\right)-K^{t_s-\tau }\zeta \left(u\right))|\hfill \\ \hfill & \le |A_n^{\tau }{P}_{n}L{P}_n(\zeta \left({\phi }^{t_s-\tau }\left(u\right)\right)-\zeta \left(u\right))|+|A_n^{\tau }{P}_{n}L(\zeta \left(u\right)-\zeta \left({\phi }^{t_s-\tau }\left(u\right)\right))|\hfill \end{array}$$

where we used the fact that $P_n\zeta =\zeta$ and that $K^t$ and L commute. Since the functionals $A_n^{\tau }{P}_{n}L{P}_n\zeta$ and $A_n^{\tau }{P}_{n}L\zeta$ are continuous and the flow ${\phi }^t$ is continuous in t, we have:

$$\begin{array}{c}\hfill \underset{t_s\to 0}{lim}\left|A_n^{\tau }{P}_{n}L{P}_n(\zeta \left({\phi }^{t_s-\tau }\left(u\right)\right)-\zeta \left(u\right))\right|=0\\ \hfill \underset{t_s\to 0}{lim}\left|A_n^{\tau }{P}_{n}L(\zeta \left(u\right)-\zeta \left({\phi }^{t_s-\tau }\left(u\right)\right))\right|=0\end{array}$$

so that:

$$\underset{t_s\to 0}{lim}\frac{1}{t_s}\left|(A_n^{t_s}-K_n^{t_s})\zeta \left(u\right)\right|=0.$$

(16)

Finally (14)–(16) yield:

$$\begin{array}{cc}\hfill & \underset{t_s\to 0}{lim}\frac{1}{t_s}\left|(log{A}_n^{t_s}-log{K}_n^{t_s})\zeta \left(u\right)\right|\hfill \\ \hfill & \le \underset{t_s\to 0}{lim}\frac{1}{t_s}\left(|(log{A}_n^{t_s}-(A_n^{t_s}-I))\zeta \left(u\right)|\right.\left.+|(log{K}_n^{t_s}-(K_n^{t_s}-I))\zeta \left(u\right)|+|(A_n^{t_s}-K_n^{t_s})\zeta \left(u\right)|\right)\hfill \\ \hfill & =0.\hfill \end{array}$$

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}_n^{t_s}=1/t_{s}log{e}^{t_s{P}_{n}L{P}_n}$ is not equal to $P_{n}L{P}_n$. This is related to the so-called system aliasing issue [25] (see also [10] for more details).

When an arbitrary basis of functionals is chosen so that $\langle W(·),w\rangle \notin \mathrm{span}\{{\zeta }_1,\cdots ,{\zeta }_n\}$, we have the following corollary.

Corollary 1.

Under the assumptions of Proposition 1, but with ${\zeta }_i(·)=\langle V_i(·),w\rangle$, $V_i:\mathcal{D}\left(V_i\right)\to \mathcal{U}$ (and $V_1=I$), the estimated nonlinear operator ${\widehat{W}}^{\left(t_s\right)}={\sum }_{i=1}^n{\tilde{\mathbf{L}}}_{i1}^{\left(t_s\right)}{V}_i$ satisfies:

$$\underset{t_s\to 0}{lim}{\widehat{W}}^{\left(t_s\right)}=\underset{\tilde{W}\in \mathrm{span}(V_1,\cdots ,V_n)}{\mathrm{argmin}}\sum _{k=1}^m{\left|〈\tilde{W}\left(u_k\right)-W\left(u_k\right),w〉\right|}^2.$$

Proof. 

It follows from (11)–(13) that:

$$\underset{t_s\to 0}{lim}\sum _{i=1}^n{\tilde{\mathbf{L}}}_{i1}^{\left(t_s\right)}{\zeta }_i(·)=P_n\langle W(·),w\rangle ,$$

which implies that:

$$\underset{t_s\to 0}{lim}〈{\widehat{W}}^{\left(t_s\right)}(·),w〉=P_n\langle W(·),w\rangle .$$

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 [10], but does not rely on the 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 [10]. For this reason, strong convergence results can be obtained in [10] 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 the exploit of the strong continuity property of the semigroup, as shown in [6], and possibly recover stronger convergence results.

4.3. 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 coefficients ${\widehat{c}}_i\triangleq {\tilde{\mathbf{L}}}_{i1}^{\left(t_s\right)}$ do not approximate the exact coefficients $c_i$. However, it follows from (12) and (13) that these coefficients satisfy the equality:

$$\sum _{i=1}^N({\widehat{c}}_i-c_i){\zeta }_i\left(u\right)=0$$

as $t_s$ goes to zero. Considering several weighting functions $w^{\left(j\right)}$, with $j=1,\cdots ,J$, we can use the proposed identification method with J different sets of basis functionals ${\zeta }_i^{\left(j\right)}=\langle W_i(·),w^{\left(j\right)}\rangle$, yielding several sets of values ${\widehat{c}}_i^{\left(j\right)}$ that satisfy the equations:

$$\sum _{i=1}^N({\widehat{c}}_i^{\left(j\right)}-c_i){\zeta }_i^{\left(j\right)}\left(u\right)=0j=1,\cdots ,J.$$

Considering the above set of equations at $u=u_k$, for $k=1,\cdots ,m$, we obtain the matrix equality:

$${\Theta }_{full}\left(\begin{array}{c}{c}_1\\ ⋮\\ c_N\end{array}\right)=\mathbf{b}$$

(17)

with:

$${\Theta }_{full}=\left(\begin{array}{c}{\Theta }_1^{\left(1\right)}\\ ⋮\\ {\Theta }_1^{\left(J\right)}\end{array}\right),\mathbf{b}=\left(\begin{array}{c}{\Theta }_1^{\left(1\right)}\left(\begin{array}{c}{\widehat{c}}_1^{\left(1\right)}\\ ⋮\\ {\widehat{c}}_N^{\left(1\right)}\end{array}\right)\\ ⋮\\ {\Theta }_1^{\left(J\right)}\left(\begin{array}{c}{\widehat{c}}_1^{\left(J\right)}\\ ⋮\\ {\widehat{c}}_N^{\left(J\right)}\end{array}\right)\end{array}\right),$$

and where ${\Theta }_1^{\left(j\right)}$ is the data matrix (4) obtained with basis functionals ${\zeta }_i^{\left(j\right)}$. Provided that J is large enough and the choice of weighting functions $w^{\left(j\right)}$ is appropriate, the matrix ${\Theta }_{full}$ can be full rank so that (17) admits a unique solution ${\Theta }_{full}^{†}\mathbf{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 ${\widehat{c}}_i^{\left(j\right)}$ can be averaged to improve the accuracy of the results. Inconsistent results can also be discarded by comparing the different sets of coefficients.

4.4. 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.

4.4.1. Nonlinear Partial Differential Equation

We aim at identifying the coefficients of the dynamics:

$$\dot{u}=-2u-0.5(1+u)\frac{\partial u}{\partial x}+(1-0.2u)\frac{{\partial }^{2}u}{\partial x^2}+0.1\frac{{\partial }^{3}u}{\partial x^3},$$

(18)

$u\in L^2[0,5]$, with homogeneous Dirichlet boundary conditions $u\left(0\right)=u\left(5\right)=0$. The PDE is used to generate $m=50$ data pairs, taken from 25 trajectories (with random, arbitrarily chosen initial conditions of the form $x(x-5)cos(a\pi x/5+b\pi )$, $a,b\in [0,1]$, which satisfy the boundary conditions). The sampling time is $t_s=0.3$. The lifting identification method is used with $n=12$ basis functionals (9), with the nonlinear operators:

$$W_i\left(u\right)\in \left\{u^j\frac{{\partial }^{k}u}{\partial x^k},j\in \{0,1,2\},k\in \{0,1,2,3\}\right\}$$

and the weighting function $w\left(x\right)=exp(-1/(1-{(x/L)}^2))$. Figure 2 shows that the coefficients $c_i$ estimates had low error.

4.4.2. Nonlinear Diffusive Dynamics on a Graphon

The graphon $G:{[0,1]}^2\to [0,1]$ is used to describe the limit of a sequence of dense graphs, and can be interpreted as the infinite-dimensional version of an adjacency matrix. Here, we consider a nonlinear diffusive dynamics on the graphon $G(x,y)=1-0.4x-0.1y-0.2xy-0.3y^2$, which is described by the integro-differential equation:

$$\dot{u}(t,x)=-0.5u(t,x)+1.5u{(t,x)}^2-u{(t,x)}^3+{\int }_0^{1}G(x,y)(u(t,y)-u(t,x))dy.$$

(19)

We generate $m=60$ data pairs, taken from 30 trajectories (with random initial conditions of the form $0.1acos(b\pi x+b\pi )$, $a,b\in [0,1]$). The sampling time is $t_s=0.5$. The lifting identification method is used with $n=10$ basis functionals (9), with the operators:

$$\begin{array}{cc}\hfill & W_1\left(u\right)=1,W_2\left(u\right)=u,W_3\left(u\right)=u^2,W_4\left(u\right)=u^3,\hfill \\ \hfill & W_5\left(u\right)\left(x\right)={\int }_0^{1}u\left(y\right)-u\left(x\right)dy,W_6\left(u\right)\left(x\right)={\int }_0^{1}x(u\left(y\right)-u\left(x\right))dy,\hfill \\ \hfill & W_7\left(u\right)\left(x\right)={\int }_0^{1}y(u\left(y\right)-u\left(x\right))dy,W_8\left(u\right)\left(x\right)={\int }_0^{1}xy(u\left(y\right)-u\left(x\right))dy,\hfill \\ \hfill & W_9\left(u\right)\left(x\right)={\int }_0^1{x}^2(u\left(y\right)-u\left(x\right))dy,W_{10}\left(u\right)\left(x\right)={\int }_0^1{y}^2(u\left(y\right)-u\left(x\right))dy.\hfill \end{array}$$

In particular, we can compute the basis functionals ${\zeta }_i(·)=\langle W_i(·),w\rangle$, with $i=5,\cdots ,8$ and we obtain:

$$\begin{array}{cc}\hfill & {\zeta }_5\left(u\right)={\int }_0^1{\int }_0^{1}w\left(x\right)(u\left(y\right)-u\left(x\right))dxdy=C_1\langle 1,u\left(x\right)\rangle +C_2\langle w\left(x\right),u\left(x\right)\rangle \hfill \\ \hfill & {\zeta }_6\left(u\right)={\int }_0^1{\int }_0^{1}w\left(x\right)x(u\left(y\right)-u\left(x\right))dxdy=C_3\langle 1,u\left(x\right)\rangle +C_2\langle xw\left(x\right),u\left(x\right)\rangle \hfill \\ \hfill & {\zeta }_7\left(u\right)={\int }_0^1{\int }_0^{1}w\left(x\right)y(u\left(y\right)-u\left(x\right))dxdy=C_1\langle x,u\left(x\right)\rangle +C_4\langle w\left(x\right),u\left(x\right)\rangle \hfill \\ \hfill & {\zeta }_8\left(u\right)={\int }_0^1{\int }_0^{1}w\left(x\right)xy(u\left(y\right)-u\left(x\right))dxdy=C_3\langle x,u\left(x\right)\rangle +C_4\langle xw\left(x\right),u\left(x\right)\rangle \hfill \end{array}$$

with $C_1={\int }_0^{1}w\left(x\right)dx$, $C_2=-1$, $C_3={\int }_0^{1}xw\left(x\right)dx$, and $C_4=-1/2$. Then it is easy to see that the functionals ${\zeta }_i$, with $i=5,\cdots ,8$, are linearly dependent for any weight function w. We therefore follow the procedure described in Section 4.3, using the weighting functions $w^{\left(j\right)}\left(x\right)=x^j$, with $j=1,\cdots ,4$. As shown in Figure 3, the coefficients are correctly estimated and, in particular, the graphon is identified.

4.4.3. Numerical Performance

In this section, we compare the lifting identification method with a simple direct identification method. For this latter method, the time derivative ${\dot{u}}_k$ at $u_k$ is estimated through (forward) finite differences:

$${\dot{u}}_k=\frac{{\phi }^{t_s}\left(u_k\right)-u_k}{t_s}$$

or equivalently:

$$\langle {\dot{u}}_k,w\rangle =\frac{\langle {\phi }^{t_s}\left(u_k\right),w\rangle -\langle u_k,w\rangle }{t_s}=\frac{{\zeta }_1\left({\phi }^{t_s}\left(u_k\right)\right)-{\zeta }_1\left(u_k\right)}{t_s}.$$

The estimated coefficients ${\widehat{c}}_i$ are obtained through least squares regression of $\langle \dot{u},w\rangle$ over basis functionals of the form (9), i.e.,

$$\left(\begin{array}{c}{\widehat{c}}_1\\ ⋮\\ {\widehat{c}}_n\end{array}\right)={\Theta }_1^{†}\left(\begin{array}{c}\langle {\dot{u}}_1,w\rangle \\ ⋮\\ \langle {\dot{u}}_m,w\rangle \end{array}\right)$$

where ${\Theta }_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^{\left(1\right)}\left(x\right)=exp(-1/(1-{(x/L)}^2))$ and $w^{\left(2\right)}\left(x\right)=exp(-0.5/(1-{(x/L)}^2))$ are also used (see Remark 6). For each method, we compute two sets of coefficients ${\widehat{c}}_i^{\left(1\right)}$ and ${\widehat{c}}_i^{\left(2\right)}$ with the two weighting functions and take the mean ${\widehat{c}}_i=({\widehat{c}}_i^{\left(1\right)}+{\widehat{c}}_i^{\left(2\right)})/2$. We finally compute the root mean square error (RMSE):

$$RMSE=\frac{\sqrt{{\sum }_{i=1}^n{(c_i-{\widehat{c}}_i)}^2}}{\sqrt{n}},$$

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 added Gaussian noise to the data with mean at zero and standard deviation $\sigma ·\mathrm{std}\left(\mathrm{data}\right)$, where $\mathrm{std}\left(\mathrm{data}\right)$ stands for the standard deviation of the data. As shown in Figure 5, the lifting method still outperforms the direct method for low noise levels (i.e., $0.001\%$ and $0.01\%$) but is not robust to larger noise levels, 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, the 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., [14]).

5. Conclusions

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 analyses. 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 [6] in semigroup theory. This could provide some insight into the results obtained with the generalized EDMD method. The possibility to consider several weight 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.