Dimensionality Reduction of Complex Metastable Systems via Kernel Embeddings of Transition Manifolds

TL;DR

This paper introduces a kernel-based machine learning method to identify low-dimensional structures in complex stochastic systems, improving robustness and efficiency over previous approaches.

Contribution

The authors develop a kernel embedding technique for transition manifolds, enhancing the mathematical framework for reaction coordinate computation in high-dimensional systems.

Findings

Kernel embeddings preserve manifold structure under mild conditions.

The new algorithm outperforms previous parametrization methods in robustness.

Distortion bounds are derived for the kernel embedding process.

Abstract

We present a novel kernel-based machine learning algorithm for identifying the low-dimensional geometry of the effective dynamics of high-dimensional multiscale stochastic systems. Recently, the authors developed a mathematical framework for the computation of optimal reaction coordinates of such systems that is based on learning a parametrization of a low-dimensional transition manifold in a certain function space. In this article, we enhance this approach by embedding and learning this transition manifold in a reproducing kernel Hilbert space, exploiting the favorable properties of kernel embeddings. Under mild assumptions on the kernel, the manifold structure is shown to be preserved under the embedding, and distortion bounds can be derived. This leads to a more robust and more efficient algorithm compared to previous parametrization approaches.

Equations

$$\tilde{\mathbb{M}}:=\{p_x^{\tau }\mid x\in \mathbb{X}\}\subset L^1,$$

$$d{X}_t=-\nabla V(X_t)dt+\sqrt{2/\beta }d{W}_t,$$

$$\underset{f\in \mathbb{M}}{\min }\Vert f-p_x^{\tau }{\Vert }_{L_{1/\rho }^2}\le \epsilon \text{for all }x\in \mathbb{X}.$$

$$\Vert f{\Vert }_{L_{1/\rho }^2}={\left({\int }_{\mathbb{X}}f(x{)}^2\frac{1}{\rho (x)}dx\right)}^{1/2}.$$

$$d{X}_t$$

$$d{Y}_t$$

$$\big{(}\mathcal{P}^{t}u\big{)}(y)=\int_{\mathbb{X}}u(x)\hskip 1.00006ptp^{t}_{x}(y)\hskip 1.00006ptdx.$$

$$\underset{i\to \infty }{\lim }|{\theta }_i^t|=0\text{and}\underset{t\to \infty }{\lim }|{\theta }_i^t|=0.$$

$${\mathcal{P}}^{t}u=\sum _{i=0}^{\infty }{\theta }_i^t{\alpha }_i{\psi }_i\approx \sum _{i=0}^d{\theta }_i^t{\alpha }_i{\psi }_i$$

$$\big{(}\Pi_{\xi}(u)\big{)}(y)=\mathbb{E}_{\rho}\big{[}u(\mathbf{x})~{}\big{|}~{}\xi(\mathbf{x})=\xi(y)\big{]},$$

$$\mathcal{P}^{t}_{\xi}u=\Pi_{\xi}\big{(}\mathcal{P}^{t}(\Pi_{\xi}u)\big{)},$$

$${\mathcal{P}}^{t}u\approx {\mathcal{P}}_{\xi }^{t}u,$$

$${\Pi }_{\xi }{\psi }_i\approx {\psi }_i,i=0,\dots ,d,$$

$$\big{\|}\psi_{i}-\tilde{\psi}_{i}\circ\xi\big{\|}_{\infty}\approx 0.$$

$$Q(f)=\underset{g\in \mathbb{M}}{\mathrm{argmin}}\Vert f-g{\Vert }_{L_{1/\rho }^2}.$$

$$\xi(x):=\big{(}\gamma\circ{Q}\big{)}(p^{\tau}_{x})$$

$$\big{\|}\psi_{i}-\tilde{\psi}_{i}\circ\xi\big{\|}_{\infty}\leq\frac{\varepsilon}{|\theta_{i}^{\tau}|}.$$

$$\frac{\left|{\theta }_d^t-{\theta }_{d+1}^t\right|}{\left|{\theta }_{d+1}^t-{\theta }_{d+2}^t\right|}\gg 1\text{for suitably large }t>0.$$

$$C{\Vert \mathcal{F}(x)-\mathcal{F}(y)\Vert }_2^{\alpha }\ge {\Vert x-y\Vert }_{\mathbb{V}},$$

$$\mathcal{F}(f)$$

$$\eta (x)$$

$$\mathcal{E}(x):=\mathcal{F}(p_x^t)=\int \eta (x^{\prime })p_x^t(x^{\prime })d{x}^{\prime }.$$

$$\underset{v\in \mathcal{F}(\mathbb{M})}{\inf }\Vert v-\mathcal{E}(x){\Vert }_{\infty }\le \Vert \eta {\Vert }_{\infty }\epsilon \text{for all }x\in \mathbb{X}.$$

$$\underset{v\in \mathcal{F}(\mathbb{M})}{\inf }\Vert v-\underset{=\mathcal{F}(p_x^t)}{\underbrace{\mathcal{E}(x)}}{\Vert }_{\infty }$$

$$={\Vert {\int }_{\mathbb{X}}\eta (x^{\prime })\left(p_{x^{\ast }}^t(x^{\prime })-p_x^t(x^{\prime })\right)d{x}^{\prime }\Vert }_{\infty }$$

$$\le \Vert \eta {\Vert }_{\infty }\underset{\le \Vert p_{x^{\ast }}^t-p_x^t{\Vert }_{L_{1/\rho }^2}}{\underbrace{\Vert p_{x^{\ast }}^t-p_x^t{\Vert }_{L^1}}}\le \Vert \eta {\Vert }_{\infty }\epsilon ,$$

$$\mathcal{F}(p_x^t)=\int \eta (x^{\prime })p_x^t(x^{\prime })d{x}^{\prime }$$

$$\xi(x):=\tilde{\gamma}\big{(}\widehat{\mathcal{E}}(x)\big{)},$$

$$\xi (x_i):=\tilde{\gamma }(z_i),i=1,\dots ,N.$$

$$\big{\|}\mathcal{E}(x_{i})-\mathcal{E}(x_{j})\big{\|}_{2}^{2}=\big{\langle}{\mathcal{E}}(x_{i}),{\mathcal{E}}(x_{i})\big{\rangle}-2\big{\langle}{\mathcal{E}}(x_{i}),{\mathcal{E}}(x_{j})\big{\rangle}+\big{\langle}{\mathcal{E}}(x_{j}),{\mathcal{E}}(x_{j})\big{\rangle}.$$

$$\big{\langle}{\mathcal{E}}(x_{i}),{\mathcal{E}}(x_{j})\big{\rangle}=\iint\big{\langle}\eta(y_{i}),\eta(y_{j})\big{\rangle}\hskip 1.00006ptp^{t}_{x_{i}}(y_{i})\hskip 1.00006ptp^{t}_{x_{j}}(y_{j})\hskip 1.00006ptdy_{i}\hskip 1.00006ptdy_{j},$$

$$\big{\langle}\widehat{\mathcal{E}}(x_{i}),\widehat{\mathcal{E}}(x_{j})\big{\rangle}=\frac{1}{M^{2}}\sum_{l_{1},l_{2}=1}^{M}\big{\langle}\eta(y_{i}^{(l_{1})}),\eta(y_{j}^{(l_{2})})\big{\rangle}.$$

Full text

1 Introduction

Many of the dynamical processes investigated in the sciences today are characterized by the existence of phenomena on multiple, interconnected time scales that determine the long-term behavior of the process. Examples include the inherently multiscale dynamics of atmospheric vortex- and current formation which needs to be considered for effective weather prediction [24, 29], or the vast difference in time scales on which bounded atomic interactions, side-chain interactions, and the resulting formation of structural motifs occur in biomolecules [19, 12, 10]. An effective approach to analyzing these systems is often the identification of a low-dimensional observable of the system that captures the interesting behavior on the longest time scale. However, the computerized identification of such observables from simulation data poses a significant computational challenge, especially for high-dimensional systems.

Recently, the authors have developed a novel mathematical framework for identifying such essential observables for the slowest time scale of a system [6]. The method—called the transition manifold approach—was primarily motivated by molecular dynamics, where the dynamics is typically described by a thermostated Hamiltonian system or diffusive motion in molecular dynamics landscapes. In these systems, local minima of the potential energy landscape induce metastable behavior, which is the phenomenon that on long time scales, the dynamics is characterized by rare transitions between certain sets that happen roughly along interconnecting transition pathways [36, 42, 17]. The sought-after essential observables should thus resolve these transition events, and are called reaction coordinates in this context [46, 4], a notion that we will adopt here. Despite of its origins, the transition manifold approach is also applicable to other classes of reducible systems (which will also be demonstrated in this article).

At the heart of this approach is the insight that good reaction coordinates can be found by parametrizing a certain transition manifold $\mathbb{M}$ in the function space $L^{1}$. This manifold has strong connections to the aforementioned transition pathway [5], but is not equivalent. Its defining property is that, for times $\tau$ that fall between the fastest and slowest time scales, the transition density functions with relaxation time $\tau$ concentrate around $\mathbb{M}$. The algorithmic strategy to parametrize $\mathbb{M}$ can then be summarized as

Choose test points in the dynamically relevant regions of the state space. 2. 2.

Sample the transition densities for each test point by Monte Carlo simulation. 3. 3.

Embed the transition densities into a Euclidean space by a generic embedding function. 4. 4.

Parametrize the embedded transition densities with the help of established manifold learning techniques.

