Efficient Approximation of Molecular Kinetics using Random Fourier Features

Slow kinetic processes of molecular systems can be analyzed by computing dominant eigenpairs of the Koopman operator or its generator. In this context, the Variational Approach to Markov Processes (VAMP) provides a rigorous way of discerning the quality of different approximate models. Kernel methods have been shown to provide accurate and robust estimates for slow kinetic processes, but are sensitive to hyper-parameter selection, and require the solution of large-scale generalized eigenvalue problems, which can easily become computationally demanding for large data sizes. In this contribution, we employ a stochastic approximation of the kernel based on random Fourier features (RFFs), to derive a small-scale dual eigenvalue problem which can easily be solved. We provide an interpretation of this procedure in terms of a finite randomly generated basis set. By combining the RFF approach and model selection by means of the VAMP score, we show that kernel parameters can be efficiently tuned, and accurate estimates of slow molecular kinetics can be obtained for several benchmarking systems, such as deca alanine and the NTL9 protein.


I. INTRODUCTION
The automated extraction of essential information about thermodynamic and kinetic properties of molecular systems from large-scale computer simulations remains one of the major open problems in computational physics and chemistry to this day.The challenge consists in processing very long time series of high-dimensional data in such a way that relevant features of the molecular system are determined automatically, with limited user input or intervention.The precise definition of what is regarded as relevant depends on the context, but in many applications, the determination of slow relaxation processes and their associated metastable sets is among the most important features to be detected.This applies, for example, to protein folding, ligand binding, or protein-protein association [1][2][3] .
Let us highlight two milestones that have shaped the field in the last two decades: The first was the inception of Markov State Models (MSMs) in the late 1990s and early 2000s [4][5][6] , which allowed for the construction of easily interpretable discrete kinetic models, from which a wealth of information could be extracted.The second milestone was the variational approach to conformational dynamics (VAC) for reversible systems 7,8 , and its extension called variational approach to Markov processes (VAMP) 9 , which provided a rigorously defined scalar objective function (called VAMP score) to assess the quality of kinetic models.Combined, these developments have led to a standard workflow for the construction of MSMs, usually consisting of a linear dimensionality reduction using time-lagged independent component analysis (TICA) 10 , geometric clustering, as well as MSM estimation and validation 11,12 .It was also shown that MSMs and the VAC are closely related to the Koopman operator framework 13,14 , and in particular the extended dynamic mode decomposition (EDMD) algorithm 15,16 , which have primarily emerged in the analysis of fluid dynamics and related fields.The Koopman framework is therefore a unifying theme for the construction of kinetic models.Extensions of the EDMD algorithm can also be used to approximate the generator of dynamical systems (e.g., the Kolmogorov backward operator for stochastic differential equations) as shown in 17 .
Despite the remarkable success of MSMs, the construction process remains somewhat unsatisfactory, at least from a theoretical perspective.The initial TICA step presupposes that there is a kinetically informative linear subspace of the initial set of descriptors used for MSM building.In order to move beyond this requirement, algorithms based on deep neural networks (DNNs) were introduced in recent years 18,19 , allowing to learn in one step the complete non-linear transfor-mation from raw atomic coordinates to low-dimensional features that optimize the VAMP score.
Both the MSM construction process (including state definition) and neural networks require a non-linear optimization step, which is typically solved by means of stochastic algorithms.A significant number of hyper-parameters need to be tuned in both cases, such as network architecture, network depth and height for DNNs, or number of cluster centers, number of iterations, convergence threshold, etc. for clustering-based MSMs.
A different avenue that has been considered is to employ reproducing kernels, which have long been known in machine learning 20,21 as a powerful non-linear model class with rich theoretical properties.Kernel methods have also been demonstrated in many applications to be fairly robust if the available data is only small-to medium-sized.Kernel-based versions of EDMD were suggested in [22][23][24] , and applied to high-dimensional molecular simulations in 25 .Extracting slow timescales from kernel representations requires the solution of a large-scale generalized eigenvalue problem, where the size of the matrix depends on the number of snapshots.Even for moderate data size, the solution of this problem can be challenging.Moreover, although kernel methods usually require only a few hyper-parameters (such as the bandwidth of the hugely popular Gaussian kernel), these need to be tuned accurately, requiring to solve the above-mentioned eigenvalue problem many times, which can be prohibitively time-consuming.
In this work, we show that kernel methods can be made tractable for the analysis of molecular kinetics by employing Random Fourier Features (RFFs) 26 , a well-known low-rank approximation technique for kernel matrices, to convert the original large-scale eigenvalue problem into a much smaller one, which can be easily solved.We derive the associated small-scale matrix problem both for the Koopman operator and its generator, and show that this approach can be interpreted as applying EDMD with a specific randomly chosen basis set.We demonstrate that the VAMP score can be effectively used for hyper-parameter tuning of the kernel method.The combination of these ideas leads to an efficient approximation algorithm, which is robust with respect to the data size and the stochastic Fourier feature selection.Therefore, our method presents a promising new direction for the automated analysis of molecular kinetics.
The remainder of this paper is structured as follows: In Section II, we review elements of the Koopman formalism, the EDMD approximation algorithm, reproducing kernel Hilbert spaces (RKHSs), and Koopman modeling using RKHSs.Then, in Section III, we present our main contribution, an approximation algorithm for the Koopman operator or generator based on Random Fourier Features.Further computational details, especially regarding model selection, are dis-cussed in Section IV, while applications to benchmarking problems in molecular dynamics are shown in Section V.

