Kernel methods for detecting coherent structures in dynamical data

TL;DR

This paper connects kernel-based dimensionality reduction techniques with dynamical systems analysis, introducing a novel approach to identify coherent structures in complex data using kernel CCA and a new method called coherent mode decomposition.

Contribution

It establishes a theoretical link between kernel CCA and transfer operators, and introduces coherent mode decomposition as a generalization of DMD for detecting coherent structures.

Findings

Kernel CCA can be interpreted via transfer operators.

Coherent sets can be computed efficiently with kernel CCA.

The proposed methods are demonstrated on fluid dynamics and molecular data.

Abstract

We illustrate relationships between classical kernel-based dimensionality reduction techniques and eigendecompositions of empirical estimates of reproducing kernel Hilbert space (RKHS) operators associated with dynamical systems. In particular, we show that kernel canonical correlation analysis (CCA) can be interpreted in terms of kernel transfer operators and that it can be obtained by optimizing the variational approach for Markov processes (VAMP) score. As a result, we show that coherent sets of particle trajectories can be computed by kernel CCA. We demonstrate the efficiency of this approach with several examples, namely the well-known Bickley jet, ocean drifter data, and a molecular dynamics problem with a time-dependent potential. Finally, we propose a straightforward generalization of dynamic mode decomposition (DMD) called coherent mode decomposition (CMD). Our results provide a generic machine learning approach to the computation of coherent sets with an objective score that can be used for cross-validation and the comparison of different methods.

Equations

$$\left(\mathcal{P}p\right)(y)$$

$$\left(\mathcal{K}f\right)(x)$$

$$\left(\mathcal{T}u\right)(y)=\frac{1}{\pi (y)}\int p_{\tau }(y\mid x)\pi (x)u(x)\mathrm{d}x.$$

$$\psi (Y)\otimes \varphi (X):f\mapsto \psi (Y){⟨\varphi (X),f⟩}_{{\mathbb{H}}_X}$$

$${\mathcal{C}}_{XX}$$

$${\mathcal{C}}_{YX}$$

$${⟨g,{\mathcal{C}}_{YX}f⟩}_{{\mathbb{H}}_Y}=\mathrm{Cov}[g(Y),f(X)]$$

$$\begin{array}{rl}{\hat{\mathcal{C}}}_{XX}& =\frac{1}{n}\sum _{i=1}^n\varphi (x_i)\otimes \varphi (x_i)=\frac{1}{n}\Phi {\Phi }^{\top },\\ {\hat{\mathcal{C}}}_{YX}& =\frac{1}{n}\sum _{i=1}^n\psi (y_i)\otimes \varphi (x_i)=\frac{1}{n}\Psi {\Phi }^{\top },\end{array}$$

$$G_{XX}$$

$$G_{YY}$$

$$\begin{array}{rl}{\mathcal{P}}_k& =({\mathcal{C}}_{XX}+\epsilon \mathcal{I}{)}^{-1}{\mathcal{C}}_{YX},\\ {\mathcal{K}}_k& =({\mathcal{C}}_{XX}+\epsilon \mathcal{I}{)}^{-1}{\mathcal{C}}_{XY}.\end{array}$$

$$\displaystyle\begin{split}\widehat{\mathcal{P}_{k}}&=(\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle XX}}+\varepsilon\mathcal{I})^{-1}\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle YX}}\\ &=\Psi\hskip 1.00006pt\big{(}G_{\mathit{\scriptscriptstyle XY}}^{-1}\hskip 1.00006pt(G_{\mathit{\scriptscriptstyle XX}}+n\varepsilon I)^{-1}\hskip 1.00006ptG_{\mathit{\scriptscriptstyle XY}}\big{)}\hskip 1.00006pt\Phi^{\top},\\ \widehat{\mathcal{K}_{k}}&=(\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle XX}}+\varepsilon\mathcal{I})^{-1}\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle XY}}\\ &=\Phi\hskip 1.00006pt(G_{\mathit{\scriptscriptstyle XX}}+n\varepsilon I)^{-1}\hskip 1.00006pt\Psi^{\top}.\end{split}$$

$$\hat{\mathcal{S}}=\Psi B{\Phi }^{\top }=\sum _{i,j=1}^n{b}_{ij}\psi (y_i)\otimes \varphi (x_j).$$

$$\underset{\begin{array}{c}f\in {\mathbb{H}}_X\\ g\in {\mathbb{H}}_Y\end{array}}{\sup }{⟨g,{\mathcal{C}}_{YX}f⟩}_{{\mathbb{H}}_Y}\text{s.t.}\{\begin{array}{l}{⟨f,{\mathcal{C}}_{XX}f⟩}_{{\mathbb{H}}_X}=1,\\ {⟨g,{\mathcal{C}}_{YY}g⟩}_{{\mathbb{H}}_Y}=1,\end{array}$$

$$\{\begin{array}{l}{\mathcal{C}}_{YX}f=\rho {\mathcal{C}}_{YY}g,\\ {\mathcal{C}}_{XY}g=\rho {\mathcal{C}}_{XX}f.\end{array}$$

$$\quad\big{(}\mathcal{C}_{\mathit{\scriptscriptstyle XX}}+\varepsilon\mathcal{I}\big{)}^{-1}\mathcal{C}_{\mathit{\scriptscriptstyle XY}}\big{(}\mathcal{C}_{\mathit{\scriptscriptstyle YY}}+\varepsilon\mathcal{I}\big{)}^{-1}\mathcal{C}_{\mathit{\scriptscriptstyle YX}}f=\rho^{2}f.$$

$$\Phi B{\Phi }^{\top }\hat{f}={\rho }^2\hat{f},$$

$$\big{(}\Phi\Phi^{\top}+n\varepsilon\mathcal{I}\big{)}^{-1}\Phi\Psi^{\top}\big{(}\Psi\Psi^{\top}+n\varepsilon\mathcal{I}\big{)}^{-1}\Psi\Phi^{\top}\widehat{f}=\rho^{2}\widehat{f}.$$

$$\{\begin{array}{l}(G_{YY}+n\epsilon I{)}^{-1}{G}_{XX}v=\rho w,\\ (G_{XX}+n\epsilon I{)}^{-1}{G}_{YY}w=\rho v,\end{array}$$

$$\left[\begin{array}{cc}0& G_{YY}\\ G_{XX}& 0\end{array}\right]\left[\begin{array}{c}v\\ w\end{array}\right]=\rho \left[\begin{array}{cc}(G_{XX}+n\epsilon I)& 0\\ 0& (G_{YY}+n\epsilon I)\end{array}\right]\left[\begin{array}{c}v\\ w\end{array}\right].$$

$$\big{(}\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle XX}}+\varepsilon\mathcal{I}\big{)}^{\!\nicefrac{{-1}}{{2}}}\hskip 1.00006pt\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle XY}}\hskip 1.00006pt\big{(}\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle YY}}+\varepsilon\mathcal{I}\big{)}^{-1}\hskip 1.00006pt\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle YX}}\hskip 1.00006pt\big{(}\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle XX}}+\varepsilon\mathcal{I}\big{)}^{\!\nicefrac{{-1}}{{2}}}\,\widetilde{v}=\rho^{2}\hskip 1.00006pt\widetilde{v}.$$

$$\big{(}\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle YY}}+\varepsilon\mathcal{I}\big{)}^{\nicefrac{{-1}}{{2}}}\hskip 1.00006pt\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle YX}}\hskip 1.00006pt\big{(}\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle XX}}+\varepsilon\mathcal{I}\big{)}^{\!\nicefrac{{-1}}{{2}}}$$

$$S(X,Y)=\sum _i^k{\rho }^r(X,Y)$$

$$(\mathcal{F}f)(y)=\int p_{\tau }(y\mid x)f(x)\mu (x)\mathrm{d}x,$$

$$(\mathcal{A}f)(y)$$

$$({\mathcal{A}}^{\ast }g)(x)$$

$${⟨{\mathcal{C}}_{YY}\mathcal{A}f,g⟩}_{{\mathbb{H}}_Y}$$

$$=\iint \frac{p(y\mid x)}{\nu (y)}f(x)\mu (x)\mathrm{d}xg(y)\nu (y)\mathrm{d}y$$

$$=\iint p(y\mid x)f(x)g(y)\mu (x)\mathrm{d}x\mathrm{d}y$$

$$=\iint p(x,y)f(x)g(y)\mathrm{d}x\mathrm{d}y$$

$$={\mathbb{E}}_{XY}[f(X)g(Y)]$$

$$={⟨{\mathcal{C}}_{YX}f,g⟩}_{{\mathbb{H}}_Y}.\text{qed}$$

$${\mathcal{K}}_k{\mathcal{A}}_{k}f={\rho }^{2}f.$$

$${⟨{\mathcal{A}}_{k}f,g⟩}_{\nu }={⟨{\mathcal{C}}_{YX}f,g⟩}_{{\mathbb{H}}_Y}={⟨f,{\mathcal{C}}_{XY}g⟩}_{{\mathbb{H}}_X}={⟨f,{\mathcal{K}}_{k}g⟩}_{\mu }.$$

$$G_{XX}={\mathbf{X}}^{\top }\mathbf{X}\text{and}{G}_{YY}={\mathbf{Y}}^{\top }\mathbf{Y},$$

Full text