The result is a reaction coordinate evaluated in the test points. This reaction coordinate has been shown to be as expressive as the dominant eigenfunctions of the transfer operator associated with the system [6], which can be considered an “optimal” reaction coordinate [20, 11, 30]. One decisive advantage of our method, however, is the ability to compute the reaction coordinate locally (by choosing the test points), whereas with conventional methods, the inherently global computation of transfer operator eigenfunctions quickly becomes infeasible due to the curse of dimensionality. Kernel-based methods for the computation of eigenfunctions alleviate this problem to some extent [43, 25]. Nevertheless, the number of dominant eigenfunctions critically depends on the number of metastable states, which can be significantly larger than the natural dimension of the reaction coordinate [6].

However, the algorithm proposed originally had several shortcomings related to the choice of the embedding function. First, in order to ensure the preservation of the manifold’s topology under the embedding, the dimension of $\mathbb{M}$ had to be known in advance. Second, the particular way of choosing the embedding functions allowed no control over the distortion of $\mathbb{M}$, which may render the parametrization problem numerically ill-conditioned.

The goal of this article is to overcome the aforementioned problems by kernelizing the transition manifold embedding. That is, we present a method to implicitly embed the transition manifold into a reproducing kernel Hilbert space (RKHS) with a proper kernel. The RKHS is—depending on the kernel—a high- or even infinite-dimensional function space with the crucial property that inner products between points embedded into it can be computed by cheap kernel evaluations, without ever explicitly computing the embedding [48, 41]. In machine learning, this so-called kernel trick is often used to derive nonlinear versions of originally linear algorithms, by interpreting the RKHS-embedding of a data set as a high-dimensional, nonlinear transformation, and (implicitly) applying the linear algorithm to the transformed data. This approach has been successfully applied to methods such as principal component analysis (PCA) [40], canonical correlation analysis (CCA) [31], and time-lagged independent component analysis (TICA) [43], to name but a few.

Due to their popularity, the metric properties of the kernel embedding are well-studied [45, 21, 47, 22, 33]. In particular, for characteristic kernels, the RKHS is “large” in an appropriate sense, and geometrical information is well-preserved under the embedding. For our application, this means that distances between points on the transition manifold $\mathbb{M}$ are approximately preserved, and thus the distortion of $\mathbb{M}$ under the embedding can be bounded. This will guarantee that the final manifold learning problem is well-posed. Moreover, if the transformation induced by the kernel embedding is able to approximately linearize the transition manifold, there is hope that efficient linear manifold learning methods can be used to parametrize the embedded transition manifold.

The main contributions of this work are as follows:

We develop a kernel-based algorithm to approximate transition manifolds and compare it with the Euclidean embedding counterpart. 2. 2.

We derive measures for the distortion of the embedding and associated error bounds. 3. 3.

We illustrate the efficiency of the proposed approach using academic and molecular dynamics examples.

In Section 2, we will formalize the definition of transition manifolds and derive conditions under which systems possess such manifolds. Section 3 introduces kernels and the induced RKHSs. Furthermore, we show that the algorithm to compute transition manifolds numerically can be written purely in terms of kernel evaluations and derive measures for the distortion of the manifold caused by the embedding into the RKHS. Numerical results illustrating the benefits of the proposed kernel-based methods are presented in Section 4 and a conclusion and future work in Section 5.

2 Reaction coordinates based on transition manifolds

In what follows, let $\{X_{t}\}_{t\geq 0}$ (abbreviated as $X_{t}$) be a reversible, thus ergodic, stochastic process on a compact state space $\mathbb{X}\subset\mathbb{R}^{n}$ and $\rho$ its unique invariant density. That is, if $X_{0}\sim\rho$, then $X_{t}\sim\rho$ for all $t\geq 0$. For $x\in\mathbb{X}$ and $t\geq 0$, let $p^{t}_{x}\colon\mathbb{X}\rightarrow\mathbb{R}$ denote the transition density function of the system, i.e., $p^{t}_{x}$ describes the probability density at time $t$, after having started in point $x$ at time [math]. For fixed $t$ and $x$, we will often consider $p^{t}_{x}$ as a point in the space $L^{1}(\mathbb{X})$ (the space of absolutely integrable functions over $\mathbb{X}$ with respect to the Lebesgue measure), and later also other related function spaces. For the sake of clarity, we will from now on omit the argument of $L^{p}$ when referring to functions over $\mathbb{X}$.

2.1 Reducibility of dynamical systems

We assume the state space dimension $n$ to be large. The main objective of this work is the identification of good low-dimensional reaction coordinates (RCs) or order parameters of the system. An $r$-dimensional RC is a smooth map $\xi\colon\mathbb{X}\rightarrow\mathbb{Y}$ from the full state space $\mathbb{X}$ to a lower-dimensional space $\mathbb{Y}\subset\mathbb{R}^{r}$, $r\ll n$. Loosely speaking, we call such an RC good if on long enough time scales the projected process $\xi(X_{t})$ is approximately Markovian and the dominant spectral properties of the operator describing its density evolution of $\xi(X_{t})$ resemble those of $X_{t}$. This ensures that important long-time statistical properties such as equilibration times are preserved under projection onto the RC. We will now introduce a method to find such RCs. The explanation that the found RCs are indeed amenable to a rigorous quantitative measure of quality will be given in Section 2.2.

The method is based on a novel framework, introduced by some of the authors in [6], that ties the existence of good RCs to certain geometrical properties of the family of transition densities. This framework is called the transition manifold framework. To motivate the main idea, we first introduce the fuzzy transition manifold:

Definition 2.1**.**

Let $\tau>0$ be fixed. The set of functions

[TABLE]

is called the fuzzy transition manifold of the system.

The main idea behind the framework is now based on the following observation: If the system is in a certain way separable into a slowly equilibrating $r$-dimensional part and a quickly equilibrating $(n-r)$-dimensional part, and the lag time $\tau$ is chosen large enough so that the latter is essentially equilibrated but the former is not, then $\widetilde{\mathbb{M}}$ concentrates around an $r$-dimensional manifold in $L^{1}$. Further, any parametrization of this manifold is a good RC. Let us illustrate this with the aid of a simple example.

Example 2.2**.**

Consider the process $X_{t}$ to be described by overdamped Langevin dynamics

[TABLE]

with the energy potential $V$, the inverse temperature $\beta$, and Brownian motion $W_{t}$. The potential depicted in Figure 1 possesses two metastable states, located around the local energy minima. Any realization of this system that started in one of the the metastable sets will likely remain in that set for a long time. Moreover, a realization that started outside of the wells will with high probability quickly leave the transition region and move into one of the two metastable sets. The probability of whether the trajectory will be trapped in the left or right well depends almost exclusively on the horizontal coordinate of the starting point $x$, or, in other words, the progress of $x$ along the transition path. Thus, for times $\tau$ that allow typical trajectories to move into one of the wells, the transition densities $p^{\tau}_{x}$ also depend only on the horizontal but not the vertical coordinate of $x$. This means that the fuzzy transition manifold $\widetilde{\mathbb{M}}$ tightly concentrates around a one-dimensional manifold in $L^{1}$. Also, any parametrization of this manifold corresponds to a parametrization of the horizontal coordinate, and thus to a good RC.  $\triangle$

The concept of (fuzzy) transition manifolds can be made rigorous by the following definition:

Definition 2.3**.**

The process $X_{t}$ is called $(\varepsilon,r)$-reducible if there exists an $r$-dimensional manifold of functions $\mathbb{M}\subset\widetilde{\mathbb{M}}$ such that for a fixed lag time $\tau$, it holds that

[TABLE]

Any such $\mathbb{M}$ is called a transition manifold of the system.

Two remarks are in order:

Note that in the above definition, the $L^{2}_{1/\rho}$ norm is used to measure distances, where $L^{2}_{1/\rho}$ is the space of (equivalence classes of) functions that are square-integrable with respect to the measure induced by the function $\frac{1}{\rho}$, and thus for $f\in L^{2}_{1/\rho}$,

[TABLE]

Closeness with respect to the $L^{2}_{1/\rho}$-norm instead of the $L^{1}$-norm is indeed a strict requirement here, as measuring the quality of a given RC will require a Hilbert space, see Section 2.2. Note that under appropriate assumptions on the system, it holds that $p^{t}_{x}\in L^{2}_{1/\rho}$ for all $x\in\mathbb{X}$. This will be shown in Lemma 3.13 and implies $\widetilde{\mathbb{M}}\subset L^{2}_{1/\rho}$, which together with the requirement $\mathbb{M}\subset\widetilde{\mathbb{M}}$ makes (2) well-defined. 2. 2.

The original definition of $(\varepsilon,r)$-reducibility (see [6, Definition 4.4]), is marginally different from the definition above: Instead of $\mathbb{M}\subset L^{2}_{1/\rho}$, we here require $\mathbb{M}\subset\widetilde{\mathbb{M}}\subset L^{2}_{1/\rho}$. The introduction of this slightly stronger technical requirement allows us to later control a certain embedding error, see Proposition 2.8. Note that the proofs in [6] regarding the optimality of the final reaction coordinate are not affected by this change.

In what follows, we always assume that the process $X_{t}$ is $(\varepsilon,r)$-reducible with small $\varepsilon$ and $r\ll n$.

Remark 2.4**.**

In addition to metastable systems, many other relevant problems possess such transition manifolds. For instance, one can show that, under mild conditions on $f$ and $g$, systems with explicit time scale separation, given by

[TABLE]

also possess a transition manifold since $p^{t}_{(x_{0},y_{0})}$ essentially depends only on $x_{0}$ for $0<\varepsilon\ll 1$ and $t\gg\varepsilon$.

2.2 A measure for the quality of reaction coordinates

We will now present a measure for evaluating the quality of reaction coordinates that is based on transfer operators, first derived in [6]. The Perron–Frobenius operator $\mathcal{P}^{t}\colon L^{1}\rightarrow L^{1}$ associated with the process $X_{t}$ is defined by

[TABLE]

This operator can be seen as the push-forward of arbitrary starting densities, i.e., if $X_{0}\sim u$, then $X_{t}\sim\mathcal{P}^{t}u$.