II. THEORETICAL BASICS
We will briefly introduce the stochastic Koopman operator and its generator as well as numerical methods to estimate these operators and their spectral properties from data.

A. Stochastic Dynamical Systems
Computer simulations are routinely used to explore the thermodynamics and kinetics of molecular systems.Many available implementations of molecular dynamics (MD) can be described as a stochastic process {X t } t≥0 , where t is the time and X t ∈ X ⊂ R d the state of the system at time t.In particular, a widely used model in this context is given by stochastic differential equations where F : R d → R d and G : R d → R d×d are vector-and matrix-valued fields, respectively, called the drift and diffusion of the SDE, and W t is d-dimensional Brownian motion.The covariance matrix of the diffusion is denoted by As far as the thermodynamics is concerned, the primary goal of molecular simulations is to sample the Boltzmann distribution where V is the potential energy at position x and β −1 = k B T is the inverse temperature.Most simulation protocols are built to achieve this goal, i.e., the Boltzmann distribution is invariant, and long time averages converge to spatial averages with respect to µ (ergodicity).A standard example is Brownian or overdamped Langevin dynamics with friction parameter γ.

B. Koopman Operator and Generator
When it comes to inferring kinetic information about X t , significant progress has been made in recent years by considering an operator-based approach to the statistics of the stochastic process X t .For a fixed lag time t ≥ 0, the Koopman operator 13 K t acts on a functions φ : X → R by taking the conditional expectation: In other words, evaluation of the function K t φ at x yields the expectation of φ after starting the dynamics at x, and evolving until time t.The linear operators K t satisfy the semigroup equation K s+t = K s K t , which implies the existence of another linear operator L , called Koopman generator, such that the following linear differential equation holds for the function φ t = K t φ : For SDE dynamics of the form (1), stochastic calculus shows that L is a second-order differential operator, and we obtain the Kolmogorov equation where we used the colon notation for the Frobenius inner product between matrices, i.e., A : B = ∑ i, j A(i, j)B(i, j).For overdamped Langevin dynamics, the generator is given by We will summarily refer to the Koopman operators or the generator as dynamical operators for the dynamical system X t .If the dynamics are reversible with respect to the invariant distribution µ, then the generator L is symmetric with respect to the inner product As a consequence, all eigenvalues of the generator L are real-valued and non-positive, i.e., for all solutions of the equation we have κ i ≤ 0. An object of prime interest for molecular kinetics are low-lying eigenvalues 0 ≤ −κ i 1, as these encode slow motions corresponding to metastable behavior encountered in many examples of molecular systems 6,27 .The connection between close-to-zero eigenvalues and metastability is also apparent in the fact that eigenfunctions of the generator are also eigenfunctions of the Koopman operator at any t, with eigenvalue λ i (t) = e −κ i t .If −κ i 1, then the corresponding eigenvalues λ i (t) decay slowly, as expected for metastable dynamics.