While coherent sets of particles are common in dynamical systems, they are notoriously challenging to identify. In this article, we leverage the combination of a suite of methods designed to approximate the eigenfunctions of transfer operators with kernel embeddings in order to design an algorithm for detecting coherent structures in Langrangian data. It turns out that the resulting method is a well-known technique to analyze relationships between multidimensional variables, namely kernel canonical correlation analysis. Our algorithm successfully identifies coherent structures in several diverse examples, including oceanic currents and a molecular dynamics problem with a moving potential. Furthermore, we show that a natural extension of our algorithm leads to a coherent mode decomposition, a counterpart to dynamic mode decomposition.

I Introduction

Representing and learning effective, low-dimensional manifolds for complex, high-dimensional data is one of the cornerstones of machine learning and complex systems theory. In particular, for dynamical processes such as turbulent flows or molecular dynamics, it is known that much of their essential, long-time dynamics can be captured by linear models acting on low-dimensional manifolds, resulting in simpler, interpretable models, and potentially large savings in process simulation, prediction, and control. Due to their simplicity, linear methods to find low-dimensional subspaces are widely used, including principal component analysis (PCA) Hotelling33:PCA , canonical correlation analysis (CCA) Hotelling_Biometrika36_CCA , independent component analysis (ICA) Hyvarinen00:ICA , time-lagged independent component analysis (TICA) MS94 ; PPGDN13 , time-lagged canonical correlation analysis (TCCA) WuNo17 , and dynamic mode decomposition (DMD) SS08 .

Since the sought manifold is usually nonlinear in the direct state representation of the system, it is important to generalize the aforementioned methods to work in nonlinear feature spaces. Two particularly important classes of learning methods that go beyond applying linear algorithms to user-defined feature functions are neural network approachesMPWN:vampnets ; LiEtAl_Chaos17_EDMD_DL ; OttoRowley_LinearlyRecurrentAutoencoder and kernel methods such as kernel PCA Scholkopf98:KPCA , kernel CCA MRB01:CCA , kernel ICA Bach03:KICA , kernel TICA HZHM03:kTICA , and kernel EDMD WRK15 . The basic idea of kernel methods is to represent data by elements in reproducing kernel Hilbert spaces associated with positive definite kernel functions.

The novel contribution of this work is to derive kernel methods for the identification of coherent structures in high-dimensional dynamical data. We do this by establishing deep mathematical connections between kernel methods that have been proposed in machine learning and dynamical systems theory. Our main results are:

We show that kernel CCA, when applied to dynamical data, admits a natural interpretation in terms of kernel transfer operators and that the resulting eigenvalue problems are directly linked to methods for the computation of coherent sets. Importantly, kernel CCA predates recent methods for coherent set identification. 2. 2.

We show kernel CCA is optimal in the variational approach for Markov processes (VAMP) WuNo17 . Therefore, kernel CCA optimally approximates the transfer operator singular values and functions within kernel methods, and VAMP is a suitable optimization method for identifying coherent sets. 3. 3.

We propose a new method called coherent mode decomposition, which can be seen as a combination of CCA and DMD.

Establishing similar connections between machine learning and dynamical systems theory have previously led to advances in different applications: By means of the variational approach of conformation dynamics (VAC) NoNu13 , optimal estimators for the leading eigenfunctions of reversible time-homogeneous transfer operators have been made, which are important to identify metastable sets and rare events SchuetteFischerHuisingaDeuflhard_JCompPhys151_146 ; Bovier06:metastability . This insight has led to the introduction of the TICA method as a way to identify slow collective variables to molecular dynamics PPGDN13 ; SP13 —a key step in the modeling of rare events in molecules. Kernel embeddings of conditional probability distributionsSHSF09 ; MFSS16 have been related to Perron–Frobenius and Koopman transfer operators Ko31 ; LaMa94 and their eigenvalue decomposition in Ref. KSM17, . In a similar way, optimization of the VAMP score can be used to derive CCA as an optimal linear algorithm to approximate the singular functions of transfer operators WuNo17 . Doing the same with neural networks as function approximators leads to VAMPnets MPWN:vampnets . Here we show that a similar connection can be made with kernel methods and transfer operator singular functions and demonstrate that the kernel CCA algorithm approximates these functions. By exploiting that the singular functions are simultaneously the eigenfunctions of the forward-backward dynamics, we can extend this framework to the identification of so-called coherent sets—a generalization of metastable sets to nonautonomous and aperiodic systems FJ18:coherent . Coherent sets are regions of the state space that are not dispersed over a specific time interval. That is, if we let the system evolve, elements of a coherent set will, with a high probability, stay close together, whereas other regions of the state space might be distorted entirely. A large number of publications investigate the numerical approximation of coherent sets with other methods, e.g., Refs. FrSaMo10:coherent, ; FrJu15, ; WRR15, ; HKTH16:coherent, ; BK17:coherent, ; HSD18, ; FJ18:coherent, , see Ref. AP15:review, for an overview of approaches for Lagrangian data. We will not address the problem of possibly sparse or incomplete data. Our goal is to illustrate relationships with established kernel-based approaches and to show that existing methods—developed independently and with different applications in mind, predating many algorithms for the computation of finite-time coherent sets—can be directly applied to detect coherent sets in Lagrangian data.

A high-level overview of data-driven approaches for the approximation of transfer operators that are relevant for our considerations and relationships with the methods proposed below are shown in Figure 1. For the derivations of some of these methods the system is assumed to be reversible and there are other subtle differences, which will not be discussed here. A comparison of methods for time-homogeneous systems can be found in Ref. KNKWKSN18, and extensions to time-inhomogeneous systems in Ref. KWNS18:noneq, . Moreover, it was shown that the role played by the eigenfunctions in the time-homogeneous setting (for the detection of metastable sets) is assumed by the left- and right singular functions in the time-inhomogeneous setting (for the detection of coherent sets)KWNS18:noneq . The right singular functions encode information about the system at initial time $t$ and the left singular functions correspond to the system’s state at final time $t+\tau$. The connections between kernel CCA and the singular value decomposition of transfer operators will be described in more detail below.

The remainder of this paper is structured as follows: In Section II, we will briefly introduce transfer operators and review the notion of positive definite kernels and induced Hilbert spaces as well as nonlinear generalizations of covariance and cross-covariance matrices. We will then define empirical RKHS operators and show that diverse algorithms can be formulated as eigenvalue problems involving such operators. The relationships between kernel CCA and coherent sets will be studied in Section III. Furthermore, coherent mode decomposition will be derived. Section IV contains numerical results illustrating how to use the presented kernel-based methods for the analysis of dynamical systems. We conclude with a summary of the main results and open problems in Section V.

II Prerequisites

We briefly introduce transfer operators, reproducing kernel Hilbert spaces, and operators mapping from one such space to another one (or itself). For more details on the properties of these spaces and the introduced operators, we refer the reader to Refs. Schoe01, ; Steinwart2008:SVM, ; SC04:KernelMethods, and Refs. Baker70:XCov, ; Baker1973, ; KSM17, , respectively.

II.1 Transfer operators

Let $\{X_{t}\}_{t\geq 0}$ be a stochastic process defined on the state space $\mathbb{X}\subset\mathbb{R}^{d}$ and let $\tau$ be a fixed lag time. We assume that there exists a transition density function $p_{\tau}\colon\mathbb{X}\times\mathbb{X}\rightarrow\mathbb{R}$ such that $p_{\tau}(y\mid x)$ is the probability of $X_{t+\tau}=y$ given $X_{t}=x$. For $1\leq r\leq\infty$, let $L^{r}(\mathbb{X})$ denote the standard space of $r$-Lebesgue integrable functions on $\mathbb{X}$. Then, for a probability density $\mu$ on $\mathbb{X}$, let $L^{r}_{\mu}(\mathbb{X})$ be the spaces of $r$-integrable functions with respect to the corresponding probability measure induced by the density $\mu$; that is, $\left\lVert f\right\rVert^{r}_{L^{r}_{\mu}(\mathbb{X})}=\int|f(x)|^{r}\hskip 1.00006pt\mu(x)\hskip 1.00006pt\mathrm{d}x$.

Given a probability density $p\in L^{1}(\mathbb{X})$ and an observable $f\in L^{\infty}(\mathbb{X})$, we define the Perron–Frobenius operator $\mathcal{P}\colon L^{1}(\mathbb{X})\to L^{1}(\mathbb{X})$ and the Koopman operator $\mathcal{K}\colon L^{\infty}(\mathbb{X})\to L^{\infty}(\mathbb{X})$ by

[TABLE]

Assuming the process admits a unique equilibrium density $\pi$, i.e., $\mathcal{P}\pi=\pi$, we can define for $u=\pi(x)^{-1}\hskip 1.00006ptp(x)$ the Perron–Frobenius operator with respect to the equilibrium density $\mathcal{T}\colon L_{\pi}^{1}(\mathbb{X})\to L_{\pi}^{1}(\mathbb{X})$ as

[TABLE]

Under certain conditions, these transfer operators can be defined on $L^{r}(\mathbb{X})$ and $L^{r}_{\pi}(\mathbb{X})$ for other choices of $r$. From now on, we will always assume that they are well-defined for $r=2$ (see Refs. LaMa94, ; BaRo95, ; KKS16, for details). This is common whenever Hilbert space properties are needed in the context of transfer operators.

Remark II.1**.**