As $L^{2}_{1/\rho}\subset L^{1}$ (see [6, Remark 4.6]) we can consider $\mathcal{P}^{t}$ as an operator on the inner product space $L^{2}_{1/\rho}$, where it has particularly advantageous properties (see [2, 38, 26]). Here, it is self-adjoint due to the reversibility of $X_{t}$. Moreover, under relatively mild conditions, it does not exhibit any essential spectrum [42]. Hence, its eigenfunctions form an orthonormal basis of $L^{2}_{1/\rho}$ and the associated eigenvalues are real. Now, the significance of the dominant eigenpairs for the system’s time scales is well-known [42]. This is the primary reason for the choice of the $L^{2}_{1/\rho}$-norm in Definition 2.3.

Let $\theta_{i}^{t}$ be the eigenvalues of $\mathcal{P}^{t}$, sorted by decreasing absolute value, and $\psi_{i}$ the corresponding eigenfunctions, where $i=0,1,\dots$. It holds that $\theta_{0}=1$ is independent of $t$, isolated and the sole eigenvalue with absolute value $1$. Furthermore, $\psi_{0}=\rho$. The subsequent eigenvalues decrease monotonously to zero both for increasing index and time. That is,

[TABLE]

The associated eigenfunctions $\psi_{1},\psi_{2},\dots$ can be interpreted as sub-processes of decreasing longevity in the following sense: Let $u\in L^{2}_{1/\rho}$, with $u=\sum_{i=0}^{\infty}\alpha_{i}\psi_{i}$, $\alpha_{i}\in\mathbb{R}$, then

[TABLE]

since for the lag time $\tau>0$ as defined above, there exists an index $d\in\mathbb{N}$ such that $|\theta_{i}^{t}|\approx 0$ for all $t\geq\tau$ and all $i>d$. Hence, the major part of the information about the long-term density propagation of $X_{t}$ is encoded in the $d$ dominant eigenpairs.

The operator $\mathcal{P}^{t}$ describes the evolution of densities of the full process $X_{t}$. In order to monitor the dependence of densities on the reduced coordinate $\xi$ only, we first introduce the projection operator $\Pi_{\xi}\colon L^{1}(\mathbb{X})\rightarrow L^{1}(\mathbb{X})$,

[TABLE]

where the $\rho$-weighted expectation value is taken with respect to the random variable $\mathbf{x}$. This operator is also known as the Zwanzig projection operator from statistical physics [50]. Intuitively, $\Pi_{\xi}$ averages a function over the individual level sets of $\xi$.

The effective transfer operator $\mathcal{P}^{t}_{\xi}\colon L^{1}(\mathbb{X})\rightarrow L^{1}(\mathbb{X})$ associated with $\xi$ is then given by

[TABLE]

see [6]. We now want to preserve the statistics of the dominant long-term dynamics of $X_{t}$ under the projection onto $\xi$, i.e.,

[TABLE]

for $t\geq\tau$, where $\tau$ is some lag time that is long enough for the fast processes, associated with the non-dominant eigenpairs, to have equilibrated. A sufficient condition for (4) is

[TABLE]

that is, the dominant eigenfunctions $\psi_{i}$ must be almost constant along the level sets of $\xi$. This motivates the following definition of a good reaction coordinate:

Definition 2.5**.**

Let $(\psi_{i},\theta_{i}^{t})$ be the eigenpairs of the Perron–Frobenius operator. Let $\tau>0$ and $d\in\mathbb{N}$ such that $\theta_{i}^{t}\approx 0$ for all $i>d$ and $t\geq\tau$. We call a function $\xi\colon\mathbb{X}\rightarrow\mathbb{R}^{r}$ a good reaction coordinate if for all $i=0,\dots,d$ there exist functions $\tilde{\psi}_{i}\colon\mathbb{R}^{r}\rightarrow\mathbb{R}$ such that

[TABLE]

If condition (5) is fulfilled, we say that $\xi$ (approximately) parametrizes the dominant eigenfunctions.

For a formal evaluation of the condition (5), see [6, Corollary 3.6].

2.3 Optimal reaction coordinates

We now justify why reaction coordinates that are based on parametrizations of the transition manifold $\mathbb{M}$ indeed fulfill condition (5). Let $Q\colon L^{2}_{1/\rho}\rightarrow L^{2}_{1/\rho}$ be the nearest-point projection onto $\mathbb{M}$, i.e.,

[TABLE]

Assume further that some parametrization $\gamma\colon\mathbb{M}\rightarrow\mathbb{R}^{r}$ of $\mathbb{M}$ is known, i.e., $\gamma$ is one-to-one on $\mathbb{M}$ and its image in $\mathbb{R}^{r}$. Then the reaction coordinate $\xi\colon\mathbb{R}^{n}\rightarrow\mathbb{R}^{k}$ defined by

[TABLE]

is good in the sense of Definition 2.5 due to the following theorem:

Theorem 2.6** ([6, Corollary 3.8]).**

Let the system be $(\varepsilon,r)$-reducible and $\xi$ defined as in (6). Then for all $i=0,\ldots,d$, there exist functions $\tilde{\psi}_{i}\colon\mathbb{R}^{r}\rightarrow\mathbb{R}$ such that

[TABLE]

Let us add two remarks:

The choice of the $L^{2}_{1/\rho}$-norm in Definition 2.3 is crucial for Theorem 2.6 to hold. 2. 2.

Metastable systems typically exhibit a time scale gap after the $d$ dominant eigenvalues, i.e.,

[TABLE]

Therefore, $\tau$ can always be chosen such that $\left|\theta_{d+1}^{\tau}\right|$ is close to zero and $\left|\theta_{i}^{\tau}\right|,i=0,\ldots,d,$ is still relatively large. Consequently, the denominator in (7) is not too small, and thus the RC (6) is indeed good according to Definition 2.5.

The main task for the rest of the paper is now the numerical computation of an (approximate) parametrization $\gamma$ of $\mathbb{M}$.

2.4 Whitney embedding of the transition manifold

One approach to find a parametrization of $\mathbb{M}$, proposed by the authors in [6], is to first embed $\mathbb{M}$ into a more accessible Euclidean space and to parametrize the embedded manifold. In order to later compare it with our new method, we will briefly describe this approach here.

To construct an embedding $\mathcal{E}$ that preserves the topological structure of $\mathbb{M}$, without prior knowledge about $\mathbb{M}$, a variant of the Whitney embedding theorem can be used. It extends the classic Whitney theorem to arbitrary Banach spaces and was proven by Hunt and Kaloshin in [23].

Theorem 2.7** (Whitney embedding theorem in Banach spaces, [23]).**

Let $\mathbb{V}$ be a Banach space and let $\mathbb{K}\subset\mathbb{V}$ be a manifold of dimension $r$. Let $k>2r$ and let $\alpha_{0}=\frac{k-2d}{k(d+1)}$. Then, for all $\alpha\in(0,\alpha_{0})$, for almost every (in the sense of prevalence) bounded linear map $\mathcal{F}\colon\mathbb{V}\rightarrow\mathbb{R}^{k}$ there exists a $C>0$ such that for all $x,y\in\mathbb{K}$,

[TABLE]

where $\|\cdot\|_{2}$ denotes the Euclidean norm in $\mathbb{R}^{k}$. In particular, almost every $\mathcal{F}$ is one-to-one on $\mathbb{K}$ and its image, and \mathcal{F}^{-1}\big{|}_{\mathcal{F}(\mathbb{K})} is Hölder continuous with exponent $\alpha$.

In particular, almost every such map $\mathcal{F}$ is a homeomorphism between $\mathbb{K}$ and its image in $\mathbb{R}^{k}$, which in short is called an embedding of $\mathbb{K}$ (see e.g. [35, §18]). This means that the image $\mathcal{F}(\mathbb{M})$ will again be an $r$-dimensional manifold in $\mathbb{R}^{k}$, provided that $k>2r$. We will apply this result to the transition manifold, i.e., $\mathbb{V}=L^{2}_{1/\rho}$ and $\mathbb{K}=\mathbb{M}$, and for simplicity restrict ourselves to the lowest embedding dimension, i.e., $k=2r+1$. Any “randomly selected” continuous map $\mathcal{F}:L^{2}_{1/\rho}\rightarrow\mathbb{R}^{2r+1}$ then is an embedding of $\mathbb{M}$.

Unfortunately, there is no practical way to randomly draw from the space of continuous maps on $L^{2}_{1/\rho}$ directly. Instead of arbitrary continuous maps, we therefore restrict our considerations to maps $\mathcal{F}\colon L^{2}_{1/\rho}\rightarrow\mathbb{R}^{2r+1}$ of the form

[TABLE]

where $\sigma$ is some distribution on the (finite-dimensional) space of $(2r+1)\times d$-matrices (e.g., Gaussian matrices). The linear map $\eta\colon\mathbb{X}\rightarrow\mathbb{R}^{2r+1}$, called feature map, is bounded due to the boundedness of $\mathbb{X}$. Maps of the form (8) are therefore continuous on $L^{1}$, and thus in particular on the subspace $L^{2}_{1/\rho}$.

By drawing from the distribution $\sigma$ of the matrices $A$, we can effectively sample maps of form (8). We assume from now on that the embedding property from Theorem 2.7 not only holds for a prevalent subset of general continuous maps, but already for a prevalent subset of the maps of form (8). In other words, we assume that a randomly drawn function of form (8) with linear $\eta$ almost surely is an embedding of $\mathbb{M}$. While the validity of this assumption is far from obvious for general manifolds, there is empirical evidence that for the typically “smooth” transition manifolds $\mathbb{M}$, it is indeed fulfilled. Still, this necessary restriction to the class of linear embedding functions represents a weak point of the transition manifold method that will later be solved by using kernel embeddings instead.

