Subspace dynamic mode decomposition for stochastic Koopman analysis

TL;DR

This paper introduces subspace DMD, a novel algorithm for Koopman spectral analysis of stochastic systems that effectively handles observation noise, providing accurate spectral estimates for random dynamical systems.

Contribution

The paper proposes subspace DMD, a new method that improves Koopman analysis of noisy stochastic systems by ensuring convergence to the true spectra under standard assumptions.

Findings

Subspace DMD accurately estimates spectra in noisy stochastic systems.

The method converges to the true Koopman spectra under standard conditions.

Empirical tests demonstrate its utility across various dynamical systems.

Abstract

The analysis of nonlinear dynamical systems based on the Koopman operator is attracting attention in various applications. Dynamic mode decomposition (DMD) is a data-driven algorithm for Koopman spectral analysis, and several variants with a wide range of applications have been proposed. However, popular implementations of DMD suffer from observation noise on random dynamical systems and generate an inaccurate estimation of the spectra of the stochastic Koopman operator. In this paper, we propose subspace DMD as an algorithm for the Koopman analysis of random dynamical systems with observation noise. Subspace DMD first computes the orthogonal projection of future snapshots to the space of past snapshots and then estimates the spectra of a linear model, and its output converges to the spectra of the stochastic Koopman operator under standard assumptions. We investigate the empirical performance of subspace DMD with several dynamical systems and show its utility for the Koopman analysis of random dynamical systems.

Figures

Equations

$${\mathbf{x}}_{t+1}=\mathbf{f}({\mathbf{x}}_t),\mathbf{x}\in \mathcal{M}$$

$$\mathcal{K}g(\mathbf{x}):=g(\mathbf{f}(\mathbf{x})).$$

$${\left[\begin{array}{ccc}K{g}_1& \cdots & K{g}_n\end{array}\right]}^{\mathsf{T}}=\mathbf{K}\mathbf{g},$$

$$\phi (\mathbf{x})={\mathbf{z}}^{\mathsf{H}}\mathbf{g}(\mathbf{x}),$$

$$K\left({\mathbf{z}}^{\mathsf{H}}\mathbf{g}(\mathbf{x})\right)={\mathbf{z}}^{\mathsf{H}}\mathbf{K}\mathbf{g}(\mathbf{x})=\lambda {\mathbf{z}}^{\mathsf{H}}\mathbf{g}(\mathbf{x}).$$

$$\mathbf{g}(\mathbf{x})=\sum _{i=1}^n{\mathbf{z}}_i^{\mathsf{H}}\mathbf{g}(\mathbf{x}){\mathbf{w}}_i=\sum _{i=1}^n{\phi }_i(\mathbf{x}){\mathbf{w}}_i,$$

$$\mathbf{g}({\mathbf{x}}_t)=\sum _{i=1}^n{\lambda }_i^t{\mathbf{c}}_i,{\mathbf{c}}_i={\phi }_i({\mathbf{x}}_0){\mathbf{w}}_i,$$

$$\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}=\mathbf{f}(\mathbf{x}),\mathbf{x}\in \mathcal{M}.$$

$${\mathcal{K}}_c^{t}g(\mathbf{x}):=g(\mathbf{\varphi }(\mathbf{x},t)),$$

$${\mathcal{K}}_c=\underset{t\to 0}{\lim }\frac{{\mathcal{K}}_c^{t}g-g}{t}.$$

$${\mathbf{Y}}_0$$

$${\mathbf{Y}}_1$$

$$\underset{m\to \infty }{\lim }\frac{1}{m}\sum _{j=0}^{m-1}\varphi ({\mathbf{x}}_j)={\mathbb{E}}_{\mathcal{M}}[\varphi (\mathbf{x})]={\int }_{\mathcal{M}}\varphi (\mathbf{x})d{\mu }_{\mathcal{M}},$$

$$\mathbf{K}{\mathbf{G}}_0={\mathbf{G}}_1,{\mathbf{G}}_0={\mathbb{E}}_{\mathcal{M}}\left[\mathbf{g}{\mathbf{g}}^{\mathsf{H}}\right],{\mathbf{G}}_1={\mathbb{E}}_{\mathcal{M}}\left[(\mathbf{g}\circ \mathbf{f}){\mathbf{g}}^{\mathsf{H}}\right],$$

$$\mathbf{A}={\hat{\mathbf{G}}}_1{\hat{\mathbf{G}}}_0^{†},{\hat{\mathbf{G}}}_0=\frac{1}{m}{\mathbf{Y}}_0{\mathbf{Y}}_0^{\mathsf{H}},{\hat{\mathbf{G}}}_1=\frac{1}{m}{\mathbf{Y}}_1{\mathbf{Y}}_0^{\mathsf{H}},$$

$${\mathbf{x}}_{t+1}={\mathbf{f}}_{\Omega }({\mathbf{x}}_t,{\omega }_t),\mathbf{x}\in \mathcal{M},\omega \in \Omega$$

$${\mathcal{K}}_{\Omega }g(\mathbf{x}):={\mathbb{E}}_{\Omega }\left[g({\mathbf{f}}_{\Omega }(\mathbf{x},\omega ))\right],$$

$$\mathbf{g}({\mathbf{x}}_{t+1})={\mathbf{K}}_{\Omega }\mathbf{g}({\mathbf{x}}_t)+{\mathbf{e}}_t,$$

$${\mathbf{e}}_t:=\mathbf{g}(\mathbf{f}({\mathbf{x}}_t,{\omega }_t))-{\mathbb{E}}_{\Omega }\left[\mathbf{g}(\mathbf{f}({\mathbf{x}}_t,{\omega }_t))\right].$$

$$\mathbf{g}({\mathbf{x}}_t)={\mathbf{K}}_{\Omega }^t\mathbf{g}({\mathbf{x}}_0)+\sum _{k=0}^{t-1}{\mathbf{K}}_{\Omega }^{t-k-1}{\mathbf{e}}_k.$$

$${\mathbb{E}}_{\Omega }\left[{\mathbf{e}}_{t^{\prime }}{\mathbf{e}}_t^{\mathsf{H}}\right]=\mathbf{P}{\delta }_{t^{\prime }t}$$