For time-homogeneous systems, the associated transfer operators depend only on the lag time $\tau$. If the system is time-inhomogeneous, on the other hand, the lag time is not sufficient to parametrize the evolution of the system since it also depends on the starting time. This is described in detail in Ref. KWNS18:noneq, . The transition density and the operators thus require two parameters; however, we will omit the starting time dependence for the sake of clarity.

II.2 Reproducing kernel Hilbert spaces

Given a set $\mathbb{X}$ and a space $\mathbb{H}$ of functions $f\colon\mathbb{X}\to\mathbb{R}$, $\mathbb{H}$ is called a reproducing kernel Hilbert space (RKHS) with inner product $\left\langle\cdot,\,\cdot\right\rangle_{\mathbb{H}}$ if there exists a function $k\colon\mathbb{X}\times\mathbb{X}\to\mathbb{R}$ with the following properties:

1. (i)

$\left\langle f,\,k(x,\cdot)\right\rangle_{\mathbb{H}}=f(x)$ for all $f\in\mathbb{H}$, and 2. (ii)

$\mathbb{H}=\overline{\operatorname{span}\{k(x,\cdot)\mid x\in\mathbb{X}\}}$.

The function $k$ is called a kernel and the first property above the reproducing property. A direct consequence is that $\left\langle k(x,\cdot),\,k(x^{\prime},\cdot)\right\rangle_{\mathbb{H}}=k(x,x^{\prime})$. That is, the map $\phi\colon\mathbb{X}\rightarrow\mathbb{H}_{\mathit{\scriptscriptstyle}}$ given by $x\mapsto k(x,\cdot)$ can be regarded as a feature map associated with $x$, the so-called canonical feature map.​***Such a feature map $\phi\colon\mathbb{X}\rightarrow\mathbb{H}_{\mathit{\scriptscriptstyle}}$ admitting the property $k(x,x^{\prime})=\left\langle\phi(x),\,\phi(x^{\prime})\right\rangle_{\mathbb{H}}$ is not uniquely defined. There are other feature space representations such as, for instance, the Mercer feature space.​Mercer ; Schoe01 ; Steinwart2008:SVM As long as we are only interested in kernel evaluations, however, it does not matter which one is considered. It is thus possible to represent data by functions in the RKHS. Frequently used kernels include the polynomial kernel and the Gaussian kernel, given by $k(x,x^{\prime})=(c+x^{\top}x^{\prime})^{p}$ and $k(x,x^{\prime})=\exp(-\|x-x^{\prime}\|_{2}^{2}/2\sigma^{2})$, respectively. While the feature space associated with the polynomial kernel is finite-dimensional, the feature space associated with the Gaussian kernel is infinite-dimensional; see, e.g., Ref. Steinwart2008:SVM, . Inner products in these spaces, however, are not evaluated explicitly, but only implicitly through kernel evaluations. This is one of the main advantages of kernel-based methods Scholkopf98:KPCA ; SC04:KernelMethods . Algorithms that can be purely expressed in terms of inner product evaluations can thus be easily kernelized, resulting, as described above, in nonlinear extensions of methods such as PCA, CCA, or TICA.

II.3 Covariance operators and Gram matrices

Let $(X,Y)$ be a random variable on $\mathbb{X}\times\mathbb{Y}$, where $\mathbb{X}\subset\mathbb{R}^{d_{x}}$ and $\mathbb{Y}\subset\mathbb{R}^{d_{y}}$. The dimensions $d_{x}$ and $d_{y}$ can in principle be different. For our applications, however, the spaces $\mathbb{X}$ and $\mathbb{Y}$ are often identical. The associated marginal distributions are denoted by $\mathbb{P}_{x}(X)$ and $\mathbb{P}_{y}(Y)$, the joint distribution by $\mathbb{P}(X,Y)$, and the corresponding densities—which we assume exist—by $p_{x}(x)$, $p_{y}(y)$, and $p(x,y)$, respectively. Furthermore, let $k$ and $l$ be the kernels associated with $\mathbb{X}$ and $\mathbb{Y}$ and $\phi$ and $\psi$ the respective feature maps. We will always assume that requirements such as measurability of the kernels and feature maps as well as separability of the RKHSs are satisfied.​†††In most cases, these properties follow from mild assumptions about $\mathbb{X}$ and $\mathbb{Y}$. For an in-depth discussion of these technical details, see Ref. Steinwart2008:SVM, . The RKHSs induced by the kernels $k$ and $l$ are denoted by $\mathbb{H}_{\mathit{\scriptscriptstyle X}}$ and $\mathbb{H}_{\mathit{\scriptscriptstyle Y}}$.

We will now introduce covariance operators and cross-covariance operators Baker70:XCov ; Baker1973 on RKHSs. In what follows, we will always assume that $\mathbb{E}_{\scriptscriptstyle\mathit{X}}[k(X,X)]<\infty$ and $\mathbb{E}_{\scriptscriptstyle\mathit{Y}}[l(Y,Y)]<\infty$, which ensures that these operators are well-defined and Hilbert–Schmidt (for a comprehensive overview of kernel covariance operators and their applications, see Ref. MFSS16, and references therein). For any $f\in\mathbb{H}_{\mathit{\scriptscriptstyle X}}$, let

[TABLE]

denote the tensor product operator Reed from $\mathbb{H}_{\mathit{\scriptscriptstyle X}}$ to $\mathbb{H}_{\mathit{\scriptscriptstyle Y}}$ defined by $\phi(X)$ and $\psi(Y)$.

Definition II.2** (Covariance operators).**

The covariance operator $\mathcal{C}_{\mathit{\scriptscriptstyle XX}}\colon\mathbb{H}_{\mathit{\scriptscriptstyle X}}\to\mathbb{H}_{\mathit{\scriptscriptstyle X}}$ and the cross-covariance operator $\mathcal{C}_{\mathit{\scriptscriptstyle YX}}\colon\mathbb{H}_{\mathit{\scriptscriptstyle X}}\to\mathbb{H}_{\mathit{\scriptscriptstyle Y}}$ are defined as

[TABLE]

Kernel covariance operators satisfy

[TABLE]

for all $f\in\mathbb{H}_{\mathit{\scriptscriptstyle X}}$, $g\in\mathbb{H}_{\mathit{\scriptscriptstyle Y}}$. Defining $\phi_{c}(X)=\phi(X)-\mathbb{E}_{\scriptscriptstyle X}[\phi(X)]$ and $\psi_{c}(Y)=\psi(Y)-\mathbb{E}_{\scriptscriptstyle Y}[\psi(Y)]$, the corresponding centered counterparts of the covariance and cross-covariance operators $\mathcal{C}_{\mathit{\scriptscriptstyle XX}}$ and $\mathcal{C}_{\mathit{\scriptscriptstyle YX}}$ are defined in terms of the mean-subtracted feature maps.

As these operators can in general not be determined analytically, empirical estimates are computed from data, i.e.,

[TABLE]

where $\Phi=[\phi(x_{1}),\dots,\phi(x_{n})]$ and $\Psi=[\psi(y_{1}),\dots,\psi(y_{n})]$ and the training data $\{(x_{i},y_{i})\}_{i=1}^{n}$ is drawn i.i.d. from $\mathbb{P}(X,Y)$. Analogously, the mean-subtracted feature maps can be used to obtain empirical estimates of the centered operators. Since in practice we often cannot explicitly deal with these operators, in particular if the feature space is infinite-dimensional, we seek to reformulate algorithms in terms of Gram matrices.

Definition II.3** (Gram matrices).**

Given training data as defined above, the Gram matrices $G_{\mathit{\scriptscriptstyle XX}},G_{\mathit{\scriptscriptstyle YY}}\in\mathbb{R}^{n\times n}$ are defined as

[TABLE]

For a Gram matrix $G$, its centered version $\widetilde{G}$ is defined by $\widetilde{G}=N_{0}\hskip 1.00006ptG\hskip 1.00006ptN_{0}$, where $N_{0}=I-\frac{1}{n}\mathds{1}\mathds{1}^{\top}$ and $\mathds{1}\in\mathbb{R}^{n}$ is a vector composed of ones Bach03:KICA . Note that centered Gram matrices are not regular.

Remark II.4**.**

In what follows, if not noted otherwise, we assume that the covariance operators $\mathcal{C}_{\mathit{\scriptscriptstyle XX}}$ and $\mathcal{C}_{\mathit{\scriptscriptstyle YY}}$ and the Gram matrices $G_{\mathit{\scriptscriptstyle XX}}$ and $G_{\mathit{\scriptscriptstyle YY}}$ are properly centered for CCA.

II.4 Kernel transfer operators

