Variational approach for learning Markov processes from time series data

Introduction

Extracting dynamical models and their main characteristics from time series data is a recurring problem in many areas of science and engineering. In the particularly popular approach of Markovian models, the future evolution of the system, e.g. state $\mathbf{x}_{t+\tau}$, only depends on the current state $\mathbf{x}_{t}$, where $t$ is the time step and $\tau$ is the delay or lag time. Markovian models are easier to analyze than models with explicit memory terms. They are justified by the fact that many physical processes – including both deterministic and stochastic processes – are inherently Markovian. Even when only a subset of the variables in which the system is Markovian are observed, a variety of physics and engineering processes have been shown to be accurately modeled by Markovian models on sufficiently long observation lag times $\tau$. Examples include molecular dynamics ChoderaNoe_COSB14_MSMs ; PrinzEtAl_JCP10_MSM1 , wireless communications konrad2001markov ; ma2001composite and fluid dynamics mezic2013analysis ; froyland2016optimal .

In the past decades, a collection of closely related Markov modeling methods were developed in different fields, including Markov state models (MSMs) SchuetteFischerHuisingaDeuflhard_JCompPhys151_146 ; PrinzEtAl_JCP10_MSM1 ; BowmanPandeNoe_MSMBook , Markov transition models WuNoe_JCP15_GMTM , Ulam’s Galerkin method dellnitz2001algorithms ; bollt2013applied ; froyland2014computational , blind-source separation Molgedey_94 ; ZieheMueller_ICANN98_TDSEP , the variational approach of conformation dynamics (VAC) NoeNueske_MMS13_VariationalApproach ; NueskeEtAl_JCTC14_Variational , time-lagged independent component analysis (TICA) PerezEtAl_JCP13_TICA ; SchwantesPande_JCTC13_TICA , dynamic mode decomposition (DMD) RowleyEtAl_JFM09_DMDSpectral ; Schmid_JFM10_DMD ; TuEtAl_JCD14_ExactDMD , extended dynamic mode decomposition (EDMD) WilliamsKevrekidisRowley_JNS15_EDMD , variational Koopman models WuEtAl_JCP17_VariationalKoopman , variational diffusion maps BoninsegnaEtAl_JCTC15_VariationalDM , sparse identification of nonlinear dynamics brunton2016discovering and corresponding kernel embeddings harmeling2003kernel ; song2013kernel ; SchwantesPande_JCTC15_kTICA and tensor formulations NueskeEtAl_JCP15_Tensor ; KlusSchuette_Arxiv15_Tensor . All these models approximate the Markov dynamics at a lag time $\tau$ by a linear model in the following form:

A direct method to estimate the matrix $\mathbf{K}$ from data is to solve the linear regression problem $\mathbf{g}(\mathbf{x}_{t+\tau})\approx\mathbf{K}^{\top}\mathbf{f}(\mathbf{x}_{t})$, which facilitates the use of regularized solution methods, such as the LASSO method tibshirani1996regression . Alternatively, feature functions $\mathbf{f}$ and $\mathbf{g}$ that allow Eq. (1) to have a probabilistic interpretation (e.g. in MSMs), $\mathbf{K}$ can be estimated by a maximum-likelihood or Bayesian methods PrinzEtAl_JCP10_MSM1 ; Noe_JCP08_TSampling .

However, as yet, it is still unclear what are the optimal choices for $\mathbf{f}$ and $\mathbf{g}$ - either given a fixed dimension or a fixed amount of data. Notice that this problem cannot be solved by minimizing the regression error of Eq. (1), because a regression error of zero can be trivially achieved by choosing a completely uninformative model with $\mathbf{f}\left(\mathbf{x}\right)\equiv\mathbf{g}\left(\mathbf{x}\right)\equiv 1$ and $\mathbf{K}=1$. An approach that can be applied to deterministic systems and for stochastic systems with additive white noise is to set $\mathbf{g}\left(\mathbf{x}\right)=\mathbf{x}$, and then choose $\mathbf{f}$ as the transformation with smallest modeling error brunton2016discovering ; brunton2016koopman .

A more general approach is to optimize the dominant spectrum of the Koopman operator. At long timescales, the dynamics of the system are usually dominated by the Koopman eigenfunctions of the Koopman operator with large eigenvalues. If the dynamics obey detailed balance, those eigenvalues are real-valued, and the variational approach for reversible Markov processes can be applied that has made great progress in the field of molecular dynamics NoeNueske_MMS13_VariationalApproach ; NueskeEtAl_JCTC14_Variational . In such processes, the smallest modeling error of (1) is achieved by setting $\mathbf{f}=\mathbf{g}$ equal to the corresponding eigenfunctions. Ref. NoeNueske_MMS13_VariationalApproach describes a general approach to approximate the unknown eigenfunction from time series data of a reversible Markov process: Given a set of orthogonal candidate functions, $\mathbf{f}$, it can be shown that their time-autocorrelations are lower bounds to the corresponding Koopman eigenvalues, and are equal to them exactly if, and only if $\mathbf{f}$ are equal to the Koopman eigenfunctions. This approach provides a variational score, such as the sum of estimated eigenvalues (the Rayleigh trace), that can be optimized to approximate the eigenfunctions. If $\mathbf{f}$ is defined by a linear superposition of a given set of basis functions, then the optimal coefficients are found equivalently by either maximizing the variational score, or minimizing the regression error in the feature space as done in EDMD WilliamsKevrekidisRowley_JNS15_EDMD – see WuEtAl_JCP17_VariationalKoopman . However, the regression error cannot be used to select the form and the number of basis functions themselves, whereas the variational score can. When working with a finite dataset, however, it is important to avoid overfitting, and to this end a cross-validation method has been proposed to compute variational scores that take the statistical error into account McGibbonPande_JCP15_CrossValidation . Such cross-validated variational scores can be used to determine the size and type of the function classes and the other hyper-parameters of the dynamical model.