The dynamical embedding of a point $x\in\mathbb{X}$ is then defined by

[TABLE]

This is the Euclidean representation of the density $p^{t}_{x}$, and the set $\{\mathcal{E}(x)\mid x\in\mathbb{X}\}\subset\mathbb{R}^{2r+1}$ is the Euclidean representation of the fuzzy transition manifold. It again clusters around an $r$-dimensional manifold in $\mathbb{R}^{2r+1}$, namely the image $\mathcal{F}(\mathbb{M})$ of the transition manifold under $\mathcal{F}$:

Proposition 2.8**.**

Let the process $X_{t}$ be $(\varepsilon,r)$-reducible with transition manifold $\mathbb{M}$, and $\mathcal{F}\colon L^{2}_{1/\rho}\rightarrow\mathbb{R}^{2r+1}$ and $\mathcal{E}\colon\mathbb{R}^{n}\rightarrow\mathbb{R}^{2r+1}$ defined as in (8) and (9). Then

[TABLE]

Proof.

Let $x\in\mathbb{X}$. By the $(\varepsilon,r)$-reducibility of $X_{t}$ (Definition 2.3), and the fact that $\mathbb{M}\subset\widetilde{\mathbb{M}}$, i.e., $\mathbb{M}$ itself consists of transition densities, there exists an $x^{*}\in\mathbb{X}$ such that $p^{t}_{x^{*}}\in\mathbb{M}$ and $\|p^{t}_{x}-p^{t}_{x^{*}}\|_{L^{2}_{1/\rho}}\leq\varepsilon$. Thus we have

[TABLE]

where $\|\cdot\|_{L^{1}}\leq\|\cdot\|_{L^{2}_{1/\rho}}$ was derived in [6, Remark 4.6]. ∎

Remark 2.9**.**

Together, Theorem 2.7 and Proposition 2.8 guarantee at least a minimal degree of well-posedness of the embedding problem: The embedded manifold $\mathcal{F}(\mathbb{M})$ has the same topological structure as $\mathbb{M}$, and $\mathcal{F}(\widetilde{\mathbb{M}})$ clusters closely around it (if $\|\eta\|_{\infty}$ is small). However, guarantees on the condition of the problem cannot be made. The manifold $\mathbb{M}$ will in general be distorted by $\mathcal{F}$, to a degree that might pose problems for numerical manifold learning algorithms. This problem is illustrated in Figure 2. Such a situation typically occurs if some of the components of the embedding $\mathcal{F}$ are strongly correlated.

Additionally, the Whitney embedding theorem cannot guarantee that the fuzzy transition manifold $\widetilde{\mathbb{M}}$ will be preserved under the embedding, as analytically $\widetilde{\mathbb{M}}$ is not a manifold. Thus, $\mathcal{F}$ is in general not injective on $\widetilde{\mathbb{M}}$.

2.5 Data-driven algorithm for parametrizing the transition manifold

Due to the implicit definition of $\mathbb{M}$, the embedded transition manifold $\mathcal{F}(\mathbb{M})$ is hard to analyze directly. However, as $\mathbb{M}\subset\widetilde{\mathbb{M}}$ and $\mathcal{F}(\widetilde{\mathbb{M}})$ concentrates \big{(}\|\eta\|_{\infty}\varepsilon\big{)}-closely around $\mathcal{F}(\mathbb{M})$, one can expect that any parametrization of the dominant directions of $\mathcal{F}(\widetilde{\mathbb{M}})$ is also a good parametrization of $\mathcal{F}(\mathbb{M})$. We now explain how $\mathcal{F}(\widetilde{\mathbb{M}})$ can be sampled numerically and how this sample can be parametrized.

Let $\mathbb{X}_{N}=\{x_{1},\dots,x_{N}\}$ be a finite sample of state space points, which covers the “dynamically relevant” part of state space, i.e., the regions of $\mathbb{X}$ of substantial measure $\rho$. The exact distribution of the sampling points is not important here. If $\mathbb{X}$ is bounded or periodic, $\mathbb{X}_{N}$ could be drawn from the uniform distribution or chosen to form a regular grid. In practice, it often consists of a sample of the system’s equilibrium measure $\rho$.

The set \mathcal{F}\big{(}\{p^{t}_{x}~{}|~{}x\in\mathbb{X}_{N}\}\big{)} will serve as our sample of $\mathcal{F}(\widetilde{\mathbb{M}})$. Its elements can be computed numerically in the following way: Let $X_{\tau}(x_{0},\omega)$ denote the end point of the time-$\tau$ realization of $X_{t}$ starting in $x_{0}\in\mathbb{X}$ and outcome $\omega\in\Omega$, where $\Omega$ is the sample space underlying the process $X_{t}$. For $x\in\mathbb{X}$, $\tau>0$ fixed as in Definition 2.3 and arbitrarily chosen $\{\omega_{1},\ldots,\omega_{M}\}\subset\Omega$, let $y^{(k)}(x):=X_{\tau}(x,\omega_{k})$. In short, the $y^{(k)}(x),k=1,\ldots,n$, sample the density $p^{t}_{x}$. In practice, the $y^{(k)}(x)$ will be generated by multiple runs of a numerical SDE solver starting in $x$ with $M$ different random seeds (“bursts of simulations”).

With the samples $y^{(k)}(x)$, we approximate $\mathcal{F}(p^{t}_{x})$ by its Monte Carlo estimator:

[TABLE]

Due to Proposition 2.8, the point cloud \mathcal{F}\big{(}\{p^{t}_{x}~{}|~{}x\in\mathbb{X}_{N}\}\big{)}, and for a large enough burst size $M$ also its empirical estimator \widehat{\mathcal{E}}\big{(}\mathbb{X}_{N}\big{)}, then clusters around the $r$-dimensional manifold $\mathcal{F}(\mathbb{M})$ in $\mathbb{R}^{2r+1}$.

Parametrizing \widehat{\mathcal{E}}\big{(}\mathbb{X}_{N}\big{)}, i.e., finding the dominant nonlinear directions in this point cloud in $\mathbb{R}^{2r+1}$, now can be accomplished by a variety of classical manifold learning methods. We assume that we have a method at our disposal that is able to discover the underlying $r$-dimensional manifold within the point cloud \widehat{\mathcal{E}}\big{(}\mathbb{X}_{N}\big{)}, and assign each of the points $\{\widehat{\mathcal{E}}(x)~{}|~{}x\in\mathbb{X}_{N}\}$ a value \tilde{\gamma}\big{(}\widehat{\mathcal{E}}(x)\big{)}\in\mathbb{R}^{r} according to its position on that manifold. For examples of such algorithms see Section 3.1. Hence, $\tilde{\gamma}\colon\widehat{\mathcal{E}}(\mathbb{X}_{N})\rightarrow\mathbb{R}^{k}$ can be seen as an approximate parametrization of $\mathcal{F}(\mathbb{M})$, defined however only at the points $\widehat{\mathcal{E}}(\mathbb{X}_{N})$. Any parametrization of $\mathcal{F}(\mathbb{M})$ in turn corresponds to a parametrization of $\mathbb{M}$, due to $\mathcal{F}$ being an embedding. Finally, any parametrization of $\mathbb{M}$ corresponds to a good reaction coordinate due to Theorem 2.6. Thus, the map $\xi(x)\colon\mathbb{X}_{N}\rightarrow\mathbb{R}^{r}$,

[TABLE]

forms a good reaction coordinate. Note however that it is only defined on the sample points $\mathbb{X}_{N}$.

The strategy of computing reaction coordinates by embedding densities sampled from $\widetilde{\mathbb{M}}$ into $\mathbb{R}^{2r+1}$ by a random linear map and learning a parametrization of the embedded manifold was first presented in [6]. The following algorithm summarizes the overall procedure:

3 Kernel-based parametrization of the transition manifold

The approach described above for learning a parametrization of the transition manifold $\mathbb{M}$ by embedding it into Euclidean spaces requires a priori knowledge of the dimension of $\mathbb{M}$. Also, more importantly, $\mathbb{M}$ might be strongly distorted by the embedding $\mathcal{F}$, as described in Section 2.4. The kernel-based parametrization, which is the main novelty of this work, will address both of these shortcomings by embedding $\mathbb{M}$ into reproducing kernel Hilbert spaces.

3.1 Kernel reformulation of the embedding algorithm

Manifold learning algorithms that can be used in Algorithm 2.1 include Diffusion Maps [14], Multidimensional Scaling [49, 28], and Locally Linear Embedding [37]. These, and many others, require only a notion of distance between pairs of data points. In our case, this amounts to the Euclidean distances between embedded points, i.e., \big{\|}{\mathcal{E}}(x_{i})-{\mathcal{E}}(x_{j})\big{\|}_{2}, which can be computed by the Euclidean inner products \big{\langle}{\mathcal{E}}(x_{i}),{\mathcal{E}}(x_{j})\big{\rangle}, as

[TABLE]

Other compatible algorithms such as Principal Component Analysis are based directly on the inner products. The inner products can be written as

[TABLE]

and the empirical counterpart is

[TABLE]

However, rather than explicitly computing the inner product between the features on the right-hand side, we now assume that it can be computed implicitly by using a kernel function $k\colon\mathbb{X}\times\mathbb{X}\rightarrow\mathbb{R}$, i.e.,

[TABLE]

This assumption, called the kernel trick, is commonly used to avoid the costly computation of inner products between high-dimensional features. However, instead of defining the kernel $k$ based on previously chosen features, one typically considers kernels that implicitly define high- and possibly infinite-dimensional feature spaces. In this way, we are able to avoid the choice of the feature map $\eta$ altogether.

Kernels with this property span a so-called reproducing kernel Hilbert space:

Definition 3.1** (Reproducing kernel Hilbert space [41]).**

Let $k\colon\mathbb{X}\times\mathbb{X}\rightarrow\mathbb{R}$ be a positive definite function. A Hilbert space $\mathbb{H}$ of functions $f\colon\mathbb{X}\rightarrow\mathbb{R}$, together with the corresponding inner product $\langle\cdot,\cdot\rangle_{\mathbb{H}}$ and norm $\|\cdot\|_{\mathbb{H}}=\sqrt{\langle\cdot,\cdot\rangle_{\mathbb{H}}}$ which fulfills

1. (i)

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

$\langle f,k(x,\cdot)\rangle_{\mathbb{H}}=f(x)$* for all $f\in\mathbb{H}$*

is called the reproducing kernel Hilbert space (RKHS) associated with the kernel $k$.

Here, $\overline{A}$ denotes the completion of a set $A$ with respect to $\|\cdot\|_{\mathbb{H}}$. Requirement (ii) implies that

[TABLE]

The inner product between general functions $f,g\in\operatorname{span}\left\{k(x,\cdot)~{}|~{}x\in\mathbb{X}\right\}$ can therefore be expressed as the weighted sum of kernel evaluations: Let

[TABLE]

where the selection of points $x_{i},x^{\prime}_{j}$ depends on $f$ and $g$, respectively. Then

[TABLE]

For functions on the boundary of $\operatorname{span}\left\{k(x,\cdot)~{}|~{}x\in\mathbb{X}\right\}$, the inner product is constructed by the usual limit procedure.

The map $\eta\colon x\mapsto k(x,\cdot)$ can be regarded as a function-valued feature map (the so-called canonical feature map). However, each positive definite kernel is guaranteed to also possess a feature map of at most countable dimension:

Theorem 3.2** (Mercer’s theorem [32]).**

Let $k$ be a positive definite kernel and $\nu$ be a finite Borel measure with support $\mathbb{X}$. Define the integral operator $\mathcal{T}_{k}\colon L^{2}_{\nu}\rightarrow L^{2}_{\nu}$ by

[TABLE]

Then there is an orthonormal basis $\{\sqrt{\lambda_{i}}\hskip 1.00006pt\varphi_{i}\}$ of $\,\mathbb{H}$ consisting of eigenfunctions $\varphi_{i}$ of $\mathcal{T}_{k}$ rescaled with the square root of the corresponding nonnegative eigenvalues $\lambda_{i}$ such that

[TABLE]

The above formulation of Mercer’s theorem has been taken from [33]. The Mercer features $\eta_{i}:=\sqrt{\lambda_{i}}\hskip 1.00006pt\varphi_{i}$ thus fulfill (10) for their corresponding kernel. The usage of the same symbol $\eta$ as for the linear feature map from Section 2.4 is no coincidence, as the Mercer features will again serve the purpose to observe certain features of the full system. In what follows, $\eta(x)$ will always refer to the vector (or $\ell^{2}$ sequence) defined by the Mercer features. If not stated otherwise, $\nu$ will be the standard Lebesgue measure.

Example 3.3**.**

Examples of commonly used kernels are:

Linear kernel: $k(x,x^{\prime})=x^{\top}x^{\prime}$. One sees immediately that (10) is fulfilled by choosing $\eta_{i}(x)=x_{i}$, $i=1,\ldots,n$ (also spanning the Mercer feature space). 2. 2.

Polynomial kernel of degree $p$: $k(x,x^{\prime})=(x^{\top}x^{\prime}+1)^{p}$. It can be shown that the Mercer feature space is spanned by the monomials in $x$ up to degree $p$. 3. 3.

Gaussian kernel: $k(x,x^{\prime})=\exp\left(-\frac{1}{\sigma}\left\lVert x-x^{\prime}\right\rVert_{2}^{2}\right),$ where $\sigma>0$ is called the bandwidth of the kernel. Let $p\in\mathbb{N}$ and $\mathbf{p}=(p_{1},\dots,p_{n})$ with $p_{1}+\ldots+p_{n}=p$ be a multi-index. The Mercer features of $k$ then take the form

[TABLE]

see **[48]**, where

[TABLE]

Let $\mathcal{F}_{k}$ denote the density embedding based on the Mercer features of the kernel $k$, i.e.,

[TABLE]

and let $\mathcal{E}_{k}(x):=\mathcal{F}_{k}(p^{t}_{x})$. The amount of information about $p^{t}_{x}$ preserved by the embedding $\mathcal{F}_{k}$ depends on the choice of the kernel $k$. For the first two kernels in Example 3.3, the information preserved has a familiar stochastic interpretation (see, e.g., [33, 39, 47]):

Let $k$ be the linear kernel. Then

[TABLE]

i.e., the means of $p^{t}_{x_{1}}$ and $p^{t}_{x_{1}}$ coincide. 2. 2.

Let $k$ be the polynomial kernel of degree $p>1$. Then

[TABLE]

i.e., the first $p$ moments $\mathbf{m}_{i}$ of $p^{t}_{x_{1}}$ and $p^{t}_{x_{1}}$ coincide.

Remark 3.4**.**

In practice, comparing the first $p$ moments often is enough to sufficiently distinguish the transition densities that constitute the transition manifold. However, densities that differ only in higher moments cannot be distinguished by $\mathcal{F}_{k}$, which means that for the above two kernels, $\mathcal{F}_{k}$ is not injective on $\mathbb{M}$. Therefore $\mathcal{F}_{k}$ does not belong to the prevalent class of maps that is at the heart of the Whitney embedding theorem 2.7. We can therefore not utilize the Whitney embedding theorem to argue that the topology of $\mathbb{M}$ is preserved under $\mathcal{F}_{k}$. Instead, in Section 3.2, we will use a different argument to show that the embedding is indeed injective for the Gaussian kernel (and others).

Still, by formally using the Mercer dynamical embedding $\mathcal{E}_{k}$ in (9) (abusing notation if there are countably infinitely many such features), and using the kernel trick, we can now reformulate Algorithm 2.1 as a kernel-based method that does not require the explicit computation of any feature vector. This is summarized in Algorithm 3.1.

3.2 Condition number of the kernel embedding

We will now investigate to what extent the kernel embedding preserves the topology and geometry of the transition manifold.

3.2.1 Kernel mean embedding

We derived the kernel-based algorithm by considering the embedding $\mathcal{F}_{k}$ of the transition manifold into the image space of the Mercer features in order to highlight the similarity to the Whitey embedding based on randomly drawn features. Of course, the Mercer features never had to be computed explicitly.

However, in order to investigate the quality of this embedding procedure it is advantageous to consider a different, yet equivalent embedding map: The transition manifold can be directly embedded into the RKHS by means of the kernel mean embedding operator.

Definition 3.5**.**

Let $k$ be a positive definite kernel and $\mathbb{H}$ the associated RKHS. Let $p$ be a probability density over $\mathbb{X}$. Define the kernel mean embedding of $p$ by

[TABLE]

Note that $\mu(p)$ and $\hat{\mu}(p)$ are again elements of $\mathbb{H}$ and that for $\nu$ in (12) being the Lebesgue measure we obtain $\mu(p)=\mathcal{T}_{k}\hskip 1.00006ptp$. Further, one sees that

[TABLE]

where the inner product $\langle\cdot,\cdot\rangle$ refers to the Euclidean inner product or the inner product in $\ell^{2}(\mathbb{N}_{0})$, dependent on whether $\mathcal{F}_{k}(p)$ is finite or countably infinite. Thus, for investigating whether the embedding $\mathcal{F}_{k}$ preserves distances or inner products between densities, we can equivalently investigate the embedding $\mu$. This is advantageous as injectivity and isometry properties of the kernel mean embedding are well-studied.

3.2.2 Injectivity of the kernel mean embedding

A first important result is that $k$ can be chosen such that $\mu$ is injective. Such kernels are called characteristic [21]. In [47], several conditions for characteristic kernels are listed, including the following:

Theorem 3.6** ([47, Theorem 7]).**

The kernel $k$ is characteristic if for all $f\in L_{2},~{}f\neq 0$ it holds that

[TABLE]

Condition (15) is known as the Mercer condition, which is, for example, fulfilled by the Gaussian kernel from Example 3.3. The Mercer features of such a kernel are particularly rich.

Theorem 3.7**.**

Assume that the kernel satisfies the Mercer condition (15). Then the eigenfunctions $\{\psi_{i}\}$ of $\mathcal{T}_{k}$ form an orthonormal basis of $L^{2}(\nu)$.

For more details, see, e.g., [41, 48]. It is easy to see that for kernels fulfilling (15), $\mu$ as a map from $L^{2}$ to $\mathbb{H}$ is Lipschitz continuous:

Lemma 3.8**.**

Let $k$ be a characteristic kernel with Mercer eigenvalues $\lambda_{i}$, $i\in\mathbb{N}_{0}$. Then $\mu\colon L^{2}\rightarrow\mathbb{H}$ is Lipschitz continuous with constant

[TABLE]

Proof.

As $\mu$ is linear, it suffices to show that $\|\mu(f)\|_{\mathbb{H}}\leq c\|f\|_{2}$ for all $f\in L_{2}$. We obtain

[TABLE]

where (11) was used in the last line. By expanding $k$ into its Mercer features via (13), this becomes

[TABLE]

By Theorem 3.7, the $\varphi_{i}$ form an orthonormal basis of $L^{2}$, and thus

[TABLE]

Thus, if the kernel is characteristic, the structure of the TM and the fuzzy TM are qualitatively preserved under the embedding.

Corollary 3.9**.**

Let $k$ be a characteristic kernel and let $X_{t}$ be $(\varepsilon,r)$-reducible. Then $\mu(\mathbb{M})$ is an $r$-dimensional manifold in $\mathbb{H}$, and for all $x\in\mathbb{X}$ it holds that