We now show how transfer operators can be written in terms of covariance and cross-covariance operators—this leads to the concept of kernel transfer operators. We assume the Perron–Frobenius operator and the Koopman operator to be well-defined on $L^{2}(\mathbb{X})$ as discussed in Section II.1. Kernel transfer operators follow from the assumption that densities and observables in $L^{2}(\mathbb{X})$ can be represented as elements of the RKHS $\mathbb{H}_{\mathit{\scriptscriptstyle X}}$. Under some technical requirements, such as $\int_{\mathbb{X}}k(x,x)\hskip 1.00006pt\mathrm{d}x=\int_{\mathbb{X}}\left\lVert\phi(x)\right\rVert_{\mathbb{H}_{\mathit{\scriptscriptstyle X}}}^{2}\hskip 1.00006pt\mathrm{d}x<\infty$, the elements of $\mathbb{H}_{\mathit{\scriptscriptstyle X}}$ are included in $L^{2}(\mathbb{X})$ when they are identified with the respective equivalence class of square integrable functions. This correspondence can be derived from the theory of $L^{2}(\mathbb{X})$ integral operatorsSteinwart2008:SVM and is often used in statistical learning theory RBD10 . We may therefore assume that we can identify RKHS elements with the corresponding equivalence classes of functions in $L^{2}(\mathbb{X})$. By requiring $\mathbb{E}_{\mu}[k(X,X)]<\infty$ for a probability density $\mu(x)$, we obtain a similar statement for $L^{2}_{\mu}(\mathbb{X})$.

We refer to Ref. KSM17, for the derivation of kernel transfer operators and a description of their relationships with kernel embeddings of conditional distributions. We will omit the technical details and directly define kernel transfer operators as the RKHS analogue of the standard transfer operators defined in Section II.1. Using the same integral representations as before and defining the transfer operators on $\mathbb{H}_{\mathit{\scriptscriptstyle X}}$ instead of $L^{2}(\mathbb{X})$, we obtain the kernel Perron–Frobenius operator $\mathcal{P}_{k}\colon\mathbb{H}_{\mathit{\scriptscriptstyle X}}\rightarrow\mathbb{H}_{\mathit{\scriptscriptstyle X}}$ and the kernel Koopman operator $\mathcal{K}_{k}\colon\mathbb{H}_{\mathit{\scriptscriptstyle X}}\rightarrow\mathbb{H}_{\mathit{\scriptscriptstyle X}}$, respectively.

By defining the time-lagged process $Y_{t}=X_{t+\tau}$, we can write kernel transfer operators in terms of covariance and cross-covariance operators KSM17 . Note that $X_{t}$ and $Y_{t}$ are defined on the same state space $\mathbb{X}$; therefore, we have $\mathbb{H}_{\mathit{\scriptscriptstyle X}}=\mathbb{H}_{\mathit{\scriptscriptstyle Y}}$ and hence $\mathcal{C}_{\mathit{\scriptscriptstyle YX}}:\mathbb{H}_{\mathit{\scriptscriptstyle X}}\rightarrow\mathbb{H}_{\mathit{\scriptscriptstyle X}}$ in this special case. We obtain the important properties $\mathcal{C}_{\mathit{\scriptscriptstyle XX}}\mathcal{P}_{k}\hskip 1.00006ptg=\mathcal{C}_{\mathit{\scriptscriptstyle YX}}g$ and $\mathcal{C}_{\mathit{\scriptscriptstyle XX}}\mathcal{K}_{k}\hskip 1.00006ptg=\mathcal{C}_{\mathit{\scriptscriptstyle XY}}g$ for all $g\in\mathbb{H}_{\mathit{\scriptscriptstyle X}}$, which allows us to write

[TABLE]

Here, $(\mathcal{C}_{\mathit{\scriptscriptstyle XX}}+\varepsilon\mathcal{I})^{-1}$ is the Tikhonov-regularized inverse of $\mathcal{C}_{\mathit{\scriptscriptstyle XX}}$ with regularization parameter $\varepsilon>0$.​‡‡‡See Refs. Gr93, ; EG96, ; EHN96, for a detailed discussion of ill-posed inverse problems and the regularization of bounded linear operators on Hilbert spaces. Note the abuse of notation, since equality in the above inverse problems is only given asymptotically for $\varepsilon\to 0$ and pointwise for feasible $\mathcal{C}_{\mathit{\scriptscriptstyle YX}}g\in\mathbb{H}_{\mathit{\scriptscriptstyle X}}$. Since $\mathcal{C}_{\mathit{\scriptscriptstyle XX}}$ is a compact operator, it does not admit a globally defined bounded inverse if the RKHS is infinite-dimensional. However, $(\mathcal{C}_{\mathit{\scriptscriptstyle XX}}+\varepsilon\mathcal{I})^{-1}$ always exists and is bounded. In fact, the operators $\mathcal{P}_{k}$ and $\mathcal{K}_{k}$ as given in the regularized form above are Hilbert–Schmidt.

The above notation and regularization of inverse covariance operators is standard in the context of kernel embeddings of conditional distributions and related Bayesian learning techniques. We refer to Refs. SHSF09, ; Song2013, ; Fukumizu13:KBR, ; Fukumizu15:NBI, ; MFSS16, for detailed discussions of properties of this ill-posed inverse problem in specific applications.

By replacing the analytical covariance operators with their empirical estimates in (2), we obtain empirical estimates for kernel transfer operators KSM17 . As done with empirical covariance operators in (1), it is possible to rewrite the empirical estimates of kernel transfer operators in terms of RKHS features in $\Phi$ and $\Psi$ (see Refs. MFSS16, ; KSM17, for the derivation):

[TABLE]

In this case, $\Phi$ and $\Psi$ both contain observations in the same space $\mathbb{H}_{\mathit{\scriptscriptstyle X}}$, since $X_{t}$ and $Y_{t}$ are both defined on $\mathbb{X}$.

II.5 Empirical RKHS operators

In what follows, we will consider finite-rank RKHS operators given by a matrix which represents the action of the operator on fixed elements in the RKHSs. We will use this general setting to formulate results about the eigenvalues and eigenfunctions of empirical RKHS operators. Given a matrix $B\in\mathbb{R}^{n\times n}$, we define the bounded finite-rank operator $\widehat{\mathcal{S}}\colon\mathbb{H}_{\mathit{\scriptscriptstyle X}}\to\mathbb{H}_{\mathit{\scriptscriptstyle Y}}$ by

[TABLE]

We remark that although $\Psi$ and $\Phi$ may contain infinite-dimensional objects, we express inner products between RKHS elements in the classical matrix-vector multiplication form. That is, we interpret the embedded RKHS elements as (potentially infinite-dimensional) column vectors. This notation has become a de-facto standard in the machine learning community MFSS16 . We can write empirical estimates of covariance operators in the form of (4). If the RKHS training features in $\Phi$ and $\Psi$ are generated i.i.d. by the joint probability distribution $\mathbb{P}(X,Y)$ of random variables $X$ and $Y$, then the cross-covariance operator $\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle YX}}$ takes the general form of an empirical RKHS operator with $B=\frac{1}{n}I$. We obtain $\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle XX}}$ as another special case with identical features $\Psi=\Phi$ drawn only from $\mathbb{P}(X)$. Furthermore, the empirical estimates of the kernel Perron–Frobenius and kernel Koopman operator are special cases of $\widehat{\mathcal{S}}$ as seen in (3) with $B=G_{\mathit{\scriptscriptstyle XY}}^{-1}\hskip 1.00006pt(G_{\mathit{\scriptscriptstyle XX}}+n\varepsilon I)^{-1}\hskip 1.00006ptG_{\mathit{\scriptscriptstyle XY}}$ and $B=(G_{\mathit{\scriptscriptstyle XX}}+n\varepsilon I)^{-1}$, respectively. Note that the roles of $\Phi$ and $\Psi$ are interchanged for the empirical estimate of the Koopman operator, i.e., it is of the form $\widehat{\mathcal{S}}=\Phi B\Psi^{\top}$.

We now show how spectral decomposition techniques can be applied to empirical RKHS operators in this general setting.​§§§In general, all considered kernel transfer operators in this paper are compositions of compact and bounded operators and therefore compact. They admit series representations in terms of singular value decompositions as well as eigendecompositions in the self-adjoint caseReed . The functional analytic details and the convergence of $\widehat{\mathcal{S}}$ and its spectral properties in the infinite-data limit depend on the specific scenario and are beyond the scope of this paper. We can compute eigenvalues and corresponding eigenfunctions of $\widehat{\mathcal{S}}$ by solving auxiliary matrix eigenvalue problems. For the sake of self-containedness, we briefly reproduce the eigendecomposition result from Ref. KSM17, .

Proposition II.5**.**

Suppose $\Phi$ and $\Psi$ contain linearly independent elements. Let $\widehat{\mathcal{S}}=\Psi B\Phi^{\top}$, then

1. (i)

$\widehat{\mathcal{S}}$* has an eigenvalue $\lambda\neq 0$ with corresponding eigenfunction $\varphi=\Psi v$ if and only if $v$ is an eigenvector of $B\hskip 1.00006ptG_{\mathit{\scriptscriptstyle XY}}$ associated with $\lambda$, and, similarly,* 2. (ii)

$\widehat{\mathcal{S}}$* has an eigenvalue $\lambda\neq 0$ with corresponding eigenfunction $\varphi=\Phi\hskip 1.00006ptG_{\mathit{\scriptscriptstyle XX}}^{-1}\hskip 1.00006ptv$ if and only if $v$ is an eigenvector of $G_{\mathit{\scriptscriptstyle XY}}\hskip 1.00006ptB$.*

For the Gaussian kernel, linear independence of elements in $\Phi$ and $\Psi$ reduces to requiring that the training data contains pairwise distinct elements in $\mathbb{X}$ and $\mathbb{Y}$, respectively. For dynamical systems applications, we typically assume that $\Phi$ and $\Psi$ contain information about the system at time $t$ and at time $t+\tau$, respectively. A more detailed version of Proposition II.5 and its extension to the singular value decomposition are described in Ref. MSKS18, . Further properties of $\widehat{\mathcal{S}}$ and its decompositions will be studied in future work. Note that we generally assume that empirical estimates of RKHS operators converge in probability to their analytical counterparts in operator norm in the infinite data limit. These statistical properties and the resulting associated spectral convergence are examined in for example in Ref. RBD10, .

