Spectral Properties of Effective Dynamics from Conditional Expectations

TL;DR

This paper investigates the spectral properties of reduced effective dynamics derived from high-dimensional stochastic systems, providing bounds and numerical insights into their approximation accuracy, especially for reversible and Langevin systems.

Contribution

It introduces a new relative error bound for the spectrum of effective generators in reversible systems and analyzes spectral approximation methods for Langevin dynamics.

Findings

Eigenfunction approximation error bounds for reversible systems

KM-type spectral approximation is robust to time window choice

Effective dynamics can be meaningfully defined for underdamped Langevin systems

Abstract

The reduction of high-dimensional systems to effective models on a smaller set of variables is an essential task in many areas of science. For stochastic dynamics governed by diffusion processes, a general procedure to find effective equations is the conditioning approach. In this paper, we are interested in the spectrum of the generator of the resulting effective dynamics, and how it compares to the spectrum of the full generator. We prove a new relative error bound in terms of the eigenfunction approximation error for reversible systems. We also present numerical examples indicating that if Kramers--Moyal (KM) type approximations are used to compute the spectrum of the reduced generator, it seems largely insensitive to the time window used for the KM estimators. We analyze the implications of these observations for systems driven by underdamped Langevin dynamics, and show how meaningful effective dynamics can be defined in this setting.

Equations

$$d{X}_t$$

$$d{X}_t=-\frac{1}{\gamma }\nabla V(X_t)dt+\sqrt{2{\beta }^{-1}{\gamma }^{-1}}d{B}_t,$$

$$0<{\eta }_1\Vert v{\Vert }^2\le v^{T}a(x)v\le {\eta }_2\Vert v{\Vert }^2,$$

$$⟨\psi ,\tilde{\psi }{⟩}_{\mu }={\int }_{\Omega }\psi (x)\tilde{\psi }(x)\mathrm{d}\mu (x)={\int }_{\Omega }\psi (x)\tilde{\psi }(x)\rho (x)\mathrm{d}x.$$

$${\mathcal{K}}_t\psi (x)={\mathbb{E}}^x\left[\psi (X_t)\right],$$

$$\mathcal{L}\psi (x)=\frac{d}{dt}{\mathbb{E}}^x\left[\psi (X_t)\right]{|}_{t=0}.$$

$$\frac{d}{dt}{\mathcal{K}}_t\psi =\mathcal{L}{\mathcal{K}}_t\psi .$$

$$\mathcal{L}\psi (x)$$

$$=b(x)\cdot \nabla \psi (x)+\frac{1}{2}a(x):{\nabla }^2\psi (x).$$

$$⟨\mathcal{L}\psi ,\tilde{\psi }{⟩}_{\mu }$$

$$\mathcal{Q}(\psi ,\tilde{\psi }):=\frac{1}{2}{\int }_{\Omega }a(x)\nabla \psi (x)\cdot \nabla \tilde{\psi }(x)\mathrm{d}\mu (x).$$

$$\Vert \cdot {\Vert }_{\mathcal{Q}}=\mathcal{Q}(\cdot ,\cdot {)}^{1/2}.$$

$$-\mathcal{L}{\psi }_i={\kappa }_i{\psi }_i.$$

$${\lambda }_i(t)=e^{-{\kappa }_{i}t}.$$

$$t_i=\frac{1}{{\kappa }_i}.$$

$$\mathcal{P}\psi (z)={\mathbb{E}}^{\mu }\left[\psi (x)\mid \xi (x)=z\right]$$

$$\mathcal{P}\psi (z)=\frac{1}{\vartheta (z)}{\int }_{\Omega }\psi (x)\delta (\xi (x)-z)\mathrm{d}\mu (x).$$

$$L_{\nu }^2=\{\phi :\hat{\Omega }\to \mathbb{R},{\int }_{\hat{\Omega }}{\phi }^2(z)\mathrm{d}\nu (z)<\infty \}.$$

$${\mathcal{L}}^{\xi }$$

$${\mathcal{L}}^{\xi }\phi =\mathcal{P}\left[\mathcal{L}\xi \right]\cdot {\nabla }_z\phi +\frac{1}{2}\mathcal{P}\left[\nabla {\xi }^{T}a\nabla \xi \right]:{\nabla }_z^2\phi .$$

$$b^{\xi }(z)$$

$${\mathcal{Q}}^{\xi }(\phi ,\tilde{\phi })=\frac{1}{2}{\int }_{\hat{\Omega }}{a}^{\xi }(z){\nabla }_z\phi (z)\cdot {\nabla }_z\tilde{\phi }(z)\mathrm{d}\nu (z).$$

$$-{⟨\mathcal{L}{\hat{\psi }}_i,{\varphi }_j⟩}_{\mu }=Q({\hat{\psi }}_i,{\varphi }_j)$$

$${\kappa }_i\le {\omega }_i\le {\hat{\omega }}_i.$$

$$\mathcal{Q}({\varphi }_i,{\varphi }_j)$$

$$-{⟨\mathcal{L}{\varphi }_i,{\varphi }_j⟩}_{\mu }$$

$$d_s^{\xi }(X_t)=\xi (X_{t+s})-\xi (X_t).$$

$$b^{\xi }(z)$$

$$-{⟨\mathcal{L}{\varphi }_i,{\varphi }_j⟩}_{\mu }\approx -\frac{1}{M}\sum _{m=1}^M\left[\frac{1}{s}{d}_s^{\xi }(x_m){\nabla }_z{\varphi }_i(x_m)+\frac{1}{2s}(d_s^{\xi }(x_m)\otimes d_s^{\xi }(x_m)):{\nabla }_z^2{\varphi }_i(x_m)\right]{\varphi }_j(x_m).$$

$$E_i=\frac{{\omega }_i-{\kappa }_i}{{\omega }_i}$$

$$E_i=\frac{{\omega }_i-{\kappa }_i}{{\omega }_i}\le C\Vert {\mathcal{P}}_{\mathcal{Q}}^{\perp }{\psi }_i{\Vert }_{\mathcal{Q}}^2.$$

Full text

1 Introduction

The description of high-dimensional dynamical systems by a reduced set of variables, usually referred to as coarse graining or model reduction, is of tremendous importance across many different fields of research. Examples range from finance to climate modeling to molecular biology. From the huge body of literature on the subject, we mention in particular the Mori–Zwanzig formalism [1, 2, 3, 4, 5] as well as the framework of averaging and homogenization for systems with explicit multiscale structure [6, 7]. Within the field of molecular physics, references [8, 9, 10, 11, 12, 13] present important contributions to this line of research. Here, we focus on model reduction for stochastic differential equations (SDEs), and follow another standard approach, which is based on conditioning along level sets of the coarse graining map [14, 15, 16]. For a detailed theoretical analysis of the method in the context of SDEs, please see [15, 17, 18, 19, 20].

For a given coarse grained description of a system, a fundamental question to address is the quality of approximation of the full system by means of the reduced system, as measured by a suitable metric (which usually depends on the problem at hand). In many cases, the approximation of spectral properties of the system’s generator is useful in this context. The generator and its associated semigroup, also called Koopman semigroup, are used to describe the time evolution of expectation values of observable functions. For metastable systems, the leading generator eigenpairs provide information on slow modes in the dynamical system. Spectral approximation results for the conditioning approach have been obtained in [17, 18].

Another important problem to consider is the analysis and parameter estimation of coarse grained models based on simulation data of the full system. In recent years, a variety of methods has been developed to learn models for the Koopman semigroup off simulation data, see [21, 22, 23, 24, 25, 26, 27] and the references therein. In Ref. [28], some of the authors of the present study presented a conceptually simple framework for the data-driven approximation of the Koopman generator, called gEDMD. This framework can also be used to identify and analyze coarse grained models within the context of the conditioning approach. The gEDMD method requires knowledge of the full system parameters. If these parameters are unknown, they need to be replaced by a suitable approximation, such as Kramers–Moyal (KM) formulae. These are based on averages of finite differences at a finite offset (time window). The quality of this approximation as a function of the offset will be addressed in this paper. Even though we focus on KM estimators here, let us mention that a multitude of more advanced methods for parameter estimation of stochastic dynamics are available, please see Ref. [29] for an overview. Spectral methods have been considered in [30, 31], while particular attention to the choice of time window has been paid in [32, 7]. The Kramers–Moyal formulae being among the simplest estimators, we take them as the starting point for our study.