[TABLE]

Proof.

By Lemma 3.8, the map $\mu\colon L^{2}\rightarrow\mathbb{H}$ is continuous (and furthermore injective), and thus $\mu(\mathbb{M})$ is again an $r$-dimensional manifold. For $x\in\mathbb{X}$, consider now any $f\in L^{2}_{1/\rho}$ with $\|f-p^{t}_{x}\|_{L^{2}_{1/\rho}}\leq\varepsilon$. For $g:=\mu(f)$, we then get

[TABLE]

Remark 3.10**.**

This result should be seen as an analogue to Proposition 2.8 for the Whitney-based TM embedding. In short, for characteristic kernels, the injectivity and continuity of $\mu$ guarantee that the image of $\mathbb{M}$ under $\mu$ is again an $r$-dimensional manifold in $\mathbb{H}$, and Corollary 3.9 guarantees that the embedded fuzzy transition manifold $\mu(\widetilde{\mathbb{M}})$ still clusters closely around $\mu(\mathbb{M})$ (if $\sqrt{\lambda_{0}}$ and $\|\sqrt{\rho}\|_{\infty}$ in Corollary 3.9 are small). This again guarantees a minimal degree of well-posedness of the problem.

3.2.3 Distortion under the kernel mean embedding

Unlike the Whitney embedding, the kernel embedding now allows us to derive conditions under which the distortion of $\mathbb{M}$ is bounded. We have to show that the $L^{2}_{1/\rho}$-distance between points on $\mathbb{M}$ is not overly decreased or increased by the kernel mean embedding. To formalize this, we consider measures for the manifold’s internal distortion, following the notions of metric embedding theory [1].

We call the embedding well-conditioned if both the

[TABLE]

are small (close to one). Here, $\mu$ denotes the embedding corresponding to a characteristic kernel.

Due to the Lipschitz continuity of $\mu$ (see Lemma 3.8) and $\|\cdot\|_{L^{2}}\leq\|\sqrt{\rho}\|_{\infty}\hskip 1.00006pt\|\cdot\|_{L^{2}_{1/\rho}}$, we have

[TABLE]

thus bounding the expansion.

Contraction bound: regularity requirement.

Unfortunately, it is not possible even for characteristic kernels to derive a bound for the contraction that holds uniformly for all $p,q\in L^{2}_{1/\rho}$, as the following proposition shows. Nevertheless, we will be able to give reasonable bounds under some regularity- and dynamic assumptions, (20) and (23), respectively.

Proposition 3.11** (Unbounded inverse embedding).**

Assume the kernel embedding operator $\mu$ has absolutely bounded orthonormal eigenfunctions $\varphi_{i}$ with corresponding nonnegative eigenvalues $\lambda_{i}$ (arranged in nonincreasing order). Assume $\lim_{i\rightarrow\infty}\lambda_{i}=0$. Then there exist functions $p,q\in L^{2}_{1/\rho}$ such that

[TABLE]

for any arbitrarily small $\varepsilon>0$.

Proof.

See Appendix A. ∎

The assumptions of Proposition 3.11 are fulfilled for example for the Gaussian kernel. A similar, but non-quantitative result has been derived in [47, Theorem 19]. The idea behind its proof and the proof of Proposition 3.11 is that, if $p$ and $q$ vary only in higher eigenfunctions $\varphi_{i}$ of the embedding operator $\mu$ (see also Theorem 3.2), the $\mathbb{H}$-distance can become arbitrarily small. If, however, we can reasonably restrict our considerations to the subclass of functions whose variation in the higher $\varphi_{i}$ is small compared to the variation in the lower $\varphi_{i}$, a favorable bound can be derived. Let the expansion of $h=p-q$ be given by

[TABLE]

with the sequence $(\tilde{h}_{0},\tilde{h}_{1},\ldots)\in\ell^{2}$. Now, for any $i_{\text{max}}\in\mathbb{N}$ such that there exists an index $i\leq i_{\text{max}}$ with $\tilde{h}_{i}\neq 0$, define the factor

[TABLE]

This factor bounds the contribution of the higher Mercer eigenfunctions to $h$ by the contribution of the lower ones, hence it is a regularity bound:

[TABLE]

Thus, for an individual $h$, we can bound the distortion of the $L^{2}$-norm under $\mu$ with the help of $c(h,i_{\text{max}})$.

Lemma 3.12**.**

Let $h\in L^{2}$, $i_{\text{max}}\in\mathbb{N}$, and $c(h,i_{\text{max}})$ be defined as in (19). Then

[TABLE]

Proof.

See Appendix A. ∎

We from now on make the assumption that for every index $i_{\text{max}}$ there exists a constant $c^{*}_{i_{\text{max}}}>0$ such that

[TABLE]

for all $p^{\tau}_{x_{1}},p^{\tau}_{x_{2}}\in\mathbb{M}$. The existence and form of this constant strongly depends on the shape of the Mercer eigenfunctions, hence the kernel. However, we motivate the existence of such a global constant by the observation that higher Mercer eigenfunctions typically consist of high Fourier modes, and that these modes decay quickly under the dynamics. Therefore, high Mercer eigenfunctions should have a negligible share of the $p^{\tau}_{x}$ and the differences $p^{\tau}_{x_{1}}-p^{\tau}_{x_{2}}$. For such $p^{\tau}_{x_{1}},p^{\tau}_{x_{2}}$, we thus have

[TABLE]

Contraction bound: dynamical requirements.

Note that (21) is only an intermediate step for deriving a contraction bound, as the relevant distance measure in Definition 2.5 is the $L^{2}_{1/\rho}$-norm, for reasons detailed in Section 2.2, and (21) measures the density distance in the $L^{2}$-norm. Unfortunately, a naive estimation yields

[TABLE]

While due to ergodicity $1/\sqrt{\rho}$ is indeed defined on all of $\mathbb{X}$, it becomes large in regions of small invariant measure $\rho$, i.e., “dynamically irrelevant” regions. This would lead to a very large upper bound for the contraction. For general $h$, a more favorable estimate is indeed difficult to obtain. For us, however, $h=p^{\tau}_{x_{1}}-p^{\tau}_{x_{2}}$, and we can utilize that these “dynamically irrelevant” regions are almost never visited by the system.

To formalize this, we require one additional assumption that can be justified by the metastability of the system. One defining characteristic of metastable systems is the phenomenon that essentially any trajectory moves nearly instantaneously111Here, “nearly instantaneously” is to be understood in the sense that there is an “attraction time” $\tau_{a}\ll\tau$ such that starting from essentially any initial condition the system will enter one of the metastable sets within time $\tau_{a}$ with overwhelming probability. Thus, choosing $\tau$ as an intermediate time that is larger than the non-metastable time scales of local fluctuations, is essential. However, as we discuss in Remark 3.15, it is also imperative not to choose $\tau$ too large. into one of the metastable sets before continuing. With $A_{i},\ldots,A_{d}\subset\mathbb{X}$ denoting these sets, we can thus assume that the probability density $p^{\tau}_{x}(y)$ to move from $x$ to $y$ in time $\tau$ depends almost only on the probabilities to (instantaneously) move to the sets $A_{i}$ from $x$ (denoted by $c_{i}(x)$, thus $\sum_{i=1}^{d}c_{i}(x)\leq 1$) and the probabilities to then move from $A_{i}$ to $y$ in time $\tau$ (denoted by $p^{\tau}_{A_{i}}(y)$), i.e.,

[TABLE]

To be more precise, we require for all $x,y\in\mathbb{X}$ that

[TABLE]

which is equivalent to

[TABLE]

for some positive function $\delta\colon\mathbb{X}\times\mathbb{X}\rightarrow\mathbb{R}$ with $\|\delta\|_{\infty}\leq\delta^{*}\ll 1$. The positivity of $\delta$ comes from the fact that there is a miniscule, but positive probability (density) to move from $x$ to $y$ without first equilibrating inside a metastable set, thus $p^{\tau}_{x}(y)>\sum_{i=1}^{d}c_{i}(x)p^{\tau}_{A_{i}}(y)$.

With this, we can bound the invariant density $\rho$ from below as follows:

[TABLE]

The $b_{i}$ can be seen as the equilibrium probability mass almost instantaneously attracted to $A_{i}$. For every important metastable set this will not be too small.

As a first step, (24) allows us to bound the $L^{2}_{1/\rho}$-norm of $p^{\tau}_{x}$ by their $L^{1}$-norm:

Lemma 3.13**.**

Let the assumption (23) hold and $b_{i},i=1,\ldots,d,$ be defined as in (24). Then for any $x\in\mathbb{X}$ it holds that

[TABLE]

Proof.

See Appendix A. ∎

This shows that indeed $p^{\tau}_{x}\in L^{2}_{1/\rho}$, as required by Definition 2.3. Further, Hölder’s inequality gives

[TABLE]

where $|\mathbb{X}|$ denotes the Lebesgue measure of the state space. This now also allows us to bound the $L^{2}_{1/\rho}$-norm of the $p^{\tau}_{x}$ by their $L^{2}$-norm:

Lemma 3.14**.**

Let the assumption (23) hold and $b_{i},i=1,\ldots,d,$ be defined as in (24). Then for any $x_{1},x_{2}\in\mathbb{X}$ it holds that

[TABLE]

Proof.

See Appendix A. ∎

Of course, due to the squared norm on the left-hand side, this is not a Lipschitz bound. However, recall that our main motivation for deriving a bound for the contraction is to show that large distances in $L^{2}_{1/\rho}$ are not overly compressed under the embedding into $\mathbb{H}$, as illustrated in Figure 2. We therefore abstain from deriving such a bound for very small distances in $L^{2}_{1/\rho}$ and only estimate the contraction of pairs of densities $p^{\tau}_{x_{1}},p^{\tau}_{x_{2}}$ with