While this approach is extremely powerful for stationary and data and reversible Markov processes, almost all real-world dynamical processes and time-series thereof are irreversible and often even non-stationary. In this paper, we introduce a variational approach for Markov processes (VAMP) that can be employed to optimize parameters and hyper-parameters of arbitrary Markov processes. VAMP is based on the singular value decomposition of the Koopman operator, which overcomes the limited usefulness of the eigenvalue decomposition of time-irreversible and non-stationary processes. We first show that the approximation error of the Koopman operator deduced from the linear model (1) can be minimized by setting $\mathbf{f}$ and $\mathbf{g}$ to be the top left and right singular functions of the Koopman operator. Then, by using the variational description of singular components, a class of variational scores, VAMP-$r$ for $r=1,2,\ldots$, are proposed to measure the similarity between the estimated singular functions and the true ones. Maximization of any of these variational scores leads to optimal model parameters and is algorithmically identical to Canonical Correlation Analysis (CCA) between the featurized time-lagged pair of variables $\mathbf{x}_{t}$ and $\mathbf{x}_{t+\tau}$. This approach can also be employed to learn the feature transformations by nonlinear function approximators, such as deep neural networks. Furthermore, we establish a relationship between the VAMP-2 score and the approximation error of the dynamical model with respect to the true Koopman operator. We show that this approximation error can be practically computed up to a constant, and define its negative as the VAMP-E score. Finally, we demonstrate that optimizing the VAMP-E score in a cross-validation framework leads to an optimal choice of hyperparameters.

Theory

The Koopman operator $\mathcal{K}_{\tau}$ of a Markov process is a linear operator defined by

where the singular value $\sigma_{i}>0$ is the square root of the $i$th largest eigenvalue of $\mathcal{K}_{\tau}^{*}\mathcal{K}_{\tau}$ or $\mathcal{K}_{\tau}\mathcal{K}_{\tau}^{*}$, the left and right singular function $\psi_{i},\phi_{i}$ are the $i$th eigenfunctions of $\mathcal{K}_{\tau}^{*}\mathcal{K}_{\tau}$ and $\mathcal{K}_{\tau}\mathcal{K}_{\tau}^{*}$ with

Consider a one-dimensional dynamical system

evolving in the state space $$, where$u_{t}$is a standard Gaussian white noise zero mean and unit variance (see Appendix K.1 for details on the numerical simulations and analysis). This system has two metastable states with the boundary close to$x=0$as shown in Fig. 1a, and the singular components are summarized in Figs. 1c and 1d. As shown in the figures, the sign structures of the second left and right singular functions clearly indicate the metastable states, and the third and forth singular functions provide more detailed information on the dynamics. An accurate estimate of the transition density can be obtained by combining the first four singular components and the corresponding relative approximation error of the Koopman operator is only$6.6\%$(see Figs. 1b and 1e). In addition, we utilize the finite-rank approximate Koopman operators to predict the time evolution of the distribution of$x_{t}$for$t=1,\ldots,256$with the initial state$x_{0}=12$, and a small error can also be achieved when the rank is only$4$ as displayed in Fig. 2, where

There are other formalisms to describe Markovian dynamics, for example, the Markov propagator or the weighted Markov propagator, also called transfer operator SchuetteFischerHuisingaDeuflhard_JCompPhys151_146 . These propagators are commonly used for modeling physical processes such as molecular dynamics, and describe the evolution of probability densities instead of observables. We show in Appendix B that all conclusions in this paper can be equivalently established by interpreting $(\sigma_{i},\rho_{1}\phi_{i},\rho_{0}\psi_{i})$ as the singular components of the Markov propagator.

2 Variational principle for Markov processes

In order to allow the optimal model (4) to be estimated from data, we develop a variational principle for the approximation of singular values and singular functions of Markov processes.

According to the Rayleigh variational principle of singular values, the first singular component maximizes the generalized Rayleigh quotient of $\mathcal{K}_{\tau}$ as

and the maximal value of the generalized Rayleigh quotient is equal to the first singular value $\sigma_{1}=\left\langle\psi_{1},\,\mathcal{K}_{\tau}\phi_{1}\right\rangle_{\rho_{0}}$. For the $i$th singular component with $i>1$, we have

and the maximal value is equal to $\sigma_{i}=\left\langle\psi_{i},\,\mathcal{K}_{\tau}\phi_{i}\right\rangle_{\rho_{0}}$. These insights can be summarized by the following variational theorem for seeking all top $k$ singular components simultaneously:

VAMP variational principle. The $k$ dominant singular components of a Koopman operator are the solution of the following maximization problem:

where $r\geq 1$ can be any positive integer. The maximal value is achieved by the singular functions $f_{i}=\psi_{i}$ and $g_{i}=\phi_{i}$ and

is called the VAMP-r score of $\mathbf{f}$ and $\mathbf{g}$.

3 Comparison with related analysis approaches

The SVD of the Koopman operator is equivalent to the eigenvalue decomposition when the Markov process is time-reversible and stationary with $\rho_{0}=\rho_{1}$, and therefore the variational principle presented here is a generalization of that developed for reversible conformation dynamics NoeNueske_MMS13_VariationalApproach ; NueskeEtAl_JCTC14_Variational . Specifically, VAMP-1 maximizes the Rayleigh trace, i.e. the sum of the estimated eigenvalues NoeNueske_MMS13_VariationalApproach ; McGibbonPande_JCP15_CrossValidation , and VAMP-2 maximizes the kinetic variance introduced in NoeClementi_JCTC15_KineticMap . See Appendix D for a detailed derivation of the reversible variational principle from the VAMP variational principle. For irreversible Markov processes, the singular functions can provide low-dimensional embeddings of kinetic distances between states like eigenfunctions of reversible processes Paul2018identification . Furthermore, the coherent sets of nonstationary Markov processes, which are the generalization of metastable states, can be identified from dominant singular functions koltai2018optimal .

The dynamics of an irreversible Markov process can also be analyzed through solving the eigenvalue problem $\mathcal{K}_{\tau}g=\lambda g$ (see, e.g., WilliamsKevrekidisRowley_JNS15_EDMD ; williams2014kernel ; KlusKoltaiSchuette_ApproximationKoopman ; KlusSchuette_Arxiv15_Tensor ), and the eigenfunctions form an invariant subspace of the Koopman operator for multiple lag times since the eigenvalue problem satisfies

However, as far as we know, there is no variational principle for approximate eigenfunctions of irreversible Markov processes, and it is difficult to evaluate errors of projections of Koopman operators to the invariant subspaces. The SVD based analysis approach overcomes the above problems and yields the optimal finite-rank approximate models. The major limitation of this approach comes from the fact that the singular functions are dependent on the choice of the lag time and the optimality of model (4) holds only for a fixed $\tau$. The optimization and error analysis of Koopman models for multiple lag times will be studied in our future work.

Estimation algorithms