The third focus of this study is model reduction for systems driven by underdamped Langevin dynamics, which is a widely used model, especially in molecular and biological physics. As the momentum variables of these dynamics often play just an auxiliary role, an interesting question to address is how to define a reduced dynamics that only involves the position state variables. As the conditioning approach does not provide meaningful answers in this case, finding meaningful effective equations remains an open problem in this setting [33, 34, 35].

In this paper, we report theoretical and numerical results on the issues raised above. The contributions of this study are as follows:

• •

Concerning the first problem, we prove a new relative error bound for the approximation of generator eigenvalues by the coarse grained generator, if the dynamics is reversible (Proposition 6.3). This bound shows that a small projection error of the full eigenfunctions with respect to the energy norm is required for a small eigenvalue error. We also derive conditions to ensure that the spectrum of the reduced generator is discrete in the first place (Proposition 6.2).

• •

Concerning the second issue, we present numerical examples indicating that if KM estimators are used within the gEDMD algorithm for reversible systems, on a good set of reaction coordinates, then the resulting eigenvalue estimates seem to be fairly insensitive to the offset used for the KM estimators (Sections 4.2 and 4.3, Conjecture 4.1).

• •

Thirdly, we suggest that if the observations of the second part can be confirmed theoretically, it is possible to use KM estimators at large offsets to define meaningful effective equations for underdamped dynamics (Corollary 5.1). The reason is that the statistics of the underdamped process approach those of an overdamped process after a suitable re-scaling of time. We provide successful illustrations of this idea using a toy example and molecular dynamics simulation data of the alanine dipeptide (Sections 5.2 and 5.3).

The rest of this paper is organized as follows: in Section 2, we recap what is needed of the theory of stochastic differential equations, their generators, the conditioning approach, and the data-driven approximation of Koopman generators. In Section 3, we present and illustrate our spectral approximation result. The technical details of the proofs are deferred to Section 6. In Section 4, we present numerical results on the spectrum of gEDMD models based on KM estimators for reversible systems, and conjecture that the observed behaviour can be expected in general. We analyze the implications of this hypothesis for systems driven by underdamped Langevin dynamics in Section 5, and provide additional numerical results for this setting. Conclusions and outlook follow in Section 7.

2 Concepts

2.1 SDEs and Generators

In this paper, we consider a reversible Markov process $X_{t}$ attaining values in a domain $\Omega\subset\mathbb{R}^{d}$. The process is governed by the stochastic differential equation

[TABLE]

Here, $B_{t}$ denotes $d$-dimensional Brownian motion, the function $b:\mathbb{R}^{d}\mapsto\mathbb{R}^{d}$ is called the drift, and $\sigma:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times d}$ is called the diffusion. We use the notation $a(x)=\sigma(x)\sigma(x)^{T}$ for the covariance matrix of the diffusion. A standard example for dynamics of type Eq. (1) are the overdamped Langevin dynamics

[TABLE]

where $V:\,\Omega\rightarrow\mathbb{R}$ is a scalar function called the potential, while $\beta,\,\gamma$ are constants corresponding to the inverse temperature and the friction in physics applications.

We assume that $X_{t}$ is ergodic with respect to a unique invariant measure $\mu$ with Boltzmann density $\rho\propto\exp(-F(x))$, where $F$ is called a generalized potential. Also, we assume the diffusion to satisfy a so-called uniform ellipticity condition

[TABLE]

for constants $\eta_{1},\eta_{2}>0$. The invariant measure $\mu$ gives rise to the Hilbert space $L_{\mu}^{2}$ of all square integrable functions with respect to that measure. We can think of these functions as physical observables. The inner product on $L^{2}_{\mu}$ is given by

[TABLE]

For a fixed time window $t\geq 0$, and an observable function $\psi\in L^{2}_{\mu}$, the Koopman operator $\mathcal{K}_{t}$ describes the evolution of the expectation value of $\psi$ by means of the dynamics (1):

[TABLE]

where $\mathbb{E}^{x}[\cdot]$ denotes expectation given that the dynamics starts deterministically at $x$. The infinitesimal generator $\mathcal{L}$ of the Markov process $X_{t}$ is then defined as a formal time-derivative of this expectation value:

[TABLE]

It follows that, by describing the system in terms of the expectations $\mathcal{K}_{t}\psi$, the non-linear dynamical system (1) turns into a linear, but infinite-dimensional system with differential equation

[TABLE]

For this reason, the Koopman operators $\mathcal{K}_{t}$ and their generator $\mathcal{L}$ have been studied extensively in past decades. A technical subtlety arising from this infinite-dimensional description is that the time-derivative (5) is not well-defined for all $\psi\in L^{2}_{\mu}$. For smooth and compactly supported functions however, stochastic calculus shows that $\mathcal{L}$ is well-defined, and acts as a second order differential operator:

[TABLE]

Here, $\nabla^{2}\psi$ is the Hessian matrix of the function $\psi$, and the colon $:$ denotes the Frobenius inner product between matrices, i.e., $A:B=\sum_{i,j}A_{ij}B_{ij}$. For the same class of functions, the generator is symmetric with respect to the inner product (4), and satisfies the important equality

[TABLE]

which requires only first order derivatives. The negative of the right-hand side of Eq. (8) is called the quadratic form

[TABLE]

Tools from functional analysis [36, 37] can be used to define (9) on a larger set of functions $\mathbb{V}_{\mu}$, usually called the form domain. Because of (3), $Q$ defines an inner product on the form domain, and $\mathbb{V}_{\mu}$ in fact turns into a Hilbert space with energy norm

[TABLE]

The reason to introduce all these concepts is that the energy norm will serve as error measure for the main results of this study.

2.2 Spectral Decomposition

We are particularly interested in eigenvalues and eigenfunctions of the (negative) generator $-\mathcal{L}$. Because of (8), the spectrum of $-\mathcal{L}$ must be part of the non-negative real axis. We will further assume that there is a complete set of eigenfunctions (i.e., they form a basis of $L_{\mu}^{2}$) corresponding to discrete eigenvalues. In other words, there are functions $\psi_{0},\psi_{1},\psi_{2},\ldots$ and non-negative numbers $\kappa_{0}<\kappa_{1}<\kappa_{2}<\ldots$ such that

[TABLE]

Conditions for the existence of a completely discrete spectrum are discussed in section 6. The assumption that all eigenvalues $\kappa_{i}$ are distinct is for simplicity only, especially with regards to the analysis in section 6. It follows again from (8) that $\kappa_{0}=0$, and $\psi_{0}\equiv 1$ is the constant function.

The physical significance of these eigenpairs is that by [38, Ch 2, Thm 2.4.], the eigenfunctions $\psi_{i}$ are also eigenfunctions of the Koopman operators $\mathcal{K}_{\tau}$ for all $t\geq 0$, corresponding to eigenvalues

[TABLE]

Due to the exponential decay of all $\lambda_{i}(t)$, it is common to refer to the $\kappa_{i}$ as rates, and to their reciprocals as implied timescales

[TABLE]

In many applications, including molecular dynamics, we expect to find a number $K$ of dominant rates $0<\kappa_{1}<\ldots<\kappa_{K}\ll\kappa_{K+1}$ separated from all others. These dominant spectral components are of particular interest as they are related to metastability, that is, the existence of long-lived macrostates such that transitions between those states are rare events [39, 22, 40].

2.3 Dimensionality Reduction