[TABLE]

for some constant $C>0$. That $C$ is reasonably large is discussed in Remark 3.15 below. For such differences, we can then relate the $L^{2}$ to the $L^{2}_{1/\rho}$-norm, i.e.,

[TABLE]

Together with Lemma 3.12 and assumption (20), this gives

[TABLE]

which is our contraction bound.

Remark 3.15**.**

For the distortion of the transition manifold under the embedding the essential property is that the global “spanning structure” of the manifold is well preserved. In other words, the embedded $p_{x}^{\tau},p_{y}^{\tau}$ should be well-separated for $x,y\in\mathbb{X}$ from different metastable sets $A_{i}$. Since the embedding is continuous, the transition paths connecting them will be preserved as well.

As $p^{\tau}_{A_{i}}\to\rho$ as $\tau\to\infty$ for every $i$, it is important that $\tau$ is not too large, such that the transition manifold is a meaningful object. In a other words, $\tau$ should be such that $p^{\tau}_{A_{i}}$ and $p^{\tau}_{A_{j}}$ are sufficiently distinct for $i\neq j$. Thus, we require that $\smash{\|p^{\tau}_{A_{i}}-p^{\tau}_{A_{j}}\|_{L^{2}_{1/\rho}}}$ is sufficiently large for $i\neq j$, hence $C$ can be chosen as a constant such that $1/C$ is reasonably small.

Remark 3.16**.**

The bounds (18) and (28) guarantee the well-posedness of the overall embedding and parametrization problem as the relevant expansion and contraction of the transition manifold cannot become arbitrarily large. We will support this statement with numerical evidence (see Section 4.2) showing that the distortion is, in practice, indeed small. It should be noted, however, that we do not expect the analytically derived bounds to perform well as quantitative error estimates, as many of the estimates that led to them are rather rough.

4 Illustrative examples and applications

4.1 Reaction coordinate of the Müller–Brown potential

As a first illustrating example, we compute the reaction coordinate of the two-dimensional Müller–Brown potential [34] via the new kernel-based Algorithm 3.1. The potential energy surface (see Figure 3 (a)) possesses three local minima, where the two bottom minima are separated only by a rather shallow energy barrier. Correspondingly, the system’s characteristic long-term behavior is determined by the rare transitions between the minima. These transitions happen predominantly along the potential’s minimum energy pathway (MEP), which is shown as white dashed line and was computed using the zero temperature string method [15, 16].

For two sets $A,B\subset\mathbb{X}$ and a starting point $x\in\mathbb{X}$, the committor function $q_{AB}(x)$ is defined as the probability that the process hits set $A$ before hitting set $B$, provided it started in $x$ at time zero. For a precise definition see [42]. For the Müller–Brown potential, the committor function associated with the top left and bottom right energy minima, shown in Figure 3 (b), can be considered an optimal reaction coordinate [18]. Therefore, we use the (qualitative) comparison with the committor function as a benchmark for our reaction coordinate. Note that the computation of the committor function requires global knowledge of the metastable sets and is often not a practical option for the identification of reaction coordinates.

The governing dynamics is given by an overdamped Langevin equation (1), which we solve numerically using the Euler–Maruyama scheme. At inverse temperature $\beta=0.05$, the lag time $\tau=0.03$ falls between the slow and fast time scales and is thus defined to be the intermediate lag time. The test points $\{x_{1},\ldots,x_{N}\}$ required by Algorithm 3.1 are given by a regular $32\times 32$ grid discretization of the domain $[-1.5,1.5]\times[-0.5,2.5]$. For the embedding, the Gaussian kernel

[TABLE]

with bandwidth $\sigma=0.1$ is used and for the manifold learning task in Algorithm 3.1, the diffusion maps algorithm with bandwidth parameter $0.1$. The reaction coordinate $\xi$ for the test points is shown in Figure 3 (c). We observe remarkable resemblance to the committor function.

Remark 4.1**.**

The kernel evaluations used for the RKHS embedding of densities and the kernel evaluations used in the diffusion maps algorithm should not be mixed up as they serve different purposes. The former is used to embed the state space densities into $\mathbb{H}$, while the latter is used to approximate the Laplace–Beltrami operator on the manifold in $\mathbb{H}$ that is to learn (this is the principle on which the diffusion maps algorithm is based). Even though the Gaussian kernel is a popular choice due to its favorable characteristics, one has great freedom in choosing a kernel for the RKHS-embedding, whereas in the classical diffusion maps algorithm, predominantly the Gaussian kernel is used, so the repeated use of the Gaussian kernel does not constitute a connection. Moreover, the fact that in this example identical bandwidth parameters were used was a mere coincidence. We do not see a way to unify these kernel evaluations, neither on a conceptual nor algorithmic level.

4.2 Distortion under the Whitney and kernel embeddings

We now demonstrate the distortion of the transition manifold under the embedding via the conventional Algorithm 2.1 (Whitney embedding) and our new Algorithm 3.1 (kernel-based embedding). To this end, we consider the two-dimensional potential depicted in Figure 4 (a). Good reaction coordinates of a diffusion process in this potential should parametrize the “horseshoe-like” MEP shown as a white dashed line. Such a reaction coordinate is depicted in Figure 4 (b).

As test points $x_{i}$, $N=200$ uniformly distributed random points in the region $\mathbb{X}=[-2,2]^{2}$ are drawn. Per point, $M=100$ short trajectories of length $\tau=2$ are computed to sample the $p^{\tau}_{x_{i}}$.

Whitney embedding

For the Whitney embedding, the expected manifold dimension $r=1$ is assumed to be known in advance. To demonstrate the different effects of “good” and “bad” embedding functions, two $2r+1$-dimensional linear observables $\eta:\mathbb{R}^{2}\rightarrow\mathbb{R}^{3}$ were chosen:

[TABLE]

and the corresponding embedding functions $\mathcal{F}_{g},\mathcal{F}_{b}$ constructed via (8). The coefficients of $A_{g}$ of the observable function $\eta_{g}$ were chosen randomly via the Matlab command

[TABLE]

which resulted in the matrix

[TABLE]

Under the embedding $\mathcal{F}_{g}$, all parts of the transition manifold can indeed be distinguished very well, see Figure 5 (a). This is due to the fact that two distinct points of the transition pathway of the potential, for example on the two opposite “branches”, are never mapped to the same point under $\eta_{g}$.

On the other hand, the matrix $A_{b}$ of the “bad” observable function $\eta_{b}$ was intentionally constructed to consist of three row vectors that are pairwise almost linearly dependent, and that essentially ignore the $x_{1}$-component of state space points:

[TABLE]

This way, points on the two opposite branches of the transition pathway but with the same $x_{2}$-coordinate are mapped to almost the same point in $\mathbb{R}^{3}$. The result is an embedding of the transition manifold in which the two branches can hardly be distinguished, see Figure 5 (b). This makes the numerical identification of the manifold structure extremely difficult.

Note that the judgment of quality of the embedding function has to be performed manually after the embedding, as it is impossible to reliably choose good embedding functions without detailed a priori knowledge of the global structure of the transition manifold or transition pathway. While in our experience, randomly chosen coefficients typically result in “good-enough” embedding functions, this uncertainty in the numerical algorithm should be seen as one of the main reasons to use the more consistent kernel embeddings instead.

Kernel embedding

For the kernel embedding, we again utilize the Gaussian kernel (29) with bandwidth $\sigma=10^{-3}$. The result of the kernel embedding of the test points is an approximation of the kernel distance matrix

[TABLE]

(see Algorithm 3.1), which cannot be visualized directly. We thus apply the Multidimensionl scaling (MDS) algorithm to $D$, in order to visualize the level of similarity between the embedded densities.

Given a distance matrix $D$, MDS generates points $z_{i}\in\mathbb{R}^{k}$ in a Euclidean space of a chosen dimension $k\in\mathbb{N}$ such that the pairwise distance between the $z_{i}$ corresponds to the distances in $D$. For an overview of different MDS methods, see for example [49]. We here use the implementation of classical MDS given by the cmdscale method in Matlab.

The MDS representation of the kernel distance matrix for $k=2$ is shown in Figure 6 (a). The curved structure of the MEP is immediately visible. Moreover, it is also possible to visualize the corresponding $L^{2}_{1/\rho}$ and $L^{2}$ distance matrices via MDS, i.e., the matrices

[TABLE]

The results are shown in Figure 6 (b) & (c).

The MDS representation of $D$ is structurally very similar to $D_{L^{2}_{1/\rho}}$ and $D_{L^{2}}$. This suggests that the $L^{2}_{1/\rho}$ and $L^{2}$ distances are preserved very well under $\mu$, up to a constant factor. To confirm this, we now compute the empirical maximum distortion of the $L^{2}_{1/\rho}$ metric based on the given test points, i.e., $\mathscr{D}_{N}(\mu):=\mathscr{C}_{N}(\mu)\hskip 1.00006pt\mathscr{E}_{N}(\mu)$ where

[TABLE]

For large enough $N$, we expect $\mathscr{D}_{N}(\mu)$ to be a good estimator for the true distortion $\mathscr{D}(\mu)$.

The blue graph in Figure 7 shows the dependence of the empirical distortion on the kernel parameter $\sigma$. Here the minimum is $\mathscr{D}_{N}(\mu)=36.7$ at $\sigma\approx 10^{-3}$. Interpreting $\mathscr{D}_{N}(\mu)$ as the condition number of the kernel-based embedding problem, the problem can be described as reasonably well-conditioned.

Analogously, we can define the empirical maximum distortion of the Whitney embedding as $\mathscr{D}_{N}(\mathcal{F}):=\mathscr{C}_{N}(\mathcal{F})\hskip 1.00006pt\mathscr{E}_{N}(\mathcal{F})$, where