II.6 Applications of RKHS operators

Decompositions of RKHS operators have diverse applications, which we will only touch upon here. We will consider a specific problem—namely, kernel CCA—in Section III.

1. (a)

By sampling points from the uniform distribution, the Mercer feature mapMercer ; Schoe01 ; Steinwart2008:SVM with respect to the Lebesgue measure on $\mathbb{X}$ can be approximated by computing eigenfunctions of $\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle XX}}$—i.e., $B=\frac{1}{n}I$ and the auxiliary matrix eigenvalue problem is $\frac{1}{n}G_{\mathit{\scriptscriptstyle XX}}\hskip 1.00006ptv=\lambda\hskip 1.00006ptv$—as shown in Ref. MSKS18, . This can be easily extended to other measures. 2. (b)

Similarly, given an arbitrary data set $\{x_{i}\}_{i=1}^{n}$, kernel PCA computes the eigenvectors corresponding to the largest eigenvalues of the centered Gram matrix $G_{\mathit{\scriptscriptstyle XX}}$ and defines these eigenvectors as the data points projected onto the respective principal components. It is well-known that kernel PCA can also be defined in terms of the centered covariance operator $\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle XX}}$. A detailed connection of the spectrum of the Gram matrix and the covariance operator is given in Ref. STWCK02, . Up to scaling, the eigenfunctions evaluated in the data points correspond to the principal components. 3. (c)

Given training data $x_{i}\sim p_{x}$ and $y_{i}=\Theta^{\tau}(x_{i})$, where $\Theta$ denotes the flow associated with the dynamical system and $\tau$ the lag time—that is, if $x_{i}$ is the state of the system at time $t$, then $y_{i}$ is the state of the system at time $t+\tau$—, we define $\Phi$ and $\Psi$ as above. Eigenvalues and eigenfunctions of kernel transfer operators can be computed by solving a standard matrix eigenvalue problem (see Proposition II.5). Eigendecompositions of these operators result in metastable sets. For more details and real-world examples, see Refs. KSM17, ; KBSS18, . The main goal of this paper is the extension of the aforementioned methods to compute coherent sets instead of metastable sets.

III Kernel CCA and coherent sets

This section contains the main results of our paper. We derive kernel CCA MRB01:CCA for finite and infinite-dimensional feature spaces from the viewpoint of dynamical systems, and show that kernel CCA can be used to approximate coherent sets in dynamical data. Furthermore, we derive the new coherent mode decomposition method.

Given two multidimensional random variables $X$ and $Y\!$, standard CCA finds two sets of basis vectors such that the correlations between the projections of $X$ and $Y$ onto these basis vectors are maximized Borga01:CCA . The new bases can be found by computing the dominant eigenvalues and corresponding eigenvectors of a matrix composed of covariance and cross-covariance matrices. Just like kernel PCA is a nonlinear extension of PCA, kernel CCA is a generalization of CCA. The goal of kernel CCA is to find two nonlinear mappings $f(X)$ and $g(Y)$, where $f\in\mathbb{H}_{\mathit{\scriptscriptstyle X}}$ and $g\in\mathbb{H}_{\mathit{\scriptscriptstyle Y}}$, such that their correlation is maximized Fukumizu07:KCCA . That is, instead of matrices, kernel CCA is now formulated in terms of covariance and cross-covariance operators. More precisely, the kernel CCA problem can be written as

[TABLE]

and the solution is given by the eigenfunctions corresponding to the largest eigenvalue of the problem

[TABLE]

Further eigenfunctions corresponding to subsequent eigenvalues can be taken into account as in the standard setting described above. In practice, the eigenfunctions are estimated from finite samples. The empirical estimates of $f$ and $g$ are denoted by $\widehat{f}$ and $\widehat{g}$, respectively.

Example III.1**.**

In order to illustrate kernel CCA, let us analyze a synthetic data set similar to the one described in Ref. Fukumizu07:KCCA, using a Gaussian kernel with bandwidth $\sigma=0.3$. Algorithms to solve the CCA problem will be described below. The results are shown in Figure 2. Classical CCA would not be able to capture the nonlinear relationship between $X$ and $Y$.  $\blacktriangle$

III.1 RKHS operator formulation

Since the inverses of the covariance operators in general do not exist, the regularized versions $(\mathcal{C}_{\mathit{\scriptscriptstyle XX}}+\varepsilon\mathcal{I})^{-1}$ and $(\mathcal{C}_{\mathit{\scriptscriptstyle YY}}+\varepsilon\mathcal{I})^{-1}$ (cf. Section II.4) are also typically used in the context of CCAFukumizu07:KCCA . Solving the first equation in (5) for $g$ and inserting it into the second equation, this results in

[TABLE]

Comparing this with the aforementioned transfer operator representations (2), $(\mathcal{C}_{\mathit{\scriptscriptstyle XX}}+\varepsilon\mathcal{I})^{-1}\mathcal{C}_{\mathit{\scriptscriptstyle XY}}$ can be interpreted as an approximation of the kernel Koopman operator, and $(\mathcal{C}_{\mathit{\scriptscriptstyle YY}}+\varepsilon\mathcal{I})^{-1}\mathcal{C}_{\mathit{\scriptscriptstyle YX}}$ as a kernel Koopman operator where now the roles of $X$ and $Y$ are reversed or as a reweighted Perron–Frobenius operator. The composition of these operators corresponds to a push-forward and subsequent pull-back of a density $f$. Eigenfunctions of the operator whose associated eigenvalues are close to one thus remain nearly unchanged under the forward-backward dynamics. This is closely related to the notion of coherence as introduced in Refs. FrSaMo10:coherent, ; Froyland13:coherent, and will be discussed in Section III.4.

Lemma III.2**.**

Replacing the covariance and cross-covariance operators by their empirical estimates, the eigenvalue problem (6) can be written as

[TABLE]

with $B=(G_{\mathit{\scriptscriptstyle XX}}+n\varepsilon I)^{-1}(G_{\mathit{\scriptscriptstyle YY}}+n\varepsilon I)^{-1}G_{\mathit{\scriptscriptstyle YY}}$.

Proof.

Inserting the definitions of the empirical covariance and cross-covariance operators yields

[TABLE]

Using $\Psi^{\top}\left(\Psi\Psi^{\top}+n\varepsilon\mathcal{I}\right)^{-1}=\left(\Psi^{\top}\Psi+n\varepsilon I\right)^{-1}\Psi^{\top}$, see Ref. MFSS16, , and a similar identity for $\Phi$ concludes the proof. ∎

That is, the empirical RKHS operator for kernel CCA is of the form $\widehat{S}=\Phi B\Phi^{\top}$. Applying Proposition II.5, we must solve the auxiliary problem

1. (i)

$(G_{\mathit{\scriptscriptstyle XX}}+n\varepsilon I)^{-1}(G_{\mathit{\scriptscriptstyle YY}}+n\varepsilon I)^{-1}G_{\mathit{\scriptscriptstyle YY}}\hskip 1.00006ptG_{\mathit{\scriptscriptstyle XX}}\hskip 1.00006ptv=\rho^{2}\hskip 1.00006ptv$, with $\widehat{f}=\Phi\hskip 1.00006ptv$, or 2. (ii)

$G_{\mathit{\scriptscriptstyle XX}}(G_{\mathit{\scriptscriptstyle XX}}+n\varepsilon I)^{-1}(G_{\mathit{\scriptscriptstyle YY}}+n\varepsilon I)^{-1}G_{\mathit{\scriptscriptstyle YY}}\hskip 1.00006ptv=\rho^{2}\hskip 1.00006ptv$, with $\widehat{f}=\Phi\hskip 1.00006pt(G_{\mathit{\scriptscriptstyle XX}}+n\varepsilon I)^{-1}\hskip 1.00006ptv$.

Since $G_{\mathit{\scriptscriptstyle XX}}$ and $(G_{\mathit{\scriptscriptstyle XX}}+n\varepsilon I)^{-1}$ as well as $G_{\mathit{\scriptscriptstyle YY}}$ and $(G_{\mathit{\scriptscriptstyle YY}}+n\varepsilon I)^{-1}$ commute, the first problem can be equivalently rewritten as $(G_{\mathit{\scriptscriptstyle XX}}+n\varepsilon I)^{-1}G_{\mathit{\scriptscriptstyle YY}}(G_{\mathit{\scriptscriptstyle YY}}+n\varepsilon I)^{-1}G_{\mathit{\scriptscriptstyle XX}}\hskip 1.00006ptv=\rho^{2}\hskip 1.00006ptv$ and the second as $(G_{\mathit{\scriptscriptstyle XX}}+n\varepsilon I)^{-1}G_{\mathit{\scriptscriptstyle XX}}G_{\mathit{\scriptscriptstyle YY}}(G_{\mathit{\scriptscriptstyle YY}}+n\varepsilon I)^{-1}\hskip 1.00006ptv=\rho^{2}\hskip 1.00006ptv$. The eigenfunction associated with the largest eigenvalue solves the CCA problem, but in order to detect coherent sets, we will need more eigenfunctions later. To obtain the function $g$ corresponding to $\rho$, we compute