C. Finite-dimensional Approximation
The basic data-driven approximation algorithm for the Koopman operator is called extended dynamic mode decomposition (EDMD) 28 , its adaptation to the generator is then called gEDMD 17 .
Both methods require choosing a finite dictionary of observable functions {φ 1 , . . ., φ n }.The projections of the dynamical operators onto the linear span of these functions can be represented by matrices: with mass and stiffness matrices given by Given a long trajectory {x 1 , x 2 , . . ., x m+1 } of the process X t , the data-based estimators for these matrices are given by Ĝ where we assumed the time step between samples to equal t in the formula for Ât .This is not required in practice if a sliding-window estimator is used, see 12,29 for details.
For reversible dynamics, the stiffness matrix A L and its estimator ÂL above can be replaced by the following energy-like expressions 17 : The latter estimator only requires first-order derivatives and retains the symmetry of A L even for finite data.
Eigenvalues of the dynamical operators can be approximated by diagonalizing the data-driven estimators Kt = Ĝ−1 Ât and L = Ĝ−1 ÂL , or equivalently, by solving one of the generalized eigenvalue problems

D. Collective Variables
For large molecular systems, the dictionary is typically defined in terms of a set of lowerdimensional descriptors z ∈ R N , N ≤ d, often called collective variables (CVs), such as intramolecular distances or angles.Formally, CVs are given by a smooth mapping ξ : X → R N .Choosing a dictionary { φi } n i=1 of functions depending only on CV space, that is φi (x) = φ i (ξ (x)), does not lead to any conceptual changes to the EDMD approach.When approximating the generator, the gEDMD method is straightforwardly adapted by computing all x-derivatives in the definition of the generator using the chain rule.This leads to a particularly elegant reformulation of the symmetric estimator (3) in terms of an effective diffusion matrix where ∇ x ξ is the d × N-dimensional Jacobian of the mapping ξ .The symmetric estimator (3)   then becomes as shown in 17 .Note that once the effective diffusion matrices have been computed, only the basis functions and their derivatives in CV-space are required.

E. Reproducing Kernel Hilbert Space
A critical modeling decision for the (g)EDMD approach is the choice of the dictionary {φ 1 , . . ., φ n }.Recent works [22][23][24][25] have shown that accurate results, even for large systems, can be obtained using methods based on reproducing kernel Hilbert spaces (RKHS) generated by a kernel function k.These approaches lead to powerful approximation spaces, yet only depend on few model parameters.
A reproducing kernel Hilbert space (RKHS) 21,30 is a Hilbert space H, with inner product •, • H , of continuous functions on X, defined by a symmetric, positive-definite two-argument function k(x, y), called the kernel.Each such kernel corresponds to a unique RKHS 30 .Given k, one can define a map from X into H, called the feature map, by Functions f in the RKHS then satisfy the reproducing property Example II.1 Arguably the most widely used kernel is the Gaussian radial basis function (RBF) analogously generates an RKHS of periodic functions 31 .
For many popular kernels, including the Gaussian RBF kernel, the associated RKHS is infinitedimensional, and densely embedded into most relevant function spaces.Kernels therefore provide powerful approximation spaces, reducing the problem of dictionary selection to the choice of a single function, the kernel, and its hyper-parameters, such as the bandwidth for the Gaussian RBF kernel.