The main topic of this study is the effect of dimensionality reduction on the generator eigenvalues $\kappa_{i}$ introduced above. Following the notation of Refs. [15, 17, 18], we consider a smooth coarse graining function $\xi$, which maps the state space $\Omega\subset\mathbb{R}^{d}$ onto a lower-dimensional space $\hat{\Omega}\subset\mathbb{R}^{m}$, where $m\leq d$. For a position $z\in\hat{\Omega}$ in reduced space, we denote the marginal probability distribution of the invariant measure $\mu$ by $\nu$, and assume it possesses a corresponding density function $\vartheta(z)$. The conditional expectation operator

[TABLE]

computes the stationary average of a function $\psi$ defined on $\Omega$, conditional to $\xi$ attaining a fixed value $z\in\hat{\Omega}$. Using the Dirac $\delta$-function, we can informally write the above expression as

[TABLE]

Consider the space $L^{2}_{\nu}$ of physical observables on reduced space $\hat{\Omega}$:

[TABLE]

By the concatenation $\varphi\circ\xi$, every such function can be viewed as a function on full state space $\Omega$. In fact, $L^{2}_{\nu}$ can be exactly identified as the subspace of functions in $L^{2}_{\mu}$ which depend only on the value of $z=\xi(x)$, with the conditional expectation operator acting as orthogonal projection onto this subspace [17]. Note that unless $\xi$ is constant, $L_{\nu}^{2}$ is an infinite-dimensional subspace.

This motivates consideration of the projected generator

[TABLE]

As discussed comprehensively in [17], this operator retains the shape of the generator of a reversible Markov process $Z_{t}$ on $\hat{\Omega}$, as for smooth compactly supported functions $\varphi\in L^{2}_{\nu}$, we have that

[TABLE]

In the above equation, $\mathcal{L}\xi$ is an $m$-dimensional vector, each entry containing the application of $\mathcal{L}$ to each component of $\xi$, and $\nabla\xi$ is the $d\times m$ Jacobian matrix of $\xi$. The coefficients

[TABLE]

serve as effective drift and effective diffusion for the process $Z_{t}$, respectively. It can also be shown that $Z_{t}$ is ergodic with respect to $\nu$ [17], and we can associate to $\mathcal{L}^{\xi}$ a form domain $\mathbb{V}_{\nu}$ with an effective quadratic form

[TABLE]

If the projected generator also possesses a discrete spectrum with eigenvalues $\omega_{i},\,i=0,1,\ldots$, comparison of those eigenvalues with the original ones $\kappa_{i}$ provides information about how well the effective dynamics $Z_{t}$ retain the relaxation processes of the original process. Our results on this topic are presented in sections 3 and 6.

2.4 Galerkin Approximation

The numerical approximation of the eigenfunctions $\psi_{i},\,i=1,\ldots,K$ is often achieved by Galerkin projection, i.e., orthogonal linear projection of the generator to a finite-dimensional subspace. After choosing such a space $\mathbb{W}\subset\mathbb{V}_{\mu}$, with a basis set $\{\phi_{i}\}_{i=1}^{N}$, the Galerkin approach consists of finding $\hat{\psi}_{i}\in\mathbb{W}$ such that

[TABLE]

Eq. (17) is a generalized matrix eigenvalue problem. The approximate eigenvalues $\hat{\omega}_{i}$ are also called Ritz values associated to $\mathbb{W}$. If $\mathbb{W}$ is chosen as a subspace of $\mathbb{V}_{\nu}$ (given a set of reduced variables as described in the previous section), it was shown in [17] that (17) serves both as a weak form for $\mathcal{L}$ in $\mathbb{V}_{\mu}$ and for $\mathcal{L}^{\xi}$ in $\mathbb{V}_{\nu}$. Importantly, the min-max-principle implies that for any such subspace

[TABLE]

Moreover, if simulation data $\{x_{m}\}_{m=1}^{M}$ of the full process (1), approximately sampling the invariant measure $\mu$, is available, the following empirical estimators

[TABLE]

will converge to the terms in Eq. (17) in the limit of infinite data. This method has been called generator Extended Dynamic Mode Decomposition (gEDMD) [28]. Consequently, for a subspace $\mathbb{W}\subset\mathbb{V}_{\nu}$ comprised of functions on reduced space, gEDMD will simultaneously approximate the eigenvalues of $\mathcal{L}$ and $\mathcal{L}^{\xi}$ by (18). Note that in this setting, simulation data of the full dynamics (1) can still be used in Eq. (19), it is not necessary to simulate the effective dynamics defined by $\mathcal{L}^{\xi}$ first. As the gradients in (19) are taken with respect to the full state variables $x$, we have to apply the chain rule when evaluating (19) for basis functions $\phi_{i}\in\mathbb{W}\subset\mathbb{V}_{\nu}$: $\nabla\phi_{i}(\xi(x))=\nabla_{z}\phi_{i}(z)^{T}\nabla\xi^{T}(x)$.

Note that the estimator for $\mathcal{Q}$ in (19) exploits reversibility of the process. An alternative estimator, which can also be applied to non-reversible processes, is

[TABLE]

We will require this last equation in the next section.

2.5 Kramers–Moyal Estimators

The parameters $b^{\xi}$ and $a^{\xi}$ in (15) involve integrals over non-linear manifolds in high-dimensional space, and they are rarely used in practice for this reason. For a process $\{X_{t}\}_{t\geq 0}$ and a positive offset $s>0$, define the projected first order finite difference at time $t$ by

[TABLE]

With this notation, the basic Kramers–Moyal (KM) formulae [41] for the approximation of $b^{\xi}$ and $a^{\xi}$ are given by:

[TABLE]

We note that many more sophisticated approximations can be found in the literature.

KM estimators can also be used in the context of gEDMD. If the parameters $b$ and $a$ of the original process (1) are unknown, gEDMD can be re-formulated upon replacing $b$ and $a$ by the first order finite differences at all data points in (20). If $s>0$ corresponds to a multiple of the integration time step in a discrete trajectory $\{x_{m}\}_{m=1}^{M}$, then (20) can be converted to:

[TABLE]

This estimator is consistent for $s\rightarrow 0$ and $M\rightarrow\infty$. There is a dual meaning to eigenpairs derived from this approximation: on one hand, they serve as an approximation to the eigenpairs of $\mathcal{L}^{\xi}$, at least for small offsets $s$. On the other hand, these are also approximations to the spectrum of the generator of a non-reversible coarse grained SDE with drift and diffusion given by (21), see again [28] for a detailed exposition. In Secs. 4 and 5, we present a numerical study on the effect of the offset $s$ on the spectrum of this generator. We note that similar approximations can certainly be built on more advanced estimators than (21), but the KM formulae will suffice for this study.

3 Spectral Properties of the Projected Generator

3.1 Summary of spectral properties

The first major result of this study concerns the approximation error for the dominant eigenvalues $\kappa_{1},\kappa_{2},\ldots$ by corresponding eigenvalues of $\mathcal{L}^{\xi}$. We only provide a high-level summary of these results here, while the technically more involved statements and their proofs can be found in section 6.

First, we introduce conditions to ensure the spectrum of the effective generator $\mathcal{L}^{\xi}$ is also discrete, see Proposition 6.2. We then show in Proposition 6.3 and Corollary 6.4 that the relative eigenvalue error

[TABLE]

can be bounded in terms of the energy norm of the projection residual $\mathcal{P}_{Q}^{\perp}\psi_{i}=(\mathcal{I}-\mathcal{P}_{Q})\psi_{i}$ of the corresponding eigenfunctions, with $\mathcal{P}_{Q}$ denoting the $\mathcal{Q}$-orthogonal projection onto $\mathbb{V}_{\nu}$:

[TABLE]

In other words, if the eigenfunction $\psi_{i}$ can be written as a function of the reduced variables $z$, up to a small error, then we can expect the eigenvalue $\kappa_{i}$ to be reproduced well by the effective dynamics on $z$. However, as the error is measured by the energy norm, Proposition 6.3 shows that not only the eigenfunction $\psi_{i}$, but also its first order derivatives must be approximated well by functions of $z$ alone. In the next section 3.2, we show that this is not merely an academic condition, but indeed necessary.