We introduce algorithms to estimate optimal dynamical models from time series data. We make the Ansatz to represent the feature functions $\mathbf{f}$ and $\mathbf{g}$ as linear combinations of basis functions $\boldsymbol{\chi}_{0}=(\chi_{0,1},\chi_{0,2},\ldots)^{\top}$ and $\boldsymbol{\chi}_{1}=(\chi_{1,1},\chi_{1,2},\ldots)^{\top}$:

Here, $\mathbf{U}$ and $\mathbf{V}$ are matrices of size $m\times k$ and $m^{\prime}\times k$, i.e. we are trying to approximate $k$ singular components by linearly combining $m$ and $m^{\prime}$ basis functions. For the sake of generality we have assumed that $\mathbf{f}$ and $\mathbf{g}$ are represented by different basis sets. However, in practice one can justify using a single basis set the joint set $\boldsymbol{\chi}^{\top}=(\boldsymbol{\chi}_{0}^{\top},\boldsymbol{\chi}_{1}^{\top})$ as an Ansatz for both $\mathbf{f}$ and $\mathbf{g}$. Please note that despite the linear Ansatz (15), the feature functions may be strongly nonlinear in the system’s state variables $\mathbf{x}$, thus we are not restricting the generality of the functions $\mathbf{f}$ and $\mathbf{g}$ that can be represented. In this section, we consider three problems: (i) optimizing $\mathbf{U}$ and $\mathbf{V}$, (ii) optimizing $\boldsymbol{\chi}_{0}$ and $\boldsymbol{\chi}_{1}$ and (iii) assessing the quality of the resulting dynamical model.

For convenience of notation, we denote by $\mathbf{C}_{00},\mathbf{C}_{11},\mathbf{C}_{01}$ the covariance matrices and time-lagged covariance matrices of basis functions, which can be computed from a trajectory $\{x_{1},\ldots,x_{T}\}$ by

If there are multiple trajectories, the covariance matrices can be computed in the same manner by averaging over all trajectories. Instead of the direct estimators (16-18), more elaborated estimation methods such as regularization methods tibshirani1996regression and reweighting estimators WuEtAl_JCP17_VariationalKoopman may be used.

We first propose a solution for the problem of finding the optimal parameter matrices $\mathbf{U}$ and $\mathbf{V}$ given that the basis functions $\boldsymbol{\chi}_{0}$ and $\boldsymbol{\chi}_{1}$ are known. Substituting the linear Ansatz (15) into the VAMP variational principle, shows that $\mathbf{U}$ and $\mathbf{V}$ can be computed as the solutions of the maximization problem:

is a matrix representation of VAMP-$r$ score, and $\mathbf{u}_{i}$ and $\mathbf{v}_{i}$ are the $i$th columns of $\mathbf{U}$ and $\mathbf{V}$. This problem can be solved by applying linear CCA hardoon2004canonical in the feature spaces defined by the basis sets $\boldsymbol{\chi}_{0}(\mathbf{x}_{t})$ and $\boldsymbol{\chi}_{1}(\mathbf{x}_{t+\tau})$, and the same solution will be obtained for any other choice of $r$. (See Appendices E.1 and E.2 for more detailed proof and analysis.) The resulting algorithm for finding the best linear model is a CCA in feature space, applied on time-lagged data. Hence we briefly call this algorithm feature TCCA:

Compute covariance matrices $\mathbf{C}_{00},\mathbf{C}_{01},\mathbf{C}_{11}$ via (16-18).

Compute $\mathbf{U}=\mathbf{C}_{00}^{-\frac{1}{2}}\mathbf{U}^{\prime}$ and $\mathbf{V}=\mathbf{C}_{11}^{-\frac{1}{2}}\mathbf{V}^{\prime}$.

Output the linear model (1) with $K_{ii}$, $f_{i}=\mathbf{u}_{i}^{\top}\boldsymbol{\chi}_{0}$ and $g_{i}=\mathbf{v}_{i}^{\top}\boldsymbol{\chi}_{1}$ being the estimates of the $i$th singular value, left singular function and right singular function of the Koopman operator.

is equal to the least square solution to the regression problem $\boldsymbol{\chi}_{1}\left(\mathbf{x}_{t+\tau}\right)\approx\mathbf{K}_{\chi}^{\top}\boldsymbol{\chi}_{0}\left(\mathbf{x}_{t}\right)$. Note that if we further assume that $\boldsymbol{\chi}_{0}=\boldsymbol{\chi}_{1}$, (21) is identical to the linear model of EDMD. Thus, the feature TCCA can be seen as a generalization of EDMD that can provide approximate Markov models for different basis $\boldsymbol{\chi}_{0}$ and $\boldsymbol{\chi}_{1}$. More discussion on the relationship between the two methods is provided in Appendix G.

2 Nonlinear TCCA: optimizing the basis functions

We now extend feature TCCA to a more flexible representation of the transformation functions $\mathbf{f}$ and $\mathbf{g}$ by optimizing the basis functions themselves:

Here, $\mathbf{w}$ represents a set of parameters that determines the form of the basis functions. As a simple example, consider $\mathbf{w}$ to represent the mean vectors and covariance matrices of a Gaussian basis set. However, $\boldsymbol{\chi}_{0}\left(\mathbf{x};\mathbf{w}\right)$ and $\boldsymbol{\chi}_{1}\left(\mathbf{x};\mathbf{w}\right)$ can also represent very complex and nonlinear learning structures, such as neural networks and decision trees.

Compute $\mathbf{w}^{*}=\arg\max_{\mathbf{w}}\left\|\mathbf{C}_{00}\left(\mathbf{w}\right)^{-\frac{1}{2}}\mathbf{C}_{01}\left(\mathbf{w}\right)\mathbf{C}_{11}\left(\mathbf{w}\right)^{-\frac{1}{2}}\right\|_{r}^{r}$ by gradient descent or other nonlinear optimization methods.

Approximate the Koopman singular values and singular functions using the feature TCCA algorithm with basis sets $\boldsymbol{\chi}_{0}\left(\mathbf{x};\mathbf{w}^{*}\right)$ and $\boldsymbol{\chi}_{1}\left(\mathbf{x};\mathbf{w}^{*}\right)$.

Unlike the estimated singular components generated by the feature TCCA, the estimation results of the nonlinear TCCA do generally depend on the value of $r$. (An example is given in Appendix E.3, where the VAMP scores can be analytically computed.) We suggest to set $r=2$ in applications for the direct relationship between the VAMP-2 score and the approximation error of Koopman operators and the convenience of cross-validation (see below). The details of the nonlinear TCCA, including the optimization algorithm and regularization, are beyond the scope of this paper. Appendix F.2 provides a brief description of the implementation, and related work based on kernel methods and deep networks can be found in andrew2013deep ; vampnet .