F. RKHS Representation of Dynamical Operators
Given a kernel k with RKHS H, the analogues of the matrices G, A t and A L in Section II C are given by linear operators on H, see 23,24 for a derivation: Here, Φ is again the feature map.The eigenvalue problem (4) then turns into an infinitedimensional generalized eigenvalue equation for eigenfunctions ψ i ∈ H: Natural empirical estimators for the RKHS operators in ( 6) are given by Ĝ The range of these operators is in the linear span of the feature map at the data sites, H m = span{Φ(x l )} m l=1 .Therefore, it makes sense to restrict these operators to the m-dimensional space H m , which leads to the matrix representations and ( 7) can be replaced by a matrix generalized eigenvalue problem The derivatives in the definition of K L X are taken with respect to first argument of the kernel function.It is important to note that, as expected for a kernel method, assembling the matrices (8)   requires only kernel evaluations at the data sites, along with kernel derivatives and coefficients of the SDE in the generator setting.For reversible dynamics, there is again an alternative symmetric formulation, which leads to the following Hermitian generalized eigenvalue problem instead of (9): The derivations of ( 8), ( 9) and (10) were shown in 23,24 , we repeat them for the reader's convenience in Appendix A.
The linear problems ( 9) and ( 10) represent large-scale generalized eigenvalue problems.Solving these problems efficiently for large data sets is one of the central challenges common to most applications of kernel methods.In the following section, we present an approach based on a randomized low-rank representation obtained from Fourier transformation of the kernel function.

III. LOW-RANK REPRESENTATIONS OF KERNEL EIGENVALUE PROBLEMS
In this section, we will present our main contribution: an approach to solving the generalized eigenvalue problems ( 9) and ( 10) using a randomized low-rank decomposition.

A. Random Fourier Features
Random Fourier Features (RFFs) have been introduced in 26 as a low-rank approximation framework for large kernel matrices.The basis for this approach is Bochner's theorem, stating that a translation-invariant kernel k, which is normalized to k(x, x) = 1 for all x, possesses the following stochastic representation 32 : Here, ρ is a probability measure in frequency space, called the spectral measure.This representation extends to the derivatives of the kernel function, namely where α ∈ N d is a multi-index and ω α = (ω 1 ) α (1) . . .d) .The derivatives in the above expression are again taken with respect to the first argument.For the Gaussian RBF kernel with bandwidth σ , the spectral measure is also Gaussian with bandwidth σ −1 .For the periodic Gaussian kernel, the spectral measure is discrete, supported on all wave vectors π L k for k ∈ Z d , with probabilities where I k is the modified Bessel function of the first kind, see 33 .Thus, drawing independent samples from the spectral measure is easily accomplished for many popular kernel functions.

B. Representation of Kernel Matrices
We apply the RFF representation to the generalized eigenvalue problems ( 9) and ( 10) by substituting it into the kernel matrices K X , K t X , K L X given by (8).We again denote samples in real space (in the form of a long trajectory 1 ) by {x l } m+1 l=1 .Additionally, we assume there are p samples {ω u } p u=1 in frequency space, sampled from the spectral measure ρ.Replacing expectation values by finite-sample averages, we obtain 1 Instead of using one long trajectory, it would also be possible to extract a training data set {(x l , y l )} m l=1 from many short trajectories, see 34 for details.
with RFF feature matrices That is, the rows of M correspond to the data points and the columns of M to the randomly sampled features.When approximating the generator, we use the representation of derivatives given in (11), to find with generator feature matrix M L , which can be written compactly using the Frobenius inner product as .
We will address the reversible case further below after analyzing the RFF-based approximation in more detail.Together, ( 12) and ( 14) provide low-rank approximations of the kernel matrices which can be easily assembled by evaluation of complex exponentials.