Our result improves on existing ones in two ways. First, in Ref. [18, Theorem 2], it was shown that the absolute eigenvalue error of the projected generator is small if $\smash{\|\mathcal{LP}^{\perp}\psi_{i}\|_{L_{\mu}^{2}}}$ and $\smash{\|\mathcal{P}^{\perp}\psi_{i}\|_{L_{\mu}^{2}}}$ are small. Proposition 6.3 (Corollary 6.4) complements these results in the sense that it bounds the relative error of eigenvalues (timescales), which is a more practical error measure for eigenvalues close to zero, i.e., large timescales. For a more detailed elaboration on the relationship of the projection error and the relative error of timescales, please refer to the text after Corollary 6.4. Second, our bound (37) is less restrictive than the conditions in [18], as it uses the energy norm involving only first order derivatives, while the term $\smash{\|\mathcal{LP}^{\perp}\psi_{i}\|_{L_{\mu}^{2}}}$ necessarily requires second derivatives. In fact, the bound assumed in [18, Theorem 2] implies our bound up to another multiplicative constant (thus, it is more restrictive), see Lemma 6.6.

3.2 Illustration of the Error Bound

In order to provide an illustration of the error bound, we consider a two-dimensional Ornstein–Uhlenbeck process, that is, overdamped Langevin dynamics (2) with potential (see Fig. 1 A):

[TABLE]

We set $\alpha_{x}=1,\alpha_{y}=5,\beta=\gamma=1.0$. The eigenvalues and eigenfunctions of this system are known analytically. For all integers $r,s\geq 0$, we have eigenvalues $\kappa_{r,s}=r\alpha_{x}+s\alpha_{y}$, with eigenfunctions $\psi_{r,s}=\frac{1}{\sqrt{r!s!}}H_{r}(x)H_{s}(y)$, using one-dimensional Hermite polynomials $H_{i}$. The first four non-zero eigenvalues correspond to eigenfunctions which are constant in $y$-direction. The first of these eigenfunctions, $\psi_{1,0}=x$, is shown in Fig. 1 C. We now consider a family of reaction coordinates

[TABLE]

For $m=0$, the reaction coordinate $\xi^{0}(x,y)=x$ perfectly captures the first eigenfunction $\psi_{1,0}$. For positive $m$, however, the level sets $\xi^{-1}(z)$ of the coarse graining map oscillate within a vertical strip of width 0.1 around $x=z$, see Fig. 1 A for a comparison of the level sets at $m=0$ and $m=10$. Due to these small scale oscillations, we still expect to find a relatively small projection error of the eigenfunction $\psi_{1,0}$, if measured by the $L^{2}_{\mu}$-norm, but an increasingly larger error if the energy norm (involving derivatives) is employed.

In order to estimate these projection errors, we use finite-dimensional subspaces $\mathbb{W}^{m}\subset\mathbb{V}^{m}_{\nu}$, where $\mathbb{V}^{m}_{\nu}$ is the form domain corresponding to reaction coordinate $\xi^{m}$. The subspaces are spanned by the first ten Hermite polynomials $H_{i}(z)$, which exactly capture the slowest eigenfunction for $m=0$. For each subspace, we calculate the Galerkin matrices in Eq. (17) by numerical integration, and then extract the first non-trivial eigenvalues $\hat{\omega}^{m}_{1}$ and eigenfunctions $\hat{\psi}^{m}_{1}$. We calculate the relative errors $E^{m}_{1}=\frac{\hat{\omega}^{m}_{1}-\kappa_{1,0}}{\hat{\omega}^{m}_{1}}$ and the eigenfunction approximation errors $\|\psi_{1,0}-\hat{\psi}^{m}_{1}\|^{2}$, measured using both the norms on $L^{2}_{\mu}$ and the energy norm. The results are shown in Fig. 1 B. We observe that the $L^{2}_{\mu}$-error remains almost constant as $m$ increases, while the energy norm error increases steadily, reflecting the increasingly oscillating shape of the approximate eigenfunctions $\hat{\psi}^{m}_{1}$ (see Fig. 1 D). In agreement with our error estimate, the energy norm approximation error provides a fairly tight bound for the relative eigenvalue error $E^{m}_{1}$. It should be noted that the quantities shown here only provide upper bounds for the relative error $E_{1}$ and the projection error $\delta_{1}^{2}$, but they suffice for the purpose of illustration.

In summary, this example confirms the eigenvalue error bound provided in Proposition 6.3, and it also highlights the importance of capturing the derivatives of generator eigenfunctions of interest when selecting a reaction coordinate for coarse graining.

4 Spectral Properties and Kramers–Moyal Estimators

The second part of this work is a numerical study on the effect of employing Kramers–Moyal type approximations when estimating spectral properties of projected generators. As explained in section 2.5, KM estimators can be incorporated into the gEDMD method by means of Eq. (22) if the full system parameters $b$ and $a$ are not known. Again, we stress that using gEDMD to calculate eigenvalues in this way possesses a dual meaning: it is an approximation to the eigenvalues of $\mathcal{L}^{\xi}$, but also an approximation to the eigenvalues of a non-reversible dynamics with parameters given by (21). We find that for projections onto a good set of reaction coordinates $\xi$, i.e., coordinates which are known to capture the slow dynamics of a system well, it is possible to use a surprisingly large offset parameter $s$ in Eq. (22).

4.1 Methods

In all of the following examples, we generate long realizations of the system under investigation, by employing the Euler-Maruyama method with discrete integration time step $\Delta_{t}$, for $M$ steps, such that the total simulation time equals $M\Delta_{t}$. The only exception is the molecular example in section 5.3, where a molecular dynamics code was used to generate the data, see the references given there. For the application of gEDMD on one-dimensional domains, we either use Gaussian basis functions $\phi_{i}$, or periodic Gaussians $\phi^{p}_{i}$ if the domain is periodic:

[TABLE]

with bandwidth $\rho$ and centers $z_{i}\in\hat{\Omega}$. On two-dimensional domains, we use products of the univariate functions defined above, centered on a regular grid.

With these basis functions, the Galerkin matrices for gEDMD are calculated in different ways. As a reference, we use estimators (19), which require knowledge of the full system parameters. In addition, we use the estimator (22) with a series of offsets $s>0$. We then solve the generalized eigenvalue problem (17) for each of these cases. Eigenvalue estimates thus obtained are denoted by $\hat{\omega}_{i}^{0}$ and $\hat{\omega}_{i}^{s}$, respectively. We keep track of the relative error

[TABLE]

and also monitor the reciprocals $\hat{t}_{i}^{0}=\left(\hat{\omega}_{i}^{0}\right)^{-1},\,\hat{t}_{i}^{s}=\left(\hat{\omega}_{i}^{s}\right)^{-1}$, which serve as estimates of the implied timescales (12). Additionally, we also extract estimates of the $K$ dominant eigenfunctions from each of these eigenvalue problems, and apply the PCCA method [42] to determine metastable decompositions of the domain based on these eigenfunctions. PCCA returns $K$ membership functions $\chi_{j}(z),j=1,\ldots,K$, such that $\sum_{j=1}^{K}\chi_{j}(z)=1$ at all points $z$, and $\chi_{j}(z)$ indicates the degree of membership of each point $z$ to metastable state $j$.

4.2 Lemon Slice Potential

We consider overdamped Langevin dynamics Eq. (2) in the “lemon slice” potential

[TABLE]

where $r,\varphi$ are two-dimensional polar coordinates, at inverse temperature $\beta=1.0$ and friction $\gamma=1.0$. A contour of the potential is shown in Fig. 2. Note that the second and the last term in (27) impose an infinite barrier along the negative $x$-axis and at the origin. The slow dynamics of this system correspond to transitions between the four main minima of the potential $V$, we therefore find three dominant timescales $t_{1}\approx 2.6,\,t_{2}\approx 0.95,\,t_{3}\approx 0.75$, see [28] for a previous analysis of the same example. Hence, we can select the polar angle $\varphi$ as a suitable reaction coordinate $\xi(x,y)=\varphi(x,y)$. In this case, the effective drift $b^{\xi}$ and diffusion $a^{\xi}$ can even be calculated analytically, see appendix A, they are indicated by the black lines in Figs. 3 A and B. The data set we use comprises $M=5\cdot 10^{6}$ data points at integration time step $\Delta_{t}=10^{-3}$.