Let us consider the stochastic system described in Example 1 again. We generate $10$ simulation trajectories of length $500$, and approximate the dominant singular components by the feature TCCA. Here, the basis functions are

which define a partition of the domain $$into$m=33$disjoint intervals. In other words, the approximation is performed based on an MSM with$33$ discrete states. Estimation results are given in Fig. 3a, where the discretization errors arising from indicator basis functions are clearly shown. For comparison, we also implement the nonlinear TCCA algorithm with radial basis functions

with smoothing parameter $w\geq 0$, where $c_{i}=\frac{40\cdot(i-0.5)}{m}-20$ for $i=1,\ldots,m$ are uniformly distributed in $$. Notice that the basis functions given in (25) are a specific case of the radial basis functions with$w=\infty$, and it is therefore possible to achieve better approximation by optimizing$w$. As can be seen from Fig. 3b, the nonlinear TCCA provides more accurate estimates of singular functions and singular values (see Appendix K.1 for more details). In addition, both feature TCCA and nonlinear TCCA underestimate the dominant singular values as stated by the variational principle.

The nonlinear TCCA is similar to the EDMD with dictionary learning (EDMD-DL) li2017extended , where the feature transformations are optimized by minimizing the regression error of (1). The major advantages of the nonlinear TCCA over EDMD-DL are: First, the uninformative model with zero regression error can be systematically excluded without any extra constraints on features. Second, the optimization objective is directly related to the approximation error of the Koopman operator (see Section 3.3). Some recent methods extend EDMD-DL for modeling Koopman operators of deterministic systems takeishi2017learning ; lusch2018deep ; otto2019linearly , and solve the first problem by using the prediction error between the observed $x_{t}$ and that predicted by the low-dimensional model. But they cannot be applied to stochastic Koopman operators of Markov processes directly.

3 Error analysis

According to (5), both feature TCCA and nonlinear TCCA lead to a rank $k$ approximation

to $\mathcal{K}_{\tau}$. We consider here the approximation error of (27) in a general case where $\mathbf{f}=\mathbf{U}^{\top}\boldsymbol{\chi}_{0}$ and $\mathbf{g}=\mathbf{V}^{\top}\boldsymbol{\chi}_{1}$ may not satisfy the orthonormal constraints due to statistical noise and numerical errors. After a few steps of derivation, the approximation error can be expressed as

Remarkably, this error decomposes into a unknown constant part (the square of Hilbert-Schmidt norm of $\mathcal{K}_{\tau}$), and a model-dependent part $\mathcal{R}_{E}$ that can be entirely estimated from data by its matrix representation:

$\mathcal{R}_{E}$, is thus a score that can be used alternatively to the VAMP-$r$ scores, and we call $\mathcal{R}_{E}$ VAMP-E score. It can be proved that the maximization of $\mathcal{R}_{E}$ is equivalent to maximization of $\mathcal{R}_{2}$ in feature TCCA or nonlinear TCCA. However, these scores will behave differently in terms of hyper-parameter optimization (see Sec. 4.1). Proofs and analysis are given in Appendix H.

Model validation

For a data-driven estimation of dynamical models, either using feature TCCA or nonlinear TCCA, we have to strike a balance between the modeling or discretization error and the statistical or overfitting error. The choice of number and type of basis functions is critical for both. If basis sets are very small and not flexible enough to capture singular functions, the approximation results may be inaccurate with large biases. We can improve the variational score and reduce the modeling error by larger and more flexible basis sets. But too complicated basis sets will produce unstable estimates with large statistical variances, and in particular poor predictions on data that has not been used in the estimation process – this problem is known as overfitting in the machine learning community. A popular way to achieve the balance between the statistical bias and variance are resampling methods, including bootstrap and cross-validation friedman2001elements . They iteratively fit a model in a training set, which are sampled from the data with or without replacement, and validate the model in the complementary dataset. Alternatively, there are also Bayesian hyper-parameter optimization methods. See ArlotCelisse_StatSurv10_CVReview ; snoek2012practical for an overview. Here, we will focus on cross-validation and describe how to use the VAMP scores in this and similar resampling frameworks.

Let $\boldsymbol{\theta}$ be hyper-parameters in feature TCCA or nonlinear TCCA that need to be specified. For example, $\boldsymbol{\theta}$ includes the number and functional form of basis functions used in feature TCCA, or the architecture and connectivity of a neural network used for nonlinear TCCA. Generally speaking, different values of $\boldsymbol{\theta}$ correspond to different dynamical models that we want to rank, and these models may be of completely different types. The cross-validation of $\boldsymbol{\theta}$ can be performed as follows:

For each hyper-parameter set $\boldsymbol{\theta}$:

The key to the above procedure is how to evaluate the estimated singular components for given test set. It is worth pointing out that we cannot simply define the validation score directly as the VAMP-$r$ score of estimated singular functions for the test data, because the singular functions obtained from training data are usually not orthonormal with respect to the test data.

A feasible way is to utilize the subspace variational score as proposed for reversible Markov processes in McGibbonPande_JCP15_CrossValidation . For VAMP-$r$ this score becomes:

Therefore, we can score the performance of estimated singular components on the test set by

It is worth pointing out that a validation score is proposed kurebayashi370optimal for cross-validation of kernel DMD based on the analysis of approximation error of transition densities, which has a similar form to that of VAMP-E. The theoretical and empirical comparisons between the two scores will be performed in our future work.

We consider here the choice of the basis function number $m$ for the nonlinear TCCA in Example 2. We use 5-fold cross-validation with the VAMP-E score to compare different values of $m$. While the average score computed by training sets keeps increasing with $m$, both the cross validation score and the exact VAMP-E score achieve their maximum value at $m=33$ as in Example 2 (see Fig. 4a). The optimality can also be demonstrated by comparing Fig. 3b and Fig. 4b. A much smaller basis set with $m=13$ yields large errors in the approximation of singular functions. When $m=250$, the estimation of singular functions suffers from overfitting and the estimated singular value is even larger than the true value due to the statistical noise.

2 Chapman-Kolmogorov test for choice of lag times

In order to address this problem, the Chapman-Kolmogorov test can be used, which is common in building Markov state models PrinzEtAl_JCP10_MSM1 . Let us consider the covariance