C. Solution of the Generalized Eigenvalue Problems
Let us now use the approximations ( 12) and ( 14) to solve the generalized eigenvalue problems (9), which read (note that the normalization 1 p can be omitted): If p ≤ m, we can employ a standard trick to express this equivalently as a lower-dimensional eigenvalue problem.Defining v i = M H w i , we can verify directly that The last equality shows that all non-zero eigenvalues of ( 15) can be equivalently computed by the lower-dimensional dual problem: This p-dimensional generalized eigenvalue equation is the central problem to be solved within the context of RFF-based kernel approximation.If p is much smaller than m, (16) already leads to computational savings compared to (9).
In fact, it is not even necessary to assemble the matrices in ( 16), because the equation can be solved by operating only on the RFF feature matrices M, M t , M L .To this end, we apply a whitening transformation, which is widely used within the context of (g)EDMD.Using a truncated singular value decomposition of M, M ≈ UΣW H , and retaining r ≤ p components of the SVD, we obtain a linear transformation T = WΣ −1 , which eliminates the matrix on the right-hand sides of ( 16): Applying the same transformation to the left-hand sides of ( 16) leads to the reduced matrices It is then sufficient to diagonalize these reduced matrices in order to solve (16).If the truncation rank of M is significantly smaller than p, additional computational savings can be harnessed this way.Algorithm 1 summarizes our solution method in compact form.

D. Interpretation as (g)EDMD Problem
The transformation of the dual problem (16) to its reduced form is completely analogous to what is called the whitening transformation for (g)EDMD 29 .This is not a coincidence: consider a finite dictionary defined by Using the notation of Section II C, the empirical mass and stiffness matrices for this basis set are given by Ĝ which are precisely the constituent matrices in ( 16) 2 .We conclude that the solution of the transposed problem is equivalent to the eigenvalue estimation based on (g)EDMD using a random dictionary of complex plane waves.With this interpretation at hand, we can also derive a convenient symmetric formulation for reversible dynamics, upon replacing the estimator for ÂL by its symmetric version (3): To compute the reduced matrix, we just have to apply the transformation T from Section III C from left and right, i.e., resulting again in a symmetric estimator.

Algorithm 1 RFF-based Spectral Approximation of the Operator or Generator
Input: data matrix X = [x 1 , . . ., x m+1 ] ∈ R d×(m+1) , kernel function k with spectral measure ρ, number of features p, truncation rule for singular values.
3: Non-reversible case: form matrix M t or M L as in ( 13) or ( 14).
4: Compute SVD of M, choose rank r according to truncation rule: M ≈ UΣW H .
5: Form reduced matrix R t or R L according to (17) or (18).
6: Compute eigenpairs of reduced problem R t u i = λi (t)u i or R L u i = κi u i .
7: Transform to original RFF basis: in the reversible case.Note that this bound assumes a dense and state-dependent diffusion field a, it reduces to linear dependence on d for diagonal diffusions, such as overdamped Langevin dynamics.Finally, the diagonalization of the reduced matrix amounts to O(r 3 ) floating point operations.Comparing this to full kernel method, we see that the most drastic computational savings will result from being able to choose p small compared to m, as it reduces a factor O(m 3 )

IV. METHODS
In this section, we will discuss aspects of the practical application of our spectral RFF method and also introduce the systems studied in the results section below.For the numerical examples,

B. Spectral Analysis and Clustering
Solution of the eigenvalue problem for the reduced matrix R provides estimates of the leading eigenvalues λi (t) or κi , respectively, and their associated eigenfunctions ψi .The first point to note is that ( 16) is a complex eigenvalue problem, so the eigenvectors can be scaled by arbitrary complex numbers.Before analyzing the eigenvectors, we simply apply a grid search over complex phase factors to determine the one that minimizes the imaginary part of all eigenvector entries and use the rescaled eigenvectors for further analysis.
Second, we transform the eigenvalues of the Koopman operator into implied timescales defined as , which have units of time, and can be interpreted as physical relaxation timescales for the i-th eigenmode of the Koopman operator.Since timescales are independent of the lag time in the absence of projection and estimation error 36 , we monitor convergence of implied timescales with increasing lag time t in order to select an optimal lag time 11 .We also compute the timescales when approximating the generator, in this case they are simply given by In order to analyze metastable sets, we evaluate the top eigenfunctions ψi at all data points, and apply the spectral clustering method PCCA 37 to assign each data point x l to each metastable set q with a degree of membership χ q (x l ), which is a number between zero and one.We then identify data point x l as part of metastable set q if χ q (x l ) ≥ 0.6, the remaining data points are treated as transition states.