1. (i)

$\widehat{g}=\frac{1}{\rho}\Psi(G_{\mathit{\scriptscriptstyle YY}}+n\varepsilon I)^{-1}G_{\mathit{\scriptscriptstyle XX}}v$, or 2. (ii)

$\widehat{g}=\frac{1}{\rho}\Psi(G_{\mathit{\scriptscriptstyle YY}}+n\varepsilon I)^{-1}G_{\mathit{\scriptscriptstyle XX}}(G_{\mathit{\scriptscriptstyle XX}}+n\hskip 1.00006pt\varepsilon I)^{-1}v$.

Algorithm III.3**.**

The CCA problem can be solved as follows:

Choose a kernel $k$ and regularization $\varepsilon$.

Compute the centered gram matrices $G_{\mathit{\scriptscriptstyle XX}}$ and $G_{\mathit{\scriptscriptstyle YY}}$.

Solve $G_{\mathit{\scriptscriptstyle XX}}(G_{\mathit{\scriptscriptstyle XX}}+n\varepsilon I)^{-1}(G_{\mathit{\scriptscriptstyle YY}}+n\varepsilon I)^{-1}G_{\mathit{\scriptscriptstyle YY}}\hskip 1.00006ptv=\rho^{2}\hskip 1.00006ptv$.

The corresponding eigenfunction $\widehat{f}$ evaluated at all data points $x_{1},\dots,x_{n}$, denoted by $\widehat{f}_{X}$, is then approximately given by the vector $v$. We can evaluate the eigenfunctions at any other point as described above, but we will mainly use the eigenfunction evaluations at the sampled data points for clustering into coherent sets.

Algorithm III.3 is based on the second problem formulation, i.e., item (ii) above. However, the first variant can be used in the same way. Alternatively, we can rewrite it as an eigenvalue problem of the form

[TABLE]

and, consequently,

[TABLE]

Other formulations can be derived in a similar fashion. The advantage is that no matrices have to be inverted. However, the size of the eigenvalue problem doubles, which might be problematic if the number data points $n$ is large.

Remark III.4**.**

In order to apply the algorithms, we first need to choose a kernel and then tune its parameters, e.g., the bandwidth $\sigma$ of the Gaussian kernel, and also the regularization parameter $\varepsilon$. If the bandwidth is too small, this leads to overfitting and to oversmoothing if it is too large. Cross-validation techniques can be used to select suitable hyperparameters. The kernel itself determines the complexity of the function space in which the eigenfunctions are sought (see Section III.4). Additionally, the results depend on the lag time $\tau$. Sets that are coherent for a given lag time are not necessarily coherent for a different lag time since these sets might be dispersed again.

The generalized eigenvalue problem (7) is almost identical to the one derived in Ref. Bach03:KICA, , with the difference that regularization is applied in a slightly different way. That is, the direct eigendecomposition of RKHS operators as proposed in Ref. KSM17, results, as expected, in variants of kernel CCA. The statistical convergence of kernel CCA, showing that finite sample estimators converge to the corresponding population counterparts, has been established in Ref. Fukumizu07:KCCA, . Kernel CCA can be extended to more than two variables or views of the data as described in Refs. Bach03:KICA, ; SC04:KernelMethods, , which might also have relevant applications in the dynamical systems context.

III.2 Finite-dimensional feature space

If the state spaces of the kernels $k$ and $l$ are finite-dimensional, we can directly solve the eigenvalue problem (5) or (6). Assuming the feature space of the kernel $k$ is $r_{x}$-dimensional and spanned by the basis functions $\{\phi_{1},\dots,\phi_{r_{x}}\}$, we define $\phi\colon\mathbb{X}\to\mathbb{R}^{r_{x}}$ by $\phi(x)=[\phi_{1}(x),\dots,\phi_{r_{x}}(x)]^{\top}$. That is, we are now using an explicit feature space representation. This induces a kernel by defining $k(x,x^{\prime})=\left\langle\phi(x),\,\phi(x^{\prime})\right\rangle$.​¶¶¶For the Mercer feature space representationMercer ; Schoe01 the functions form an orthogonal basis, but orthogonality is not required here. We could, for instance, select a set of radial basis functions, monomials, or trigonometric functions. Analogously, we define a vector-valued function $\psi\colon\mathbb{Y}\to\mathbb{R}^{r_{y}}$, with $\psi(y)=[\psi_{1}(y),\dots,\psi_{r_{y}}(y)]^{\top}$, where $r_{y}$ is the dimension of the feature space of the kernel $l$. Any function in the respective RKHS can be written as $f=\alpha^{\top}\phi$ and $g=\beta^{\top}\psi$, where $\alpha\in\mathbb{R}^{r_{x}}$ and $\beta\in\mathbb{R}^{r_{y}}$ are coefficient vectors.

Given training data $\{(x_{i},y_{i})\}_{i=1}^{n}$ drawn from the joint probability distribution, we obtain $\Phi\in\mathbb{R}^{r_{x}\times n}$ and $\Psi\in\mathbb{R}^{r_{y}\times n}$ and can compute the centered covariance and cross-covariance matrices $\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle XX}}$, $\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle XY}}$, and $\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle YY}}$ explicitly.

Algorithm III.5**.**

Given explicit feature maps, we obtain the following CCA algorithm:

Select basis functions $\phi$ and $\psi$ and regularization $\varepsilon$.

Compute (cross-)covariance matrices $\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle XX}}$, $\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle XY}}$, $\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle YX}}$, and $\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle YX}}$.

*Solve the eigenvalue problem *

\big{(}\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle XX}}+\varepsilon\mathcal{I}\big{)}^{-1}\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle XY}}\big{(}\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle YY}}+\varepsilon\mathcal{I}\big{)}^{-1}\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle YX}}v=\rho^{2}\hskip 1.00006ptv.

The eigenfunctions are then given by $\widehat{f}(x)=\left\langle v,\,\phi(x)\right\rangle$. Expressions for $\widehat{g}$ can be derived analogously.

The difference between the Gram matrix approach described in Section III.1 and the algorithm proposed here is that the size of the eigenvalue problem associated with the former depends on the number of data points and permits the dimension of the feature space to be infinite-dimensional, whereas the eigenvalue problem associated with the latter depends on the dimension of the feature space but not on the size of the training data set. This is equivalent to the distinction between extended dynamic mode decomposition (EDMD) WKR15 and kernel EDMD WRK15 (or the variational approach NoNu13 and kernel TICA SP15 , where the system is typically assumed to be reversible; see Ref. KSM17, for a detailed comparison) with the small difference that often the Moore–Penrose pseudoinverse Penrose is used for EDMD in lieu of the Tikhonov-regularized inverse.

III.3 Relationships with VAMP

Defining v=\big{(}\widehat{\mathcal{C}}_{\mathit{\scriptscriptstyle XX}}+\varepsilon\mathcal{I}\big{)}^{\!\nicefrac{{-1}}{{2}}}\hskip 1.00006pt\widetilde{v}, the eigenvalue problem in Algorithm III.5 becomes

[TABLE]

The transformed eigenvectors $\widetilde{v}$ are thus equivalent to the right singular vectors of the matrix

[TABLE]

and the values $\rho$ are given by the singular values, which we assume to be sorted in nonincreasing order.

This form is the kernel version of the TCCA method that has been first derived as a way to approximate transfer operator singular functions using VAMP, see Ref. WuNo17, . VAMP is an optimization principle that defines a score function whose optimum leads to specific data-driven algorithms. The VAMP-$r$ score is defined as

[TABLE]

where $\rho^{r}(X,Y)$ are the singular value estimates obtained from an SVD of, e.g., (8), and $r$ is a positive integer.

When using $r=1$ and kernel feature functions, we obtain the kernel CCA algorithms III.3 or III.5. When doing the same on an explicit basis set of feature functions, we obtain TCCA, i.e., time-lagged CCA WuNo17 ; No18 . However, since with VAMP a score (or loss) function is available, we have now turned the identification of coherent sets into a generic machine learning problem. For example, training neural networks with VAMP results in VAMPnets, a deep learning method to low-rank approximattion of the transfer operators and the identification of metastable or coherent sets.

On the other hand, the fact that kernel CCA results from maximizing the VAMP score within a kernel approach shows that we can use the VAMP score in the context of cross-validation in order to optimally determine hyperparameters such as the kernel function. Furthermore we can explore other choices as $r=1$. For example, the choice $r=2$ has an interesting interpretation in terms of kinetic maps, which are embeddings of the dominant eigenspace or singular space of a transfer operator where Euclidean distances are related to timescales of transitions NoeClementi_JCTC15_KineticMap .

III.4 Relationships between kernel CCA and transfer operators

We have seen in Section III.1 that the resulting eigenvalue problem (6) involves expressions resembling kernel transfer operators. The goal now is illustrate how this eigenvalue problem is related to the operators derived in Ref. BK17:coherent, for detecting coherent sets. We first introduce a forward operator $\mathcal{F}\colon L_{\mu}^{2}(\mathbb{X})\to L^{2}(\mathbb{Y})$ by

[TABLE]