Since applying gEDMD with a positive offset in (22) corresponds to approximating the generator of an SDE with coefficients (21), and we also have analytical expressions for the exact effective parameters $b^{\xi}$ and $a^{\xi}$, we first provide a comparison between these analytical parameters and histogrammed estimates of Eq. (21). We find in Figs. 3 A and B, that the exact drift and diffusion are recovered well for small $s$, as expected, while very different results are obtained as $s$ increases. In particular, the effective diffusion is no longer constant as a function of $\varphi$ for large $s$.

Next, we compare the dominant spectra of the generators corresponding to these different dynamics. To this end, we employ gEDMD with fifteen Gaussian basis functions (25), centered at equal distance between $z_{i}=-2.8$ and $z_{i}=2.8$, each of width $\rho=0.1$, as described in section 4.1. In Fig. 3 C, we show the first three timescales $\hat{t}_{i}^{s}$ compared to their reference values, and also the relative errors $E_{i}^{s}$, see (26). As $s$ increases by about two orders of magnitude, the relative errors remain within a ten to twenty percent margin around the references, which is generally acceptable from a practical point of view. As shown in Fig. 3 D, the membership functions $\chi_{j},j=1,\ldots,4$ generated by the PCCA method for $s=0.1$, barely differ from the ones computed using full system parameters. We therefore conclude that the dominant spectrum is approximately retained by all models up to $s=0.1$. In other words, the structure of the spectrum is unchanged as long as the offset is clearly smaller than the slow timescales, even though the KM estimates for drift and diffusion are very different from the analytical coefficients $b^{\xi}$ and $a^{\xi}$.

4.3 Prototypical Molecular Potential

In order to confirm the observations made in the previous section, we study a more complex example. The system is designed to mimic a small molecule consisting of five atoms. The three-dimensional Euclidean coordinates of all five atoms thus define the system’s fifteen-dimensional state space. We imagine these five atoms to be linked by bonds like a chain. The dynamical model is again the overdamped Langevin dynamics (2), with a potential energy comprised of typical molecular interactions, namely harmonic bond, bond angle, and dihedral potentials. The precise parameters can be found in appendix B. Note, however, that the system is only qualitatively similar to a small molecule, since parameter values and the units of time and energy do not correspond to physical values. We also do not include any solvent molecules or velocities. The data set comprises $M=2\cdot 10^{6}$ data points at integration time step $\Delta_{t}=5\cdot 10^{-3}$.

The system parameters are tuned in such a way that the two dihedral angles $\phi_{1},\,\phi_{2}$, spanned by those five atoms, capture the slow dynamics of the system, corresponding to transitions between six symmetrically arranged minima of the effective free energy in the $\phi_{1}$-$\phi_{2}$-plane, see Fig. 4 A. We run the reference gEDMD analysis using 100 two-dimensional Gaussian basis functions centered on a ten-by-ten grid, with $\rho=0.05$. Five dominant implied timescales between $t_{1}\approx 25$ and $t_{5}\approx 10$ (black lines in Fig. 4 C) are determined, and the six energy minima are recovered as metastable states by a PCCA analysis, as shown in Fig. 4 B.

We repeat the same experiment as for the previous example: using the same basis set of 100 Gaussians, we apply gEDMD with KM estimators (22), for a series of offsets between $s=5\cdot 10^{-3}$ and $s=1.0$. The corresponding timescale estimates $\hat{t}^{s}_{i}$ and the relative errors $E^{s}_{i}$ are also shown in Fig. 4 C. The relative errors are generally larger than for the previous example, but still within an acceptable range of ten to forty percent. Applying the PCCA analysis to the model extracted at $s=1.0$, we are still able to reproduce the correct metastable decomposition of reaction coordinate space. In summary, we can still conclude that the dominant spectrum is retained as long as the offset is at least an order of magnitude smaller than the fastest interesting timescale.

4.4 Summary of Observations

The previous two examples have indicated that, at least for reversible systems, and for a good selection of reaction coordinates, it is possible to use the KM estimators within the framework of gEDMD without disrupting the overall structure of the dominant spectrum. We state this observation as a conjecture, which will be investigated and made quantitative in future work:

Conjecture 4.1**.**

Let $\xi$ be a good set of reaction coordinates for a reversible process $X_{t}$, which capture the slow part of the dynamics well (e.g. small projection errors $\delta_{i}$ as in Proposition 6.3). Let $s$ be an offset such that $0<s\ll t_{K}$. Then it is possible to recover the dominant spectrum by applying gEDMD with KM estimators, as in (22), at least in a qualitative sense.

Note that Conjecture 4.1 is equivalent to saying that the effective dynamics with parameters (21) will retain the dominant spectrum well in the setting described above. Also, the quality of a gEDMD approximation depends on the basis set. Conjecture 4.1 should be understood in the sense of using gEDMD with a powerful basis set, such that the basis set error does not play a major role.

5 Underdamped Langevin Dynamics

The third topic of this paper is to study the implications of the first two results for systems driven by underdamped Langevin dynamics. Let us recall that the theoretical analysis presented in sections 3 and 6 hinges on the reversible setting. Also, the numerical results on KM estimators shown in section 4 were using the reversible overdamped Langevin process (2) (OL process from now on). However, a popular dynamical model, especially in molecular physics, are underdamped Langevin dynamics (UL process) in position and momentum space $(q,\,p)$, where $q\in\Omega$ is the position of the system and $p\in\mathbb{R}^{d}$ is its momentum. The equations of motion are

[TABLE]

Here $V,\gamma,\beta$ have the same physical meaning of energy, friction, and inverse temperature, as in (2) above. The invariant density of the process (28 - 29) factors as

[TABLE]

The underdamped process is not reversible, so the error theory developed in this paper does not apply. However, it is well known that for large enough friction, if we observe the position coordinate $q_{t}$ only at every $s$-th step, for $s$ large, this time-rescaled process behaves like an OL process. This phenomenon is called the overdamped limit, and gives rise to the following idea: if we select a reaction coordinate $\xi=\xi(q)$, which does not depend on momentum, then $\xi(q_{ks}),\,k=1,2,\ldots$, observed at large offset $s$, behaves like a projected OL process. If we employ KM estimators or related approaches at the same offset within the context of gEDMD, we effectively compute models for the OL dynamics. Our error theory does apply to the reversible OL process, and if our Conjecture 4.1 is correct, we still get the dominant spectrum of these OL dynamics right by using KM formulae at a large offset. This idea will be illustrated by numerical examples in the following subsections.

5.1 Projection and Re-scaling of the Underdamped Process

We start by discussing the projection of the underdamped Langevin process by a map $\xi=\xi(q)$ depending only on the position coordinates. Oftentimes, one is not interested in quantities that explicitly depend on the momenta, which renders this a realistic setting (see [35] for an approach to model reduction which includes the momenta). Unfortunately, the coefficients of the effective dynamics (13) are identically zero in this case, see also [43, 34]. To see this, it is readily checked using the definition of the generator that $\mathcal{L}\xi=p\cdot\nabla_{q}\xi$ and $\nabla\xi^{T}a\nabla\xi=0$. Using the definition of the parameters (15), the factorization (30) of the invariant measure, and the fact that $\mu_{P}(p)$ is the density of a multivariate normal distribution with mean zero, we find that

[TABLE]

How can one define a suitable effective dynamics in this case? As already mentioned above, if the friction $\gamma$ is sufficiently large, the positional component $q_{t}$ behaves like a reversible overdamped dynamics (2) on long time scales. More precisely, if we observe the positions $q_{ks},\,k=1,2,\ldots$ for an offset $s\gg\frac{1}{\gamma}$, then the statistics of this process will be approximately the same as those of an OL process $X_{t}$ in position space, observed at the same offset. A particular pair of statistics to observe is given by the KM formulae (21). As a consequence, if the reaction coordinate $\xi$ captures the slow dynamics of the OL process $X_{t}$ well, and if in addition our Conjecture 4.1 is correct, then we can use the KM estimators on the underdamped data to build a good model (or a suitable effective dynamics) for the overdamped process $X_{t}$.