C. Systems and Simulation Setups
We study four different examples in order to illustrate different aspects of the performance of the proposed method.The first is overdamped Langevin dynamics (2) in a two-dimensional model potential, which is a minor modification of the Lemon-Slice potential 38 , given in polar coordinates as see Figure 1 for a visualization.We can easily generate large amounts of simulation data by applying the standard Euler-Maruyama integrator, the elementary simulation time step is set to ∆ t = 10 −3 .Friction γ and inverse temperature β are both set to one.Reference implied timescales are obtained from a Markov state model based on a k-means discretization of the two-dimensional state space, using 50 discrete states and m = 10 5 data points.
The second system is molecular dynamics simulation data of the alanine dipeptide in explicit water, at temperature T = 300 K, employing the Langevin thermostat.We generated a total of one million data points at 1 ps time spacing.Reference implied timescales were computed by a Markov state model in the well-known two-dimensional space of backbone dihedral angles φ and ψ, using a 30 × 30 box discretization and all one million data points.
The third example is molecular dynamics simulation data of the deca alanine peptide in explicit water at temperature T = 300 K, see 39 for a description of the simulation setup.The data comprises a total of four million time steps at 1 ps time spacing.As reference values, we refer to the MSM analysis presented in 40 , which is built on a 500-state k-means discretization after applying timelagged independent component analysis (TICA) 10 in the space of the peptide's sixteen backbone dihedral angles.
Finally, we re-analyze the 39-residue protein NTL9, which was simulated on the Anton supercomputer by D.E.Shaw Research, see 41 for details.The data comprises around 1.2 million frames corresponding to a time spacing of ∆ t = 2 ns.As a reference, we apply linear TICA on the set of 666 minimal heavy atom inter-residue distances.We also rank these distance according to the fraction of simulation time where a contact is formed, that is, the distance does not exceed 0.35 nm, see 40 , and calculate TICA models using only the first [20, 100, 200, 300, 400, 500, 666] of these distances.
V. RESULTS

A. Lemon-Slice Potential
We start by analyzing the Lemon-Slice potential (19), shown in Figure 1 A, which features four metastable states corresponding to the wells of the potential energy.We use this system as a first test case for approximating the generator L , due to the availability of all relevant quantities in analytical form.
Figure 1 B shows the VAMP score as a function of the bandwidth for a selection of data sizes m and feature numbers p.We see that while m = 1000 data points are not sufficient, for m = 5000 data points, the score stabilizes for a range of bandwidths σ ∈ [0.3, 0.8], and already attains optimal values for very small p.We confirm this finding by plotting the first three non-trivial eigenvalues of the generator, for m = 5000 and p = 50, as a function of the bandwidth in Figure 1 C. Indeed, accurate estimates for eigenvalues can be obtained in the suggested regime of bandwidths.After calculating the eigenvectors of the reduced matrix for σ = 0.4, and applying the PCCA analysis, we see that the metastable structure is perfectly recovered, see Figure 1 D. Thus, the VAMP score reliably guides the hyper-parameter search towards a very efficient approximation of the leading spectral components of this system.