between observables $f$ and $g$ of lag time $n\tau$, which can be estimated from data as

where $\rho_{0}(n\tau)$ is the empirical distribution of the simulation data excluding $\{x_{t}|t>T-n\tau\}$. If our methods provide an ideal Markov model of lag time $\tau$, the Koopman operator $\mathcal{K}_{n\tau}$ can be approximated by $\hat{\mathcal{K}}_{\tau}^{n}$, and the covariance can also be predicted as

Therefore, the lag time $\tau$ can be selected according to the following criteria in applications: (i) The lag time is smaller than the timescale that we are interested in. (ii) The equation

holds approximately for multiple observables $f,g$ and lag times $n\tau$. In this paper, we simply set $f,g$ to be the estimated leading singular functions since they dominate the dynamics of the Markov process.

Numerical examples

Let’s consider a stochastic double-gyre system defined by:

where $\mathbf{W}_{t,1}$ and $\mathbf{W}_{t,2}$ are two independent standard Wiener processes. The dynamics are defined on the domain $\times$ with reflecting boundary. For $\varepsilon=0$, it can be seen from the flow field depicted in Fig. 5a that there is no transport between the left half and the right half of the domain and both subdomains are invariant sets with measure $\frac{1}{2}$ froyland2009almost ; froyland2014almost . For $\varepsilon>0$, there is a small amount of transport due to diffusion and the subdomains are almost invariant. Here we used the parameters $A=0.25$, $\epsilon=0.1$, and lag time $\tau=2$ in analysis and simulations. The first two nontrivial singular components are shown in Fig. 5c, where the two almost invariant sets are clearly visible in $\psi_{2},\phi_{2}$ and $\psi_{3},\phi_{3}$ are associated with the rotational kinetics within the almost invariant sets.

We generate $10$ trajectories of length $4$ with step size $0.02$, and perform modeling by nonlinear TCCA with basis functions

where $\mathbf{c}_{1},\ldots,\mathbf{c}_{m}$ are cluster centers given by k-means algorithm, and the smoothing parameter $w$ is determined via maximizing the VAMP-2 score given in (24) (see Appendix K.2 for more details of numerical computations). The size of basis set $m=37$ is selected by the VAMP-E based cross-validation proposed in 4.1 with $5$ folds (see Fig 5b), and it can be observed from Figs. 5c, 5d and 5e that the leading singular components are accurately estimated. In contrast, as shown in Figs. 5f and 5g, a much small value of $m$ leads to significant approximation errors of singular components, while for a much larger value, the estimates are obviously influenced by statistical noise. Fig. 6 illustrates that the Koopman operator estimated by nonlinear TCCA can successfully predict the long-time evolution of the distribution of the state. The Chapman-Kolmogorov test results are displayed in Fig. 7, which confirm that $\tau=2$ is a suitable choice of the lag time.

2 Stochastic Lorenz system

As the last example, we investigate the stochastic Lorenz system which obeys the following stochastic differential equation:

with parameters $s=10$, $r=28$ and $b=8/3$. The deterministic Lorenz system with $\epsilon=0$ is known to exhibit chaotic behavior sparrow1982lorenz with a strange attractor characterized by two lobes as illustrated in Fig. 8a. We generate $20$ trajectories of length $25$ with $\epsilon=0.3$ by using the Euler–Maruyama scheme with step size $0.005$, and one of them is shown in Fig. 8b. As stated in chekroun2011stochastic , all the trajectories move around the deterministic attractor with small random perturbations and switch between the two lobes.

The leading singular components computed from the simulation data by the nonlinear TCCA are summarized in Fig 8c, where the lag time $\tau=0.75$ is determined via the Chapman-Kolmogorov test (see Fig. 8d), $\boldsymbol{\chi}_{0}=\boldsymbol{\chi}_{1}$ consist of $m$ normalized radial basis functions similar to those used in Section 5.1, and the selection of $m$ is also implemented by $5$-fold cross-validation. According to the patterns of the singular functions, the stochastic Lorenz system can be coarse-grained into a simplified model which transitions between four macrostates corresponding to inner and outer basins of the two attractor lobes. In particular, the sign-boundary of $\psi_{1}$ closely matches that between the almost invariant sets of the Lorenz flow froyland2009almost .

Next, we map the simulation data to a higher dimensional space via the nonlinear transformation $\boldsymbol{\eta}_{t}=\boldsymbol{\eta}(x_{t},y_{t},z_{t})$ defined by

Fig. 9a plots the transformed points of the illustrative trajectory in Fig. 8b. We utilize the nonlinear TCCA to compute the singular components in the space of $\boldsymbol{\eta}_{t}=(\eta_{t}^{1},\ldots,\eta_{t}^{6})$ by assuming that the available observable is $\boldsymbol{\eta}_{t}$ instead of $(x_{t},y_{t},z_{t})$, show in Fig. 9b the projections of the singular functions back on the three-dimensional space

Conclusion

The linearized coarse-grained models of Markov systems are commonly used in a broad range of fields, such as power systems, fluid mechanics and molecular dynamics. Although the models were developed independently in different communities, the VAMP proposed in this paper provides a general framework for analysis of them, and the modeling accuracy can be quantitatively evaluated by the VAMP-$r$ and VAMP-E scores. Moreover, a set of data-driven methods, including feature TCCA, nonlinear TCCA and VAMP-E based cross-validation, are developed to achieve optimal modeling for given finite model dimensions and finite data sets.

The major challenge in real-world applications of VAMP is how to overcome the curse of dimensionality and solve the variational problem effectively and efficiently for high-dimensional systems. One feasible way of addressing this challenge is to approximate singular components by deep neural networks, which yields the concept of VAMPnet vampnet . The optimal models can therefore be obtained by deep learning techniques. Another possible way is to utilize tensor decomposition based approximation approaches. Some tensor analysis methods have been presented based on the reversible variational principle and EDMD NueskeEtAl_JCP15_Tensor ; KlusSchuette_Arxiv15_Tensor ; Klus_Arxiv16_Tensor , and it is worth studying more general variational tensor method within the framework of VAMP in future.

One drawback of the methods developed in this paper is that the resulting models are possibly not valid probabilistic models with nonnegative transition densities if only the operator error is considered, and the probability-preserving modeling method requires further investigations. Moreover, the applications of VAMP to detection of metastable states DeuflhardWeber_LAA05_PCCA+ , coherent sets froyland2014almost and dominant cycles conrad2016finding will also be explored in next steps.