[TABLE]

For a given embedding $\mathcal{F}$, this distortion can again be computed numerically. For the “good” embedding $\mathcal{F}_{g}$, we obtain $\mathscr{D}_{N}(\mathcal{F}_{g})\approx 6\cdot 10^{2}$, while for the “bad” embedding $\mathcal{F}_{b}$, we obtain $\mathscr{D}_{N}(\mathcal{F}_{b})\approx 7\cdot 10^{3}$. The kernel embedding is therefore much better conditioned than both Whitney embeddings.

Remark 4.2**.**

Analogously, we can also define and compute the maximum distortion $\mathscr{D}_{N}(\mu)$ of the $L^{2}$-metric (red graph in Figure 7). Here, for $\sigma\approx 10^{-1}$ we obtain $\mathscr{D}_{N}(\mu)=2.7$, i.e., the embedding becomes nearly isometric. This is not surprising as it has been shown in [47] that for radial kernels $k_{\sigma}(x,y)=\sigma^{-d}g(\sigma^{-1}\|x-y\|)$ where $g$ is bounded, continuous, and positive definite, it holds that

[TABLE]

The Gaussian kernel belongs to this class of kernels. We thus expect that by increasing the sample number $M$ of the transition densities and further decreasing $\sigma$, the distortion can be reduced further. However, recall that for our application, only the distortion of the $L^{2}_{1/\rho}$ distance is relevant.

4.3 Alanine dipeptide

We now demonstrate the applicability of Algorithm 3.1 to realistic, high-dimensional molecular systems by computing reaction coordinates of the Alanine dipeptide. The peptide, depicted in Figure 8 (a), consists of 22 atoms, the state space $\mathbb{X}$ thus has the dimension $n=66$.

It is well-known that the essential long-term behavior of this system is governed by the metastable transitions between four local minima of the potential energy surface (PES) [13, 44]. These minima are clearly visible when projecting the PES onto two specific backbone dihedral angles $(\varphi,\psi)$ that we call essential from now on (see Figure 8 (b)). The transition between the metastable states happens along minimal energy pathways that we aim to reveal with our reaction coordinate. Note however that no information about the existence of the two essential dihedral angles was used in our experiments, and we perform all of the analysis on the full 66-dimensional data.

The simulations were performed using the Gromacs molecular dynamics software [3]. We consider the molecule in explicit aqueous solution at temperature $400\,$K (the water molecules are discarded prior to further analysis). To generate the test points $x_{i}$, $N=1000$ snapshots from a long, equilibrated trajectory were subsampled. This guarantees that the $x_{i}$ cover the dynamically relevant regions of $\mathbb{X}$, i.e., the metastable sets and transition pathways. The values of the dihedral angles $\varphi$ and $\psi$ of the test points are shown in Figure 10 (the $x$- and $y$-coordinates of the points). We see that the metastable sets and transition pathways from Figure 8 (b) are adequately covered. Note however that the projection onto the $(\varphi,\psi)$-space here serves only illustrative purposes; we continue to work with the test points in the full 66-dimensional space.

The intermediate lag time $\tau=20\,\mathrm{ps}$ falls between the slow and fast time scales. For strategies to estimate $\tau$ prior to simulation, see [5]. For each test point $x_{i}$, $M=100$ simulations of length $\tau$ were performed, which took $40\,\mathrm{h}$ on a 96 core cluster. The resulting point clouds $\{y_{i}^{(l)}$, $l=1,\ldots,M\}$ are samplings of the densities $p^{\tau}_{x_{i}}$.

To compute the kernel distance matrix $D$ from the simulation data, the Gaussian kernel (29) with bandwidth $\sigma=0.1$ was chosen. For the plug-in manifold learning algorithm that is applied to $D$, the diffusion maps algorithm with bandwidth parameter $0.01$ was used. The analysis of the simulation data was performed in Matlab and took $4$ minutes on a 4 core laptop.

Figure 9 (a) shows the leading spectrum of the diffusion map matrix that was computed based on $D$. The first diffusion map eigenvalue is always one, and the associated eigenvector carries no structural information. Therefore, a spectral gap after the third sub-dominant eigenvalue indicates that the underlying transition manifold is intrinsically three-dimensional. The associated three subdominant eigenvectors now are the final reaction coordinate. For each of the 1000 test points, the values of the three eigenvectors are shown in Figure 9 (b). This can be seen as the embedding of the test points into the reaction coordinate space. Here we observe four clusters of points, and three connecting paths. In Figure 10, the values of the dihedral angles $\varphi$ and $\psi$ at the test points are compared to the values of the three components of the computed reaction coordinates, shown in color. We see that areas of almost constant color correspond to the four metastable sets from Figure 8. Thus, the computed reaction coordinate is able to identify the metastable sets and resolve transitions between them.

5 Conclusion and future work

In this work, we have analyzed the embedding of manifolds that lie in certain function spaces into reproducing kernel Hilbert spaces. Moreover, we have proposed efficient numerical algorithms for learning parametrizations of these embedded manifolds from data. The question is motivated by the recent insight that parametrizations of the so-called transition manifold, a manifold consisting of the transition density functions of a stochastic system, are strongly linked to reduced coordinates for that system. The method can thus be used for coarse graining a given system.

Compared to previous approaches based on random embeddings into a Euclidean space, the new kernel-based approach eliminates the need to know the transition manifold dimension a priori. Furthermore, if a universal kernel is used, the topological structure of the transition manifold is guaranteed to be preserved under the embedding. We have derived bounds for the geometric distortion of the transition manifold under the RKHS embedding, which can be interpreted as the condition of the overall coarse graining procedure. Correspondingly, the numerical algorithm was demonstrated to be very robust, especially when compared to random embeddings, and, in realistic applications, we obtained very favorable results regarding algorithmic distortion bounds.

There are several new avenues to use the broader theory of kernel embeddings to characterize the kernel embedding of transition manifolds. First, we plan to improve the theoretic distortion bounds derived in Section 3.2 by considering different established interpretations of the metric defined by $d(p,q)=\|\mu(p)-\mu(q)\|_{\mathbb{H}}$. For an overview, see [47].

Recently, the spectral theory of transfer operators was extended to reproducing kernel Hilbert spaces in [27]. The usefulness of this new theory for the data-driven conformation analysis of molecular systems was demonstrated in [25]. As the transition manifold can be defined via the transfer operator222The fuzzy transition manifold is the image of all Dirac densities under the transfer operator, i.e., $\widetilde{\mathbb{M}}=\{\mathcal{T}^{t}\delta_{x}\mid x\in\mathbb{X}\}$., it seems natural to attempt to relate the embedded transition manifold to the kernel transfer operators and corresponding embedded transfer operators defined in [27].

Finally, as illustrated in [7, 8, 9], RKHSs can act as linearizing spaces in the sense that performing linear analysis in the RKHS can capture strong nonlinearities in the original system. A typical example is the problem of linear separability in data classification: A data set which is not linearly separable might be easily separated when mapped into a nonlinear feature space. In our current context, this means that efficient linear manifold learning methods might be suitable to parametrize the embedded manifold, if the kernel is chosen appropriately. We will investigate whether a corresponding theory can be developed.

Acknowledgements

The authors would like to thank the anonymous reviewers for constructive comments and suggestions that helped to improve the paper.

This research has been partially funded by Deutsche Forschungsgemeinschaft (DFG) through grant CRC 1114, projects A01, B03, and B06.

Appendix A Proof of the distortion bounds

Proof of Proposition 3.11.

First, note that $L^{2}_{1/\rho}\subset L^{2}$, as for $p\in L^{2}_{1/\rho}$

[TABLE]

Let $(\lambda_{i},\varphi_{i})$ be the eigenpairs of the integral operator $\mathcal{T}_{k}$, ordered in decreasing order of $\lambda_{i}$. For arbitrary $p\in L^{2}_{1/\rho}$ consider the decomposition into the basis \big{\{}\varphi_{i}\big{\}}_{i\in\mathbb{N}} of $L^{2}$:

[TABLE]

Select $i_{\text{max}}\in\mathbb{N}$ such that there exists an index $i\leq i_{\text{max}}$ with $\tilde{h}_{i}\neq 0$, and such that \lambda_{i}<\big{(}\varepsilon/\|\sqrt{\rho}\|_{\infty}\big{)}^{2} for all $i\geq i_{\text{max}}$, and define

[TABLE]

Then,

[TABLE]

Further, using that \big{\{}\sqrt{\lambda_{i}}\hskip 1.00006pt\varphi_{i}\big{\}}_{i\in\mathbb{N}} forms an orthonormal basis of $\mathbb{H}$, and that $\mu$ is a linear operator, we get

[TABLE]

Thus we get

[TABLE]

and with (30) finally

[TABLE]

Proof of Lemma 3.12.

As \big{\{}\sqrt{\lambda_{i}}\hskip 1.00006pt\varphi_{i}\big{\}}_{i\in\mathbb{N}_{0}} forms an orthonormal basis of $\mathbb{H}$, we obtain

[TABLE]

Further, \big{\{}\varphi_{i}\big{\}}_{i\in\mathbb{N}_{0}} forms an orthonormal basis of $L_{2}$, and so

[TABLE]

Thus,

[TABLE]

Proof of Lemma 3.13.

With assumption (23) and (24), we can write the left-hand side of (25) as

[TABLE]

As for all $x\in\mathbb{X}$ it holds $c_{i}(x)\geq 0$ and $\sum_{i=1}^{d}c_{i}(x)\leq 1$, we obtain

[TABLE]

With this, we can estimate the integral $(\star)$ as

[TABLE]

Proof of Lemma 3.14.

The proof is completely analogous to the proof of Lemma 3.13, while in the estimate corresponding to (31) we use that

[TABLE]