B. Alanine Dipeptide
For alanine dipeptide, we approximate the Koopman operator with lag time t, using the periodic Gaussian kernel on the well-known two-dimensional reaction coordinate space given by the dihedral angles φ and ψ. Figure 2 A shows the familiar free energy landscape of the system in the reaction coordinate space.For all RFF-based approximations, we downsample the data set to m = 20 000 points, corresponding to a time spacing of 50 ps.We apply the same model selection protocol as before, and find that a range of kernel bandwidths σ ∈ [0.4,1.0] can be stably identified as optimal in terms of the VAMP score.We observe in Figure 2 B that again a small number of Fourier features p = 50 is sufficient to arrive at an optimal VAMP score.We also confirm that the optimal parameter regime is stable across different lag times.Note that the VAMP score necessarily decreases with increasing lag time, as all eigenvalues but λ 0 (t) decay exponentially with t.
As shown in Figure 2 C, the optimal range of bandwidths indeed allows for accurate estimation of the three leading implied timescales of the system.Further, we also verify that the metastable decomposition of the dihedral space into four states is correctly recovered using the RFF model, see Figure 2 D. The dramatic rank reduction for the kernel matrix achieved by using a small number of Fourier features allows us to easily compute eigenvalues and test scores for a broad range of parameters, using m = 20000 data points in real space.Solving the same problem by means of standard methods would be significantly more costly if the same data size was used.

C. Deca Alanine
We use the deca alanine example to demonstrate that, by means of low-rank kernel methods, generator models can be efficiently learned on a reaction coordinate space that is more than just one-or two-dimensional, and that dynamical hypotheses can be tested this way.We project the system onto the space of its sixteen interior backbone torsion angles.The generator of the full molecular dynamics is not readily available in closed form.It is well-known, however, that if observed only on a subset of position space, many systems driven by thermostatted molecular dynamics behave like the reversible overdamped Langevin dynamics at long timescales 42 .For that reason, we set the full space diffusion matrix to 1 β times the identity and employ our reversible estimator.The effective diffusion (5) then turns into where ∇ξ is the Jacobian of the mapping from Euclidean coordinates to the peptide's backbone dihedral angles, which can be evaluated analytically.Note that this approach can be seen as a dynamical hypothesis: we choose to model the position space dynamics as an overdamped Langevin process, and test if its projection onto the dihedral angle space can recover all relevant features of the full MD simulation.
We find in Figure 3 A that a kernel bandwidth σ ∈ [2.0, 8.0] leads to an optimal VAMP score for the projected generator.Using about m = 1000 data points in real space and p = 300 Fourier features is sufficient to obtain stable results, rendering the calculation highly efficient.The optimality of this regime of bandwidths in confirmed by plotting the first three non-trivial eigenvalues of the generator, as shown in Figure 3 B. We note that these eigenvalues differ from those of the reference MSM by a constant factor.This is due to the re-scaling of time inherent in assuming overdamped Langevin dynamics in position space, rendering the effective dynamics too fast.However, we also see in Figure 3 that all slow eigenvalues are indeed re-scaled by the same factor.Moreover, the PCCA analysis of the slowest eigenvectors obtained for the optimal RFF model also recovers the metastable states of the deca alanine system, corresponding to the folded and unfolded states as well as two intermediates.The effective overdamped Langevin model learned by the kernel method thus retains all major features of the original molecular dynamics, up to the re-scaling of time.

D. NTL9
We use the NTL9 data set to demonstrate once again the ability of the proposed method to efficiently and robustly analyze large data sets.We downsample the available trajectory data to around m = 15.000data points at 200 ns time spacing.As mentioned before, m = 15.000 is still a significant data size for kernel methods.We project the MD simulations onto the 666-dimensional space of minimal heavy atom inter-residue distances.It is known from previous publications that these coordinates capture the folding process of NTL9 very well 40,43 .Indeed, just using linear TICA on the top-ranked 300 distances allows for an accurate estimation of the slowest implied timescale t 1 , associated to the folding process, see the top green line in Figure 4 B.