The following corollary provides a formal derivation of this argument, again in a qualitative sense, thus connecting the results of sections 3 and 4.

Corollary 5.1**.**

Let $X_{t}\in\Omega\subset\mathbb{R}^{d}$ denote the OL process (2). Moreover, assume $\xi$ is good reaction coordinate for $X_{t}$ (e.g. small projection errors $\delta_{i}$ in Proposition 6.3). Let $(q_{t},p_{t})$ denote the UL process on $\Omega\times\mathbb{R}^{d}$ with the same parameters as $X_{t}$.

(i) For $\gamma$ sufficiently large, the statistics of the positional component $q_{t}$ of the underdamped process, and of the overdamped process $X_{t}$, approximately agree at offsets much larger than $\frac{1}{\gamma}$.

(ii) If in addition Conjecture 4.1 is true, then application of the Kramers–Moyal estimator (22) at offsets $t_{K}>s>\frac{1}{\gamma}$ to $q_{t}$, allows to recover the dominant spectrum by means of gEDMD.

Proof.

(i) Setting $\epsilon=\frac{1}{\gamma}$, and re-scaling the UL process (28 - 29) by $(q^{c}_{t},p^{c}_{t})=(q_{c\gamma t},p_{c\gamma t})=(q_{ct/\epsilon},p_{ct/\epsilon})$, the re-scaled equations of motion are

[TABLE]

For sufficiently large $\gamma$ and $c\geq 1$, Theorem 18.1 of [6] implies that the law of $q_{t}^{c}$ is close to that of

[TABLE]

We note that $Q^{c}_{t}$ is a time re-scaling of $X_{t}$ via $Q^{c}_{t}=X_{c\gamma t}$, hence we have for functions $f=f(q)\in L^{2}_{\mu_{Q}}(\Omega)$:

[TABLE]

(ii) In particular, the approximate equality (31) applies directly to estimator (22) for $s>\frac{1}{\gamma}$. If Conjecture 4.1 is true, the resulting gEDMD models will qualitatively retain the leading timescales of the process $X_{t}$. ∎

5.2 Langevin Toy Model

To illustrate the ideas outlined in the previous section, we first consider another two-dimensional toy potential

[TABLE]

shown in Fig. 5. For $\gamma=10,\,\beta=0.4$, we generate data of both the UL and OL processes for this potential, each data set comprising $M=10^{7}$ points at integration time step $\Delta_{t}=10^{-3}$. For both dynamical models, the slowest transition in this energy landscape is the crossing of the barrier around $x=0$, thus $\xi(x,y)=x$ is a suitable reaction coordinate in both instances. The associated implied timescale in the OL model is $t_{1}\approx 7.3$. The effective drift in the overdamped case simply corresponds to the $x$-derivative of the $x$-dependent part the potential, while the effective diffusion remains constant, see the black dashed lines in Figs. 6 A and B.

We first illustrate the relationship between the overdamped limit and estimates of drift and diffusion by means of KM formulae. In Figures 6 A and B, we show histogrammed estimates of the KM expressions (21), using both the OL and UL simulation data. We see that for a small offset $s=10^{-3}$, UL estimates are almost zero, while the OL data lead to estimates close to the analytical values. Both findings are as expected. For a much larger offset $s=0.5$, however, both estimates are significantly different from what we find for small offsets, and most importantly, the overdamped and underdamped estimates agree well. This in line with our qualitative argument in Corollary 5.1 (i).

We apply gEDMD to both data sets, using the KM estimators as in (22), for a series of offsets between $s=10^{-3}$ and $s=1.0$. The basis set is comprised of fifteen Gaussian basis functions (25), centered uniformly between $x=-2.5$ and $x=2.5$, each of width $\rho=0.1$. From all of these models, we extract the slowest implied timescale $\hat{t}_{1}^{s}$ and the relative errors $E^{s}_{1}$, and present the results in Fig. 6 C. As expected, OL estimates are highly accurate at small offsets. For large offsets, these estimates remain within the same error margin that we have seen before. Estimates based on the UL data, however, are far off for small offsets (as the coefficients of the corresponding dynamics are essentially zero), but once the offset exceeds the critical value $s>\frac{1}{\gamma}=0.1$, they are about as accurate as those based on OL data. We also verify in Fig. 6 D that for both data sets, the two metastable states of the effective dynamics can be correctly identified by applying PCCA to the eigenvectors of the gEDMD model at a large offset $s=0.5$.

This example shows that upon increasing the offset it is possible to find a sweet spot where $s$ is larger than the critical relaxation time $\frac{1}{\gamma}$, but smaller than the slowest timescale $t_{1}$, such that meaningful effective dynamics for the UL process along $x$ can be defined in this regime.

5.3 Alanine Dipeptide

Our final numerical example is a more complex data set, namely molecular dynamics simulations of the alanine dipeptide. The data we use is the same as in Refs. [44, 45], consisting of $M=10^{6}$ points at $1\mathrm{ps}$ time spacing. The dynamical model is the UL process as in section 5.2, using the AMBER 99 molecular force field [46] as potential $V$. Friction is set to $\gamma=0.1\mathrm{ps}^{-1}$, and $\beta$ is derived from the temperature $T=300\mathrm{K}$ via $\beta^{-1}=k_{B}T$, using the Boltzmann constant $k_{B}$. As is well-known from numerous previous studies, the slow dynamics of alanine dipeptide can be represented well in the space of backbone dihedral angles $\phi,\,\psi$. Figure 7 A shows the effective free energy landscape in the space of these reaction coordinates. Three major minima can be identified, which correspond to three metastable states of the full dynamics. The corresponding transition timescales $t_{1}\approx 1\,\mathrm{ns}$ and $t_{2}\approx 0.1\,\mathrm{ns}$ are indicated by the black lines in Fig. 7 B.

For gEDMD, we employ a similar basis set as in section 4.3, comprised of 100 periodic Gaussians, centered on a ten-by-ten grid between $-2.5$ and $2.5$, with bandwidth $\rho=0.05$. We compute gEDMD models for a range of offsets between $s=1\mathrm{ps}$ and $s=50\mathrm{ps}$. The resulting first two implied timescales $\hat{t}_{1}^{s},\,\hat{t}_{2}^{s}$ and the corresponding relative errors $E_{i}^{s}$ are shown in Fig. 7 B. These results confirm the findings of the previous examples, as both leading timescales are roughly reproduced for all offsets considered. The relative error $E_{2}^{s}$ for the second timescale is generally larger than what we observed in the previous examples, but it is still acceptable in the context of this example. An interesting observation is that for small offsets, we are able to fully recover both slow processes by a PCCA analysis. The corresponding state decomposition is shown in Fig. 7 C. For a larger offset which is comparable to the timescale $t_{2}$, the faster transition within the left part of the plane appears to be blurred out. However, applying PCCA with only two metastable states still recovers the slowest process, as indicated by the decomposition in Fig. 7 D.

6 Precise statements on spectral properties and their proofs

In this section, we provide detailed proofs of the spectral approximation results outlined in section 3.

6.1 Form Domain

We consider an open domain $\Omega\subset\mathbb{R}^{d}$, and make the assumption of uniform ellipticity (3). The generator $\mathcal{L}$ in (7) can be defined initially on the set of smooth functions or smooth and compactly supported functions. The form domain $\mathbb{V}_{\mu}$ can then be obtained as the closure of this initial domain with respect to the Dirichlet norm

[TABLE]

We note that on domains with a boundary, the choice of initial domain has an impact on the boundary conditions. In addition, we can restrict all function spaces to the orthogonal complement of the constant one function without explicitly changing the notation.

The assumption of uniform ellipticity implies that $\mathbb{V}_{\mu}$ is also a Hilbert space if equipped with the energy norm

[TABLE]