where $\mu$ is some reference density of interest. Furthermore, let $\nu=\mathcal{F}\mathds{1}$ be the image density obtained by mapping the indicator function on $\mathbb{X}$ forward in time. Normalizing $\mathcal{F}$ with respect to $\nu$, we obtain a new operator $\mathcal{A}\colon L_{\mu}^{2}(\mathbb{X})\to L_{\nu}^{2}(\mathbb{Y})$ and its adjoint $\mathcal{A}^{*}\colon L_{\nu}^{2}(\mathbb{Y})\to L_{\mu}^{2}(\mathbb{X})$, with

[TABLE]

It holds that $\left\langle\mathcal{A}f,\,g\right\rangle_{\nu}=\left\langle f,\,\mathcal{A}^{*}g\right\rangle_{\mu}$. Consequently, $\mathcal{A}$ plays the role of a reweighted Perron–Frobenius operator, whereas $\mathcal{A}^{*}$ can be interpreted as an analogue of the Koopman operator (note that $\mathcal{A}$ and $\mathcal{A}^{*}$ are defined on reweighted $L^{2}$-spaces). A more detailed derivation can be found in Ref. BK17:coherent, , where the operator $\mathcal{A}^{*}\mathcal{A}$ (or a trajectory-averaged version thereof) is used to detect coherent sets. We want to show that this is, up to regularization, equivalent to the operator in (6).

Proposition III.6**.**

Assuming that $\mathcal{A}f\in\mathbb{H}_{\mathit{\scriptscriptstyle Y}}$ for all $f\in\mathbb{H}_{\mathit{\scriptscriptstyle X}}$, it holds that $\mathcal{C}_{\mathit{\scriptscriptstyle YY}}\mathcal{A}f=\mathcal{C}_{\mathit{\scriptscriptstyle YX}}f$.

Proof.

The proof is almost identical to the proof for the standard Perron–Frobenius operator (see Ref. KSM17, ). For all $g\in\mathbb{H}_{\mathit{\scriptscriptstyle Y}}$, we obtain

[TABLE]

We define the RKHS approximation of the operator $\mathcal{A}$ by $\mathcal{A}_{k}=(\mathcal{C}_{\mathit{\scriptscriptstyle YY}}+\varepsilon\mathcal{I})^{-1}\mathcal{C}_{\mathit{\scriptscriptstyle YX}}$. Note that the operator technically depends not only on $k$ but also on $l$, which we omit for brevity. In practice, we typically use the same kernel for $\mathbb{X}$ and $\mathbb{Y}$. As a result, the eigenvalue problem (6) can now be written as

[TABLE]

The adjointness property for $\varepsilon=0$, i.e., assuming that the inverse exists without regularization,​∥∥∥Conditions for the existence of the inverse can be found, for instance, in Ref. Song2013, and in Section III.2. can be verified as follows:

[TABLE]

We have thus shown that the eigenvalue problem for the computation of coherent sets and the CCA eigenvalue problem are equivalent, provided that the RKHS is an invariant subspace of $\mathcal{T}_{k}$. Although this is in general not the case—depending on the kernel the RKHS might be low-dimensional (e.g., for a polynomial kernel), but could also be infinite-dimensional and isometrically isomorphic to $L^{2}$ (e.g., for the Gaussian kernel)—, we can use the kernel-based formulation as an approximation and solve it numerically to obtain coherent sets. This is the mathematical justification for the claim that CCA detects coherent sets, which will be corroborated by numerical results in Section IV.

III.5 Coherent mode decomposition

Borrowing ideas from dynamic mode decomposition (DMD) Schmid10 ; TRLBK14 , we now introduce a method that approximates eigenfunctions or eigenmodes of the forward-backward dynamics using linear basis functions and refer to it as coherent mode decomposition (CMD)—a mixture of CCA and DMD.​******In fact, the method described below is closer to TICA than DMD, but other variants can be derived in the same fashion, using different combinations of covariance and cross-covariance operators. The relationships between DMD and TICA (including their extensions) and transfer operators are delineated in Refs. KNKWKSN18, ; KSM17, . DMD is often used for finding coherent structures in fluid flows, dimensionality reduction, and also prediction and control; see Ref. KBBP16, for an exhaustive analysis and potential applications.

Let us assume we have high-dimensional time-series data but only relatively few snapshots. That is, $\mathbf{X},\mathbf{Y}\in\mathbb{R}^{d\times n}$ with $d\gg n$, where $\mathbf{X}=[x_{1},\dots,x_{n}]$ and $\mathbf{Y}=[y_{1},\dots,y_{n}]$. This is, for instance, the case for fluid dynamics applications where the, e.g., two- or three-dimensional domain is discretized using (un)structured grids. It is important to note that this analysis is now not based on Lagrangian data as before, where we tracked the positions of particles or drifters over time, but on the Eulerian frame of reference.

Using Algorithm III.5 with $\phi(x)=x$ and $\psi(y)=y$ is infeasible here since the resulting covariance and cross-covariance matrices would be prohibitively large; thus, we apply the kernel-based counterpart. The linear kernel $k\colon\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$ is defined by $k(x,x^{\prime})=\phi(x)^{\top}\phi(x^{\prime})=x^{\top}x^{\prime}$ and the Gram matrices are simply given by

[TABLE]

where $G_{\mathit{\scriptscriptstyle XX}},G_{\mathit{\scriptscriptstyle YY}}\in\mathbb{R}^{n\times n}$.

Algorithm III.7**.**

Coherent mode decomposition.

Choose regularization $\varepsilon$.

Compute Gram matrices $G_{\mathit{\scriptscriptstyle XX}}$ and $G_{\mathit{\scriptscriptstyle YY}}$.

*Solve the eigenvalue problem *

$(G_{\mathit{\scriptscriptstyle XX}}+n\varepsilon I)^{-1}(G_{\mathit{\scriptscriptstyle YY}}+n\varepsilon I)^{-1}G_{\mathit{\scriptscriptstyle YY}}\hskip 1.00006ptG_{\mathit{\scriptscriptstyle XX}}\hskip 1.00006ptv=\rho^{2}\hskip 1.00006ptv$.

The eigenfunction $\widehat{f}$ evaluated in an arbitrary point $x\in\mathbb{R}^{d}$ is then given by

[TABLE]

where we define the coherent mode $\xi$ corresponding to the eigenvalue $\rho$ by $\xi=\mathbf{X}v$. That is, $\xi$ contains the coefficients for the basis functions $\phi$. Analogously, we obtain

[TABLE]

where $w=\frac{1}{\rho}(G_{\mathit{\scriptscriptstyle YY}}+n\varepsilon I)^{-1}G_{\mathit{\scriptscriptstyle XX}}v$ and $\eta=\mathbf{Y}w$.

As mentioned above, DMD (as a special case of EDMD WKR15 ) typically uses the pseudoinverse to compute matrix representations of the corresponding operators. Nonetheless, a Tikhonov-regularized variant is described in Ref. EMBK17:DMD, .

IV Numerical results

As we have shown above, many dimensionality reduction techniques or methods to analyze high-dimensional data can be regarded as eigendecompositions of certain empirical RKHS operators. We now seek to illustrate how kernel CCA results in coherent sets and potential applications of the coherent mode decomposition.

IV.1 Coherent sets

We will first apply the method to a well-known benchmark problem, namely the Bickley jet, and then to ocean data and a molecular dynamics problem.

IV.1.1 Bickley jet

Let us consider a perturbed Bickley jet, which is an approximation of an idealized stratospheric flow Rypina07:coherent and a typical benchmark problem for detecting coherent sets (see, e.g., Refs. HKTH16:coherent, ; BK17:coherent, ; HSD18, ; FJ18:coherent, ). The flow is illustrated in Figure 3. For a detailed description of the model and its parameters, we refer to Ref. BK17:coherent, . Here, the state space is defined to be periodic in the $x_{1}$-direction with period $20$. In order to demonstrate the notion of coherence, we arbitrarily color one circular set yellow and one red and observe their evolution. The yellow set is dispersed quickly by the flow; the red set, on the other hand, moves around but barely changes shape. The red set is hence called coherent.

We generate 10000 uniformly distributed test points $x_{i}$ in $\mathbb{X}=[0,20]\times[-3,3]$ and then simulate their progression in time. For the computation of the coherent sets, we use only the start and end points of each trajectory, i.e., we define $y_{i}=\Theta^{\tau}(x_{i})$, where $\Theta^{\tau}$ denotes the flow associated with the dynamical system. We set $\tau=40$. From the vectors $x_{i}$ and $y_{i}$, we then compute the Gram matrices $G_{\mathit{\scriptscriptstyle XX}}$ and $G_{\mathit{\scriptscriptstyle YY}}$ using the same Gaussian kernel. Here, we define the bandwidth to be $\sigma=1$ and the regularization parameter to be $\varepsilon=10^{-7}$.

A few dominant eigenfunctions are shown in Figure 4 (a)–(d). The first eigenfunction distinguishes between the top and bottom “half” and the second one between the middle part and the rest. The subsequent eigenfunctions pick up combinations of the vortices. Applying $k$-means with $k=9$ to the first $9$ eigenfunctions results in the coherent sets shown in Figure 4 (e). This is consistent with the results presented in Ref. BK17:coherent, as shown in Figure 4 (f), where we apply space-time diffusion maps to the trajectory data (comprising 40 snapshots). While the results are qualitatively the same although kernel CCA uses only two snapshots, the coherent sets computed by our approach are less noisy, which might be due to the smoothing effects of the Gaussian kernel.