Appendix

for every measurable set $A$, and define the matrix of scalar products:

for $\mathbf{a}=(a_{1},a_{2},\ldots,a_{m})^{\top}$, $\mathbf{b}=(b_{1},b_{2},\ldots,b_{n})^{\top}$ and $\mathbf{g}=(g_{1},g_{2},\ldots)^{\top}$. In addition, $\mathcal{N}(\cdot|c,\sigma^{2})$ denotes the probability density function of the normal distribution with mean $c$ and variance $\sigma^{2}$.

Appendix A Analysis of Koopman operators

where $\mathcal{P}_{t}$ denotes the Markov propagator defined in (63). We can then conclude that the estimates of $\mathbf{C}_{00},\mathbf{C}_{11},\mathbf{C}_{01}$ given by (16-18) are unbiased and consistent as $S\to\infty$.

In more general cases where trajectories $\{\mathbf{x}_{t}^{1}\}_{t=1}^{T_{1}},\ldots,\{\mathbf{x}_{t}^{S}\}_{t=1}^{T_{S}}$ are generated with different initial conditions and different lengths, the similar conclusions can be obtained by defining $\rho_{0},\rho_{1}$ as the averages of marginal distributions of $\{\mathbf{x}_{t}^{s}|1\leq t\leq T_{s}-\tau,1\leq s\leq S\}$ and $\{\mathbf{x}_{t}^{s}|1+\tau\leq t\leq T_{s},1\leq s\leq S\}$ respectively.

A.2 Proof of Theorem 2.1

Because $\mathcal{K}_{\tau}$ is a Hilbert-Schmidt operator from $\mathcal{L}_{\rho_{1}}^{2}$ to $\mathcal{L}_{\rho_{0}}^{2}$, there exists the following SVD of $\mathcal{K}_{\tau}$:

Due to the orthonormality of right singular functions, the projection of any function $g\in\mathcal{L}_{\rho_{1}}^{2}$ onto the space spanned by $\{\phi_{1},\ldots,\phi_{k}\}$ can be written as $\sum_{i=1}^{k}\left\langle g,\phi_{i}\right\rangle_{\rho_{1}}\phi_{i}$. Then $\hat{\mathcal{K}}_{\tau}$ defined by (5) is the approximate Koopman operator deduced from model (4), and it is the best rank $k$ approximation to $\mathcal{K}_{\tau}$ in Hilbert-Schmidt norm according to the generalized Eckart-Young Theorem (see Theorem 4.4.7 in hsing2015theoretical ).

Since the adjoint operator $\mathcal{K}_{\tau}^{*}$ of $\mathcal{K}_{\tau}$ satisfies

A.3 Transition densities deduced from Koopman operators

The Koopman operator can also be written as

if the transition density is given, which implies that

Then the transition density deduced from the approximate Koopman operator $\hat{\mathcal{K}}_{\tau}$ defined by (5) is

i.e., the operator error between $\hat{\mathcal{K}}_{\tau}$ and $\mathcal{K}_{\tau}$ can be represented by the error between $\hat{p}_{\tau}$ and $p_{\tau}$.

It is worth pointing out that the approximate transition density in (54) satisfies the normalization constraint with

but $\hat{p}_{\tau}(\mathbf{x},\mathbf{y})$ is possibly negative for some $\mathbf{x},\mathbf{y}$. Thus, the approximate Koopman operators and transition densities are not guaranteed to yield valid probabilistic models, although they can still be utilized to quantitative analysis of Markov processes.

A.4 Sufficient conditions for Theorem 2.1

We show here $\mathcal{L}_{\rho_{0}}^{2},\mathcal{L}_{\rho_{1}}^{2}$ are separable Hilbert spaces and $\mathcal{K}_{\tau}:\mathcal{L}_{\rho_{1}}^{2}\mapsto\mathcal{L}_{\rho_{0}}^{2}$ is Hilbert-Schmidt if one of the following conditions is satisfied:

The state space of the Markov process is a finite set.

The proof is trivial by considering $\mathcal{K}_{\tau}$ is a linear operator between finite-dimensional spaces, and thus omitted.

A.5 Koopman operators of deterministic systems

For the completeness of paper, we prove here the following proposition by contradiction: The Koopman operator $\mathcal{K}_{\tau}$ of the deterministic system $\mathbf{x}_{t+\tau}=F(\mathbf{x}_{t})$ defined by

is not a compact operator from $\mathcal{L}_{\rho_{1}}^{2}$ to $\mathcal{L}_{\rho_{0}}^{2}$ if $\mathcal{L}_{\rho_{1}}^{2}$ is infinite-dimensional.

Assume that $\mathcal{K}_{\tau}$ is compact. Then, the SVD (50) of $\mathcal{K}_{\tau}$ exists with $\sigma_{i}\to 0$ as $i\to\infty$, and there is $j$ so that $0\leq\sigma_{j}<1$. This implies $\left\langle\mathcal{K}_{\tau}\psi_{j},\mathcal{K}_{\tau}\psi_{j}\right\rangle_{\rho_{0}}=\sigma_{j}^{2}<1$. However, according to the definition of the Koopman operator, $\left\langle\mathcal{K}_{\tau}\psi_{j},\mathcal{K}_{\tau}\psi_{j}\right\rangle_{\rho_{0}}=\left\langle\psi_{j},\psi_{j}\right\rangle_{\rho_{1}}=1$, which leads to a contradiction. We can conclude that $\mathcal{K}_{\tau}$ is not compact and hence not Hilbert-Schmidt.

Appendix B Markov propagators

The Markov propagator $\mathcal{P}_{\tau}$ is defined by

Where the following normalizations were used:

The SVD of $\mathcal{P}_{\tau}$ can be written as

Appendix C Proof of the variational principle

Notice that $\mathbf{f}$ and $\mathbf{g}$ can be expressed as

the optimization problem can be equivalently written as

under the constraint $\mathbf{D}_{0}^{\top}\mathbf{D}_{0}=\mathbf{I},\mathbf{D}_{1}^{\top}\mathbf{D}_{1}=\mathbf{I}$. The variational principle can then be proven by considering

when the first $k$ rows of $\mathbf{D}_{0}$ and $\mathbf{D}_{1}$ are identity matrix.

Appendix D Variational principle of reversible Markov processes