In many cases of practical interest, the energy and Dirichlet norms are equivalent [36]. We state this as an assumption and generally use the energy norm on $\mathbb{V}_{\mu}$ in what follows:

Assumption 1: The Dirichlet norm (33) and energy norm (34) are equivalent. .

6.2 Solution Operator and Discrete Spectrum

The solution operator $\mathcal{T}:\,\mathbb{V}_{\mu}\rightarrow\mathbb{V}_{\mu}$ associated to $\mathcal{Q}$ is defined by

[TABLE]

Assumption 2: The solution operator is compact on $\mathbb{V}_{\mu}$ with norm $\|\cdot\|_{\mathcal{Q}}$.

As a consequence of this assumption, the generator possesses a complete set of eigenfunctions $\psi_{i}\in\mathbb{V}_{\mu}$ with eigenvalues $\kappa_{i},\,i=1,2,\ldots$, which are given as reciprocals of the eigenvalues of $\mathcal{T}$. There is a number of well-known settings where Assumption 2 can be shown to hold. These include bounded Lipschitz domains with Dirichlet, Neumann or periodic boundary conditions, as well as overdamped Langevin dynamics with a potential satisfying suitable growth conditions on the potential $V$ [36].

6.3 Coarse-Grained Generator and its Spectrum

Analogously to section 6.1, the projected generator $\mathcal{L}^{\xi}$ can be defined initially on the set of smooth or smooth and compactly supported functions in $L^{2}_{\nu}$, which is an infinite-dimensional subspace of $L^{2}_{\mu}$ [17]. Using the effective quadratic form $\mathcal{Q}^{\xi}$ (16), the effective form domain $\mathbb{V}_{\nu}$ can be defined again by completion of the initial domain with respect to the corresponding Dirichlet norm $\|\phi\|^{2}_{1,\xi}=\|\phi\|^{2}_{L^{2}_{\nu}}+\mathcal{Q}^{\xi}(\phi,\phi)$. Due to the relations [17]

[TABLE]

Assumption 1 carries over to the effective Dirichlet norm and the effective energy norm. In order to proceed from this point, some care needs to be taken with regard to the coarse graining map:

Assumption 3: The coarse graining map $\xi$ is such that the effective form domain $\mathbb{V}_{\nu}$ is a subspace of $\mathbb{V}_{\mu}$.

Remark 6.1**.**

Assumption 3 will hold in many cases of practical interest, which include the projection of a periodic domain onto a lower-dimensional periodic domain by a function $\xi$ which respects the periodic boundary conditions. As a negative example, however, consider an SDE on a rectangle in two dimensional space, with absorbing boundary conditions. Choose the coarse graining function as the projection onto the first coordinate axis. The form domain and the effective form domain will be given as first order Sobolev spaces with zero boundary conditions, but the effective form domain is not contained in the full form domain as its elements do not vanish on parts of the full boundary which are parallel to the first coordinate axis.

If Assumption 3 holds, it makes sense to define the $Q$-orthogonal projection from $\mathbb{V}_{\mu}$ onto $\mathbb{V}_{\nu}$, which we denote by $\mathcal{P}_{\mathcal{Q}}$. As a first main result, we show that assumptions 1, 2 and 3 are sufficient to ensure that the spectrum of $\mathcal{L}^{\xi}$ is also discrete:

Proposition 6.2**.**

If assumptions 1, 2, and 3 hold, the spectrum of $\mathcal{L}^{\xi}$ is discrete.

Proof.

By assumption 1, the effective solution operator $\mathcal{T}^{\xi}$ is uniquely defined on $\mathbb{V}_{\nu}$ by

[TABLE]

It is readily checked that $\mathcal{T}^{\xi}=\mathcal{P}_{\mathcal{Q}}\mathcal{T}\mathcal{P}_{Q}$. Since $\mathcal{P}_{\mathcal{Q}}$ is bounded and $\mathcal{T}$ is compact by assumption 2, so is $\mathcal{T}^{\xi}$. Hence its spectrum (and correspondingly that of $\mathcal{L}^{\xi}$) is discrete. ∎

6.4 Approximation Result

Next, we derive a bound on the eigenvalue error of $\mathcal{L}^{\xi}$ in terms of the energy-norm approximation error of the dominant eigenfunctions of $\mathcal{L}$. The idea is to apply classical Galerkin error estimates to a sequence of finite-dimensional subspaces $\mathbb{W}^{h}$ in $\mathbb{V}_{\nu}$, and exploit that these provide approximations to the (nonzero) eigenvalues of both $\mathcal{L}$ and $\mathcal{L}^{\xi}$.

Proposition 6.3**.**

Denote the projection error of eigenfunction $\psi_{i}$ with respect to the energy-norm by $\delta_{i}=\|\mathcal{P}^{\perp}_{\mathcal{Q}}\psi_{i}\|_{\mathcal{Q}}$. The relative error between the $i$-th eigenvalues of $\mathcal{L}^{\xi}$ and $\mathcal{L}$ is bounded by

[TABLE]

Proof.

We consider a sequence $\mathbb{W}^{h}$ of finite-dimensional subspaces in the reduced space $\mathbb{V}_{\nu}$. For every $\mathbb{W}^{h}$, the $\mathcal{Q}$-orthogonal projection from $\mathbb{V}_{\mu}$ onto $\mathbb{W}^{h}$ is labeled by $\mathcal{P}^{h}_{\mathcal{Q}}$, and the corresponding $\mathcal{Q}^{\xi}$-orthogonal projection in $\mathbb{V}_{\nu}$ is called $\mathcal{P}^{h}_{\mathcal{Q},\xi}$. We assume that all $\mathbb{W}^{h}$ contain the projections $\mathcal{P}_{\mathcal{Q}}\psi_{j}$ for $j=1\,\ldots,i$, and satisfy the following approximability condition in $\mathbb{V}_{\nu}$ (which holds in any separable Hilbert space):

[TABLE]

Note that the spaces $\mathbb{W}^{h}$ are just of auxiliary nature; their sole purpose being to reconcile the theory of [47]—that uses finite-dimensional approximation spaces—with our infinite-dimensional $\mathbb{V}_{\nu}$. There is no need to specify them in detail, and $h$ only serves as a formal parameter here. Using (18), and by applying [47, Theorem 3.2] to the approximation of $\mathcal{L}$, we find that for any of the approximation spaces $\mathbb{W}^{h}$:

[TABLE]

where $\hat{\omega}_{i}$ are the Ritz values associated to $\mathbb{W}^{h}$, see section 2.4. The last equality holds because $\mathcal{P}_{Q}\psi_{i}\in\mathbb{W}^{h}$ by assumption. It remains to study the pre-factor in the limit of $h\rightarrow 0$. We conclude from (38) that $\hat{\omega}_{i}\rightarrow\omega_{i}$ as $h\rightarrow 0$, which already yields the first term in the pre-factor. Regarding the second term, we first observe that $\mathcal{P}^{h}_{\mathcal{Q}}=\mathcal{P}^{h}_{\mathcal{Q},\xi}\mathcal{P}_{\mathcal{Q}}$, since for any $\psi\in\mathbb{V}_{\mu}$ and $\phi\in\mathbb{W}^{h}$:

[TABLE]

The third equality is due to (36), while the last equality is due to the definition of $\mathcal{P}^{h}_{\mathcal{Q},\xi}$. With this result, and using pointwise convergence (38) of $\mathcal{P}^{h}_{\mathcal{Q},\xi}$ to the identity in $\mathbb{V}_{\nu}$, we conclude for any $\psi\in\mathbb{V}_{\mu}$:

[TABLE]

Combing this pointwise convergence of $\mathcal{I}-\mathcal{P}^{h}_{\mathcal{Q}}$ towards $\mathcal{P}_{Q}^{\perp}$ in $\mathbb{V}_{\mu}$, and compactness of the solution operator, we have

[TABLE]

with respect to the operator norm. This establishes the second term in the pre-factor and hence the proposition. ∎

To make the bound from the previous result more accessible, we will bound the term involving the operator $\mathcal{T}$.