VS
A: VAMP Score m=1000, p=100 m=1000, p=300 m=1000, p=500 m=5000, p=100 m=5000, p=300 m=5000, p=500 the VAMP score and implied timescales for different bandwidths, keeping the number of Fourier features fixed at p = 300.We see in Figures 4 A-B that both the optimal value of the VAMP score as well as the timescale estimates remain stable across random selections of the input distances (magenta lines in panel A, red lines in panel B).Convergence of timescales is therefore mainly dependent on the number of Fourier features p, which can be automatically tuned using the VAMP score, and not on the selection of input coordinates.Moreover, we also find that under the random distance selection, the optimal regime of kernel bandwidths changes as the dimension of the input coordinate space is significantly smaller.This is reliably reflected by the behavior of the VAMP score, see again Figure 4 A. magnitude.
The main promise of our method is the ability to analyze large data sets using a very generic, and in fact randomly generated basis set.The success on the method depends on the number of required Fourier features.Future research will need to test if this number indeed remains small even for larger systems, which would make the method scalable beyond the benchmarking examples studied here.It will also be interesting to see if kernel functions with more specific properties can be used, for instance, kernels incorporating more detailed symmetries than only translational invariance, or kernels acting on mixed types of domain (periodic and non-periodic).Finally, combining our approach with re-weigthing techniques would allow application of the method to en-hanced sampling simulation data.
Let us break down the computational effort required for Algorithm 1: The SVD of M requires O(mp 2 ) operations.Assembling the reduced matrix costs O(rpm) operations for the Koopman operator.In the case of generator approximation, we obtain an additional dependence on the dimension, due to the required contraction of derivatives of the kernel function.We can upperbound the resulting effort as O(mp(d 2 + r)) in the non-reversible case and by O(m(rpd + r 2 d 2 ))

FIG. 1 .
FIG. 1. Results for the Lemon-Slice potential.A: Contour of the potential (19).B: VAMP score as a function of the kernel bandwidth for different data sizes m and feature numbers p. C: Leading non-trivial eigenvalues of the generator for m = 5000 and p = 50, as a function of the bandwidth.Black lines indicate the Markov model reference values.D: Decomposition into four metastable states based on eigenvectors for σ = 0.4, p = 50, m = 5000.All error bars are based on twenty independent simulations.

FIG. 2 .
FIG. 2. Results for alanine dipeptide.A: Free Energy (in kJ/mol) in two-dimensional dihedral space.B: VAMP Score for selected feature sizes p and lag times t as a function of the kernel bandwidth.C: Implied timescales for t = 100 ps and p = 50 as a function of the bandwidth.D: Metastable decomposition obtained for σ = 0.6, p = 50 and t = 100 ps.

Figure 4 A
Figure4A shows that using Gaussian bandwidth parameters σ ∈[10, 50] can be identified as optimal in terms of the VAMP score, again across different lag times.In Figure4B, we first compare the slowest timescale t 1 obtained for σ = 15 as a function of the lag time, for different feature numbers p (blue lines).We see that already for p ≥ 300, satisfactory convergence of t 1 with t can be obtained, reducing the computational cost by about a factor two thousand compared to the full kernel eigenvalue problem.Interestingly, timescale convergence improves with p in about the same way as the performance of TICA improves as more and more of the ranked distances are added.However, we show that the kernel method is more robust with respect to the selection of input coordinates, by randomly selecting d = 50 out of the 666 distances, and recomputing both

FIG. 3 .
FIG. 3. Results for the Koopman generator on backbone dihedral angle space of the deca alanine peptide.A: VAMP score for selected data sizes m and feature sizes p as a function of the kernel bandwidth.B: Eigenvalues for m = 1000 and p = 300 as a function of the bandwidth.The reference MSM results are shown in black.Re-scaling MSM eigenvalues by the average ratio between optimal RFF and MSM eigenvalues leads to the magenta lines.C-F: Representative structures for each of the four PCCA states based on the RFF model at m = 1000, p = 300, σ = 4.0.

FIG. 5 .
FIG.5.Top Row: Representations of the two PCCA states obtained for σ = 15, p = 300, t = 2 µs.For each PCCA state, we show the fraction of simulation time during which each residue-residue pair forms a contact, see also40 .Bottoms Row: Representative protein structure for both PCCA states.