The variational principle of reversible Markov processes can be summarized as follows: If the Markov process $\{\mathbf{x}_{t}\}$ is time-reversible with respect to stationary distribution $\mu$ and all eigenvalues of $\mathcal{K}_{\tau}$ is nonnegative, then

for $r\geq 1$ and the maximal value is achieved with $f_{i}=\psi_{i}$, where $\psi_{i}$ denotes the eigenfunction with the $i$th largest eigenvalue $\lambda_{i}$. The proof is trivial by using variational principle of general Markov processes and considering that the eigendecomposition of $\mathcal{K}_{\tau}$ is equivalent to its SVD if $\{\mathbf{x}_{t}\}$ is time-reversible and $\rho_{0}=\rho_{1}=\mu$.

Appendix E Analysis of estimation algorithms

We show in this appendix that the feature TCCA algorithm described in Section 3.1 solves the optimization problem (19).

Let $\mathbf{U}^{\prime}=\mathbf{C}_{00}^{\frac{1}{2}}\mathbf{U}=(\mathbf{u}_{1}^{\prime},\ldots,\mathbf{u}_{k}^{\prime})$ and $\mathbf{V}^{\prime}=\mathbf{C}_{00}^{\frac{1}{2}}\mathbf{V}=(\mathbf{v}_{1}^{\prime},\ldots,\mathbf{v}_{k}^{\prime})$, (19) can be equivalently expressed as

According to the Cauchy-Schwarz inequality and the conclusion in Section I.3.C of marshall1979inequalities , we have

under the constraints $\mathbf{U}^{\prime\top}\mathbf{U}^{\prime}=\mathbf{I},\mathbf{V}^{\prime\top}\mathbf{V}^{\prime}=\mathbf{I}$, where $s_{i}$ is the $i$th largest singular value of $\mathbf{C}_{00}^{-\frac{1}{2}}\mathbf{C}_{01}\mathbf{C}_{11}^{-\frac{1}{2}}$. Considering the equalities hold in the above when $\mathbf{U}^{\prime},\mathbf{V}^{\prime}$ are the first $k$ left and right singular vectors of $\mathbf{C}_{00}^{-\frac{1}{2}}\mathbf{C}_{01}\mathbf{C}_{11}^{-\frac{1}{2}}$, we get

and the correctness of the feature TCCA can then be proved.

E.2 Feature TCCA of projected Koopman operators

E.3 An example of nonlinear TCCA

where $u_{t}$ is Gaussian white noise with mean zero and variance $1$. By setting

to be the stationary distribution and basis functions

with parameter $w\in[0.01,1]$, we can obtain

The maximal VAMP-$r$ score for a given $w$ can then be analytically computed by

according to (24). We evaluate $\mathcal{R}_{r}(w)$ at $9901$ equally spaced points of $w$ in the interval $[0.01,1]$ for $r=1,2$, and the maximal values of $\mathcal{R}_{1},\mathcal{R}_{2}$ are achieved at $w=0.3157$ and $w=0.7069$ respectively.

Appendix F Implementation of estimation algorithms

For convenience of notation, here we define

In this paper, we utilize principal component analysis (PCA) to explicitly reduce correlations between basis functions as follows: First, we compute the empirical means of basis functions and the covariance matrices of mean-centered basis functions:

Next, perform the truncated eigen decomposition of the covariance matrices as

where the diagonal of matrices $\mathbf{S}_{0,d},\mathbf{S}_{1,d}$ contain all positive eigenvalues that are larger than $\epsilon_{0}$ and absolute values of all negative eigenvalues ($\epsilon_{0}=10^{-10}$ in our applications). Last, the new basis functions are given by

Then the feature TCCA algorithm with de-correlation of basis functions can be summarized as:

Compute covariance matrices $\mathbf{C}_{00},\mathbf{C}_{01},\mathbf{C}_{11}$ by

Perform the truncated SVD $\mathbf{C}_{00}^{-\frac{1}{2}}\mathbf{C}_{01}\mathbf{C}_{11}^{-\frac{1}{2}}=\mathbf{U}_{k}^{\prime}\hat{\boldsymbol{\Sigma}}_{k}\mathbf{V}_{k}^{\prime\top}$.

Notice that the estimated $\mathbf{C}_{00}$, $\mathbf{C}_{01}$ and $\mathbf{C}_{11}$ in the above algorithm satisfy

where $\mathbf{C}\succeq 0$ means $\mathbf{C}$ is a positive semi-definite matrix. According to the Schur complement lemma, we have

where $\mathbf{I}$ denotes an identity matrix of appropriate size. So the estimated $\sigma_{1}\leq 1$.

which implies that $1$ is the largest singular value of $\mathbf{C}_{00}^{-\frac{1}{2}}\mathbf{C}_{01}\mathbf{C}_{11}^{-\frac{1}{2}}$.

F.2 Parameter optimization in nonlinear TCCA

can be solved by direct search as in our examples (see Appendix K.1). But for a high-dimensional parameter vector $\mathbf{w}$, it is more efficient to perform the optimization by the gradient descent method in the form of

where $\eta$ is the step size. When $r=2$, the gradient of $\mathcal{R}_{r}$ with respect to an element $w_{i}$ in $\mathbf{w}$ can be written as

where $\mathbf{X},\mathbf{Y}$ have the same definitions as in Appendix F.1. If the data size is too large, we can approximate the gradient based on a random subset of data in each iteration, and update $\mathbf{w}$ in a stochastic gradient descent manner andrew2013deep ; vampnet .

Like feature TCCA, the nonlinear TCCA also suffers from the numerical singularity when $\mathbf{C}_{00}$ or $\mathbf{C}_{11}$ is not full rank. This problem can be addressed by the de-correlation of basis functions when performing direct search. For the gradient descent method (or stochastic gradient descent method), we can replace the objective function $\mathcal{R}_{r}(\mathbf{w})$ by a regularized one

where $\epsilon>0$ is a hyperparameter and can be selected by the cross-validation.

Appendix G Relationship between VAMP and EDMD

The proof of (21) is trivial. Here, we only show that the eigenvalue problem of $\hat{\mathcal{K}}_{\tau}$ given by the feature TCCA is equivalent to that of matrix $\mathbf{K}_{\chi}$ as

under the assumption that $\boldsymbol{\chi}_{0}=\boldsymbol{\chi}_{1}=\boldsymbol{\chi}$ and $\mathbf{C}_{00}$ is invertible, which is consistent with the spectral approximation theory in EDMD. First, if $g$ and $\lambda$ satisfy $\mathcal{K}_{\tau}g=\lambda g$, there must exist vector $\mathbf{b}$ so that $g=\mathbf{b}^{\top}\boldsymbol{\chi}$. Then

