Optimal Data-Driven Estimation of Generalized Markov State Models for Non-Equilibrium Dynamics

: There are multiple ways in which a stochastic system can be out of statistical equilibrium. It might be subject to time-varying forcing; or be in a transient phase on its way towards equilibrium; it might even be in equilibrium without us noticing it, due to insufﬁcient observations; and it even might be a system failing to admit an equilibrium distribution at all. We review some of the approaches that model the effective statistical behavior of equilibrium and non-equilibrium dynamical systems, and show that both cases can be considered under the uniﬁed framework of optimal low-rank approximation of so-called transfer operators. Particular attention is given to the connection between these methods, Markov state models, and the concept of metastability, further to the estimation of such reduced order models from ﬁnite simulation data. All these topics bear an important role in, e.g., molecular dynamics, where Markov state models are often and successfully utilized, and which is the main motivating application in this paper. We illustrate our considerations by numerical examples.


Introduction
The term "equilibrium" is not used uniformly throughout the literature. So, to start off with, an equilibrium process in this paper means a reversible one in the stochastic sense, see Table 1.
Metastable molecular systems under non-equilibrium conditions caused by external fields have attracted increasing interest recently. For instance, new experimental techniques like atomic force microscopy or simulation studies regarding the potential effects of electromagnetic radiation on the human body tissue have been extensively investigated in the literature. Specifically adapted molecular dynamics (MD) simulations have proved particularly useful for understanding the response of biomolecular conformations to external fields. Despite this significance, reliable tools for the quantitative description of non-equilibrium phenomena like the conformational dynamics of a molecular system under external forcing are still lacking.
For MD simulations in equilibrium such specific and reliable tools have been developed: Markov State Models (MSMs) allow for an accurate description of the transitions between the main conformations of the molecular system under investigation. MSMs for equilibrium MD have been well developed over the past decade in theory [SS13, PWS + 11], applications (see the recent book [BPN14] for an overview), and software implementations [STSM + 12, HSH + 17]. They now form a set of standard tools. The principal idea of equilibrium MSMs is to approximate the MD system (in continuous state or phase space) by a reduced Markovian dynamics over a finite number of (macro-)states (i.e., in discrete state space). These (macro-)states represent the dominant metastable sets of the system, i.e., sets in which typical MD trajectories stay substantially longer than the system needs for a transition to another such set [SS13,SNL + 11]. In equilibrium MD, these metastable sets are the main conformations of the molecular system under consideration which, often enough, are given by the main wells in its energy landscape. It has been shown that for many (bio)molecular systems the Markovian dynamics given by an MSM allows very close approximation of the longest relaxation processes of the underlying molecular system under equilibrium conditions [SNS10,DSS12].
However, in the non-equilibrium setting the above tools are not guaranteed to continue working. Note that there are different possibilities to deviate from the "equilibrium" situation, and this makes the term "non-equilibrium" ambiguous. To avoid confusion, we consider one of the following cases when referring to the non-equilibrium setting (again, see Table 1 on terminology).
(i) Time-inhomogeneous dynamics, e.g., the system feels a time-dependent external force, for instance due to an electromagnetic field or force probing. (ii) Time-homogeneous non-reversible dynamics, i.e., where the governing laws of the system do not change in time, but the system does not obey detailed balance, and, additionally we might want to consider the system in a non-stationary regime. (iii) Reversible dynamics but non-stationary data, i.e., the system possesses a stationary distribution with respect to which it is in detailed balance, but the empirical distribution of the available data did not converge to this stationary distribution. Even though we consider genuinely stochastic systems here, the algorithm of section 5 can be used for deterministic systems as well-and indeed it is, see Remark 5.2 and references therein.
Note that with regard to the considered dynamics (i)-(iii) represent cases with decreasing generality. For (i), time-dependent external fields act on the system, such that the energy landscape depends on time, i.e., the main wells of the energy landscape can move in time. That is, there may no longer be time-independent metastable sets in which the dynamics stays for long periods of time before exiting. Instead, the potentially metastable sets will move in state space. Generally, moving "metastable" sets cannot be considered metastable anymore. However, the so-called coherent sets, which have been studied for non-autonomous flow fields in fluid dynamics [FSM10,Fro13], permit to get a meaning to the concept of metastability [KCS16]. For (iii), the full theory of equilibrium Markov state modeling is at one's disposal, but one needs to estimate certain required quantities from non-equilibrium data [WNP + 17]. Case (ii) seems the most elusive, due to the fact that on the one hand it could be handled by the timeinhomogeneous approach, but on the other hand it is a time-homogeneous system and some structural properties could be carried over from the reversible equilibrium case that are out of reach for a time-inhomogeneous analysis. For instance, if the dynamics shows cyclic behavior, it admits structures that are well captured by tools from the analysis of time-homogeneous dynamical systems (e.g., Floquet theory and Poincaré sections [SW15,FK17]), and a more general view as in (i) might miss them; however, cyclic behavior is not present in reversible systems, such that the tools from (iii) are doomed to failure in this respect. In order to avoid confusion, however, it should be emphasized that the three cases distinguished above do not suffice to clarify the discussion about the definition of equilibrium or non-equilibrium, e.g., see the literature on non-equilibrium steady state (NESS) systems [SS10,LLP17].
Apart from MSMs the literature on kinetic lumping schemes offers several other techniques for finding a coarse-grained descriptions of systems [YCB + 13, BMH13,KS15]. These techniques are, however, not built on the intuition of metastable behavior in state space. What we consider here can be seen in connection to optimal prediction in the sense of the Mori-Zwanzig formalism [Mor65, Zwa73, CHK00, CHK02], but we will try to choose the observables of the system such that projecting the dynamics on these keeps certain properties intact without including memory terms.
The aim of this article is to review and unify some of the theoretical and also datadriven algorithmic approaches that attempt to model the effective statistical behavior of non-equilibrium systems. To this end, a MSM, or, more precisely, a generalized MSM is sought, i.e., a possibly small matrix T k that carries the properties of the actual system that are of physical relevance. In the equilibrium case, for example, this includes the slowest timescales on which the system relaxes towards equilibrium (section 2). The difference of generalized to standard MSMs is that we do not strictly require the former to be interpretable in terms of transition probabilities between some regions of the state space (section 3), however usually there is a strong connection between the matrix entries and metastable sets. We will, however, focus on a slightly different characteristic of the approximate model, namely its "propagation error", that will allow for a straightforward generalization from equilibrium (reversible) to all our non-equilibrium cases (section 4); and even retain the physical intuition behind true MSMs through the concept of coherent sets. We will show in section 5 how these considerations can be carried over to the case when only a finite amount of simulation data is available. The above non-equilibrium cases (ii)-(iii) can be then given as specific instances of the construction (section 6). The theory is illustrated with examples throughout the text.
We note in advance that in course of the (generalized) Markov state modeling we will consider different instances of approximations to a certain linear operator T : S → S mapping some space to itself (and sometimes to a different one). On the one hand, there will be a projected operator T k : S → S, where T k = QT Q with a projection Q : S → V having a k-dimensional range V ⊂ S. On the other hand, we will consider the restriction of the projected operator T k to this k-dimensional subspace, i.e., T k : V → V, also called V-restriction of T k , which has a k ×k matrix representation (with respect to some chosen basis of V) that we will denote by T k .
A stochastic (Markov) process is called... time-homogeneous if the transition probabilities from time s to time t depend only on t − s (in analogy to the evolution of an autonomous ODE).
stationary if the distribution of the process does not change in time (such a distribution is also called invariant, cf. (2)).
reversible if it is stationary and the detailed balance condition (5) holds (reversibility means that time series are statistically indistinguishable in forward and backward time).
describing diffusion in the potential energy landscape given by W . Here, β is the nondimensionalized inverse temperature, and w t is a standard Wiener process (Brownian motion). The transition density function p t : X × X → R ≥0 of a time-homogeneous stochastic process {x t } t≥0 is defined by That is, p t (x, y) is the conditional probability density of x t = y given that x 0 = x. We also write x t ∼ p t (x 0 , ·) to indicate that x t has density p t (x 0 , ·).
With the aid of the transition density function, we will now define transfer operators, i.e., the action of the process on functions of the state. Note, however, that the transition density is in general not known explicitly, and thus we will need data-based approximations to estimate it. We assume that there is a unique stationary density µ, such that {x t } t≥0 is stationary with respect to µ; that is, it satisfies x 0 ∼ µ and ( Let now f be a probability density over X, u = f /µ a probability density with respect to µ (meaning that uµ is to be interpreted as a physical density), and g a scalar function of the state (an "observable"). We define the following transfer operators, for a given lag time τ : (a) The Perron-Frobenius operator (also called propagator), The Perron-Frobenius operator with respect to the equilibrium density (also called transfer operator, simply), evolves densities with respect to µ. (c) The Koopman operator (3) evolves observables. All our transfer operators are well-defined non-expanding operators on the following Hilbert spaces: 1 , P τ : The equilibrium density µ satisfies P τ µ = µ, that is, µ is an eigenfunction of P τ with associated eigenvalue λ 0 = 1. The definition of T τ relies on µ, we have thus P τ µ = µ translates into T τ 1 = 1, where 1 = 1 X is the constant one function on X.

Reversible equilibrium dynamics and spectral decomposition
An important structural property of many systems used to model molecular dynamics is reversibility. Reversibility means that the process is statistically indistinguishable from its time-reversed counterpart, and it can be described by the detailed balance condition The process generated by (1) is reversible and ergodic, i.e., it admits a unique positive equilibrium density, given by µ(x) ∝ exp(−βW (x)), under mild growth conditions on the potential W [MS02,MSH02]. Note that the subsequent considerations hold for all stochastic processes that satisfy reversibility and ergodicity with respect to a unique positive invariant density and are not limited to the class of dynamical systems given by (1). See [SS13] for a discussion of a variety of stochastic dynamical systems that have been considered in this context. Furthermore, if p t (·, ·) is a continuous function in both its arguments for t > 0, then all the transfer operators above are compact, which we also assume from now on. This implies that they have a discrete eigen-and singular spectrum (the latter meaning it has a discrete set of singular values). For instance, the process generated by (1) has has continuous transition density function under mild growth and regularity assumptions on the potential W .
As a result of the detailed balance condition, the Koopman operator K τ and the Perron-Frobenius operator with respect to the equilibrium density T τ become identical and we obtain i.e., all the transfer operators become self-adjoint on the respective Hilbert spaces from above. Here, ·, · ν denotes the natural scalar products on the weighted space L 2 ν , i.e., f, g ν = X f (x)g(x)ν(x) dx; the associated norm is denoted by · ν . Due to the self-adjointness, the eigenvalues λ τ i of P τ and T τ are real-valued and the eigenfunctions form an orthogonal basis with respect to ·, · 1/µ and ·, · µ , respectively.
Ergodicity implies that the dominant eigenvalue λ 1 is the only eigenvalue with absolute value 1 and we can thus order the eigenvalues so that The eigenfunction of T τ corresponding to λ 1 = 1 is the constant function φ 1 = 1 X . Let φ i be the normalized eigenfunctions of T τ , i.e. φ i , φ j µ = δ ij , where δ ij denotes the Kronecker-delta. Then any function f ∈ L 2 µ can be written in terms of the eigenfunctions For more details, we refer to [KNK + 18] and references therein. For some k ∈ N, we call the k dominant eigenvalues λ τ 1 , . . . , λ τ k of T τ the dominant spectrum of T τ , i.e., λ dom (T τ ) = {λ τ 1 , . . . , λ τ k }. Usually, k is chosen in such a way that there is a spectral gap after λ τ k , i.e. 1 − λ τ k λ τ k − λ τ k+1 . The (implied) time scales on which the associated dominant eigenfunctions decay are given by If {T t } t≥0 is a semigroup of operators (which is the case for every time-homogeneous process, as, e.g., the transfer operator associated with (1)), then there are Assuming there is a spectral gap, the dominant time scales satisfy ∞ = t 1 > . . . ≥ t k t k+1 . These are the time scales of the slow dynamical processes, also called rare events, which are of primary interest in applications. The other, fast processes are regarded as fluctuations around the relative equilibria (or metastable states) between which the relevant slow processes travel.

Markov state models for reversible systems in equilibrium
In the following, we will fix a lag time τ > 0, and drop the superscript τ from the transfer operators for clarity of notation.

Preliminaries on equilibrium Markov state models
Generally, in the equilibrium case, a generalized MSM (GMSM) is any matrix T k ∈ R n k ×n k , n k ≥ k, that approximates the k dominant time scales of T , i.e., its dominant eigenvalues; It is natural to ask for some structural properties of T to be reproduced by T k , such as: • T is a positive operator ←→ all entries of T k are non-negative; • T is probability-preserving ←→ each column sum of T k is 1. These two properties together make T k to a stochastic matrix, and in this case T k is usually called a MSM. We shall use the term Generalized MSM for a matrix T k that violates these requirements but still approximates the dominant spectral components of the underlying operator. Another structural property that one would usually ask for is to have apart from the time scales/eigenvalues also some approximation of the associated eigenvectors of T , as these are the dynamic observables related to the slow dynamics. This is incorporated in the general approach, which we discuss next.
The question is now how to obtain a GMSM T k for a given T ? To connect these objects, a natural and popular approach is to obtain the reduced model T k via projection. To this end, let Q : L 2 µ → V ⊂ L 2 µ be a projection onto a n k -dimensional subspace V. The GMSM is then defined by the projected transfer operator and T k can now be taken as the matrix representation of the V-restriction of the projected operator T k : V → V with respect to a chosen basis of V. Is there a "best" choice for the projection? If we also ask for perfect approximation of the time scales, i.e., λ dom (T k ) = λ dom (T ), the requirement of parsimony-such that the model size is minimal, i.e., n k = k-leaves us with a unique choice for V, namely the space spanned by the dominant (normalized) eigenfunctions φ i of T , i = 1, . . . , k. This follows from the so-called variational principle (or Rayleigh-Ritz method) [NN13, NKPH + 14]. In fact, it makes a stronger claim: every projection to a k-dimensional space V yields a GMSM T k : V → V which underestimates the dominant time scales, i.e., λ i (T k ) ≤ λ i (T ), i = 1, . . . , k; and equality holds only for the projections on the eigenspaces.
Note that the discussion about the time scales (equivalently, the eigenvalues) involves only the range of the projection, the space V. However, there are multiple ways to project on the space V. It turns out, that the µ-orthogonal projection given by is superior to all of them, if we consider a stronger condition than simply reproducing the dominant time scales. This condition is the requirement of minimal propagation error, and it will be central to our generalization of GMSMs for non-equilibrium, or even time-inhomogeneous systems. Let us define the best k-dimensional approximation T k to T , i.e., the best projection Q, as the rank-k operator satisfying where · denotes the induced operator norm 2 for operators mapping L 2 µ to itself. Equivalently, this can be viewed as a result stating that T k is the k-dimensional approximation of T yielding the smallest (worst-case) error in density propagation: where x * = arg min x h(x) means that x * minimizes the function h, possibly subject to constraints that are listed under arg min.
To summarize, the best GMSM (9) in terms of (11) (or, equivalently, (12)) is given by the projection (10). This follows from the self-adjointness of T and the Eckard-Young theorem; details can be found in [WN17] and in Appendix A. Caution is needed however, when interpreting T k f as the propagation of a given probability density f . The projection to the dominant eigenspace in general does not respect positivity (i.e., , thus T k f loses its probabilistic meaning. This is the price to pay for the perfectly reproduced dominant time scales. We can retain a physical interpretation of a MSM if we accept that the dominant time scales will be slightly off, as we discuss in the next section.

Metastable sets
There is theoretical evidence [SS13] that the more pronounced the metastable behavior of system is (in the sense that the size of the time scale gap t 1 ≥ . . . ≥ t k t k+1 is large), the more constant the dominant eigenfunctions φ i are on the metastable sets M 1 , . . . , M d , given the lag time with respect to which the transfer operator T = T τ is taken satisfies τ t k+1 . Assuming such a situation, the eigenfunctions of T can approximately be combined from the characteristic functions over the metastable sets, i.e., with the abbreviation where the c ij are components of the linear combination, such that the φ i are orthonormal. Using the "approximate eigenfunctions" φ i defined in (13), the modified projection where P µ [ · | x ∈ M] denotes the probability measure that arises if x ∈ M has distribution µ (restricted to M). That is, T k has the transition probabilities between the metastable sets as entries, giving a direct physical interpretation of the MSM. Note, however, that for this approximation to reproduce the dominant time scales well, i.e., to have t i ≈ t i , i = 1, . . . , k, we need a strong separation of time scales in the sense that t k t k+1 has to hold, and the lag time τ needs to be chosen sufficiently large [SNS10].

Markov state models for time-inhomogeneous systems
As all our non-equilibrium cases will be special instances of the most general, timeinhomogeneous case, we consider this next.

Minimal propagation error by projections
Conceptual changes. The above approach to Markov state modeling is relying on the existence of an stationary distribution and reversibility. In the case of a timeinhomogeneous system there will not be, in general, any stationary distribution µ. Additionally, we are lacking physical meaning, since it is unclear with respect to which ensemble the dynamical fluctuations should be described. From a mathematical perspective there is a problem as well, since the construction relies on the reversibility of the underlying system, which gives the self-adjointness of the operator T with respect to the weighted scalar product ·, · µ . Time-inhomogeneous systems are not reversible in general. Additionally to these structural properties, we might need to depart from some conceptional ones as well. As time-inhomogeneity usually stems from an external forcing that might not be present or known for all times, we need a description of the system on a finite time interval. This disrupts the concept of dominant time scales as they are considered in equilibrium systems, because there it relies on self-similarity of observing an eigenmode over and over for arbitrary large times. It also forces us to re-visit the concept of metastability for two reasons. First, many definitions of metastability rely on statistics under the assumption that we observe the system for infinitely long times. Second, as an external forcing may theoretically arbitrarily distort the energy landscape, it is a priori unclear what could be a metastable set.
As a remedy, we aim at another property when trying to reproduce the effective behavior of the full system by a reduced model; this will be minimizing the propagation error, as in (12). Remarkably, this will also allow for a physical interpretation through so-called coherent sets; analogous to metastable sets in the equilibrium case.
A prototypical time-inhomogeneous system can be given by where the potential W now depends explicitly on time t. In this case, a lag time τ is not sufficient to parametrize the statistical evolution of the system, because we need to know when we start the evolution. Thus, transition density functions need two time parameters, e.g., p s,t (x, ·) denotes the distribution of x t conditional to x s = x. Similarly, the transfer operators P, T , K are parametrized by two times as well, e.g., P s,t propagates probability densities from initial time s to final time t (alternatively, from initial time s for lag time τ = t − s). To simplify notation, we will settle for some initial and final times, and drop these two time parameters, as they stay fixed.
Adapted transfer operators. Let us observe the system from initial time t 0 to final time t 1 , such that its distribution at initial time is given by µ 0 . Then, if P denotes the propagator of the system from t 0 to t 1 , then we can express the final distribution at time t 1 by µ 1 = Pµ 0 . As the transfer operator in equilibrium case was naturally mapping L 2 µ to itself (because µ was invariant), here it is natural to consider the transfer operator mapping densities (functions) with respect to µ 0 to densities with respect to µ 1 . Thus, we define the transfer operator T : L 2 µ 0 → L 2 µ 1 by which is the non-equilibrium analogue to (4). This operator naturally retains some properties of the equilibrium transfer operator [Den17]: • T 1 = 1, encoding the property that µ 0 is mapped to µ 1 by the propagator P.
• Its adjoint is the Koopman operator K : An optimal non-stationary GMSM. As already mentioned above, it is not straightforward how to address the problem of Markov state modeling in this time-inhomogeneous case via descriptions involving time scales or metastability. Instead, our strategy will be to search for a rank-k projection T k of the transfer operator T with minimal propagation error, to be described below.
The main point is now that due to the non-stationarity the domain L 2 µ 0 (where T maps from) and range L 2 µ 1 (where T maps to) of the transfer operator T are different spaces, hence it is natural to choose different rank-k subspaces as domain and range of T k too. In fact, it is necessary to choose domain and range differently, since f ∈ L 2 µ 0 has a different meaning than f ∈ L 2 µ 1 . Thus, we will search for projectors Q 0 : has essentially optimal propagation error. In quantitative terms, we seek to solve the optimization problem where · denotes the induced operator norm of operators mapping L 2 µ 0 to L 2 µ 1 . As an implication of the Eckart-Young theorem [HE15,Theorem 4.4.7], the solution of (19) can explicitly be given through singular value decomposition of T ; yielding the variational approach for Markov processes (VAMP) [WN17]. More precisely, the k largest singular values σ 1 ≥ . . . ≥ σ k of T have right and left singular vectors solves (19), see Appendix A.

Coherent sets
Similarly to the reversible equilibrium case with pronounced metastability in section 3.2, it is also possible in the time-inhomogeneous case to give our GMSM (18) from section 4.1 a physical interpretation-under some circumstances.
In the reversible equilibrium situation, recall from (13) (14), we modify the projections Q 0 , Q 1 from (20) to Q 0 : L 2 µ 0 → V 0 , Q 1 : L 2 µ 1 → V 1 by using φ i and ψ i instead of φ i and ψ i , respectively, and define the modified GMSM by T k = Q 1 T Q 0 . An analogous computation to (15) yields for the matrix representation T k of the restriction T k : In other words, the entries of T k contain the transition probabilities from the sets M 0,i (at initial time) into the sets M 1,j (at final time). Thus, T k has the physical interpretation of a MSM, with the only difference to the reversible stationary situation being that the "metastable" sets at initial and final time are different. This can be seen as a natural reaction to the fact that in the time-inhomogeneous case the dynamical environment (e.g., the potential energy landscape governing the dynamics of a molecule) can change in time.  This indicates that a rank-3 GMSM is sufficient to approximate the system, and that we have three show coherent sets at initial time, and left singular vectors the associated coherent sets at final time. 242 We can identify the three wells as three coherent sets. Figure 4 shows that they are coherent indeed: the respective parts of the initial ensemble µ 0 is to a large extent mapped onto the corresponding part of the final ensemble µ 1 , cf. Figure 2 and (23). Computing the MSM from the transition probabilities between the coherent sets as in (22)

Example: diffusion in shifting triple-well potential
Let us consider the diffusion (1) in the time-dependent potential landscape with β = 5 and on the time interval [t 0 , t 1 ] = [0, 10]; cf. Figure 3 (left). Taking the initial distribution µ 0 ∝ exp(−βW (0, ·)), we build the transfer operator (17), and consider its dominant singular values: This indicates that a rank-3 GMSM is sufficient to approximate the system, and that we have three coherent sets. We observe the characteristic almost constant behavior (21) of the left and right singular vectors over the respective coherent sets; Figure 3 (middle and right). Recall that right singular vectors show coherent sets at initial time, and left singular vectors the associated coherent sets at final time. We can identify the three wells as three coherent sets. Figure 4 shows that they are coherent indeed: the respective parts of the initial ensemble µ 0 is to a large extent mapped onto the corresponding part of the final ensemble µ 1 , cf. Figure 2 and (23). Computing the MSM from the transition probabilities between the coherent sets as in (22) gives the stochastic matrix The singular values of T k as mapping from the µ 0 -weighted R 3 to the µ 1 -weighted R 3 are σ 1 = 1.000, σ 2 = 0.733, σ 3 = 0.534 ; they are in good agreement with the true singular values of T . We repeat the computation with a different initial distribution µ 0 , where only the left and right well are initially populated, as shown in Figure 5. The largest singular values of T , σ 1 = 1.000, σ 2 = 0.643, σ 3 = 0.030, σ 4 ≈ 0 , already show that there are only two coherent sets, as the third singular value is significantly smaller than the second one. The left well forms one coherent set, and the union of the middle and right ones form the second coherent set.

Data-based approximation
Setting and auxiliary objects. We would like to estimate the GMSM (18) from trajectory data. In the time-inhomogeneous setting, let us assume that we have m data points x 1 , . . . , x m at time t 0 , and their (random) images y 1 , . . . , y m at time t 1 . We can think of the empirical distribution of the x i and y i being estimates of µ 0 and µ 1 , respectively. Let us further define two sets of basis functions χ 0,1 , . . . , χ 0,n and χ 1,1 , . . . , χ 1,n , which we would like to use to approximate the GMSM. If we would like to estimate the first k dominant modes, the least requirement is n ≥ k; in general we have n k. The vector-valued functions are basis functions at initial and final times, respectively. One can take χ 0 and χ 1 to have different lengths too, we just chose them to have the same lengths for convenience. Now we can define the data matrices The following n × n correlation matrices C 00 , C 01 , C 11 will be needed later: C 00,ij = χ 0,i , χ 0,j µ 0 , C 01,ij = χ 1,i , T χ 0,j µ 1 , C 11,ij = χ 1,i , χ 1,j µ 1 .
Their Monte Carlo estimates from the trajectory data are given by products of the data-matrices, as Note that the approximations in (24) become exact if we take µ 0 , µ 1 to be the empirical distributions µ 0 = 1 m i δ(· − x i ) and µ 1 = 1 m i δ(· − y i ), where δ(·) denotes the Dirac delta. We assume that C 00 , C 11 , just as their data-based approximations in (24) are invertible. If they are not, all the occurrences of their inverses below need to be replaced by their Moore-Penrose pseudoinverses. Alternatively, one can also discard basis functions that yield redundant information, until C 00 , C 11 are invertible. Further strategies to deal with the situation where the correlation matrices are singular or illconditioned can be found in [WNP + 17].
Projection on the basis functions. To find the best GMSM representable with the bases χ 0 and χ 1 , we would like to solve (19) under the constraint that the ranges of Q 0 and Q 1 are in W 0 := span(χ 0 ) and W 1 := span(χ 1 ), respectively. To the knowledge of the authors it is unknown whether this problem has an explicitly computable solution, because it involves a non-trivial interaction of W 0 , W 1 and T .
Instead, we will proceed in two steps. First, we compute the projected transfer operator T n = Π 1 T Π 0 , where Π 0 and Π 1 are the µ 0 -and µ 1 -orthogonal projections on W 0 and W 1 , respectively. Second, we reduce T n to its best rank-k approximation T k (best in the sense of density propagation).
Thus, the restriction T n to W 0 → W 1 is simply the µ 1 -orthogonal projection of T on W 1 , giving the characterization It is straightforward to compute that with respect to the bases χ 0 and χ 1 the matrix representation T n of T n : W 0 → W 1 is given by see [WN17].
Best low-rank approximation. To find the best rank-k projection of T n , let us now switch to the basesχ 0 = C −1/2 00 χ 0 andχ 1 = C −1/2 11 χ 1 . We can switch between representations with respect the these bases by kχ0,k , wherec = C 1/2 00 c , and similarly for χ 1 andχ 1 . Again, a direct calculation shows thatχ 0 andχ 1 build orthonormal bases, i.e., χ 0,i ,χ 0,j µ 0 = δ ij and χ 1,i ,χ 1,j µ 1 = δ ij . This has the advantage, that for any operator S n : W 0 → W 1 having matrix representation S n with respect to the basesχ 0 andχ 1 we have where · 2 denotes the spectral norm of a matrix (i.e., the matrix norm induced by the Euclidean vector norm). The matrix representation of T n in the new bases is However, finding now the best rank-k approximation T k of T n amounts, written in these new bases, to Again, by the Eckart-Young theorem [HE15,Theorem 4.4.7], the solution to this problem is given byT whereŨ ,Ṽ ∈ R n×k are the matrices with columns being the right and left singular vectors ofT n to the largest k singular values σ 1 ≥ . . . ≥ σ k , andΣ is the diagonal matrix with these singular values on its diagonal. Thus, the best GMSM in terms of propagation error is given with respect to the bases χ 0 and χ 1 by The resulting algorithm to estimate the optimal GMSM is now identical to the timelagged canonical correlation algorithm (TCCA) that results from VAMP and is described in [WN17].
Algorithm 1 TCCA algorithm to estimate a rank-k GMSM.
2. Estimate the correlation matrices C 00 , C 01 , C 11 from data, as in (24). Remark 5.1 (Reversible system with equilibrium data): If the system in consideration is reversible, the data samples its equilibrium distribution, i.e., µ 0 = µ 1 = µ, and also χ 0 = χ 1 , then C 00 = C 11 , and by the self-adjointness of T from (6) we have C 01 = C T 01 . Thus,T n in (28) is a symmetric matrix, and as such, its singular value and eigenvalue decompositions coincide. Hence, the construction for the best GMSM in this section (disregarding the projection on the basis functions) coincides with the one from section 3. This is not surprising, as both give the best model in terms of propagation error.

Build the projectionT
Remark 5.2 (Other data-based methods): The approximation (26) of the transfer operator has natural connections to other data-based approximation methods. It can be seen as a problem-adapted generalization of the so-called Extended Dynamic Mode Decomposition (EDMD) [WKR15,KKS16]. Strictly speaking, however, EDMD uses an orthogonal projection with respect to the distribution µ 0 of the initial data {x i }, and so approximation (32) below is equivalent to it. EDMD has been shown in [KNK + 18] to be strongly connected to other established analytic tools for (molecular) dynamical data, such as time-lagged independent component analysis (TICA) [PHPG + 13, SP13], blind source separation [MS94], and the variational approach to conformation analysis [NN13].

Time-homogeneous systems and non-stationary data
In this final section we illustrate how the above methods can be used to construct a GMSM for and assess properties of a stationary system, even if the simulation data at our disposal does not sample the stationary distribution of the system. In the first example we reconstruct the equilibrium distribution of a reversible system-hence we are able to build an equilibrium GMSM. In the second example we approximate a nonreversible stationary system (i.e., detailed balance does not hold) by a (G)MSM, again from non-stationary data.
Of course, all the examples presented so far can also be computed by the data-based algorithm of section 5.

Equilibrium MSM from non-equilibrium data
When working with simulation data, we need to take into account that this data might not be in equilibrium. Then, obviously, the empirical distribution does not reflect the stationary distribution of the system. In general, any empirical statistical analysis (e.g., counting transitions between a priori known metastable states) will be biased in such a case.
Let us consider a reversible system with equilibrium distribution µ, and let the available trajectory data be µ ref -distributed. Then, it is natural to describe the system by its transfer operator T ref : Note that µ corr := µ/µ ref is the stationary distribution of this transfer operator, hence we can retrieve the equilibrium distribution of the system by correcting the reference distribution, µ = µ corr µ ref .
In the data-based context, we choose the same basis χ 0 = χ 1 for initial and final times, since the system is time-homogeneous. In complete analogy to (26) above, the µ ref - We will now apply this procedure to the double-well system from section 3.3 with initial points x 1 , . . . , x m distributed as shown in Figure 6 (gray histogram). We chose the number of points to be m = 10 5 , the basis functions χ 0,i to be indicator functions of subintervals of an equipartition of the interval [−2, 2] into n = 100 subintervals, and the lag time τ = 10. In a preprocessing step we discard all basis functions that do not have any of the points x i in their support, thus obtaining a non-singular C 00 , and use the remaining 77 to compute T n,ref ∈ R 77×77 . We obtain λ 2 = 0.894 giving a time scale It is now simple to reconstruct the approximation T n of T , the transfer operator with respect to the equilibrium density. Let D corr denote the diagonal matrix with the elements of µ corr as diagonal entries. Then, T n = D −1 corr T n,ref D corr approximates the matrix representation of T n with respect to our basis of step functions.
Remark 6.1 (Koopman reweighting): One can make use of the knowledge that the system that one estimates is reversible, even though due to the finite sample size m this is not necessarily valid for T ref,n . In [WNP + 17], the authors add for each sample pair (x i , y i ) also the pair (x i+m = y i , y i+m = x i ) to the sample set, thus numerically forcing the estimate to be reversible. In practice, one defines the diagonal matrix W with diagonal χ T µ corr , builds the reweighted correlation matricesC 00 = 1 2 (χ 0 W χ T 0 + χ 1 W χ T 1 ) andC 01 = 1 2 (χ 1 W χ T 0 + χ 0 W χ T 1 ), and uses them instead of C 00 , C 01 .
6.2. A non-reversible system with non-stationary data Reversible dynamics gives rise to self-adjoint transfer operators, and their theory of Markov state modeling is well developed. However, transfer operators of non-reversible systems are not self-adjoint, hence their spectrum is in general not purely real-valued. Thus, the definition of time scales, and in general the approximation by GMSMs is not fully evolved. Complex eigenvalues indicate cyclic behavior of the process. As this topic is beyond the scope of this paper, we refer the reader to [DJ99,DCWS16,FK17] and to [CBS15,KS15] for Markov state modeling with cycles. We will consider a non-reversible system here, and show that restricting its behavior to the dominant singular modes of its transfer operator is able to reproduce its dominant long-time behavior, and even allows for a good, few-state MSM. Note that the best rankk GMSM (18) maps to the k-dimensional subspace V 1 of left singular vectors, thus its eigenvectors also fall into this subspace.
The system in consideration consists of two driving "forces", one is a reversible part F r (x) = −∇W (x) coming from the potential and the other is a circular driving given by where β = 2 is the inverse temperature, as in (1). The dynamics now is governed by the SDE dx t = (F r + F c )(x t ) dt + 2β −1 dw t . It is a diffusion in a 7-well potential (the wells are positioned uniformly on the unit circle) with an additional clockwise driving that is strongest along the unit circle and decreases exponentially in the radial distance from this circle. For our data-based analysis we simulate a trajectory of this system of length 500 and sample it every 0.01 time instances to obtain an initial set of 5 · 10 4 points. Every point herein is taken as initial condition of 100 independent simulations of the SDE for lag time τ = 1, thus obtaining 5 · 10 6 point pairs (x i , y i ). We observe in Figure 7 (left) that the empirical distribution of the x i did not yet converge to the invariant distribution of the system, which would populate every well evenly.
To approximate the transfer operator we use χ = χ 0 = χ 1 consisting of the characteristic functions of a uniform 40 × 40 partition of [−2, 2] × [−2, 2], and restrict this basis set to those 683 partition elements that contain at least one x i and y j . The associated projected transfer operator, T n from (26) is then used to computeT n from (28), and its singular values σ 1 = 1.000, σ 2 = 0.872, σ 3 = 0.588, . . . , σ 7 = 0.383, σ 8 = 0.052, indicating a gap after seven singular values. Thus, we assemble a rank-7 GMSM T k via (30). This GMSM maps L 2 µ 0 to L 2 µ 1 , thus to make sense of its eigenmodes, we need to transform its range to densities with respect to µ 0 instead of µ 1 . As a density u with respect to µ 1 is made by µ 1 u µ 0 to a density with respect to µ 0 , T k = C −1 00 C 11 T n rescales the GMSM to map L 2 µ 0 to itself. 4 We are also interested in the system's invariant distribution. As in section 6.1, we can correct the reference distribution µ ref = µ 0 by the first eigenfunction µ corr of T k to yield the invariant distribution µ = µ corr µ ref , cf. Figure 7 (middle). The dominant eigenvalues of T k are λ k,1 = 0.998 + 0.000i, λ k,2/3 = 0.803 ± 0.261i, λ k,4/5 = 0.511 ± 0.230i, λ k,6/7 = 0.378 ± 0.077i , Note that λ k,1 < 1. This is due to our restriction of the computation to certain partition elements, as specified above. This set of partition elements is not closed under the process dynamics, and this "leakage of probability mass" (about 0.2%) is reflected by the dominant eigenvalue. All eigenvalues of T k are within 0.5% error from the dominant eigenvalues of the transfer operator T n with respect to the stationary distribution (projected on the same basis set, and computed with higher accuracy), which is a surprisingly good agreement.
The 8-th eigenvalue of T n is smaller in magnitude than 0.03. As indicated by this spectral gap, we may obtain a few-state MSM T k here as well. To this end we need to find "metastable sets" (although in the case of this cyclically driven system the term metastability is ambiguous) on which we can project the system's behavior. Let v i = (v i,1 , . . . , v i,n ) T denote the i-th eigenvector of T k . As in the reversible case, where eigenvectors are close to constant on metastable sets, we will seek also here for regions that are characterized by almost constant behavior of the eigenvectors. More precisely, if the p-th and q-th partition elements belong to the same metastable set, then we expect v i,p ≈ v i,q for i = 1, . . . , 7. Thus, we embed the p-th partition element into C 7 ≡ R 14 (i.e., a complex number is represented by two coordinates: its real and imaginary parts) by p → (v 1,p , . . . , v 7,p ) T , and cluster the hence arising point cloud into 7 clusters by the k-means clustering algorithm. 5 The result is shown in Figure 7 (right). Taking these sets we can assemble the MSM T k ∈ R 7×7 via (15). We obtain a MSM that maps a Markov state (i.e., a cluster) with probability 0.62 to itself, with probability 0.29 to the clockwise next cluster, and with probability 0.06 to the second next cluster in clockwise direction. The probability to jump one cluster in the counterclockwise direction is below 0.001. The eigenvalues of T k , λ 1 = 0.998 + 0.000i, λ 2/3 = 0.800 ± 0.260i, λ 4/5 = 0.507 ± 0.227i, λ 6/7 = 0.374 + 0.076i , are also close to those of T n (below 1% error), justifying this MSM. optimal rank-k approximation A k of A in the sense that where · denotes the induced operator norm, is given by where σ i , ψ i , φ i are the singular values (in non-increasing order), left and right normalized singular vectors of A, respectively. The optimum is unique iff σ k > σ k+1 .
Proof. Let A k be defined as in (33).
Let now B k be any rank-k operator from H 0 to H 1 . Then, there exist linear functionals c i : H 0 → R and vectors v i ∈ H 1 , i = 1, . . . , k, such that For every i, since c i has one-dimensional range, its kernel has co-dimension 1, thus the intersection of the kernels of all the c i has co-dimension at most k. Thus, any (k + 1)dimensional space has a non-zero element w with c i (w) = 0 for i = 1, . . . , k. By this, we can find scalars γ 1 , . . . , γ k+1 such that k+1 i=1 γ 2 i = 1 and w = γ 1 φ 1 + . . . + γ k+1 φ k+1 satisfies c i (w) = 0 for i = 1, . . . , k. By construction w 0 holds. It follows that This with (34) proves the claim.
As a corollary, if A : H → H is a self-adjoint operator, then its eigenvalue and singular value decompositions coincide, giving ψ i = φ i , and thus A k in (33) is the projection on the dominant eigenmodes.