$${\hat{\mathbf{G}}}_0$$

$${\hat{\mathbf{G}}}_1$$

$${\mathbf{G}}_0$$

$${\mathbf{G}}_1$$

$$={\int }_{\mathcal{M}\times \Omega }\mathbf{g}({\mathbf{f}}_{\Omega }(\mathbf{x},\omega ))\mathbf{g}(\mathbf{x}{)}^{\mathsf{H}}d{\mu }_{\mathcal{M}}d{\mu }_{\Omega }$$

$$\mathbf{h}({\mathbf{x}}_t,s_t):=\mathbf{g}({\mathbf{x}}_t)+\mathbf{w}(s_t),\mathbf{x}\in \mathcal{M},s\in S,$$

$$\begin{array}{c}{\mathbb{E}}_S\left[{\mathbf{w}}_t\right]=0,{\mathbb{E}}_S\left[{\mathbf{w}}_{t^{\prime }}{\mathbf{w}}_t^{\mathsf{H}}\right]=\mathbf{Q}{\delta }_{t^{\prime }t},\\ {\mathbb{E}}_{\Omega ,S}\left[{\mathbf{e}}_{t^{\prime }}{\mathbf{w}}_t^{\mathsf{H}}\right]=\mathbf{R}{\delta }_{t^{\prime }t},\end{array}$$

$${\mathbb{E}}_{\Omega }\left[\mathbf{g}({\mathbf{x}}_t)\right]$$

$${\mathbb{E}}_{\Omega }\left[\mathbf{g}({\mathbf{x}}_{t^{\prime }})\mathbf{g}({\mathbf{x}}_t{)}^{\mathsf{H}}\right]$$

$${\mathbb{E}}_{\mathcal{M}}\left[{\mathbf{G}}_{t,t}\right]$$