Corollary 6.4**.**

In the setting of Proposition 6.3, the following error bound holds:

[TABLE]

Proof.

As already mentioned, it follows from (35) that

[TABLE]

Since the $\psi_{i}$ are $\mathcal{Q}$-orthogonal, it follows that $\|\mathcal{T}\|_{\mathcal{Q}}=\frac{1}{\kappa_{1}}$. Further, since $\mathcal{P}_{Q}$ is a $\mathcal{Q}$-orthogonal projection, so is $\mathcal{I}-\mathcal{P}_{Q}$, and thus both have $\mathcal{Q}$-norm at most one. It follows that $\|(\mathcal{I}-\mathcal{P}_{Q})\mathcal{T}\|_{\mathcal{Q}}\leq\frac{1}{\kappa_{1}}$, and hence the claim. ∎

Remark 6.5**.**

Due to the equivalence of the norms $\|\cdot\|_{\mathcal{Q}}$ and $\|\cdot\|_{1}$, the bound in Proposition 6.3 also applies to the latter norm. Let $\mathcal{P}_{1}$ denote the orthogonal projection onto $\mathbb{V}_{\nu}$ with respect to the Dirichlet norm. Then

[TABLE]

6.5 Comments

The bound from Corollary 6.4 allows some interpretation. Note that in (39) the ratio of the eigenvalues $\kappa_{1},\kappa_{i}$ of $\mathcal{L}$ and the relative difference of eigenvalues $\kappa_{i}$ and $\omega_{j}$ of $\mathcal{L}$ and $\mathcal{L}^{\xi}$ (respectively) play a role. Let us fix $i$ and assume that the first $i$ timescales are comparable, meaning that $\kappa_{i}/\kappa_{1}$ is a moderate number. Let us further assume that the squared projection errors $\delta_{\ell}$ are all much smaller than the relative eigenvalue gaps, i.e.,

[TABLE]

Then, inductively from $j=1$ to $j=i$ it follows that $(\omega_{j}-\kappa_{j})/\omega_{j}\approx\mathcal{O}(\delta_{j}^{2})$, since we find that $\omega_{j}\approx\kappa_{j}$, thus the relative differences between $\kappa_{\ell}$ and $\omega_{j}$ are large for $j\neq\ell$, which makes the term in the square brackets on the right-hand side of (39) moderate. Hence, in this situation the relative error of the $i$-th timescale’s approximation is effectively governed by the projection error $\delta_{i}^{2}$.

We also note that the bound assumed in [18, Theorem 2] implies our bound up to another multiplicative constant:

Lemma 6.6**.**

Let us consider the diffusion (1), satisfying (3), on a bounded domain with smooth boundary and reflecting boundary conditions. If $\|\mathcal{L}\mathcal{P}^{\perp}\psi_{i}\|_{L_{\mu}^{2}}\leq\delta_{i}$, then there is a $C>0$ such that the assumptions of Proposition 6.3 are satisfied with

[TABLE]

Proof.

We use the weighted Sobolev spaces $H^{k}_{\mu}$ [48]. Since the spatial domain is bounded and $\mu$ is smooth, it is a weight in Muckenhoupt class, cf. [48, Eq. (1.2)]. The regularity conditions [48, Eqs. (2.2)–(2.4)] are satisfied by assumption. Now, the result follows directly from the weighted Agmon–Douglis–Nierenberg estimate [48, Theorem 2.4] giving $\|u\|_{H^{2}_{\mu}}\leq C\|\mathcal{L}u\|_{L^{2}_{\mu}}$ for some $C>0$ independent of $u$. Using assumption 1 and (3), we have

[TABLE]

concluding the proof. ∎

7 Conclusions

We have investigated the approximation of high-dimensional diffusion processes by effective dynamics defined on the lower-dimensional space of reduced variables. For reversible diffusions, a new relative error bound for the approximation of eigenvalues of the infinitesimal generator by the eigenvalues of a reduced generator was proved. Our bound shows that a small projection error of the corresponding generator eigenfunctions, measured by the energy norm, is sufficient for small eigenvalue errors.

In addition, we have presented numerical examples regarding the data-driven estimation of the eigenvalues of projected generators by means of the gEDMD method. If the full system parameters are unknown, they need to be approximated, for example using Kramers–Moyal formulae. We have presented numerical examples that for reversible systems, and good reaction coordinates, the resulting spectral estimates seem to remain stable across a long range of time windows in the KM-estimators.

Finally, we have suggested a strategy to define meaningful effective equations for underdamped Langevin dynamics on a subset of its position space. Exploiting the overdamped limit, and using KM estimators at large time windows, we can effectively model a projected overdamped process on the same domain. Numerical examples have confirmed the feasibility of this approach.

Future work will focus on providing a theoretical foundation for the observations stated in Conjecture 4.1. Also, the relation between the positional coordinate of an (underdamped) Langevin process on long timescales, and the corresponding overdamped Langevin equation (cf. Corollary 5.1 (i)), needs to be analyzed in more detail.

Funding and Acknowledgments

This research was funded by Rice University Academy of Fellows (F.N.), by Deutsche Forschungsgemeinschaft (DFG) through CRC 1114 ”Scaling Cascades in Complex Systems”, Project Number 235221301, project A01 (P.K.), by the Einstein Foundation Berlin (C.C.), by the National Science Foundation (Nos. CHE-1265929, CHE- 1738990, and PHY-1427654) and the Welch Foundation (No. C-1570) (F.N., L.B., C.C.).

We would like to thank Ralf Banisch, Stefan Klus, Thomas Swinburne, Tony Lelièvre and Frank Noé. We are also indebted to the Institute for Pure and Applied Mathematics (IPAM) for hosting the long program ”Complex High-Dimensional Energy Landscapes”.

Appendix A Lemon Slice Potential: Effective Drift and Diffusion Coefficients

Here, we calculate the effective drift and diffusion if the lemon slice potential Eq. (27) is projected onto the polar angle $\xi(x,y)=\varphi(x,y)$. We start by expressing the generator $\mathcal{L}$ in polar coordinates. The Laplacian operator in polar coordinates is

[TABLE]

Moreover, it follows from the chain rule and the definition of $\xi$ that for any function $f$,

[TABLE]

The generator in polar coordinates then becomes

[TABLE]

Applying the generator to the reaction coordinate $\xi$, only one of the terms above is non-zero, resulting in

[TABLE]

In order to evaluate the first term in Eq. (15), we note that the Jacobian determinant of $\xi$ is $J(x,y)=\frac{1}{r^{2}}$, and that the stationary distribution factors, canceling the $\varphi$-dependent term. We are left with

[TABLE]

with the definitions $C_{1}:=\int_{0}^{\infty}\frac{1}{r}\exp(-10(r-1)^{2}-\frac{1}{r})\,\mathrm{d}r$ and $C_{2}:=\int_{0}^{\infty}r\exp(-10(r-1)^{2}-\frac{1}{r})\,\mathrm{d}r$. For the diffusion, we obtain

[TABLE]

Inserting this into the second term in (15), and using the same arguments as before, we end up with

[TABLE]

Appendix B Parameters of Prototypical Molecular Example

The system we considered in section 4.3 is defined as follows: denote the three-dimensional position vectors of atoms one through five by $\mathbf{r}_{i}$, $1\leq i\leq 5$. We denote the pairwise distances between atoms $i,\,j$ by $d_{i,j}$, the bond angle formed by atoms $i,j,k$ by $\theta_{i,j,k}$, and the dihedral angle formed by atoms $i,j,k,l$ by $\phi_{i,j,k,l}$. The potential energy is then defined as a superposition of harmonic bond terms, harmonic bond angle terms, and dihedral angle terms:

[TABLE]

The constants appearing above are set to

[TABLE]

The overdamped Langevin dynamics (2) with potential (41) are simulated at inverse temperature $\beta=15$, for $10^{7}$ discrete steps at integration time step $\Delta_{t}=10^{-3}$. For the purpose of our analysis, we retain only every fifth step of this trajectory.