Choosing a finite-dimensional feature space explicitly, as described in Section III.2, by selecting a set of radial basis functions whose centers are given by a regular grid leads to comparable results. Currently, only start and end points of trajectories are considered. As a result, points that drift apart and then reunite at time $\tau$ would constitute coherent sets. Applying kernel CCA to less well-behaved systems might require more sophisticated kernels that take entire trajectories into account, e.g., by employing averaging techniques as suggested in Ref. BK17:coherent, .

IV.1.2 Ocean data

Ocean currents are driven by winds and tides, as well as differences in salinity. There are five major gyres as illustrated in Figure 5 (a), which has been reproduced with permission of the National Ocean Service (NOAA).​††††††NOAA. What is a gyre? https://oceanservice.noaa.gov/facts/gyre.html Our goal now is to detect these gyres from virtual buoy trajectories. In order to generate Lagrangian data, we use the OceanParcels toolbox‡‡‡‡‡‡OceanParcels project: http://oceanparcels.org/ (see Ref. OceanParcels17, for details) and data from the GlobCurrent repository,​******GlobCurrent data repository: http://www.globcurrent.org/ provided by the European Space Agency. More precisely, our drifter computations are based on the Eulerian total current at significant wave height from the sum of geostrophic and Ekman current components, starting on the 1st of January 2016 and ending on the 31st of December 2016 with 3-hourly updates.

We place 15000 uniformly distributed virtual drifters in the oceans and let the flow evolve for one year, which thus constitutes the lag time $\tau$. Let $x_{i}$ denote the initial positions and $y_{i}$ the new positions of the drifters after one year. The domain is $\mathbb{X}=[-180^{\circ},180^{\circ}]\times[-80^{\circ},80^{\circ}]$, where the first dimension corresponds to the longitudes and the second to the latitudes. For the coherent set analysis, we select a Gaussian kernel $k(x,x^{\prime})=\exp\left(-\frac{d(x,x^{\prime})^{2}}{2\sigma^{2}}\right)$ with bandwidth $\sigma=30$, where $d(x,x^{\prime})$ is the distance between the points $x$ and $x^{\prime}$ in kilometers computed with the aid of the haversine formula. The regularization parameter $\varepsilon$ is set to $10^{-4}$. The first two dominant eigenfunctions computed using kernel CCA are shown in Figure 5 (b) and (c) and a $k$-means clustering of the six dominant eigenfunctions in Figure 5 (d). CCA correctly detects the main gyres—the splitting of the South Atlantic Gyre and the Indian Ocean Gyre might be encoded in eigenfunctions associated with smaller eigenvalues—and the Antartic Circumpolar Current. The clusters, however, depend strongly on the lag time $\tau$. In order to illustrate the flow properties, typical trajectories are shown in Figure 5 (e). The trajectories belonging to different coherent sets remain mostly separated, although weak mixing can be seen, for instance, at the borders between the red and purple and red and green clusters.

IV.1.3 Time-dependent energy potential

As a last example, we will analyze a molecular-dynamics inspired problem, namely diffusion in a time-dependent two-dimensional energy landscape, given by the stochastic differential equation

[TABLE]

with

[TABLE]

The parameter $\beta$ is the dimensionless inverse (absolute) temperature, $W_{t}$ a standard Wiener process, and $s$ specifies the number of wells. This is a generalization of a potential defined in Ref. BKKBDS18, , whose wells now move periodically towards and away from the center and which furthermore slowly rotates. We set $s=5$. The resulting potential for $t=0$ is shown in Figure 6 (a). Particles will typically quickly equilibrate in radial direction towards the closest well and stay in this well, which moves over time. Particles trapped in one well will remain coherent for a relatively long time. The probability of escaping and moving to another one depends on the inverse temperature: The higher $\beta$, the less likely are transitions between wells.

We generate 1000 uniformly distributed test points in $\mathbb{X}=[-2.5,2.5]\times[-2.5,2.5]$ and integrate the system with the aid of the Euler–Maruyama method and the step size $h=10^{-3}$ from $t=0$ to $t=10$. As before, we use only the start and end points of the trajectories and a Gaussian kernel (here, $\sigma=1$ and $\varepsilon=10^{-6}$) for the coherent set analysis.

Due to the centering of the Gram matrices, the eigenvalue $\lambda=1$ vanishes and—depending on the parameter $\beta$—four eigenvalues close to one remain as illustrated in Figure 6 (b). Figure 6 (c) shows a clustering of the dominant four eigenfunctions for $\beta=3$ based on PCCA+ Roeblitz2013 , resulting in the expected five coherent sets. The clustering at $t=10$ (see Figure 6 (d)) illustrates that the computed sets indeed remain largely coherent.

Standard methods for the computation of metastable sets such as Ulam’s method, EDMD, or their variants are in general not suitable for non-equilibrium dynamics; see also Ref. KWNS18:noneq, and Section III.3.

IV.2 Coherent mode decomposition

In order to illustrate the coherent mode decomposition outlined in Algorithm III.7, we consider the classical von Kármán vortex street and generate data using a simple Python implementation.​*†††Palabos project: http://wiki.palabos.org/numerics:codes It is important to note that here we take into account the full trajectory data $\{z_{0},\dots,z_{n}\}$, where $z_{i}$ is the state at time $t=20\hskip 1.00006pti$, and define $X=[z_{0},\dots,z_{n-1}]$ and $Y=[z_{1},\dots,z_{n}]$, whereas we generated uniformly distributed data for the coherent set analysis in the previous subsection and furthermore used only the start and end points of the trajectories. We set $n=100$ and $\varepsilon=0.1$. Some snapshots of the system are shown in Figure 7 (a)–(d). Applying CMD results in the modes depicted in Figure 7 (e)–(h), where the color bar is the same as in Figure 4. As described above, we obtain two modes, denoted by $\xi$ and $\eta$, for each eigenvalue $\rho$, where $\eta$ can be interpreted as the time-lagged counterpart of $\xi$.

For this standard DMD benchmark problem, which we chose for illustration purposes, CMD and (regularized) DMD lead to modes that look highly similar. The interpretations, however, are different. While the DMD modes, which correspond to eigenvectors, are objects that are mapped onto scalar multiples of themselves, the CMD modes, which correspond to singular vectors, encode information about how coherent structures at time $t$ are transported by the flow to time $t+\tau$. In fact, the DMD eigenvalues associated with the DMD modes resembling the CMD modes shown in Figure 7 are negative and close to $-1$, implying periodicity. Analogously, the CMD modes $\eta$ are akin to $-\xi$, which also implies periodic motion. Further applications of CMD pertaining to, for instance, more complicated fluid flows or also non-sequential data, will be investigated in future work.

V Conclusion

We demonstrated that several kernel-based dimensionality reduction techniques can be interpreted as eigendecompositions of empirical estimates of certain RKHS operators. Moreover, we showed that applying CCA to Lagrangian data results in coherent sets and illustrated the efficiency of the methods using several examples ranging from fluid to molecular dynamics. This approach worked out of the box, although taking into account entire trajectories might improve the results even further, which would then necessitate dedicated kernels. In this work, we analyzed only low-dimensional benchmark problems. Nevertheless, the kernel-based algorithms can be easily applied to more complex problems and also non-vectorial domains such as graphs or strings.

As a byproduct of the coherent set analysis, we derived a method called CMD that is a hybrid of CCA and DMD (or TICA). This method can, for instance, be applied to high-dimensional fluid flow or video data. For specific problems, CMD and DMD—unsurprisingly, given the close proximity—result in highly similar modes. An interesting topic for future research would be to systematically analyze the relationships between these methods. Furthermore, as with the transfer operators and embedded transfer operators as well as their kernel-based estimates KSM17 , there are again several different combinations and variants of the proposed algorithms.

Another open problem is the influence of different regularization techniques on the numerical results. How does Tikhonov regularization compare to approaches based on pseudoinverses or other spectral filtering methods? And how do we choose the kernel and the regularization parameters in an optimal way, preferably without cross-validation? Additionally, future work includes analyzing the properties of the empirical estimate $\widehat{S}$. Can we show convergence in the infinite-data limit? Which operators can be approximated by $\widehat{S}$ and can we derive error bounds for the resulting eigenvalues and eigenfunctions?

We expect the results in this paper to be a starting point for further theoretical research into how RKHS operators in the context of dynamical systems could be approximated and, furthermore, how they connect to statistical learning theory. Additionally, the methods proposed here might be combined with classical modifications of CCA in order to improve the numerical performance. The experiments here were performed using Matlab, and the methods have been partially reimplemented in Python and are available at https://github.com/sklus/d3s/.

Acknowledgements

We would like to thank Péter Koltai for the Bickley jet implementation as well as helpful discussions related to coherent sets, Ingmar Schuster for pointing out similarities between CCA and kernel transfer operators, and the reviewers for many helpful suggestions for improvements. We gratefully acknowledge funding from Deutsche Forschungsgemeinschaft (DFG) through grant CRC 1114 (Scaling Cascades in Complex Systems, project ID: 235221301, projects A04 and B06), Germany’s Excellence Strategy (MATH+: The Berlin Mathematics Research Center, EXC-2046/1, project ID: 390685689, projects AA1-2, AA1-6, and EF1-2), and European Commission through ERC CoG 772230 “ScaleCell”.