$${\mathbf{G}}_{t^{\prime },t}=\{\begin{array}{ll}{\mathbf{K}}_{\Omega }^{t^{\prime }-t}{\mathbf{G}}_{t,t},& t^{\prime }\ge t,\\ {\mathbf{G}}_{t^{\prime },t^{\prime }}({\mathbf{K}}_{\Omega }^{t-t^{\prime }}{)}^{\mathsf{H}},& t^{\prime }<t.\end{array}$$

$${\mathbf{G}}_{\tau }={\mathbf{K}}_{\Omega }^{\tau }\mathbf{G}.$$

Full text

I Introduction

Operator-theoretic approaches for the analysis of dynamical systems, which rely on the Perron–Frobenius operator Lasota and Mackey (1994) or its adjoint, the Koopman operator Koopman (1931), are attracting attention for use in mathematical and engineering applications. The Koopman operator is an infinite-dimensional linear operator that acts on a space of observation functions (observables), and the analysis based on it has been intensively studied recently (see, e.g., Mezić (2005); Budišić et al. (2012); Mezić (2013)). Several methods for conducting the Koopman spectral analysis have been proposed, such as the generalized Laplace analysis Budišić et al. (2012) and the Ulam–Galerkin method Froyland et al. (2013).

Dynamic mode decomposition (DMD) Rowley et al. (2009); Schmid (2010) is a data-driven method that can be utilized for Koopman spectral analysis. It has been applied to a wide range of scientific and engineering subjects including fluid mechanics Schmid et al. (2011), power system analysis Susuki and Mezić (2014), medical care Bourantas et al. (2014), epidemiology Proctor and Eckhoff (2015), robotic control Berger et al. (2015), neuroscience Brunton et al. (2016a), image processing Kutz et al. (2016a), nonlinear system identification Mauroy and Goncalves (2016), finance Mann and Kutz (2016), and chaotic systems Brunton et al. (2017). DMD computes a set of modes along with the corresponding frequencies and decay rates, given a sequence of measurements from the target dynamics. Those modes coincide with the ones obtained by the Koopman spectral analysis under certain conditions, which we briefly review in Section II.

In practice, popular implementations of DMD (e.g., Schmid (2010); Tu et al. (2014)) suffer from observation noise. Several researchers have addressed this issue; Duke *et al. * Duke et al. (2012) and Pan *et al. * Pan et al. (2015) conducted error analyses on the DMD algorithms, and Dawson *et al. * Dawson et al. (2016) and Hemati *et al. * Hemati et al. (2017) proposed reformulating DMD as a total-least-squares problem to treat the observation noise explicitly. Moreover, there is a line of research on the low-rank approximation of dynamics, including optimized DMD Chen et al. (2012), optimal mode decomposition Wynn et al. (2013), sparsity-promoting DMD Jovanović et al. (2014), and the closed-form solution for a low-rank constrained problem Héas and Herzet (2017). In addition, Takeishi *et al. * Takeishi et al. (2017) suggested a Bayesian formulation of DMD to incorporate uncertainties. Those studies provide clear perspectives on the treatment of the observation noise. Note that, however, they focus on deterministic dynamical systems, i.e., they do not explicitly deal with process noise, which limits their applicability to situations where the underlying dynamics contain random effects.

In fact, the Koopman analysis can also be applied to dynamical systems with process noise via the stochastic Koopman operator Mezić (2005). The spectra of the stochastic Koopman operator may convey information on the process noise; Bagheri *et al. * Bagheri (2014) investigated the effects of weak noise on the spectra of the Koopman operator for oscillating flows. The DMD algorithms are applicable even to stochastic systems Williams et al. (2015), unless observation noise is present. However, the existing variants of DMD do not explicitly consider both observation and process noise, and, in fact, most of them cannot compute the spectra of the stochastic Koopman operator accurately from noisy observations, which is partly demonstrated in Section IV using numerical examples.

In this paper, we present an algorithm based on the stochastic Koopman operator for decomposing nonlinear random dynamical systems from noisy observations. The proposed algorithm is referred to as subspace DMD because it has a strong connection to the subspace system identification methods developed in control theory. Subspace DMD is aware of both the observation noise and process noise at the same time, and we show its validity with numerical examples.

The remainder of this paper is organized as follows. In Section II, we introduce the fundamental concepts required to understand the purpose and procedures of Koopman analysis and DMD. Section III includes the main results of this paper, the algorithm of subspace DMD. In Section IV, we introduce numerical examples to show the empirical performance of subspace DMD. In Section V, we mention the important elements of DMD that are not fully addressed in this paper. This paper ends with the conclusions in Section VI.

II Background

We briefly review the method for decomposing nonlinear dynamical systems based on the spectra and the invariant subspace of the Koopman operator, along with a corresponding numerical procedure. Regarding the theory of Koopman spectral analysis, readers can consult studies such as Mezić (2005); Budišić et al. (2012); Mezić (2013) for more details.

II.1 Koopman spectral analysis on Koopman invariant subspace

Consider a discrete-time dynamical system

[TABLE]

with a map $\bm{f}:\mathcal{M}\to\mathcal{M}$ and time index $t\in\mathbb{T}=\{0\}\cup\mathbb{N}$, where $(\mathcal{M},\Sigma_{\mathcal{M}},\mu_{\mathcal{M}})$ is a probability space associated with a phase space $\mathcal{M}$. Instead of trajectories in the phase space, we analyze the evolution of an observable $g:\mathcal{M}\to\mathbb{C}$ in a function space $\mathcal{G}\subset L^{2}(\mathcal{M},\mu_{\mathcal{M}})$. Koopman operator $\mathcal{K}:\mathcal{G}\to\mathcal{G}$ is an infinite-dimensional linear operator defined as

[TABLE]

Here, suppose that there exists a subspace $G\subset\mathcal{G}$ that is invariant to $\mathcal{K}$, i.e., $\exists G~{}\text{s.t.}~{}\mathcal{K}g\in G,\forall g\in G$ 111The theory of decomposition based on the spectra of the Koopman operator does not necessarily require the existence of the invariant subspace. We introduce the notion of the invariant subspace here merely for clarifying the connection between the Koopman spectral analysis and dynamic mode decomposition.. Let us consider the restriction of $\mathcal{K}$ to $G$ and denote it by $K$. If $G$ is finite-dimensional, then $K$ also becomes a finite-dimensional linear operator. Suppose we have a set of observables $\{g_{1},\dots,g_{n}\}$ ($n<\infty$) that spans $G$, and let $\bm{K}$ be the representation of $K$ with regard to $\{g_{1},\dots,g_{n}\}$, i.e.,

[TABLE]

where $\bm{g}=\begin{bmatrix}g_{1}&\cdots&g_{n}\end{bmatrix}^{\mathsf{T}}$. Now let $\varphi$ be the eigenfunction of $K$ corresponding to an eigenvalue $\lambda$. Then, $\varphi$ with regard to $\{g_{1},\dots,g_{n}\}$ is expressed as

[TABLE]

where $\bm{z}$ is the left-eigenvector of $\bm{K}$ corresponding to eigenvalue $\lambda$, since

[TABLE]

Let $\bm{w}_{i}$ and $\bm{z}_{i}$ respectively be the right- and the left-eigenvector of $\bm{K}$ corresponding to an eigenvalue $\lambda_{i}$ for $i=1,\dots,n$. In the sequel, without loss of generality, we assume that $\bm{w}$ and $\bm{z}$ are normalized so that $\bm{w}_{i^{\prime}}^{\mathsf{H}}\bm{z}_{i}=\delta_{i^{\prime}i}$ ($\delta_{i^{\prime}i}$ is [math] if $i^{\prime}=i$ and $1$ otherwise). We assume that all the eigenvalues of $K$ are distinct, i.e., their multiplicities are one. Then, any values of $\bm{g}$ are expressed as

[TABLE]

where $\varphi_{i}$ is the eigenfunction of $K$ corresponding to eigenvalue $\lambda_{i}$. Applying $K$ on both sides of Eq. (5) repeatedly starting at $\bm{x}=\bm{x}_{0}$, we obtain the modal decomposition of the values of the observables, i.e.,

[TABLE]

where coefficient $\bm{c}_{i}$ (or $\bm{w}_{i}$) is referred to as a Koopman mode.

The same sort of discussion is also possible for a continuous-time dynamical system

[TABLE]

Instead of the Koopman operator, the Koopman semigroup $\{\mathcal{K}_{c}^{t}\}_{t\in\mathbb{R}^{+}}$ on this dynamical system is defined as

[TABLE]

where $\bm{\phi}(\bm{x},t)$ is the flow map that takes $\bm{x}$ as the initial state and returns the state after a time interval of length $t$. The infinitesimal generator of the Koopman semigroup is given as

[TABLE]

Let $\lambda_{c}$ be an eigenvalue of $\mathcal{K}_{c}$. This “continuous-time” eigenvalue can be computed from the corresponding “discrete-time” eigenvalue $\lambda$ by $\lambda_{c}=\ln(\lambda)/\Delta t$, where $\Delta t$ is the temporal interval in discrete-time dynamical system (1).

II.2 Dynamic mode decomposition

DMD Rowley et al. (2009); Schmid (2010); Kutz et al. (2016b) is a decomposition method for numerical datasets whose output converges to the modal decomposition via the Koopman operator under some conditions. Suppose that $\{g_{1},\dots,g_{n}\}$ spans $G$ and that we have data matrices generated with $\bm{g}=\begin{bmatrix}g_{1}&\cdots&g_{n}\end{bmatrix}^{\mathsf{T}}$:

[TABLE]

The popular algorithm of DMD Schmid (2010); Tu et al. (2014) leverages a compact singular value decomposition (SVD) to avoid a direct eigendecomposition of a large matrix, and its procedure is shown in Algorithm 1.

Algorithm 1** (DMD Schmid (2010); Tu et al. (2014)).**

Build a pair of data matrices $(\bm{Y}_{0},\bm{Y}_{1})$ as in Eq. (10). 2. 2.

Compute the compact SVD as $\bm{Y}_{0}=\bm{U}_{r}\bm{S}_{r}\bm{V}_{r}^{\mathsf{H}}$ with $\bm{U}_{r}\in\mathbb{C}^{n\times r}$, $\bm{S}_{r}\in\mathbb{C}^{r\times r}$ and $\bm{V}_{r}\in\mathbb{C}^{m\times r}$, where $r=\operatorname{rank}(\bm{Y}_{0})$. 3. 3.

Define matrix $\tilde{\bm{A}}=\bm{U}_{r}^{\mathsf{H}}\bm{Y}_{1}\bm{V}_{r}\bm{S}_{r}^{-1}$. 4. 4.

Compute the eigenvalues $\lambda$ and eigenvectors $\tilde{\bm{w}}$ of $\tilde{\bm{A}}$. 5. 5.

Return dynamic modes $\bm{w}=\lambda^{-1}\bm{Y}_{1}\bm{V}_{r}\bm{S}_{r}^{-1}\tilde{\bm{w}}$ and the corresponding eigenvalues $\lambda$.

The convergence of an output of Algorithm 1 in the large sample limit can be shown with the assumption of ergodicity as follows.

Assumption 1**.**

The time average of a measurable function $\phi:\mathcal{M}\to\mathbb{C}$ converges to its space average, i.e.,

[TABLE]

for almost all $\bm{x}_{0}\in\mathcal{M}$.

Proposition 1**.**

Suppose Assumption 1 holds. If all the modes are sufficiently excited in the data (i.e., $r=\dim(G)$) and all the nonzero eigenvalues of $\bm{A}=\bm{Y}_{1}\bm{Y}_{0}^{\dagger}$ are distinct, then the dynamic modes calculated by Algorithm 1 converge to the eigenvectors corresponding to the non-zero eigenvalues of $\bm{K}$ in $m\to\infty$ with probability one.

Proof.

Taking the inner product of both sides of Eq. (3) with $\bm{g}$, we have

[TABLE]

and thus the minimum-norm solution for $\bm{K}$ is given as $\bm{K}=\bm{G}_{1}\bm{G}_{0}^{\dagger}$, where $\bm{G}_{0}^{\dagger}$ is the Moore–Penrose pseudoinverse of $\bm{G}_{0}$. In contrast, from the definition of $\bm{A}$,

[TABLE]

and by Assumption 1, empirical matrices $\hat{\bm{G}}_{0}$ and $\hat{\bm{G}}_{1}$ converge to $\bm{G}_{0}$ and $\bm{G}_{1}$, respectively, in $m\to\infty$ with probability one. Moreover, because $r=\dim(G)$, the rank of $\hat{\bm{G}}_{0}$ is always $\dim(G)$ and thus $\hat{\bm{G}}_{0}^{\dagger}$ converges to $\bm{\bm{G}}_{0}^{\dagger}$ Rakočević (1997). Further, because Algorithm 1 returns the eigenvectors corresponding to all the non-zero eigenvalues of $\bm{A}$ Tu et al. (2014), if those non-zero eigenvalues are distinct (i.e., their multiplicity is one), the outputs of Algorithm 1 are continuous with respect to $\bm{A}$. Therefore, the dynamic modes converge to the eigenvectors corresponding to the non-zero eigenvalues of $\bm{K}$ with probability one. ∎

Remark 1*.*

We have shown the convergence utilizing the assumption of ergodicity. However, the algorithm defined by Tu et al. Tu et al. (2014) does not require the sequential sampling in Eq. (10) as long as the corresponding columns of the data matrices are sampled with a fixed temporal interval. One can also prove the convergence (in probability) for this case using the law of large numbers if one assumes the snapshots in $\bm{Y}_{0}$ are independently sampled from $\mu_{\mathcal{M}}$.

Remark 2*.*

For an output of DMD to coincide with the modal decompositions via the Koopman operator, data must be generated from observables that span a subspace invariant to the Koopman operator. In this paper, we assume that the data in hand are intrinsically generated with such observables $\bm{g}$, as most of the existing DMD variants do. Some ways to design such observables manually are reviewed in Section V.

III Stochastic Koopman analysis with noisy observations

In the previous section, we considered systems with no stochastic elements. However, stochasticity often comprises an essential part of a variety of physical phenomena and sensing. In this section, we introduce the notions of process noise on dynamics and observation noise on observables, and discuss Koopman analysis and DMD for stochastic noisy systems.

III.1 Process noise on dynamics

Instead of deterministic dynamical system (1), consider a random dynamical system (RDS) Arnold (1998)

[TABLE]

with a measure-preserving base flow $\vartheta:\Omega\to\Omega$, where $(\Omega,\Sigma_{\Omega},\mu_{\Omega})$ is a probability space of process noise. We assume that $\omega_{t}$ is independent from $\bm{x}_{0},\dots,\bm{x}_{t}$. A one-step evolution of observables $g$ with regard to the RDS can be characterized by stochastic Koopman operator $\mathcal{K}_{\Omega}$ Mezić (2005), defined as

[TABLE]

where $\mathbb{E}_{\Omega}[\cdot]$ denotes expectation in sample space $\Omega$. Note that deterministic Koopman operator $\mathcal{K}$ can be regarded as a special case of $\mathcal{K}_{\Omega}$. Now let $K_{\Omega}$ be the restriction of $\mathcal{K}_{\Omega}$ to its invariant subspace $G$, suppose that a set of observables $\{g_{1},\dots,g_{n}\}$ spans $G$, and let $\bm{g}=\begin{bmatrix}g_{1}&\dots&g_{n}\end{bmatrix}^{\mathsf{T}}$. In addition, let $\bm{K}_{\Omega}\in\mathbb{C}^{n\times n}$ be the representation of $K_{\Omega}$ with regard to the components of $\bm{g}$. Then, we have

[TABLE]

where

[TABLE]

Given $\bm{x}_{0}$, the solution of (13) then becomes

[TABLE]

The modal decomposition of $\bm{g}$ via $K_{\Omega}$ can be obtained likewise, as shown in Section II.1. Regarding the characteristics of the spectra of the stochastic Koopman operator, Bagheri Bagheri (2014) elaborated on the effects of weak noise in a limit cycle, and Williams *et al. * Williams et al. (2015) applied a variant of DMD to the data obtained from a stochastic differential equation.

The standard DMD (Algorithm 1) is also applicable to the RDS and $\mathcal{K}_{\Omega}$ if there is no observation noise and the ergodicity (Assumption 1) also holds for the RDS. This can be shown in a manner similar to the one in Proposition 1, except for the definition of $\bm{G}_{0}$ and $\bm{G}_{1}$, as follows. Let $\bm{Y}_{0}$ and $\bm{Y}_{1}$ be the data matrices generated from RDS $\bm{f}_{\Omega}$ and observable $\bm{g}$ as in Eq. (10), and let us assume the whiteness on process noise.

Assumption 2**.**

Process noise $\omega$ is independently and identically distributed in time, i.e., for all $t^{\prime},t\in\mathbb{T}$,

[TABLE]

for some $\bm{P}\in\mathbb{C}^{n\times n}$.

Then, from the law of large numbers and the assumption of ergodicity, the empirical matrices

[TABLE]

respectively converge to

[TABLE]

with probability one. One can use this convergence property to show the applicability of Algorithm 1 for the RDS, as in the proof of Proposition 1.

III.2 Observation noise on observables

In addition to the process noise, let us take the observation noise into account. Consider a new (noisy) observable $\bm{h}:\mathcal{M}\times S\to\mathbb{C}^{n}$:

[TABLE]

where $\bm{w}:S\to\mathbb{C}^{n}$ is a random variable on a probability space $(S,\Sigma_{S},\mu_{S})$ of the observation noise. Hereafter, we denote $\bm{w}(s_{t})$ by $\bm{w}_{t}$ for notational simplicity. Now assume that $s$ is independent from $\bm{x}$ and that $\bm{w}$ is a white noise.

Assumption 3**.**

Observation noise $\bm{w}$ is zero-mean, has time-invariant finite variance, and is temporally uncorrelated, i.e., for all $t^{\prime},t\in\mathbb{T}$,

[TABLE]

for some $\bm{Q},\bm{R}\in\mathbb{C}^{n\times n}$.

Note that under the presence of observation noise, an output of DMD (Algorithm 1) no longer converges to the spectra of the Koopman operator. An output of total-least-squares DMD Dawson et al. (2016); Hemati et al. (2017) is unbiased even for noisy observations as long as the dynamics are deterministic, but it is biased as a realization of $K_{\Omega}$ for the RDS. These inconsistencies in the existing methods are partly revealed in the numerical examples in Section IV.

III.3 Statistics of noisy observables on RDS

We would like to develop a DMD algorithm for stochastic Koopman analysis that is always aware of both the process noise and observation noise. To this end, we summarize the statistics of noisy observable $\bm{h}$ on RDS $\bm{f}_{\Omega}$. Proofs of the lemmas in this section are deferred to Appendix for clarity of presentation.

First, assume that $\bm{g}$ is quasi-stationary (see Ljung (1999)), i.e.,

Assumption 4**.**

For almost all $\bm{x}_{0}\in\mathcal{M}$ and all $t^{\prime},t\in\mathbb{T}$,

[TABLE]

for some $\bm{G}\in\mathbb{C}^{n\times n}$.

Then, the second-order moment of $\bm{g}$, $\bm{G}_{t^{\prime},t}$, satisfies the following properties.

Lemma 1**.**

$\bm{G}_{t^{\prime},t}$* is expressed as*

[TABLE]

Corollary 1**.**

Denote $\mathbb{E}_{\mathcal{M}}\left[\bm{G}_{t+\tau,t}\right]$ by $\bm{G}_{\tau}$. Then,

[TABLE]

Now let us define $\bm{H}_{t^{\prime},t}\coloneqq\mathbb{E}_{\Omega,S}\left[\bm{h}(\bm{x}_{t^{\prime}})\bm{h}(\bm{x}_{t})^{\mathsf{H}}\right]$, where we have dropped argument $s$ of $\bm{h}$ for ease of notation. Then, $\bm{H}_{t^{\prime},t}$ satisfies the following properties.

Lemma 2**.**

$\bm{H}_{t^{\prime},t}$* is expressed as*

[TABLE]

Corollary 2**.**

Denote $\mathbb{E}_{\mathcal{M}}\left[\bm{H}_{t+\tau,t}\right]$ by $\bm{H}_{\tau}$. Then,

[TABLE]

III.4 Subspace DMD

Finally, we introduce a numerical method to compute an instance of the stochastic Koopman operator given noisy observations, namely, subspace DMD. Analogously to Eq. (10), let us define the data matrix as a concatenation of $m$ observations starting at time $t$, i.e.,

[TABLE]

Then, using a data quadruple $(\bm{Y}_{0},\bm{Y}_{1},\bm{Y}_{2},\bm{Y}_{3})$, we can obtain a calculation for $K_{\Omega}$ using the following theorem.

Theorem 1**.**

Define $\bm{Y}_{p},\bm{Y}_{f}\in\mathbb{C}^{2n\times m}$ by

[TABLE]

and let $\bm{O}=\bm{Y}_{f}\mathbb{P}_{\bm{Y}_{p}^{\mathsf{H}}}\in\mathbb{C}^{2n\times m}$ be the orthogonal projection of rows of $\bm{Y}_{f}$ onto the row space of $\bm{Y}_{p}$. Further, consider a compact SVD

[TABLE]

with $\bm{U}_{q}\in\mathbb{C}^{2n\times q}$, $\bm{S}_{q}\in\mathbb{C}^{q\times q}$, and $\bm{V}_{q}\in\mathbb{C}^{m\times q}$, where $q=\operatorname{rank}(\bm{O})$. Moreover, let $\bm{U}_{q1}$ be the first $n$ rows and $\bm{U}_{q2}$ be the last $n$ rows of $\bm{U}_{q}$. If $\operatorname{rank}(\bm{Y}_{p})=2n$ and $\operatorname{rank}(\bm{K}_{\Omega}\bm{G}+\bm{R})=n$, then in $m\to\infty$,

[TABLE]

with probability one.

Proof.

Let $\hat{\bm{H}}$ be the empirical matrix such that

[TABLE]

In $m\to\infty$, $\hat{\bm{H}}_{t+\tau,t}$ converges to $\bm{H}_{\tau}$ with probability one for all $t\in\mathbb{T}$ and $\tau\geq 0$, because, from the law of large numbers and the assumption of ergodicity,

[TABLE]

Because we have assumed $\operatorname{rank}(\bm{Y}_{p})=2n$, in $m\to\infty$,

[TABLE]

with probability one, where $\bm{D}=\bm{K}_{\Omega}\bm{G}+\bm{R}\in\mathbb{C}^{n\times n}$ and

[TABLE]

Because we have assumed $\operatorname{rank}(\bm{D})=n$, the rank of both $\bm{O}_{1}$ and $\bm{O}_{2}$ is $n$. Hence, in $m\to\infty$, $q$ also becomes $n$. Remember that by compact SVD (19), we have the decomposition of $\bm{O}$ into two rank-$n$ matrices, i.e.,

[TABLE]

Therefore, from Eq. (21), in $m\to\infty$, we have

[TABLE]

with probability one, where $\bm{T}\in\mathbb{C}^{n\times n}$ is an arbitrary unitary matrix. Consequently, $\bm{U}_{q1}$ and $\bm{U}_{q2}$ become $\bm{T}$ and $\bm{K}_{\Omega}\bm{T}$ respectively, and Eq. (20) holds. ∎

Based on Theorem 1, we present a subspace DMD algorithm as follows.

Algorithm 2** (Subspace DMD).**

Build matrices $\bm{Y}_{p}$ and $\bm{Y}_{f}$ like Eq. (18) from a data quadruple $(\bm{Y}_{0},\bm{Y}_{1},\bm{Y}_{2},\bm{Y}_{3})$. 2. 2.

Compute orthogonal projection $\bm{O}=\bm{Y}_{f}\mathbb{P}_{\bm{Y}_{p}^{\mathsf{H}}}$. 3. 3.

Compute compact SVD $\bm{O}=\bm{U}_{q}\bm{S}_{q}\bm{V}_{q}^{\mathsf{H}}$ and define $\bm{U}_{q1}$ and $\bm{U}_{q2}$ by the first and the last $n$ rows of $\bm{U}_{q}$, respectively. 4. 4.

Compute compact SVD $\bm{U}_{q1}=\bm{U}\bm{S}\bm{V}^{\mathsf{H}}$ and define $\tilde{\bm{A}}=\bm{U}^{\mathsf{H}}\bm{U}_{q2}\bm{V}\bm{S}^{-1}$. 5. 5.

Compute the eigenvalues $\lambda$ and eigenvectors $\tilde{\bm{w}}$ of $\tilde{\bm{A}}$. 6. 6.

Return dynamic modes $\bm{w}=\lambda^{-1}\bm{U}_{q2}\bm{V}\bm{S}^{-1}\tilde{\bm{w}}$ and corresponding eigenvalues $\lambda$.

Remark 3*.*

With subspace DMD, we can naturally conduct a low-rank approximation of dynamics by replacing the compact SVD in Step 3 with a truncated SVD. In contrast, in Algorithm 1 and total-least-squares DMD Dawson et al. (2016); Hemati et al. (2017), the low-rank approximation is achieved via the truncated proper orthogonal decomposition (POD). Note that there is also a line of research on the low-rank approximation of DMD, such as Chen et al. (2012); Wynn et al. (2013); Jovanović et al. (2014); Héas and Herzet (2017).

Remark 4*.*

Again note that we suppose that data are obtained with observable $\bm{g}$ that spans a subspace invariant to the Koopman operator, as in Algorithm 1 and other popular variants of DMD. See Section V for ways to design such observables.

IV Numerical Examples

We present numerical examples for the application of subspace DMD to several types of dynamical systems to show its empirical performance. When describing target dynamical systems in the following examples, we denote Gaussian white process noise by $\bm{e}$ and Gaussian white observation noise by $\bm{w}$. The standard deviation of the process noise is referred to as $\sigma_{p}$ and that of the observation noise as $\sigma_{o}$. Moreover, we denote the number of snapshots fed into algorithms by $m$ and the dimensionality of the data by $n$.

IV.1 Oscillation perturbed by noise

The stochastic Stuart–Landau equation on a complex-valued function $z(t)=r(t)\exp(i\theta(t))$ is defined as

[TABLE]

where $e(t)$ is Gaussian white noise with unit variance, and $i$ denotes the imaginary unit. The solution of this equation evolves on a limit cycle at $|z|=\sqrt{\mu}$ in the absence of process noise (i.e., $\sigma_{p}=0$). Bagheri Bagheri (2014) has analyzed the effects of weak process noise on the Koopman eigenvalues of the limit cycle of the Stuart–Landau equation, which can be summarized as follows; the continuous-time eigenvalues lie on the imaginary axis if process noise is absent because the data are completely periodic, but in contrast, when perturbation (phase diffusion, as shown in Figure 1a) is present owing to the process noise, a line of the eigenvalues is “bent” as shown in Figure 1c. Hence, by investigating the distribution of the eigenvalues, one can anticipate the presence and magnitude of the phase diffusion. To this end, we must eliminate the effects of observation noise if any, which produces an extra bias on the eigenvalues, leaving the effects of process noise.

Following the scheme in Dawson et al. (2016), we generated data using the following discretized Stuart–Landau equation in polar-coordinates with process noise:

[TABLE]

and a set of noisy trigonometric observables:

[TABLE]

where the magnitude of the observation noise was fixed to $\sigma_{o}=0.05$. We estimated the continuous-time eigenvalues using subspace DMD (Algorithm 2), standard DMD (Algorithm 1), total-least-squares DMD (tls-DMD) Dawson et al. (2016); Hemati et al. (2017), and noise-corrected DMD (nc-DMD) Dawson et al. (2016).

Figure 1b shows the eigenvalues without any process noise (i.e., $\sigma_{p}=0$), and Figure 1c shows the ones with process noise of $\sigma_{p}=0.5$. In both plots, we also show the “clean” eigenvalues computed with the data without the observation noise. When the process noise is present (in Figure 1c), while the eigenvalues estimated by tls-DMD differ from the clean ones (as reported in Dawson et al. (2016)), subspace DMD successfully estimates them. Note that the estimation by nc-DMD also coincides with the clean eigenvalues, but nc-DMD needs a precise estimation of magnitude of observation noise, which is often difficult to obtain. Subspace DMD can eliminate the effects of observation noise without such information while keeping the effects of the process noise.

IV.2 Noisy damping modes

Let us consider the stochastic Burger’s equation

[TABLE]

with $k>0$ and Gaussian space-time white noise $e(x,t)$. In fact, the eigenvalues of the Koopman operator on Burger’s equation (without process noise, i.e., $\sigma_{p}=0$) can be analytically obtained via the Cole–Hopf transformation and they correspond to the decaying modes of the solution Budišić et al. (2012); Kutz et al. (2016c). When process noise is present ($\sigma_{p}>0$), the solution of Burger’s equation becomes “rough,” but its global appearance remains similar to the case of no process noise, as shown in Figure 2a.

We approximated the solution of the stochastic Burger’s equation with $k=0.01$ and $\sigma_{p}=0.01$ using Crank–Nicolson–Maruyama method (see, e.g., Hausenblas (2003)) with initial condition $u(x,0)=\sin(2\pi x)$ and Dirichlet boundary condition $u(0,t)=u(1,t)=0$, setting the ranges by $x\in\left[0,1\right]$ and $t\in\left[0,1\right]$ and the discretization step sizes by $\Delta x=1\times 10^{-2}$ and $\Delta t=5\times 10^{-5}$. Based on the approximated solution, we finally generated data with observation noise

[TABLE]

where the magnitude of the observation noise was set $\sigma_{o}=0.001$. The estimated eigenvalues are plotted in Figure 2b. While the eigenvalues obtained by tls-DMD lie approximately on the imaginary axis, the estimation by subspace DMD agrees well with the eigenvalues computed with data that contain no observation noise. Again note that, though the estimation by nc-DMD also aligns with the clean eigenvalues, it requires a precise estimation of observation noise magnitude.

IV.3 Quantitative investigation of effects of noises

Let us investigate the performance of subspace DMD quantitatively using a simple linear system. We generated data using a linear time-invariant system

[TABLE]

whose Koopman eigenvalues obviously contain $\lambda$ and $\bar{\lambda}$. Moreover, we used the identity observable with observation noise

[TABLE]

We fixed the standard deviation of the process noise to $\sigma_{p}=0.1$, the eigenvalue to $\lambda=ri$ with $r=1.0$ or $0.9$, and the initial state to $\bm{x}_{0}=\begin{bmatrix}1&1\end{bmatrix}^{\mathsf{T}}$. Hence, this system exhibits oscillation perpetuated by the process noise when $r=1.0$ and is damped while being excited by the process noise when $r=0.9$. We applied subspace DMD, standard DMD, tls-DMD, and optimized DMD (opt-DMD) Chen et al. (2012) to multiple datasets generated with different random seeds. In those experiments, we have found that opt-DMD is unstable when $r=1.0$ and it does not output much reasonable results because it needs to compute exponentials of eigenvalues. Hence, the results of opt-DMD when $r=1.0$ are not plotted in Figures 3, 5, and 5.

In Figure 3, we show the 95% confidence intervals of the estimated eigenvalues for 1,000 random trials with observation noise of magnitude $\sigma_{o}=0.1$ and $m=$ 1,000 snapshots. When $r=1.0$, the estimations by subspace DMD, standard DMD, and tls-DMD scatter around the true value, while the results of standard DMD deviate a little more than the others. When $r=0.9$, the estimations by standard DMD, tls-DMD, and opt-DMD deviate from the true value; only the outputs of subspace DMD distribute around the true value.

In Figure 5, the relative errors $\varepsilon=|\tilde{\lambda}-\lambda|/|\lambda|$ of estimated eigenvalues $\tilde{\lambda}$ are plotted against different magnitudes of the observation noise $\sigma_{o}$, with the number of snapshots fed into the algorithms being fixed by $m=$ 1,000. When $r=1.0$, subspace DMD, standard DMD, and tls-DMD work almost equally well. When $r=0.9$, while the errors of standard DMD and opt-DMD rapidly grow and tls-DMD generates a regular bias, subspace DMD produces almost no bias and is tolerant to the observation noise.

In Figure 5, relative errors $\varepsilon$ are plotted against different $m$ with fixed $\sigma_{o}=0.1$. When $r=1.0$, subspace DMD, standard DMD, and tls-DMD converge when $m$ becomes large. When $r=0.9$, only subspace DMD converges, which is expected from Theorem 1.

IV.4 Low-rank high-dimensional data

We have shown the convergence of subspace DMD in the large sample limit in Theorem 1, but in practice, DMD is often applied in a high-dimensional setting, where the number of snapshots $m$ is much less than dimensionality of data $n$. To simulate such circumstances, we generated 500-dimensional data using a linear time-invariant system:

[TABLE]

where $\bm{L}\in\mathbb{R}^{500\times 2}$ satisfies $\bm{L}^{\mathsf{T}}\bm{L}=\bm{I}$, with $\sigma_{p}=0.1$ and $\sigma_{o}=0.1$. We prepared two datasets with different sizes, $m=50$ and $m=200$, and applied the following DMD variants: subspace DMD, standard DMD, tls-DMD, optimal low-rank DMD (lr-DMD) Héas and Herzet (2017), and optimal mode decomposition (OMD) Wynn et al. (2013). For each method, we introduced the way to obtain a low-rank solution; only the first two POD modes were used in the standard DMD and tls-DMD, the rank parameter was set to two in lr-DMD and OMD, and only the first two columns of $\bm{U}_{q}$ were used in subspace DMD.

In Figure 6a, we show the $95$% confidence intervals and averages of the estimated eigenvalues for 1,000 random trials with $m=50$. While the estimations by standard DMD and tls-DMD deviate far from the true value because of the process noise (and observation noise), the estimations by subspace DMD, lr-DMD, and OMD distribute around the true value. In particular, the variance of the estimation by subspace DMD is smaller than that of lr-DMD and OMD. Figure 6b shows the results with $m=200$; in this case, the distributions of the estimations by subspace DMD, lr-DMD, and OMD almost coincide.

IV.5 Application: cylinder wake

As an example of an application, we applied subspace DMD to a two-dimensional flow past a circular cylinder with Reynolds number $\mathrm{Re}=100$. We generated data using a solver based on the fast immersed boundary method with the multi-domain technique Taira and Colonius (2007); Colonius and Taira (2008) with four nested grids, each of which contains $450\times 200$ points. The diameter of the cylinder corresponds to $50$ points in the finest grid. The solver uses the third-order Runge–Kutta method with time-step $\Delta t=0.02$. We collected $400$ snapshots of the vorticity fields with intervals of size $10\Delta t$ from the limit cycle characterized by the Kármán vortex street. An example of the snapshots (without observation noise) and the time-variation of vorticity at two locations (A and B) are shown in Figure 7. We applied subspace DMD, standard DMD, and tls-DMD to the data contaminated with observation noise of $\sigma_{o}=0.1$. Every method was run with a low-rank approximation of $r=15$ because the first $15$ POD modes contained about $99.9$% of the energy of the original data.

In Figure 8a, the eigenvalues estimated with the noisy dataset and the noise-free dataset are plotted; subspace DMD and tls-DMD generate smaller biases than standard DMD does. The eigenvalues are numbered from one to seven in Figure 8a, according to their frequency (the magnitude of the imaginary part). In the remainder of Figure 8, we show the dynamic modes corresponding to eigenvalue 1 ($\sim 10i$) in the upper row and eigenvalue 4 ($\sim 40i$) in the lower row. We confirm that no adversarial effect is present in the dynamic modes computed by subspace DMD. Note that, in this cylinder wake experiment, no process noise (except for small errors due to the numerical integration) is involved. Subspace DMD is also applicable to classical (but frequent) situations of data analysis like this, where almost no process noise is present.

V Discussion

This section contains additional discussion on related methods and important properties of Koopman spectral analysis and DMD algorithms that have not been addressed sufficiently in this paper. These issues may imply future directions of research.

V.1 Relation to subspace system identification and extension to controlled systems

Subspace DMD has a strong connection to the methods called subspace system identification (see, e.g., Van Overschee and De Moor (1996); Katayama (2005)) in their computational methodology. Subspace system identification is a series of methods mainly for the identification of linear time-invariant systems, whereas in this paper, we present a similar methodology for nonlinear dynamical systems involving the observables and the stochastic Koopman operator.

Subspace system identification has been studied from the viewpoint of control theory and admits the presence of input signals distinguished from the process noise. Therefore, an extension of subspace DMD to controlled dynamical systems would be straightforward following the methodologies developed in the research of subspace system identification. Also, one may take a closed-loop controlled system into consideration. In the context of DMD, Proctor *et al. * Proctor et al. (2016) have discussed a variant of DMD for data obtained from the controlled systems.

V.2 Construction of Koopman invariant subspace

The key point of DMD as a numerical realization of the Koopman analysis lies in preparing a set of observables that spans a subspace invariant to $\mathcal{K}$. Several researchers have worked on this issue; Williams *et al. * Williams et al. (2015) proposed using a user-defined dictionary of observables to adopt DMD to highly nonlinear systems, and Brunton et al. Brunton et al. (2016b) utilized an identification technique based on a sparse regression Brunton et al. (2016c) to identify the dynamic-specific observables to be used. In addition, Kawahara Kawahara (2016) defined the Koopman analysis for observables in reproducing kernel Hilbert spaces to build a theory of DMD based on the reproducing kernels.

Another option, especially for deterministic systems, is to use delay coordinates, i.e., stacking observations at neighboring timestamps in each column of the data matrices. In general, a Krylov-like sequence of observables $\{g,\mathcal{K}g,\mathcal{K}^{2}g,\dots\}$ rapidly becomes almost linearly dependent, and thus can be used to obtain a subspace that is approximately invariant to $\mathcal{K}$. Based on the delayed measurements, we obtain a data matrix as a Hankel matrix. The use of delay coordinates for DMD was first discussed by Tu *et al. * Tu et al. (2014), and Brunton *et al. * Brunton et al. (2017) mentioned DMD based on Hankel matrices, referring to the well-known Taken’s theorem Takens (1981). Susuki and Mezić Susuki and Mezić (2015) defined an approximation of the Koopman analysis using Prony’s method, which also uses Hankel matrices. Arbabi and Mezić Arbabi and Mezić (2016) showed the convergence of DMD on Hankel matrices build with an observable contained in a Koopman invariant subspace. However, since delay coordinates with a linear monomial cannot span a Koopman invariant subspace of nonlinear systems exactly, one should use a combination of the nonlinear observables and the delay coordinates.

VI Conclusion

In this work, we developed subspace DMD, an algorithm for stochastic Koopman analysis with noisy observations. We have shown that the output of the proposed algorithm converges to the spectra of the stochastic Koopman operator in the large sample limit even if both process noise and observation noise are present. Moreover, we have shown its empirical performance with the numerical examples on different types of random dynamical systems. We also discussed the possible future directions of research, such as an extension to controlled dynamical systems and development of a way to obtain Koopman invariant subspaces.

Acknowledgements.

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

Appendix A Proofs of lemmas in Section III

A.1 Proof of Lemma 1

From Eq. (15), for the case of $t^{\prime}>t$, we have

[TABLE]

where the last equality holds because $\bm{e}$ is zero-mean and because of Assumption 2. For the case of $t^{\prime}>t$, from the definition of $\bm{G}_{t^{\prime},t}$ and the above equation, we have

[TABLE]

A.2 Proof of Lemma 2

From Eqs. (15) and (16), for the case of $t^{\prime}>t$, we have

[TABLE]

where the third and the fourth equalities are from Assumption 3. When $t^{\prime}=t$, from Assumption 3, we have

[TABLE]