Second, if $\mathbf{K}_{\chi}\mathbf{b}=\lambda\mathbf{b}$,

Appendix H Analysis of the VAMP-E score

Considering $\{\phi_{i}\}$ is an orthonormal basis of $\mathcal{L}_{\rho_{1}}^{2}$, we have

H.2 Relationship between VAMP-2 and VAMP-E

We first show that the feature TCCA algorithm maximizes VAMP-E. Notice that

where $\left\|\cdot\right\|_{F}$ denotes the Frobenius norm and $\mathbf{U}^{\prime}=\mathbf{C}_{00}^{\frac{1}{2}}\mathbf{U}$, $\mathbf{V}^{\prime}=\mathbf{C}_{11}^{\frac{1}{2}}\mathbf{V}$. It can be seen that the feature TCCA algorithm maximizes the first term on the right-hand side of (121) and therefore maximizes VAMP-E.

For the optimal model generated by the nonlinear TCCA, the first term on the right-hand side of (121) is equal to zero and the second term is maximized as a function of $\mathbf{w}$. Thus, the nonlinear TCCA also maximizes VAMP-E.

In addition, for $\mathbf{K},\mathbf{U},\mathbf{V}$ provided by both feature TCCA and nonlinear TCCA,

Appendix I Subspace variational principle

The variational principle proposed in Section 2.2 can be further extended to singular subspaces of the Koopman operator as follows:

The approximate Koopman operator in the form of (27) can also be written as

Notice that substituting $\mathbf{f}=\mathbf{U}^{\top}\boldsymbol{\chi}_{0},\mathbf{g}=\mathbf{V}^{\top}\boldsymbol{\chi}_{1}$ into (128) yields (37).

Appendix K Details of numerical examples

For convenience of analysis and computation, we partition the state space $$into$2000$bins$S_{1},\ldots,S_{2000}$ uniformly, and discretize the one-dimensional dynamical system described in Example 1 as

where $s_{i}$ is the center of the bin $S_{i}$, and the local distribution of $x_{t}$ within any bin is always uniform distribution. All numerical computations and simulations in Examples 1, 2 and 3 are based on (131), and the initial state $x_{0}$ is distributed according to the stationary distribution $\rho_{0}=\rho_{1}=\mu$.

In Example 1, the the stationary distribution and singular components of the Koopman operator are analytically computed by the feature TCCA with basis functions $\chi_{0,i}(x)=\chi_{1,i}(x)=1_{x\in S_{i}}$ as follows:

Compute $\mathbf{U}=[U_{ij}]=\mathbf{C}_{00}^{-\frac{1}{2}}\mathbf{U}^{\prime}$ and $\mathbf{V}=[V_{ij}]=\mathbf{C}_{11}^{-\frac{1}{2}}\mathbf{V}^{\prime}$.

Output the stationary distribution $\mu(x)=\sum_{i}50\pi_{i}\cdot 1_{x\in S_{i}}$ and singular components $(\sigma_{i},\psi_{i}(x),\phi_{i}(x))=(\sigma_{i},\sum_{j}U_{ji}\cdot 1_{x\in S_{j}},\sum_{j}V_{ji}\cdot 1_{x\in S_{j}})$.

The transition density of the projected Koopman operator $\hat{\mathcal{K}}_{\tau}=\sum_{i=1}^{k}\sigma_{i}\left\langle\cdot,\phi_{i}\right\rangle_{\rho_{1}}\psi_{i}$ is obtained by

(see Appendix A.3) and the corresponding approximate transition matrix is

the long-time transition density in Fig. 2 is given by

and the cumulative error of $\hat{p}_{n\tau}(x,y)$ is

In Examples 2 and 3, the smoothing parameter $w$ are optimized by the golden-section search algorithm press2007numerical as follows for nonlinear TCCA:

Let $a=-6$, $b=6$, $c=0.618a+0.382b$, $d=0.382a+0.618b$.

Compute $\mathcal{R}_{2}(\exp a)$, $\mathcal{R}_{2}(\exp b)$, $\mathcal{R}_{2}(\exp c)$ and $\mathcal{R}_{2}(\exp d)$, where $\mathcal{R}_{2}(w)=\left\|\mathbf{C}_{00}\left(w\right)^{-\frac{1}{2}}\mathbf{C}_{01}\left(w\right)\mathbf{C}_{11}\left(w\right)^{-\frac{1}{2}}\right\|_{F}^{2}$ and $\|\cdot\|_{F}$ denotes the Frobenius norm.

If $\max\{\mathcal{R}_{2}(\exp a),\mathcal{R}_{2}(\exp b),\mathcal{R}_{2}(\exp c)\}>\max\{\mathcal{R}_{2}(\exp b),\mathcal{R}_{2}(\exp c),\mathcal{R}_{2}(\exp d)\}$, let $(a,b,c,d):=(a,d,0.618a+0.382d,c)$. Otherwise, let $(a,b,c,d):=(c,b,d,0.618b+0.382c)$.

If $|a-b|<10^{-3}$, output $\log w\in\{a,b,c,d\}$ with the largest value of $\mathcal{R}_{2}(w)$. Otherwise, go back to Step 2.

Furthermore, $w$ is computed in the same way when perform nonlinear TCCA in Sections 5.1 and 5.2.

K.2 Double-gyre system

For the double-gyre system in Section 5.1, we first perform the temporal discretization by the Euler–Maruyama scheme as

where $\mathbf{x}_{t}=(x_{t},y_{t})^{\top}$ and $\Delta=0.02$ is the step size. Then perform the spatial discretization as

Here $S_{1},\ldots,S_{1250}$ are $50\times 25$ bins which form a uniform partition of the state space $\times$ and $(s_{i,x},s_{i,y})$ represents the center of $S_{i}$. Simulation data and the “true” singular components are all computed by using (138) with the initial distribution of $(x_{0},y_{0})$ being the stationary one.

In Fig. 6, the transition density of lag time $n\tau$ is computed from the estimated singular components $(\mathbf{K},\mathbf{U}^{\top}\boldsymbol{\chi}_{0},\mathbf{V}^{\top}\boldsymbol{\chi}_{1})$ as

is the approximate transition matrix, and $\boldsymbol{\rho}_{1}=[\boldsymbol{\rho}_{1i}]$ with

References