From Biased to Unbiased Dynamics: An Infinitesimal Generator Approach
Authors: Timothée Devergne
Vladimir R Kostic
Michele Parrinello
Massimiliano Pontil
Venue: NeurIPS 2024 poster (2024)
Keywords: Molecular dynamics
stochastic differential equations
Primary area: machine_learning_for_physical_sciences
TL;DR: we prove that learning the true dynamics  rom biased simulations is possible via the Infinitesimal Generator and develop the method that does it efficently

Abstract: We investigate learning the eigenfunctions of evolution operators for time-reversal invariant stochastic processes, a prime example being the Langevin equation used in molecular dynamics. Many physical or chemical processes described by this equation involve transitions between metastable states separated by high potential barriers that can hardly be crossed during a simulation. To overcome this bottleneck, data are collected via biased simulations that explore the state space more rapidly. We propose a framework for learning from biased simulations rooted in the infinitesimal generator of the process {and the associated resolvent operator}. We contrast our approach to more common ones based on the transfer operator, showing that it can provably learn the spectral properties of the unbiased system from biased data.  In experiments, we highlight the advantages of our method over transfer operator approaches  and recent developments based on generator learning, demonstrating its effectiveness in estimating eigenfunctions and eigenvalues. Importantly, we show that even with datasets containing only a few relevant transitions due to sub-optimal biasing, our approach
recovers relevant information about the transition mechanism.

From Biased to Unbiased Dynamies: An Infinitesimal Generator Approach
Timothée Devergne
Vladimir R. Kostic
CSML & ATSIM, Istituto Italiano di Tecnologia
CSML, Istituto Italiano di Tecnologia University of Novi Sad vladimir kostic
timothee.co timothee.devergueee der.iieeeedverreiiii..iii.
Michele Parrinello
Massimiliano Pontil
ATSIM, Istituto Italiano Tecnologia michele, par.eeeeerriiellliiii..
CSML, Istituto Italiano Tecnologia AI Centre, University College London massimiliano pon.illiiiiitttte
Abstract
We investigate learning the eigenfunctions of evolution operators for time-reversal invariant stochastic processes, a prime example being the Langevin equation used in molecular dynamics. Many physical or chemical processes described by this equation involve transitions between metastable states separated by high potential barriers that can hardly be crossed during a simulation. To overcome this bottle- neck, data are collected via biased simulations that explore the state space more rapidly. We propose a framework for learning from biased simulations rooted the infinitesimal generator of the process and the associated resolvent operator. We contrast our approach to more common ones based on the transfer operator, showing that
provably learn the spectral properties of the unbiased system from biased data. experiments, we highlight the advantages of our method over transfer operator approaches and recent developments based on generator learning, demonstrating its effectiveness in estimating eigenfunctions and eigenvalues. Im- portantly, show that even with datasets containing only few relevant transitions due to sub-optimal biasing, our approach recovers relevant information about the transition mechanism.
Introduction
Dynamical systems and stochastic differential equations (SDEs) provide general mathematical framework to study natural phenomena, with broad applications science and engineering. Langevin SDEs, the main focus of this paper, widely used to simulate physical processes such as protein folding or catalytic reactions [see e.g. [4]] and references therein]. main objective is to describe the dynamics of the process, forecast its evolution from starting state, ultimately gaining insights macroscopic properties of the system.
In molecular dynamics, the motion of a molecule is sampled according to potential energy U Uxx), where the state vector æ represents the positions all the atoms. Specifically, the Langevin equation d.X, - VU (X.)dtt dt + odW, describes the stochastic behavior of the system at thermal equilibrium, where X, is the random position of the state at time t, the scalar or is a multiple of the square root the system's 's temperature, and W, is vector random variable describing thermal fluctuations (Brownian motion). Most often, the atoms evolve in metastable states that are separated by barriers which can hardly be crossed during simulation. For instance, for protein the free energy barrier between the folded and unfolded states is larger than thermal agitation, making the transition between the
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

two states a rare event. Consequently, long trajectories need to be simulated before such interesting events are observed. In fact, one needs observe many events to get the relevant thermodynamies (free energy) and kinetics (transition rates) information R4 Beyond molecular dynamics, the slow mixing behavior of many systems modeled by SDEs is major bottleneck the study of rare events, and so designing methodologies which can accelerate the process is paramount.
A general idea to overcome the above problem is to perturb the system dynamics. One important approach which has been put in place in molecular dynamics is the so-called "bias potential enhanced sampling" 25 44 IT. The main idea to add to the potential energy bias potential V. thereby lowering the barrier and allowing the system state to be explored more rapidly. To make this approach tractable in large systems, V often chosen as a function of few wisely selected variables called collective variables (CVs). For instance, if chemical reaction involves bond breaking, physical intuition suggests to choose the distance between the reactive atoms 26 29, However, for complex processes, hand-crafted CVs might be "suboptimal", meaning that some of the degrees of freedom important for the transition are not taken into account, making the biasing process inefficient.
In recent years, machine learning approaches have been employed to find the most relevant CVs 277. A key idea is to use available dynamical information to construct the CVs 7 42 For instance, if one can identify the slowest degrees of freedom of the system, one can accurately describe the transitions between metastable states. These approaches are based on learning the transfer operator of the system, which models the conditional expectation of a function (or observable) of the state at a future time, given knowledge of the state the initial time. It is learned from the behavior of dynamical correlation functions at large lag times which reflects the slow modes of the system. The leading eigenfunctions of the learned transfer operator can then be used as CVs in biased simulations. Moreover, they provide valuable insights into the transition mechanism, such as the location of the transition state ensemble [48]. Still, this approach suffers from the same shortcoming described above, namely if the system is slowly mixing, long trajectories are needed to learn the transfer operator and extract good eigenfunctions.
42
9
More recently, there has been growing interest in learning the infinitesimal generator of the process [1 I 50) 20]. which allows one overcome the difficult choice of the lag-time. The statistical learning properties of generator learning have been addressed in [21], where an approach based on the resolvent operator has been proposed in order to bypass the unbounded nature of the generator. However the key difficulty of learning from biased simulations remains an open question. In this work, we prove that the infinitesimal generator is the adequate tool to deal with dynamical information from biased data. Leveraging on the statistical learning considerations in [23 211, we introduce novel procedure to compute the leading eigenpairs of the infinitesimal generator from biased dynamics, opening the doors to numerous applications in computational chemistry and beyond.
Contributions In summary, our main contributions are: 1) We introduce a principle approach, based on the resolvent of the generator, to extract dynamical properties from biased data; 2) We present a method learn the generator from a prescribed dictionary of functions; 3) We introduce a neural network loss function for learning the dictionary. with provable learning guarantees; 4) We report experiments on popular molecular dynamics benchmarks, showing that our approach outperforms saute-of-the-ann transfer operator and recent generator learning approaches in biased simulations. Remarkably, even with datasets containing only few relevant transitions due to sub-optimal biasing, our method effectively recovers crucial information about the transition mechanism.
Paper organization In Section 2 introduce the learning problem. Section explores limitations of transfer operator approaches. In Section 4 we review a recent generator learning approach [21] and adapt it nonlinear regression with finite dictionary of functions. Section 5 presents our method for learning from biased dynamics. Finally, in Section 6 we report our experimental findings.
2 Learning dynamical systems from data
In this section, we address learning stochastic dynamical systems from data. After introducing the main objects, we review existing data-driven approaches and conclude with practical challenges. We ground the discussion in the recently developed statistical learning theory, 24),24 contributing in particular to the existence of physical priors and feasibility of data acquisition for successful learning.
2

---

Stochastic differential equations (SDEs) and evolution operators While our observations in the paper naturally extend to general forms of SDEs [see e.g. simplify the exposition, we focus on the Langevin equation, which is most relevant to our discussion of biased simulations. Specifically, we consider the overdamped Langevin equation
dX,- = dx + 2 and Xu =x,
describing dynamics in a (state) space X' S IRd, governed by the potential V : R" + R at the temperature 3- T where W2 is Rd -ddiienssionl standard Brownian motion.
The SDE T admits unique strong solution
(Xt)au that a Markov process to which we can associate the semigroup of Markov transfer operators (Ti) defined, for every t 2 0, as
[Tef((>) == [Xx, = z], z z E. f: R.
For (I) the distribution of Xx converges to the invariant measure π on X' called the Boltzmann distribution, given by π (dx) x e e=)V In such cases, one can define the semigroup on L (X), and characterize the process by the infinitesimal generator
= li (Tx -I)/t
defined on the Sobolev space (2) of functions in (xt) whose gradient are also in L (X)). The transfer operator and the generator are linked one to another by the formula Te = exp(tt). Moreover, it can be shown (see Appendix A that the generator acts on f: X - R as
Cf = - (VU. 8
- Af
(3)
which, integrating by parts, gives - RVV, Vg) d== ff((gg),,, showing that L self-adjoint. If C has only a discrete spectrum, one can solve (T) by computing the spectral
decomposition
==
EiEnNdeei @ Si,
Using 22 and the exponential relation between the transfer operator and the generator, one can write [T)](>>= = F(S(X:) = x] Eie (f,)) x E E f: X + RR
(5)
where the timescales of the process appear as the inverses of the generator eigenvalues. Consequently, the eigenpairs of the generator offer valuable insight about the transitions within the studied system.
Learning from simulations The main difference underpinning the development of learning algo- rithms for the transfer operator and the generator lies in the nature of the data used. While for the transfer operator we can only observe noisy evaluation of the output to learn a compact operator, in the case of the generator, knowing the drift and diffusion coefficients allows us to compute the output, albeit at the cost of learning for the former align with vector-valued regression in function spaces (22). whereas methods for the latter, discussed in the following section, more akin to physics-iffrrmed regression algorithms. In both settings, we learn operators defined on a function (hypothesis) space, formed by the linear combinations of preseribed set of basis functions (dictionary) 2j X - R, E [m].
unbounded differential operator. Consequently, learning methods
H hn.. = [jj(m]"jj
(ui, Um)) E RRM
(6)
The choice of the dictionary, instrumental designing successful learning algorithms, may be based on prior knowledge the process learned from data [4 30 50]. The space H is naturally equipped with the geometry induced by the norm || 11 : دو u². Moreover, every operator A: H > H. can be identified with matrix A E Rm by Ah,, = =(.)TAu. v. In the following, we will refer to A and as the same object, explicitly stating the difference when necessary.
Transfer operator learning Learning the transfer operator Tr can be simply seen as the vector-valued regression problem (22) in which the action of T.: (Z() - L (X') on the domain H C L2 (X') is estimated by an operator T;: H + H. This aims to minimize the mean square error (MSE) w.r.t the invariant distribution. Given dataset D,, := (I4, y4 2(+1)}= of time-lag > 0 consecutive states from a trajectory of the process, a common approach is to minimize the regularized empirical MSE, leading to the ridge regression (RR) estimator î, C, Cr, where the empirical covariance matrices are C - 1
12 Eietn =()))=(;;) and C, Eienn = We then estimate the eigenpairs (A,, fi) in 4 by the eigenpairs (ū., ū.) of T, as x; In (P./t) and fi:=2()TTO..
3

---

We stress that transfer operator approaches crucially relies on the definition of the time-lag from which dynamics is observed. Setting this value is a delicate task, depending on the events one wants to study. If / is chosen too small, the cross-covariance matrices will be too noisy for slowly mixing processes. On the other hand, if is too large, because the relevant phenomena occur at large time scales, very long simulation is needed to compute the covariance matrices. In order to overcome this problem biased simulations can be used, which we discuss next.
Learning from biased simulations
As discussed above, in molecular dynamics, the desired physical phenomena often cannot be observed within an affordable simulation time. To address this, one solution is to modify the potential,
U'(z) = U(x) + (xx). x x X
where we assume that the introduced perturbation (a form of bias in the data) V(x) is known. For example the bias potential V may be constructed from previous system states to promote transitions to not yet visited regions. One of the prototypical examples is metadynamics [25], where V is a sum of Gaussians built on the fly in order to reduce the barrier between metastable states. However, the bias potential alters the invariant distribution [12. making challenging to recover the unbiased dynamics from biased data. Denoting the invariant measure of the perturbed process by x' and its generator by L H (X) H H (X), our principal objective is thus to:
Gather data from simulations generated by L' to learn the spectral decomposition of the unperturbed generator L.
To tackle this problem, we note that since the eigenfunctions of the generator C are also eigenfunctions of every transfer operator To c²C we address the related problem of learning the transfer operator from perturbed dynamics. Unfortunately, there an inherent difficulty doing so. While one typically knows the perturbation in the generator, that is L' = L + (VV, (-)),, this knowledge is not easily transferred to the perturbation of the transfer operator. Indeed, recalling that T == Ti e and since the differential operator (VV, V-)) in general does share the same eigenstructure of L. one has that
== -C-(VVN))) * (T-(VVVV))
Simply put, the generator depends linearly the bias, while the transfer operator does not. One strategy to overcome the data distribution change, is to adapt the notion of the risk. To discuss this idea, recall that the invariant distribution of overdampted Langevin dynamics the Boltzmann distribution defined by the potential. Hence, we have that
T (dxx e- (w) n 7'xx)
- (T--Bu(()d) and dm
(x) dr
(v)
(x) dm f
(U-
(')(()
where the last term the Radon-Nikodye Consequently, we can express the covariance operators for the unperturbed process as weighted expectations of the perturbed data features
derivative, which exposes the distriistribttion change.
C C=Ex (X'))(X')x(XX)))
(8)
However, since the transition kernel of the process (X))r20 generated by L' is different from that of the original process, the above reasoning does not hold for the cross-covariance matrix, that
C, :=E=xo ddx' (Xo) 2 2(Xx) * Ex Ldwr (Xx) 2 =(X,)= =: C,,
Consequently, the estimator T, obtained by minimizing the reweighed risk functional R' (T)))
does not minimize the true risk since R' RT(T.) + R(T.). Despite this difference, whenever the perturbation is small or controlled and the time-lag t small enough, estimating the true transfer operator of the process from the perturbed dynamics via reweighed operators has been systematically used as the state-of-the art approach in the field of atomistic simulations 7 31 49 The (limited) success of such approaches is based on a delicate balance small enough aime and biased potential, since for small > 0 one can approximate C, by C₁ and minimize R'(T) xe R(T') over T: H. + H.
dr (Xp) Ex dge'
covariance/croos-covrrii

---

4 Infinitesimal generator learning
In this section, we address generator learning. While there has been significant progress on this topic [24 35 50 20), we follow the recent approach in [2I] for learning the generator L on an a priori fixed hypothesis space through its resolvent. Leveraging on its strong statistically guarantees, we adapt it from kernel regression to nonlinear regression over dictionary of basis functions, setting the stage for the development of our deep-learning method.
While transfer operator learning does not require any prior knowledge of the system's 's drift and diffusion, making use of this information helps learning the generator and avoids the need for setting the time lag parameter. We briefly discuss how to achieve this for over-damped Langevin processes when the constant diffusion term is known. We estimate the generator indirectly via its resolvent (n/ C)- where > 0 is a prescribed parameter. To this end, we observe that the action of the resolvent in H. can be expressed as (t,))))) = (x)xx)Tu, where Xny is the embedding of the resolvent (nI-C)) into H. given by Xxx(2) So E[ nt [x===]dd dt x X, see [1T]. We then aim to approximate Xx,(x) = G (z) by matrix G e Rm*m, Unfortunately the embedding of the resolvent is not known in close form. To overcome this, we contrast the resolvent by defining a regularized
energy kernel E (x) x H (x) - R. given by ef[f.g = E (nf(2)g(a)
which using (3) becomes
Cf gf
(9)
and, due to the identity f f = - -
also
[_ (x)))))))))))))))))))))))))))()) E (x)
ke
(10)
Since L is negative semi-definite, the above kernel induces the regularized squared energy norm eg: : (X') + +o0) by Cf[f]= := @flf. fl = Egnu (x) [x)
It counteracts the resolvent and balances the transient dynamics (energy) of the process with the invariant distribution
т. In nutshell, instead of using the mean square error of fxx) := [xx,(x) - G GT2(x)) to define the
risk, we "fight fire with fire" and penalize the energy to formulate the generator regression problem
min G: H-HH
Ro(G) = RG(G) = E [1lx,())) -))))]]
(11)
Indeed, this risk overcomes the difficulty of not knowing Xn- To show this, let us define the space (2) == {f E H2)( [fl < )o} associated to the energy norm || == V and recalling that the operator G: H + H is identified with a matrix G E Rmxm via Gh,, = (Gu). detine the (injection) operator Z: Rm + H2 by Zu T(, for every RTT. Then, since HS (RTh HD) = HS (H, H2), the norm is the sum of squared H norm over the standard basis in Rm, and one obtains
Ro(G)= = [1 -) Gl
= 11) GlHiSHHA + ||TT [[(T-PP)) Px) HS(
(12)
projected problem
represcntation error p(H))
where PH is the orthogonal projector in HC(()) onto H. In learning theory p(H) is known as the approximation error of the hypothesis space H [see e.g. e.g.43) 43). While this error may vanish for infiniteddimensione spaces, when H is finite dimensional, controlling p(H)) is crucial to achieving statistical consistency. This can be accomplished by minimizing (TT), which equivalent to
min G: H-H
| - = 2) Z Zn(nI ZG G (S)))"H))
(13)
where (.)$ is the Moore-Penrose" pscudoinverse. Using the covariance matrices
(nI-C) C= (E)nn=((()((((())))0)))
and -==---- (E #,]) (.emm)) (14)
w.r.t. the invariant distribution and energy, respectively, gives the ridge regularized (RR) solution G = (W + >1) CC O. The induced RR estimator of the resolvent, G H + H. is given, for every hú E H. by h. := Z(W + -1)-1Cu = z(:)TWW +WW++0 71) -IC- and can be estimated given data from by replacing expectation and the energy in [44 with their empirical counterparts.
5

---

Unbiased learning of the infinitesimal generator from biased simulations
In this section, we present the main contributions of this work: approximating the leading eigenfunc- tions (corresponding to the slowest time scales) of the infinitesimal generator from biased dhta.While While the general pipeline for the method can be found in figure in the following, we first address re- gressing the generator on an a priori fixed hypothesis space H. Then we introduce our deep-learning method to either build suitable space H. or even directly learn the eigenfunctions.
(di(z))
(Kjddzll
Choose intecutar tystem study
Bun biss simulation with auboptimal CV
Learn representaaton
Regueex ganieaaoor
Compute esgenpairs
Figure 1: Pipeline of our method: from biased simulations to timescales and metastable states.
Unbiasing generator regression Whenever π is absolutely continuous w.r.t. T', the regularized energy kernel (9) satisfies the simple identity
Cf E"
der: [5 di
dr de'
f.geH2?
(15)
which, recalling the rightmost equation in 7 implies that when the bias V and the diffusion coefficient& 3 are known, the energy kernel can be empirically estimated through samples from x'
via (TO). Moreover, when the potential U is known too, we can use ( Now, leveraging on (15) we
directly obtain that
R.a = R2(G) En (x)x ((xx)))))) 2x)))l2 da (x) S Ky KV C (xx, - -G= =()))l2 (16)
du
where Ky = ess sup goot de (z), which recalling is finite whenever the bias V essentially bounded. Therefore, sharp contrast transfer operator learning, whenever the true embedding (z) can be estimated, one can derive principled estimators of the true generator L's dominant eigenpairs from the biased dynamics generated by L' This is established by the following proposition, the proof of which is presented in Appendix B
Theorem 1. Let D. (z()iejn) be the biased dataset generated.froe n'. Let w(x) e Vvxx and define the empirical covariances w.r.t the empirical distribution ='=n"* Eieini dm; by
C = ==(E<>>> (1)))==))((")) (.j[mm]
and W -EE [Vizi, v (wzz)) [.j[m]]
(17)
Compute the eigenpairs (V,, Uz)ie(m) of the RR estimator G, (W 11:1)-1C. and estimate the eigenpairs in 4 (Ne, f.) (n 1/v1, 1/v,, v.). If the elements of H. and their gradients are essentially bounded, and lim p(H)) then for every
there exist (m, n, 7) EN N X.R.. R.,
m2-000
such that, for every E [m]]. AA- 2 < and sin,3 (4(f), f.)) fi)) SE, with high probability.
Note that, due the form of the estimator, the normalizing constant, f ww(x)dr does not need be computed. Moreover, relying on the upper bound in (66) we can alternatively compute and without the weights and still ensure that the above result holds true.
Neural network based learning Theorem T guarantees successful estimation of the eigendecom- position of the generator in (4) whenever the bbaey representation error p(H)) in (12) is controlled. therefore natural to minimize p(H) by choosing an appropriate basis function 2; Inspired by the recent work (24), we parameterize them by neural network, and optimize them to span the leading invariant subspace of the generator.
Let - = =2( (e m X' - Rm be neural network (NN) embedding parameterized by 0 O weights with continuously differentiable activation functions, and let X, i [m], be real non-positive

---

(trainable) weights. We propose to optimize the NN to find the slowest time-scales X that solve the eigenvalue equation Czς = Xx iE [m]. Letting Zq: R + (x) be the (pa- rameterized) injection operator, given, for every u E RR by Zzu= -iemm 2fu,, and denot- ing Ag = (nI- diag(A?...... (.....,,))))+ the eigenvalue equations for the resolvent then become (nI L) Zn ZgAz. In other words, aim to find the best rank-m decomposition of resol- vent (nI - 1 = ZoAg Therefore, for some hyperparameter 0 we introduce the loss
E, (0) = 1 - -) 2 28llss HiSHH)) 11 -) (SHH)) +o
Σ (1e) (9)=)))2 81. 4.jE(m]
While the first term measures the approximation error in the energy space, it cannot be used as a loss, because the action of the resolvent is not known. To mitigate this, the second term is introduced, under the assumption that (1yT - --)- E HS (HC) (2)) (see Appendix for a discussion). The third term is optional; specifically, if the goal is not only to identify the proper invariant subspace of the generator (a = 0), but also to optimize the neural network to extract eigenfunctions as features, then this last term (a 0) encourages the orthonormality of features in L (X), an idea successfully exploited in machine learning and computational chemistry [see e.g. 24 and references therein].
Recalling (14 and denoting by Ca and Wo the covariance matrices associated to the parameterized features, after some algebra, we obtain that
(a = tr [CoA]WoA7 2CoA} + (Co - -)
(18)
In turn, this can be estimated from biased data by two independent samples 74 and 732 as
f (0)= (f)= CA 2 C'A w C -110 (19)
where C and We a are the empirical covariances given by (7) for distribution T%e while @* Eo("), kE[2]. Importantly, the computational complexity of the loss (19) is of the order O(nm²d), where d is the state dimension and the sample size, however it can be reduced to (nmd) (see Appendix C) allowing its application to learn large dictionaries for high-dimensional problems with big amounts of (biased) data.
The following result, linked to controlling of the representation error as detailed in Theorem | pro- vides theoretical guarantees for our approach. The proof and discussion are provided in Appendix B
Theorem 2. Given compact operator
(nT- 17 O, (2") de(m)) S HT ( all EE then = - -ie(( (m-
E E
EhA (0)
for all O,
(20)
where w- Egwor [w(x))] Moreover. if 0 and Xm+2 Am+11, then the equality holds if and only if (X 25+)= = (d., fi) T-..e. up to the ordering of indices and choice of eigenfinction signs for i E [m].
This theorem provides a justification for minimizing the loss in (19). which can be achieved by stochastic optimization algorithms, to obtain approximation of either the leading invariant subspace of the resolvent (17 - -)- (without orthonormality loss, i.e. a = 0), on which the estimator in Theorem can be computed, or even the individual eigenpairs (a > 0). A pseudocode of our method is provided below. The main advantage of this method is that it exploits the knowledge of the process. namely, if only the bias V and the diffusion coefficient& B are known, recalling [0)), the computation of loss relies just of the gradient of the features. On the other hand, the knowledge of the potential can also be exploited via (9). Finally, even if the neural network features are not perfectly learned, one can still resort to Theorem to compute the approximate eigendecomposition of L.
6 Experiments
In this section, we test the method described above on well-established [4 32 36 molecular dynamics benchmarks, featuring biased simulations of increasing complexity. We first start by showing the efficiency of our method on a simple one dimensional double well potential. We then proceed to the Huller-Brown potential which is 2D potential, where this time, sampling is
7

---

Algorithm From biased to unbiased dynamics via infinitesimal generator
1: Parameters 17 > 0 shift of the generator, m number of wanted cigenfunctions, K number of optimization steps, 7 > 0 and a > ( regression and NN hyperparameters
2: Inputs Dataset D, (at)) [e(] gathered from a simulation with bias potential V
3: Compute weights wr(ar))
(xx)) E [n)] be used in line 7
4: if the dictionary of function a does not already exist then
5: Initialization: randomly initialize neural networks weights of A and (z)) set ls 0
6: while ls < K do
7: Compute C2 and W j 1,2, using (77 for two independent batches f// and ft'
8: Compute loss (0) using (99) and backpropagate
9: end while
10: end if
11: Compute C and W using (17 the datatset Dr
12: Compute the eigenpairs (x),ejm) of Gm==(W++mp))
13: Output Estimated eigenpairs of L are (2x, (1)-111/ (-)es) [m]]
accelerated by a bias potential built on the fly. Finally, we study the conformational landscape of alanine dipeptide. This small molecule is a classical testing ground for rare event methods. To showcase the efficiency of our method we analyse two different sets of data both generated in a metadynamics-lie like approach and showcase the efficiency of our approach, even with small number of transitions in the training set. The codes used train the models can be found in the following repository: http://gghuub.oomDeeeeegg
imothee/Genl.enm
One dimensional double well potential We first showcase the efficiency of our method simple one dimensional toy model. We sample transitions from U + V, where U is a double well potential and V is a bias potential. The results shown Figure in the appendíx, where our method clearly outperforms transfer operator approaches and recovers the true underlying dynamics.
Muller Brown potential with metadynamics biasing Muller Brown is 2 dimensional potential presenting metastable states often studied in the context of enhanced sampling [99 37 14 9.
presents two minima, with one of them separated into two sub-basins. We thus expect two relevant eigenpairs: the slowest one corresponding to the transition between the two basins and the second slowest one describing the transition between the two sub-basins. However, at low temperature crossing the barrier occurs rarely. To expedite the rate of transition we use metadynamics and instead of having a predefined bias potential, as in the previous section, the bias is built on the fly using metadynamics 233). The results of the training procedure are presented in Figure 2 We compare the results with TeCp and a state of the art generator learning approach in (50). From this figure, we see that we managed to accurately learn the dynamical behavior of the system despite the fact that the dynamics was performed using bias potential. As expected, it is elearly outperforming transfer operator approaches. We achieve similar or slightly better results (particularly near the transition state) on the qualitative shape of the eigenfunctions. On the other hand, our method performs better than previous work on generator learning on the estimation of eigenvalues, and is the closest to the ground truth eigenvalues. This is likely to be due to the fact that the method in [50] requires the tuning of hyperparameters in the loss function, while in our case, these coefficients are trainable. It should be noted that here, the eigenfunctions were fitted with well-learned features. However, we present in the appendix results where the features are not perfectly learned, but we still manage to recover the eigenfunctions.
Alanine dipeptide with OPES biasing We next treat a more concrete molecular dynamics example with the study of the conformational change of alanine dipeptide in gas phase. It is molecule containing 22 atoms, of which 10 are heavy. For the remaining of this study, we will only take into account the positions of the heavy atoms, making it 30 dimensional system. This molecule has widely been used to test methods in enhanced sampling 99 30 421: it presents a conformational change which is rare event described by the angles Φ and v.. In the studies made on this system, the angle & has been shown to be good CV: the transition between the two states is very well described, and thus a bias potential can easily be built with this CV. On the other hand, the angle & misses most of the transition and is non optimal CV. We generated biased datasct using variation of
8

---

Ground truth
Generator EF (this work)
Zhang et al. (2022)
deepTICA
2.0
1.5
#A=1310
1.0
0.5
15
0
A=1.13
10
0.5
0.0
-1
Figure 2: Muller Brown potential. Comparison of the ground truth two first relevant eigenfunctions of the potential (first column) with this work (second column), transfer operator approach deepTICA (7) (third column) and the work of Zhang et al. (50) (fourth column). x and axis are the coordinates of the system and points are colored according to the value of the eigenfunction. The underlying potential is represented by the level lines in white. Associated eigenvalues À are also reported.
metadynamics called the fly probability enhanced sampling (OPES) [16), which allows a more extensive and faster exploration of the state space than metadynamics 255
Dataset 1: 800ns simulation, biasing the s dihedral angle, with OPES leading to few transitions
between the two states. The bias potential was built during the first 100ns of the simulation. For the remaining 700ns, the potential built during the first part was kept fixed to enhance transitions. Dataset 2: 50ns simulation, biasing on ø dihedral angle, with OPES leading to many transitions between the states. The bias potential was built during the first 20ns of the simulation. For the remaining 30ns, the potential built during the first part was kept to enhanced transitions.
Dataset mimics situations where one has only a basic prior knowledge of the system: only "suboptimal" CV used yielding to only few transitions between the metastable states within the affordable simulation time. In order to ensure translational and rotational invariant vectors, we
use Kabsch [88 algorithm, which has been used previous studies [50 4 I0) to transform positions of the atoms. The results presented in Figure Panels a) and b) display the first and second eigenfunctions learned by our method respectively. Notice that, even though only transitions present in dataset 1, the first eigenfunction separates the two metastable states, and the second identifies a faster transition in one metastable state. Panel c) showcases the good out-of-sample generalization ability of the method. It visualizes the first eigenfunction obtained above, but this time visualized on points from dataset 2 and in the plane of dihedral angles ø and 8. Interestingly, we discover that a linear relationship is present in the transition region, in agreement with recent findings in the molecular dynamics literature 19
To further improve the description of the transition and to enhance the training set without any prior knowledge of the mechanism, one could perform biased simulations using the first eigenfunction. Nonetheless, this is not the scope of this paper. To push our method further and see its capabilities when training on a good dataset, we trained on Dataset One key quantity in molecular dynamics the committor function for metastable states A and B, which is defined the probability of, starting from A, going to B before going back to A. Theory tells us that the committor is linearly related to the first eigenfunction of the generator, result going back to Kolmogorov [see[5] for a discussion]. This relation is exposed in panel d) of Figure 3 when comparing to the committor model obtained in 19 indicating the good performance of our method.
Chignolin miniprotein In this section, we report the results our method obtained on a larger scale experiment: the folding/umfolding mechanism of the chignolin miniprotein. This system has

---

First eigenfunction
Second elgenfunction
c) First eigenfunction
Committor correlation
0.75
10
0.00
-0.5
-0.50
02
-0.75
00
02
0.6 0.8
9
Figure 3: Alanine Dipeptide. Results of our method trained on Dataset 1 a) and b) first and second eigenfunctions represented on dataset 1, in the plane of the and dihedral angles. e) first eigenfunction represented on dataset 2. in the plane of the & and & dihedral angles, indicating that our method is effective even when trained from poor CVs (see for more discussion). On all three panels, points are colored according to the value of the cigenfunction. d) Comparison of our method with the committor (g) of
We first performed a jus biased simulation using the to gather transitions. Then we chose descriptors as input of the neural networks that are known to describe well the folding process (7]. Finally, we trained the method described in the current work with this trajectory and compared it with the results obtained when training on a 106us unbiased trajectory provided by D.E. Shaw research (28). The results are
extensively been studied
19
39
deep-TDD collective variable 46
37
presented in figure 4 showing very good agreement between the training on unbiased trajectory and on a biased one.
a) First eigenfunction (DE Shaw)
b) First eigenfunction (biased data)
1.75
0/ASP3N-THRB9 (HRB0)
1.50
(O8HNL-NESS
1.25
1.00
0.50
0.25
0.25
0.50
0.75
1.00
1.25
1.50
1.75
0.25
0.50
0.75
1.00
1.25
1.50
1.75
0VASP3N-GLY70)
d(ASP3N-GLYTO)
Figure 4: Our method for the chignolin miniprotein. The data points are represented in the plane of the distance between the nitrogen atom of the residue 3: ASP (ASP3N) and the oxygen atom of the residue 7: Gly (O)y70) and the distance between ASP3N and the oxygen atom of residue THR (THR8) which allow visualizing the folded and unfolded states.
7 Conclusions
We presented method learn the eigenfunctions and eigenvalues of the generator of Langevin dynamics from biased simulations, with strong theoretical guarantees. We contrasted this approach with those based on the transfer operator and recent generator learning approach based on Rayleigh quotients. In experiments, we observed that our approach is effective even when trained from sub- optimal biased simulations. In the future our method could be applied to larger-scale simulations to discover rare events such as poteen--ligand binding or catalytic processes. A main limitation of our method that, in its current form, formulated for homooeooggss bias potentials. However, the proposed framework could be naturally extended to dime-dependent biasing, broadening its applicability in computational chemistry. Furthermore, given the quality of our results on alanine dipeptide, in the future, we can use our method to compute accurate eigenfunctions from old, possibly poorly converged, metadynamics simulations, thereby gaining novel and more accurate physical information.
10

---

Acknowledgements
This work was partially funded by the European Union NextGenerationEU initiative and the Italian National Recovery and Resilience Plan (PNRR) from the Ministry of University and Research (MUR). under Project PE0000013 CUP J53C22003010006 "Future Artificial Intelligence Research (FAIR)". We thank D.E Shaw research for providing the chignolin trajectory
References
[1] A Cabannes, V. and Bach, F. (2024). The Galerkin method beats graph-based approaches for spectral algorithms. Proceedings of The 27th International Conference on Artificial Intelligence and Statistics, volume 238, pages 45 -499. PMLR.
[2] Abraham, M. J., Murtola, T., Schulz, R., Páll, S., Smith, J. C., Hess, B., and Lindahl, E. (2015). Gromacs: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX, 1-2:19-25.
[3] Bakry, D., Gentil, I., and Ledoux, M. (2014). Analysis and Geometry of Markov Diffusion Operators. Springer.
[4] Belkacemi, Z., Gkeka, P., Lelièvre, T., and Stoltz, G. (2022). Chasing collective variables using autoencoders and biased trajectories. Journal of Chemical Theory and Computation, 18(1):59-775
[5] Berezhkovskii, A. and Szabo, A. (2004). Ensemble of transition states for two-state protein folding from the eigenvectors of rate matrices. The Journal of Chemical Physics, 121(18)99188 186- 9187.
[6] Bolhuis, P. G., Dellago, C., and Chandler, D. (2000). Reaction coordinates of biomolecular isomerization. Proceedings of the National Academy of Sciences, 97(11):5877-5882
[7] Bonati, L.,, Piccini, G., and Parrinello, M. (2021). Deep learning the slow modes for rare events sampling. Proceedings of the National Academy of Sciences, 118(44)e211333311
[8] Bonati, I.., Trizio, E.,, Rizzi, A., and Parrinello, M. (2023). unified framework for machine learning collective variables for enhanced sampling simulations: mlcolvar. The Journal of Chemical Physics, 159(1).
[9] Chen, H., Roux, B., and Chipot, C. (2023). Discovering reaction pathways, slow variables, and committor probabilities with machine learning. Journal of Chemical Theory and Computation, 19(14):4414 4426.
[10] Chen, W., Tan, A. R., and Ferguson, A. L. (2018). Collective variable discovery and enhanced sampling using autoencoders: Innovations in network architecture and error function design. The Journal of Chemical Physics, 149(7):072312.
[11] Comer, J., Gumbart, J. C.,, Hénin, J., Lelièvre, T., Pohorille, A., and Chipot, C. (2015). The adaptive biasing force method: Everything you always wanted to know but were afraid ask. The Journal of Physical Chemistry B, 119(3):1 129-1151.
[12) Dama, J. F., Parrinello, M., and Voth, G. A. (2014). Well-tempered metadynamics converges asymptotically. Phys. Rev. Lett., 112:240602.
[13] Davis, C. and Kahan, W. M. (1970). The rotation of eigenvectors by perturbation. iii. SIAM Journal on Numerical Analysis, 7(1):1-46.
[14] France-Lanord, A., Vroylandt, H., Salanne, M., Rotenberg, B., Saitta, A. M., and Pietrucci, F. (2024). Data-driven path collective variables. Journa of Chemical Theory and Computation, 20(8):306930844
[15] Hou, B., Sanjari, S., Dahlin, N., Bose, S., and Vaidya, U. (2023). Sparse learning of dynamical systems in RKHS: An operator- approach. In Proceedings of the 40th International Conference Machine Learning. volume 202, pages
PMLR

---

[16] Invernizzi, M. and Parrinello, M. (2020). Rethinking metadynamics: From bias potentials to probability distributions. The Journal of Physical Chemistry Letters, 11(7):2731-2736.
[17] Invernizzi, M. and Parrinello, M. (2022). Exploration vs convergence speed in adaptive-biat enhanced sampling. Journal of Chemical Theory and Computation, 18(6):3988-3996.
(18) Kabsch, W. (1976). A solution for the best rotation to relate two sets of vectors. Acta Crystallographica Section A, 32(5):922-923.
[19] Kang. P. Trizio, E., and Parrinello, M. (2024). Computing the committor with the committor to study the transition state ensemble. Nature Computational Science, 4(6):451-460.
[20] Klus, S., Nüske, F., Peitz, S., Niemann, J.-H., Clementi, C., and Schütte, C. (2020). Data-driven approximation of the koopman generator: Model reduction, system identification, and control. Physica D: Nonlinear Phenomena, 1006:13246.
[21] Kostic, V., Lounici, K., Halconruy, H., Devergne, T., and Pontil, M. (2024a). Learning the infinitesimal generator stochastic diffusion processes. Preprint, arXiv:2405. 12940.
[22] Kostic, V., Novelli, P., Maurer, A., Ciliberto, C., Rosasco, L.,, and Pontil, M. (2022). Learning dynamical systems via Koopman operator regression in reproducing kernel Hilbert spaces. In Advances in Neural Information Processing Systems.
[23] Kostic, V. R., Lounici, K., Novelli, P., and Pontil, M. (2023). Sharp spectral rates for koopman operator learning. Thirty-seventh Conference on Neural Information Processing Systems.
[24] Kostic, V. R., Novelli, P., Grazzi, R., Lounici, K., and Pontil, M. (2024b). Deep projection networks for learning time- homooenogg dynamical systems. In The Twelfth International Conference on Learning Representations.
[25] Laio, and Parrinello, M. (2002). Escaping free-energy minima. Proceedings of the National Academy of Sciences, (2520)):2222125555
[26] Leitold, C., Mundy, C. J., Baer, M. D., Schenter, G. K., and Peters, B. (2020). Solvent reaction coordinate for an sn2 reaction. Journal of Chemical Physics, 153(2).
[27] Lelièvre, T., Pigeon, T., Stoltz, G., and Zhang, (2024). Analyzing multimodal probability measures with autoencoders. The Journal of Physical Chemistry B, 128(11):2007-2631
[28] Lindorff-Larsen, K. Piana, S., Dror, R. O., and Shaw, D. E. (2011). How fast-folding proteins fold. Science, 334(6055):517-522
[29] Magrino, T., Huet, L.,, Saitta, A. M., and Pietrucci, F. (2022). Critical assessment of data-driven versus heuristic reaction coordinates in solution chemistry. The Journal of Physical Chemistry A, 126(47)88887-800)
[30] Mardt, A., Pasquali, I.., Wu, H., and Noé, (2018). VAMPnets for deep learning of molecular kinetics. Nature Communications, 9(1).
[31] McCarty, J. and Parrinello, M. (2017). A variational conformational dynamics approach to the selection of collective variables in metadynamics. The Journal of Chemical Physics, 147(20):204109.
[32] Mori, Y., Okazaki, K.-i., Mori, T., Kim, K., and Matubayasi, N. (2020). Learning reaction coordinates via cross-entropy minimization: Application to alanine dipeptide. The Journal of Chemical Physics, 053(5):05115.
[33] Oksendal, B. (2013). Stochastic Differential Equations: an Introduction with Applications. Springer Science & Business Media.
[34] Pietrucci, F. (2017). Strategies for the exploration of free energy landscapes: Unity in diversity and challenges ahead. Reviews in Physics, 2:32 45.
[35] Pillaud- Pillaud-Vivien, I.. and Bach, F. (2023). Kernelized diffusion maps. The Thirty Sixth Annual Conference Learning Theory, pages 5236-5259. PMLR.
12

---

[36] Plainer, M.. Stark, H., Bunne. C., and Günnemann, S. (2023). Transition path sampling with Boltzmann generator- based MCMC moves. NeurIPS 2023 Generative Al and Biology (GenBio) Workshop.
[37] Ray, D., Trizio, E., and Parrinello, M. (2023). Deep learning collective variables from transition path ensemble. The Journal of Chemical Physics, 158(20):204102.
[38] Reed, M. and Simon, B. (1980). Functional Analysis (Modern Mathematical Physics, Volume
1). Academic Press.
[39] Rydzewski. Chen, M., Ghosh, K., and Valsson, O. (2022). Rewcighted manifold learning of collective variables from enhanced sampling simulations. Journal of Chemical Theory and Computation, 18(12)7117--112. PMID: 36367826.
[40] Salomon-Ferrer, R., Case, D. A., and Walker, R. C. (2013). overview of the amber biomolecular simulation package. WIREs Computational Molecular Science, 3(2):198-21))
[41] Schwantes, C. R. and Pande, V. (2015). Modeling molecular kinetics with tica and the kernel trick. Journal of Chemical Theory and Computation, 11(2):600 608.
[42] Shmilovich, K. and Ferguson, A. L.. (2023). Girsanov reweighting enhanced sampling technique (grest): On-the-fly data-driven discovery of and enhanced sampling in slow collective variables. The Journal of Physical Chemistry A, (7)349334773331
[43]
Steinwart, I. and Christmann, A. (2008). Support Vector Machines. Springer New York.
[44] Torrie, G. and Valleau, J. (1977). Nonphysical sampling distributions in monte carlo free-energy estimation: Umbrella sampling. Journal of Computational Physics, 23(2):187-199
[45] Tribello, G. A., Bonomi, M., Branduardi, D.,, Camilloni, C., and Bussi, G. (2014). Plumed New feathers for an old bird. Computer Physics Communications, (85(2):604-613.
[46] Trizio, E. and Parrinello, M. (2021). From enhanced sampling reaction profiles. The Journal of Physical Chemistry Letters, 12(35):8211-6266 PMID: 34469175.
[47] Tuckerman, M. E. (2023). Statistical Mechanics: Theory and Molecular Simulation. Oxford university press.
[48] Fijnden, (2006). Transition Path Theory, pages 453 493. Springer Berlin Heidelberg, Berlin, Heidelberg.
[49] Yang, Y. I. and Parrinello, M. (2018). Refining collective coordinates and improving free energy representation variational enhanced sampling. Journal of Chemical Theory and Computation, 14(6):2889-2994
[50] Zhang, W., Li, T., and Schüitte, C. (2022). Solving eigenvalue pdes of metastable diffusion processes using artificial neural networks. Journal of Computational Physics, 11377.
13

---

Appendix
The appendix contains additional background on stochastic differential equations, proofs of the results omitted in the main body, and more information about our learning method and the numerical experiments.
notation
meaning
notation
mcaning
[-]
set (1.2.....)
[+]+
nonnegative part of a number
X'
state space of the Markov process
(X4)
time-homogeneous Markov process
dW,
Brownian motion
inverse temperature of the system bias potential
U
potential energy of the system
V
x
Boltzmann distribution of potential U
n'
Boltzmann distribution of potential U V+
empirical version of ग density of x w.r.t n'
H'
empirical version of n²'
di dr/dxt /dm
w
exponential weights of drr / /nr' Sobolev space w.r.t on X
L2 (X))
L space on X' w.r.t the measure transfer operator with lag time t generator the true process shift parameter
Te
T.
empirical estimator of generator the biased process
C
L'
77
(nl (2)-
resolvent at basis functions
7
hypothetical Hilbert space
2
[1]---]
regularized energy kemnel
e"[-]
regularized energy norm
723 (X'
space associated to the energy norm
Z
injection operator
Ra
risk functional
Ra
empirical risk functional
C
covariance
Č
empirical covariance
W
covariance in energy space
W
empirical covariance in energy space
A
generator eigenvalue
f
generator eigenfunction
X
empirical estimator of A neural network embedding
f.
empirical estimator of
z"
0
neural network weights
Zo
neural network injection operator
neural network weights
E,
regularized loss function
E.
regularized empirical loss function
Table 1: Summary of used notations.
A
Background
In this work we consider stochastic differential equations (SDE) of the form
dxx,- = a(X,)dl + b(X,)dW, and X = x.
(21)
The special case of Langevin equation considered in the main body of the paper corresponds to a(z) = (z) and b(z) / Equation (2) describes the dynamics of the random vector X, in the state space X' S Rd, governed by the drift a Rd - Rd and the diffussion b Rd -> Rdxp coefficients, where W, is IRP oP-diimenioonal standard Brownian motion. Under the usual conditions [see e.g. 33) that a and b are globally Lispchitz and sub-linear, the SDE (2I) admits an unique strong solution X= = (X+)20 that is Markov process to which we can associate the semigroup of Markov transfer operators (Ti) defined, for every 2 0, as
[Tef](>) == E[f(X) = z], Z:
(22)
For stable processes, the distribution of X, converges to invariant measure X, such that Xy ~ π implies that X₂ ~ π for all t > O. In such cases, one can define the semigroup on L (X),, and characterize the process by the infinitesimal generator of the semi-group (Te)
T 1
C = lim
(23)
L
14

---

defined on the Sobolev space (X)) of functions in (A)) whose gradient are in (A)), too, i.c.
C: L3 (X') + (X') and dom -H (H,((X), The transfer operator and the generator are linked to each other by the formula T, exp(12).
After defining the infinitesimal generator for Markov processes by for solution processes of equations like (21). Given smooth function f E (2', R), Itö's formula [see for instance 3 p. 495] provides for E R+,
we provide its explicit form
L.
f(X,) - ((Xo)
Σ 0 O.f(X.) + Σ xx W+=1
I VS(X,) dX, + ( [=Xx](()))Xx,,,,]]
Recalling (11. we get
[(x.)= /''
+ % VV(X))6(X.))WW
(24)
Provided f "and b are smooth enough, the expectation of the last stochastic integral vanishes so that we get
[[f(x,) = =] = f(x) + /
Xo x ds
Recalling that
lim (Tef ff/t,, we get for every X,
r00
([f(x, Xo = x] fxx)
Cf(x) = lim
t-+0
Iim t-0
/
do Xo = x
- [V(())
=V((=)'a(z) -
(25)
which provides the closed formula for the IG associated with the solution process of (1 In partieular, for Langevin dynamics this reduces to (3)
Next, recalling that for a bounded linear operator A some Hilbert space H the resolvent set of the operator A defined as p(A) = (A E C A >/ is bijective). and its spectrum Sp(A) = CV C))((()) let A Sp(AA) be the isolated part of the spectra, i.c. both À and μ Sp(A) \ X are closed in Sp(A). Then, the Riesz spectral projector PA : H + H is defined by
2x /. (=I d2.
PA
(26)
where I' is any contour in the resolvent set Res (AA) with À in its interior and separating À from μ.. Indeed, we have that P2 PA and H. (m() @ Ker(P,) where Im(P,) and Ker( Px) are both invariant under A, and we have Sp( Ap(AAmcp)) = À and Sp( (A ( = Moreover, PA + PA =
and PAP, = P,PPa = 0.
Finally if A is compact operator, then the Riesz-Schauder theorem [see e.g. 38 assures that Sp(T) a discrete set having no limit points except possibly = 0. Moreover, for any nonzero À E Sp(T". then À is eigenvalue (i.e. it belongs the point spectrum) of finite multiplicity, and, hence, can deduce the spectral decomposition in the form
A
Σ A AES)(A)
(27)
where the geometric multiplicity of X, rλ = rank (A),, is bounded by the algebraic multiplicity of X. If additionally A is a normal operator, i.e. AA" A, then PA Pr is orthogonal projector for each À E Sp(AA) and PA = Σnz Wi @ th, where eh: are normalized eigenfunctions of A corresponding to À and rx is both algebraic and geometric multiplicity of X.
We conclude this section by stating the well-known Davis-Kahan perturbation bound for eigenfunc- tions of self-adjoint compact operators.
15

---

Proposition ((13)). Let A be compact self-adjoint operator a separable Hilbert space H. Given a pair (f, C x H such that [[ 1, let X be the eigenvalue of A that is closest to f and let f be its normalized eigenfinction. g min A/l Sp(A) \ {A}} 0, then sin(<a(,.)) f)) & | 1/1/0
B Unbiased generator regression
In this section, we prove Theorem relying on recently developed statistical theory of generator learning [21]. To that end, let 6(I) 2(x) H. be feature map of the RKHS space H of dimension dim(H) = m. Let Way = ofu, , the eigenvalue decomposition bi = uu/ This induces the SVD of the injection operator, Zu; = oz² for 2j = 2.)
and let
Since H C H2 (X), we have that
Cy sup > (i) 2,,001,(()00))) (x)([ū]](({[ < to0,
2
and, denoting W = W + 17y1
I[W= S m.
In addition denote the empirical version of W, as W,, == W 17y.
Now, we can apply the following propositions from 21 to our setting, recalling the notation for normalizing constant == [w=]] for which w( w(-)/w dx dz/dxt /dx'.
Proposition 2. Given > with probability in the ii.d draw of (s:) n, it holds that
P( - -WI S En)( 2 21
where
e, 2|
2 2| C(5)
4 trr(W)
3n ((6) +
and L(8) log 5 [[
(28)
n
Proposition 3. Given
O, with probability in the i.i.d. draw (x;) from n, from
holds that
flIW, 1/W w)W 1/21 s (,. 8)} 2 -1-6 1,
(29)
where
20, 3n
2cr c)(1,0),
()
L (),)0) V
(30)
and
4
tr tr(w W) || , W|
(7,
In 8
Moreover.
P P [|W1/2W w 1 (,(
1
2 1-8 8.
(31)
Proposition 4. With probability in the i.i.d. draw of (x:) from æ, it holds
P | 1/2 C)IlF < S (7, 1-8
where
(r, 4 4 m W|
C8 Cr
In
(32)
n2
We are now ready prove Theorem ! which we restate here for convenience.
Theorem 1. Let D, (x())etm) be the biased dataset generated from n'. Let w(xx) (v and define the empirical covariances w.r.t the empirical distribution म' = n - -1Eei daz by
C = (o [je(m] and W = =EE (V Viezl)
(17)
16

---

Compute the eigenpairs (Vi, ki)ke(mi m of the RR estimator Gn, (W + - - and estimate the eigenpairs 4 as (Ai, Fi) = (1) 1/1; T.)) If the elements of and their gradients are essentially bounded, and lim p (L) = 0, then for every > 0, there exist (m, n, 7) ENN N R+,
such that, for every i E [m]] A A < and sin (4(fi. fi) S., with high probability.
Proof. We first show that Ra (G(G) for big enough m, E N and small enough 7 > 0.
Observe that
= = / 1I and C =[CC//
(33)
and, hence
G = W w = (E(W,))-
due to cancellation of
Given & > 0, let m E IN be such that p(H) = 11(1
(isc,, < < Next, since
we have that
= 2G -
as y - O. Hence, let y > 0 be such that
< <//3.
Finally, using the decomposition of the risk
-- R,(G S 1-00--PPrth-------2--2-lll------ (1 - -nnn)llsssee....
remains to show that for large enough n we have |z(Gn Gnn))lnsSRR S 8//3. To that end observe that
w - Gn,) w// (Č W.
- w1/2 / W 1/2 C) - (C-) -) 1 -W W/
Thus, by multiplying the above expression by w and applying Propositions 3 and there exists n E IN such that | - G, HS( *2) S e/3.
we obtain that
Next, assuming that ||W W|| is small, for the normalization of the estimated eigenfunctions have that
10 |
I
2) Wu
o. We, oW W)u,
A (W) W W|
Am (W) 2/| w|
where we have that Am (W) > 0 due fact that (2z) are linearly independent.
Therefore, conclude the proof, we apply 21 Proposition which directly relying on Proposition
I yields the result.
At last we remark, based on the observation that W = E[Wl/w and C E(C)/w, one can readily obtain stronger version of Theorem in the general RKHS setting of [21.
C Unbiased deep learning of spectral features
In this section, provide details DNN method and prove Theorem 2 To that end, let us denote the terms in the loss as
(0)= [[ - C) - 2 Il 11n C) 1| HS(())) + ++ 2
29)) B1,)2, (34)
expected loss E
orthonormality loss Eon

---

and recall that the injection operator is (g=) RT + (x)).. It is easy to show from the definition of the adjoint that, for every f HT we have
/ = Exna -0 (c)(mn1-))())
(35)
Thus, 2g Eaww (x) (x))a.smmm is denoted by Co-
(² 29) [m)] which we denote by Wo. while 1z(nT Zo =
Theorem 2. Given compact operator (n1-2) 0, if (E S (x) (H for all E then
2r E[Ea a (0)) w² (,) - (n-A,) for all E O,
(20)
where w = En [w][) Moreover, if >>>>0> 0 and Am+ Am+ then the equality holds if and only if (A% 2/ (Ms, fi) a-a.e. up to the ordering of indices and choice of eigenfunction signs for i E [m]
Proof. Let Px: Hz(X) - (x)) be spectral projector of L corresponding to the largest eigen- values. Now, consider
-14 lisce)))
HS(H2)
Due to Eckhart- Mochhart-Yunggg Young theorem, have that for every O the best rank-m approximation of (nI--1 [ L) Pm = (n --) Pmss for > m, and holds
Ex (0)2
Σ (n x
-) Σ n1 Aa)
j=++1
As before, expanding in 8 the HS norm via the trace, we obtain =0 -2
and, hence, by Cauchy-Schwarr inequality, have that [E()) - E0 < 3ep
Σ kE(m]
e
Observing that 2f E (X),,, i.e [3[2 < oc, we conclude that, for every E O, lim (²) k4000
(0). Therefore, noting that (0) 2 E()), inequality in (20) is proven. Since the equality clearly holds for the leading eigenpairs of the generator, to prove the reverse, it suffices to recall the uniqueness result of the best rank-m estimator, which is given by (n)- .ic. (e])se(()) span the leading invariant subspace (fi)eelm) of the generator. So,
= = Σ (1- j=1
and 0, have that (f) 0, implying that (2))e(m) is orthonormal basis, and, hence Pm (nI-) i.. (2. (2])) ZoA Zz. The result follows.
To show that
E, = [0)]
we rewrite (18) to encounter the distribution change, noting that the empirical covariances are reweighted but not normaliezed by w. So, we have that
()) T [[C 2w "ECC//A To Pa([[(C - w1)
But, since
Ea = (()))
=F 4 (x))((a) Eegea V (cb)a(*))),
where 2/ and 2h are two independent with law n', the proof is completed.
18

---

train loss
100
-50
500
1000
1500
2000
2500
3000
opoch
Figure Typical behavior the loss function during training.
Ground truth
Features
Fitted EF
2.0
10
1.0
1.0
0.0
0.0
0.0
-1
-1
2.0
2.0
15
15
15
0.
05
05
-05
-05
0.5
Figure 6: Comparison on Muller Brown potential with ground truth, learned features and fitted
eigenfunctions
D Training of neural network model
D.1 Evolution of the loss with eigenfunctions
It is interesting note that during the training process, the loss reaches progressively lower plateaus. This due to the fact that the NN has found novel eigenfunction orthogonal the previously found ones, starting with the constant until a new relevant direction is found. Typical behavior is shown in Figure 5 It was obtained for the case of a double well potential (see appendix below), but the same behavior observed in all the training sessions. This nice property is handy tool in properly optimizing the loss and understanding the proper stopping time.
Then during the plateau phase, the subspace is being explored
D.2 Training with imperfect features
One of the main advantages of our method is that even with features that are not eigenfunctions of the generator, but that were trained with our method, we can recover the good eigenfunction estimates as proven in Theorem [ In Figure 66 we illustrate such situation potential: the trained features do not represent the ground truth, however, using our fitting method on the same dataset, we managed to recover eigenfunctions close to the ground truth. This is to the best of our knowledge the first time this kind of "learn and fit" method has been applied to the learning of the infinitesimal generator.
simulation of the Muller Brown
19

---

1.0
0.8
0.6
0.4
0.2
0.0
-1.0
-0.5
0.0
0.5
1.0
X
Figure 7: Simulation of the double well potential (black dashed lines) under an effective biased potential (blue dashed lines). Our method (blue points) compared to ground truth (black line) and transfer operator based approach (red points)
D.3 Activation functions and structure of the neural network
In of our experiments we used the hyperbolic tangent activation function. This choice was made because it is widely used, bounded function with continuous derivative. It thus satisfies all the criteria needed for this method. Finally, when looking for eigenpairs, instead of having a neural network with outputs, we choose to have m neural networks with one output.
D.4 Hyperparameters
Besides common hyperparameters such as learning rate, neural network architecture and activation function, our method requires only two hyperparameters: 1 and a. Other methods such as [50) do not require 7), but on the other hand requires one weight per searched eigenfunction.
E
Experiments
For all the experiments we used pytorch 1.13, and the optimizations of the models were performed using the ADAM optimizer. The version of python used is 3.9.18. All the experiments were performed on workstation with AMD@ Ryzen threadripper pro 3975wx 32-cores processor and an NVIDIA Quadro RTX 4000 GPU. In all the experiments, the datasets were randomly split into training and validation dataset. The proportion were set 80% for training and 20% for validation. The training of deep models was performed using the micolvar package [
E.1 One dimensional double well potential
In this subsection, we showcase the efficiency of our method on a simple one dimensional toy model.
The target potential we want to sample has the form (zg(x) 4( 4.5 ex.5exxp(-8rrf + x²), which has form of wells separated by a high barrier which can hardly be crossed during simulation. In order to observe more transitions between the two wells and efficiently sample the space, we lower the barrier by running simulations under the following potential: Usm 44-0.5 exp -0zz)) x*))
which thus makes a bias potential: Viias (x) = (2) (z) - -4 4(eep((ppp-p) (oxp(-80k 80x (expp(-0x²) )). In Figure 7 we compare our method based on kernel methods (infinite dimensional dictionary of functions) with the ground truth and transfer operator baselines, namely deepTTCC 7 which is
method for molecular dynamics simulations, For this experiment we have used a Gaussian kernel of lengthscale 0.1, 0.1 and regularization parameter of 10-5
F.2 Muller Brown potential
For this experiment, the dataset was generated using an in-house code implementing the Fuler- Maruyama scheme to discretize the overdamped Langevin equation. The simulation performed at
20

---

Dataset
Dataset 2
Lime(ps)
Umelps)
Figure 8: Time evolution of the & angle in both datasets. Points are colored according to the value of the first nontrivial eigenfunction
a temperature of I (arbitrary unit) and a timestep of 10-3 The simulation was run for 10" timesteps. The bias potential is built according the following equation
|
Uzxt)
(36)
center added to this list. This kind of bias potential is called metadynamics [25] and allows reducing the height of the barrier for better exploration of the space. order to have static potential, no more centers are added after 300000 timesteps. We use learning rate of 5.10 the architecture of the neural network used is multilayer perceptron with layers of size 2 (inputs), 20, 20 and The parameter 7) was chosen be 0.05.
where the centers zy are built the fly: every 500 timestep
E.3 Alanine dipeptide
Simulation details All the simulations
with GROMACS 2022.3 (2 and patched with plumed 2.10 [45] in order to perform enhanced sampling simulations. We used the Amber99-SB [40] force field. The Langevin dynamics was sampled with timestep of 21s with a damping coefficent 7% m000 at target temperture of 300K. For both dataset, in order to make proper comparison, we used the explore version of OPES 17,, with barrier parameter of 20 kJ/mol and a pace of deposition of 500 timesteps.
Neural networks training For neural networks training, assume overdamped Langevin dynamics with the same value of the friction coefficient as in the simulation, as done in other works [50]. We learning of 10-3, and the architecture of the neural network used is multilayer perceptron with layers of size 30 (inputs), 20, 20 and 1. The parameter v) was chosen to be 0.1.
E.4 Chignolin
Simulation details The simulations are with GROMACS 2022.3 (2) and patched with plumed 2.10 [45). They share the same setup used for the D.E Shaw trajectory (28) used in the main text. For the same reason, kept the simulation condition consistent with that work. All simulations were performed with an integration time step of 2 fs and sampling NVT ensemble at 340K. The deepTDD model used for biasing the simulation is the obtained in ref (46
Neural networks training For our neural networks training, we assume overdamped Langevin dynamics. We use learning rate of 5.1044, and the architecture of the neural network used is multilayer perceptron with layers size 210 (inputs), 50, 50 and 1. The parameter 7) was chosen to be 0.2.
21

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: Justification: in Section we provide experimental proofs of our claims.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results
be expected to generalize other settings.
is fine to include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes],
Justification: see the conclusion section.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification. asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low images are taken in low lighting. Or speech-to-text system might be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question:
For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
22

---

Justification: see Theorem and 1 proofs are provided in the Appendix C
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Answer: [Yes]
Justification: the experimental setting is reported in Section andd expanded in the appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make clear how to reproduce that algorithm. (b) the contribution is primarily a new model architecture, the paper should describe the architccture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reprodueibility. the case of closed-source models, may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
23

---

Answer: [Yes]
Justification: The experimental setting is reported in Section 6 and expanded in the appendix. The code is also shared on a github repository. The codes used to the simulations all open source.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.cm cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
7 //ipps.ccppblic//uu.dde ss ssinisolicyyyy for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/eetails
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: please see Appendix
E
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes].
Justification: whenever appropriate we have reported these.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emror bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).

---

It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: please see Appendix
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics http://aaaaag..cc//ba/a///ccctee
Answer: [Yes]
Justification:
We enjoy research and we respect other people work.
Guidelines:
The answer NA means that the authors have reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification:
We briefly present in the conclusions possible broader impacts of this work.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no socictal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
25

---

The conference expects that many papers will be foundational research and not tied particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system Iearns from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
his paper is proof of concept, therefore we do not use large scale data or do not need any safeguard
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes)
Justification: We used open-source codes and algorithms, such as Gromacs and Plumed which are duly cited. D.E. Shaw has been acknowledged for providing the chignolin trajectory
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms use in package should be provided. For popular datasets, paperswithcode com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
26

---

For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: The github repository is provided.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/moodl as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: not applicable.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: applicable.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
27


Reviews:
[{"review_id": "YzW4DgH9Az", "reviewer": "Reviewer_MSBv", "summary": "This paper studies an infinitesimal generator approach for learning the eigenfunctions of evolution operators for Langevin SDEs. Due to the slow mixing caused by the high potential barriers, direct learning from simulation data can be sample inefficient. Biased simulation (based on a biased potential) is used to explore the space faster, and importance weights are constructed to get an unbiased loss function; such unbiased construction is more natural and feasible in the generator approach compared to the existing transfer operator approach. The minimizer of the resulting quadratic loss function, given a dictionary of basis functions, can be obtained by ridge regression type algorithms as it is a linear method. The authors also propose a loss function to learn a good dictionary that approximates the space of eigenfunctions more accurately; this nonlinear dictionary learning improves accuracy. Experiments on molecular dynamics benchmarks demonstrate the effectiveness of the approach.", "questions": "- Is the estimator in Section 4 simply the Galerkin approximation of the inverse of $\\eta I - \\mathcal{L}$ with the basis functions in $\\mathcal{H}$? \nI found the descriptions in Section 4 complicated and not easy to digest. As an example, for equation (11), by definition, the formula will differentiate the term $\\|\\chi_\\eta(x) - \\hat{G}^Tz(x)\\|_2$ (which is not differentiable at zero) which looks strange. And I didn't understand the sentence \"we contrast the resolvent ...\" on page 5, line 182: what do the authors mean by \"contrast\" here? And the explanation of motivation of using the generator regression problem (11) rather than the mean square error also appears less clear.\n\n- I am also curious about the mechanism behind the improvement of accuracy in Figure 2. It appears to me that the approach in this paper first uses the loss function in equation (19) to find approximate eigenfunctions and then uses the generator approach in section 4 to refine the estimation of the eigenfunctions and eigenvalues. If so, an understanding of which step of the two plays the significant role will be useful. For example, if the authors apply the generator approach in section 4 to the approximate eigenvalues obtained in the work of Zhang et al. (2022), will the resulting accuracy in Figure 2 also be significantly improved?\n\n- Potential typos: Page 2, line 67: \"principle\" -> \"principled\".\nPage 3, line 125-126: the notation of $C_{\\gamma}$ is not introduced. Moreover, it seems $\\lambda_i = \\log \\mu_i / t$ rather than $\\log (\\mu_i/t)$.\nPage 8, line 282, \"approqch\"\n\n- Page 3, line 129: \"If t is chosen too small, the cross-covariance matrices will be too noisy for slowly mixing processes\" Could the authors elaborate more on this?", "limitations": "Limitations are discussed.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The introduction and motivation are exceptionally well written and demonstrate the benefits of the generator approach compared to the transfer operator approach for biased simulations. The numerical experiments, especially in Figure 1, show significant improvements in accuracy compared to previous methods.", "weaknesses": "The mathematical presentation of the technical details, namely sections 4 and 5, could be made clearer. Many notations are introduced, but the ideas seem simple; see Questions. I would appreciate an overarching description highlighting the key idea and the difference between the proposed approach and existing works using generators.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer's insightful evaluation and valuable comments. In what follows, we aim to address the highlighted weaknesses and respond to the reviewer's questions.\n\n## Weaknesses:\n\nThank you for motivating us to present the overarching summary of our approach. While we discuss this in the general reply, we would like to provide a few more relevant details here. \n- First, data-driven generator methods emerged relatively recently in the literature [1,15,20,21,46], and up to our knowledge, the only method that possesses finite-sample statistical guarantees for the proper spectral estimation is [21]. However, all of the above approaches analyzed learning from unbiased simulations, which is unfeasible in many realistic scenarios. Moreover, the method in [21] suffers from scalability issues inherent to kernel methods. On the other hand, deep learning methods for the generator were applied to biased simulations [46], but their statistical consistency was not proven.  __In sharp contrast to all of the above, our work is the first to show that one can learn eigenvalues and eigenfunctions__ of the infinitesimal generator via its resolvent __in a scalable and statistically consistent way from biased simulations__. \n- To summarize, we extend the method in [21] by __incorporating biased dynamics__, developing appropriate __representation learning__, and designing the overall scalable and statistically consistent approach.  Our method is twofold: first, we learn a basis set using Theorem 2 (representation), and then we learn the generator’s resolvent on this basis set using Theorem 1 (regression). This approach is common in many fields, notably in transfer operator approaches [24,29], but is applied for the first time ever, to the best of our knowledge, to the case of the infinitesimal generator. In doing so, we followed the core idea of [21] to formulate the problem of learning the resolvent so that the geometry of the process is efficiently exploited via the energy norm, see also reply to Q1. Practically, this means that one can avoid inverting the shifted generator by rescaling the norm in which one learns the resolvent. This allows one to work in the latent space where the state becomes represented as $(\\eta I - L)^{1/2} z$ and all the relevant inner products and norms can be computed using the features $z$ and their gradients. \n\n## Questions:\n\n- __Q1:__ In fact, the estimator is not a Galerkin projection of the resolvent onto $\\mathcal{H}$ as a subspace of $\\mathcal{L}^2_\\pi$. This would be the case if in (11) we replace the energy with the expectation w.r.t. invariant distribution. However, note that this is not feasible since we cannot compute/estimate the action of the resolvent. To solve this issue we change the domain of the problem from  $\\cal{L}^2_\\pi$ (expectation) to $\\cal{H}\\^\\eta\\_\\pi$ (energy). So, the reviewer can think of our estimator as a __(regularized) Galerkin projection of the resolvent acting on $\\cal{H}\\^\\eta\\_\\pi$ onto a subspace $\\cal{H}$.__ The reason why we do so, as shown in (12), lies in the fact that Hilbert Schmidt norm of an operator $A$ on $\\cal{H}^\\\\eta_\\\\pi$  becomes $\\rm{tr}(A^*(\\eta I-L)A)$. That is, the features $z\\in\\mathcal{L}^2_\\pi$ are transformed into features $\\tilde z:=(\\eta I - L)^{1/2}) z$, so inner products $\\langle \\widetilde z, \\widetilde z’ \\rangle = \\langle z,(\\eta I - L) z’ \\rangle$ and $\\langle \\widetilde z, (\\eta I - L)^{-1} \\widetilde z’ \\rangle = \\langle z, z’ \\rangle$ can be now estimated from data. Using this trick we neatly avoid estimating the action of the resolvent, and no differentiation in (11) is needed since integral in $\\chi_\\eta$ is essentially canceled out by differential operator $\\eta I -L$. This is what we meant by contrasting the resolvent in line 182. To further clarify our approach to the reader, we will include this discussion in the revision.  \n\n- __Q2:__ Do we correctly understand that the reviewer is referring to Figure 1 instead of 2, since there we compare to [46] (Zhang et al. 2022)? If so, let us note that, as pointed out, we can use both methods for representation learning (loss (19) or loss 39 from [46]) to identify a subspace and then apply our regression method. Based on the theoretical background given in [46], we believe that their method indeed reduces representation error in Theorem 1, making this approach sound. However, when we performed this additional test, the first eigenvalue went from 0.043 to 0.013, the second one went from 1.13 to 4.19. This indicates that even with regression the result does not significantly improve. To strengthen our point, we also did the opposite operation: we computed the Rayleigh quotient from our features without doing any regression, just as is done in the work [46] we obtain for the first eigenvalue 1.79e-4 and 0.36 for the second one, which are close to the ones found by regression and very close to the ground truth values. We are speculating that the main reason is in the efficiency of our loss and its robustness on the hyperparameter choices.  \n\n- __Q3 and Q4:__ Thank you for pointing out typos, we will correct them. Also, sorry for the confusion, we intended to say that if $t$ is chosen too small w.r.t. the mixing time, the noise in the corresponding transfer operator regression dominates the signal, meaning that the learning becomes harder. We will correct this discussion in the revision.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "wKNLT7MA3O", "reviewer": "Reviewer_txkX", "summary": "In this paper, the authors investigate the possibility of estimating the leading eigenvalues and the corresponding eigenvectors for the evolution operators of Langevin dynamics using biased simulations. To this end, they rely on strong statistical guarantees and on the use of deep learning regression to build a suitable Hiltbert space that approximates the evolution operator of the unbiased process using only biased data. They evaluate the reliability of this approach in three increasingly complex problems and show that this approach successfully recovers the slowest relaxation modes in toy models and obtains better approximate values for the eigenvalues than competing state-of-the-art methods.", "questions": "I did not understand well the physical meaning of the modes they extract with their approach. Could they give a physical explanation of them or of the order parameter associated?", "limitations": "Authors discuss correctly the limitations", "rating": 7, "confidence": 1, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The work deals with a very important problem in chemistry and physics, namely the identification of the slowest modes of a dynamic process by means of simulations. In complex cases, simulations are not long enough to observe the most prominent transitions, and enhancing sampling methods are used to facilitate the observation of rare events. The problem is that it is often difficult to identify the key collective variables to facilitate the required slow movements. In this paper, the authors propose a simple way to do this using biased simulation data based on strong theoretical guarantees.", "weaknesses": "I find it difficult to read the paper and follow it, perhaps because I am too far from the field", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer's insightful evaluation and valuable comments. In what follows, we aim to address the highlighted weaknesses and respond to the reviewer's questions.\n\n## Weaknesses:\nWe hope that our general reply brings more clarity. We commit to incorporate all the suggestions and additional discussions in the revised version of the paper. If the reviewer would like us to clarify any other specific aspects, we are happy to address them during the discussion period.\n\n## Questions:\nThank you for raising this question, in the following we try to briefly review the usefulness of generator eigenfunctions (“modes”) for interpretability and control of Langevin dynamics. We plan to add it in the revision to further clarify the rationale of our method for the general reader. From the spectral decomposition of the generator in Equation (4) and its link to the transfer operator in Equation (5) , it can be seen that each mode $f_i(x)\\langle f_i,f \\rangle$ (essentially the eigenfunction $f_i$) is related to a timescale of the process ($-\\lambda_i>0$). This means that eigenfunctions corresponding to the slow time-scales (eigenvalues close to zero ) give information on __metastable states__. Recalling that normalized eigenfunctions (in $\\mathcal{L}_\\pi^2$) have mean zero and unit variance, they tend to have constant values for ensembles $x$ in such meta-stable states (meaning at time-scale $-\\lambda_i>0$ there is no dynamics), while it changes sign when ensembles move from one meta-stable state to another. This is what the reviewer called “order parameter”. But it is more than this, because contrary to an order parameter which only discriminates two states, the eigenfunction will also give valuable information on how the transition takes place, see illustrative example in appendix E.1 and corresponding Figure 5. This feature of the eigenfunctions is of particular relevance, since the knowledge of rare transitions can be used by experimentalists to accelerate or slow down the process, see e.g. Wei Zhang, Christof Schütte, _Understanding recent deep‐learning techniques for identifying collective variables of molecular dynamics_. PAMM 2023, 23 (4).", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "txQ3k2U5vR", "reviewer": "Reviewer_BQ7h", "summary": "@Authors, I would appreciate it if you could point out any inaccuracies in the following summary since it took me a long time to understand your paper and am still not completely certain my understanding is correct.\n\nThe paper aims to address a common problem for molecular dynamics simulations, which is to obtain collective variables that can be used to bias simulations such that, e.g., rare transitions occur more frequently and important observables can be computed faster. The authors follow previous work and learn the system's generator, whose eigenvalues inform transition timescales and whose eigenfunctions can be used as collective variables to bias simulations along. The novelty is in deriving an approach that allows learning the generator from biased enhanced sampling simulations which converge faster and observe more transitions than conventional simulation. The authors evaluate their method on toy systems.", "questions": "I would appreciate any time you can take to answer some of my questions.\n1.  107 \"for the  transfer operator we can only observe a noisy evaluation of the output\": Do you mean a noisy version of the expectation that defines a transfer operator? Or that the output one can observe after simulating stochastic dynamics is \"noisy\"?\n2. I assume your underlying goal is to discover CVs as the eigenfunctions of your learned generators. Your generators can be learned from biased enhanced sampling simulations where e.g. more transitions occur. Then we can obtain the CVs from your generator. With the CVs we can run biased simulations. Why are these biased simulations more informative than the biased simulations which you ran to train your generator? Can we compute different quantities from your CV biased simulations? Is it just a matter of your CVs being better biases than the bias of the original enhanced sampling bias?", "limitations": "The paper discusses some limitations such as not being time dependent - the significance of which I cannot assess. Limitations such as the limited evaluations and missing understanding of when the method fails seem more significant to me but are not mentioned or explained why they are present.", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 2, "contribution": 4, "strengths": "1. Novel approach to accelerate obtaining scientifically important collective variables by enabling their computation from simulations biased by enhanced sampling. This approach requires developing new theory and to derive a training objective. The procedure to do so seems non trivial to me but I cannot sufficiently judge the value of the provided theory as I can only follow the broad strokes and am not well versed in the used math. \n   I hope other reviewers can provide more useful signal for this aspect.\n2. Impactful application: The paper provides useful ideas for a problem with significance for scientific applications that can have large downstream impacts.\n3. Very well written introduction and overview of the transfer operator and generator formalism.", "weaknesses": "My points contain many understanding questions and some \"criticizing\" questions. It would be great if you could also answer the understanding questions.\n1. Experiments: The experiments are carried out on three toy systems and the proposed approach visually improves upon two deep learning baselines in one of the three experiments. \n\t1. I understand that the cited works on generator learning provide the same or fewer experiments. However, from skimming the cited works it also seems to me that one of the main areas of application would be simulations on proteins. Is that incorrect and if not, why do you and others not provide experiments on proteins? Aren't there proteins with well known dynamics on which we could evaluate the methods?\n\t2. **I think this method should be tried on systems of increasing sizes until it fails. The experiments of a paper introducing a new method ideally provide information on the capabilities of a method.**  Why do you (and the rest of the generator learning papers) not do so?\n\t3. Why do you omit the deep learning baselines from the ALDP experiments?\n2. For an understanding of how your method fits into the broader landscape of approaches to determine CVs, I think it would be important to also provide comparisons with non-deep learning approaches. The relationship to DL methods is clear but to understand the general \"usefulness\" it would be nice to have delineation from classical methods in terms of approach and experiments. Are there classical approaches that can be used as baselines that would outperform all DL methods?\n3. Speed comparisons: it seems to me that the underlying tradeoff for all methods is between speed and accuracy and that accuracy alone might be less meaningful. Is there nothing to be said about the runtimes for data collection and training?\nMinor:\n1. To make the paper more broadly accessible to the wider deep learning community, I think it would be helpful to describe your final operator learning approach more procedurally on e.g. the ALDP example. What is the neural network input and output dimensionality? Why do you have separate neural networks for each operator output dimension? What is the input and output to your neural network? \n   Once you trained your neural network, how do you obtain the eigenfunctions from it.\n2. How do you obtain your plots in e.g. the ALDP figure Figure 2.? To understand this: what concretely is the eigenfunction in practice once you computed it? How is it represented? Once you have it, how does it assign different values to each molecular structure?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer's insightful evaluation and valuable comments. In what follows, we aim to address the highlighted weaknesses and respond to the reviewer's questions.\n\n## Summary & strengths:\n- Thank you for your summary, while it captures the general idea of the paper, maybe it lacks one important aspect that, to the best of our knowledge, __we are the first to show__ that covariance reweighting in transfer operator (TO) models built from biased simulations fails, while infinitesimal generator (IG) models do not suffer from this effect. Indeed __we prove that one can learn eigenvalues and eigenfunctions__ of the IG via its resolvent __in a scalable and statistically consistent__ way from a single biased trajectory.\n\n- Related to the first strength, please refer to our general reply, part _Clarity of presentation_. Besides proposing Alg. 1 (methodological contribution), we also __theoretically prove its statistical consistency__. Let us provide the context: data-driven IG methods emerged relatively recently in the literature [1,15,20,21,46], and up to our knowledge, the only method that possesses finite-sample statistical guarantees for the proper spectral estimation is [21]. However, all of the above approaches analyzed learning from _unbiased_ simulations, which is unfeasible in many realistic scenarios. Moreover, the method in [21] suffers from scalability issues inherent to kernel methods. On the other hand, a deep learning method for the IG [46] was applied to biased simulations , but its statistical consistency was an open question.  \n\n## Weaknesses:\n\n1. Noting that learning the spectral decomposition of IG is a problem of interest besides protein dynamics, e.g. crystallization, we agree with the suggested best practice. We did our best to meet such a requirement, as reported in _Experiments_ of our general reply and the attached pdf. \n  - We believe that the main reason why IG and TO methods are not initially tested on larger scale systems is the complexity of the task. Namely, recovering proper time scales and the transition state ensemble is a complex task that highly depends on the descriptors of the system and their effective dynamics. Indeed for large proteins one of the main challenges is to design appropriate descriptors so that effective Langevin dynamics can reveal the true dynamics of a protein [D]. These are highly dependent on the protein, and their design is a problem per se. Once this has been resolved, the next challenging task is to explore the rare transitions by biasing dynamics and recover proper time scales and transition state ensembles. While our method is designed to solve the latter, the former is an active field of research that currently limits the development of the full pipeline for large scale complex molecular systems.  \n  - In order to fully address the question on the capability of our method, we have performed an additional experiment. To the best of our knowledge, the most complex benchmark for all TO and IG methods is the (very long) unbiased simulation dynamics of a small protein, typically Chignolin [7,19]. So, we designed a biased simulation experiment based on our method, which is, to the best of our knowledge, the first experiment of its kind for generator learning.\n  - Concerning baselines for ALDP, since we don’t have ground-truth it would be just a qualitative comparison.\n2. Most non DL CV are heuristics based on chemical/physical intuition, they are system dependent and often suboptimal. This has stimulated intense research towards ML based alternatives. In this context, TO/IG methods are well-motivated, since eigenfunctions of these operators are “optimal” CVs for biasing simulations, see e.g. [A]. While some classical ML methods are available for this [21,23], they are intractable at large scale since at each time step one needs to compute the kernel function with respect to the current ensemble of the simulation. \n3. Our method is as expensive as standard CV based methods, but it reliably leads to more accurate results, which is important when one wants to reconstruct the free energy profile and use them iteratively to learn the generator. See also our additional experiment on Chignolin, see the global reply and the attached pdf where we highlight the computational speedup relative to unbiased bruteforce simulations.\n4. We incorporated your suggestion in the general reply on the general pipeline. Here we address your related questions. First, since the input of the DNN encoder in the representation learning part are the descriptors, the input dimension depends on the problem. On the other hand, the output (latent features) dimension is typically very small,  reflecting the number of slow time-scales to learn (typically around 3-5). Second, the choice of using separated DNNs per latent feature is made to aid learning of orthonormal features, since, depending on the architecture, shared weights may implicitly introduce prohibitive constraints. Third, eigenfunctions, defined in line 13 of Alg. 1, map descriptors of the system to a real number that, when properly scaled, has mean zero and variance 1 w.r.t. equilibrium distribution.\n5. Eigenfunctions can be plotted in the plane of relevant features to analyze the data. For example, for alanine dipeptide we have for each configuration the values of the angles $\\phi$ and $\\psi$ and the values of the eigenfunction, therefore, we may plot the points in the $\\phi-\\psi$ plane colored according to the eigenfunction value. \n\n## Questions:\n1. If we understand well, two statements are equivalent. Indeed the noise of the feature is $z(X_{t+s})  - \\mathbb{E}[z(X\\_{t+s})\\vert X_t] = z(X\\_{t+s}) -[ {\\cal T}\\_s z ] (X_t)$, and from one trajectory we only observe an instance of the r. v. $z(X_{t+s})$, and not the output of TO.  \n2. Please see our general reply.\n\nRef:\n[D] Zhang et al., Effective dynamics along given reaction coordinates and reaction rate theory. Faraday Discuss. (2016)", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We wish to thank all reviewers for their insightful evaluation of our paper. We appreciate all their comments and remarks, which we’ll incorporate in our revision. Before addressing each review in detail, we’d like to point out some general remarks that apply to all of them. 

##   Clarity of presentation 
To improve the clarity of our work in Sec. 4 and 5, and emphasize its impact, let us provide a general overarching view of our method. Dealing with rare events implies a large timescale, thus generator’s eigenvalues close to zero, whose eigenfunctions reveal transitions between meta-stable states. To solve these learning tasks via transfer operator approaches one typically needs a prohibitively large amount of data from simulations with very small time-increments. We overcome this problem by proposing to work with the infinitesimal generator which is, contrary to transfer operators, time-scale independent and, as we show, can be successfully combined with biased simulations.  
In particular, we propose to work with the resolvent of the infinitesimal generator, which shares the same eigenfunctions with the generator. However, this approach introduces a new issue: an operator needs to be inverted. To overcome this, we work in a new Hilbert space $\cal{H}\_\pi^\eta$ with inner product defined in Eq. (10). In this space, the Hilbert Schmidt norm of an operator $A$ on $\cal{L}^2_\pi$ becomes $\rm{tr}(A^*(\eta I-L)A)$, that is, features $z$ are transformed into features $(\eta I-L)^{1/2}) z$. Using this trick, __the resolvent may be regressed easily__ in a theoretically grounded way knowing only the diffusion part of $L$, see also [21] for more details on this approach using kernels. On the other hand, to obtain a method that is truly scalable to large problems, one needs a good dictionary of functions on which to regress the generator. This  __representation learning problem__ is addressed in Section 5, where we use the same weighted norm to define a novel  loss. While we show in Thm. 2 that deep neural network representations, when properly trained, converge towards the eigenfunctions of the resolvent, due to imperfect dataset or parameter choices, this convergence might be slow. However, the learned features are useful to identify a subspace of the $\cal{L}^2_\pi$ domain with small representation error, and one can efficiently regress the resolvent with the method presented in Sec. 4 to obtain the overall statistical consistency from Thm. 1.  

## General pipeline 
Following the comments of the referees, we also believe that a description of our method pipeline would be beneficial. We propose it in the form of Figure 1 in the attached pdf, and commit to include it in the revision, together with the following discussion and more detailed version of the current pseudocode (Alg. 1) in the appendix:
First, one chooses a molecular system to study, and, based on this choice, identifies a __collective variable (CV)__ on which to bias the dynamics to observe transitions from one state to the other. Often this CV is heuristic and suboptimal, leading to only a few transitions. To apply our method (and any deep-learning based method), one then needs to __choose descriptors__ of the system to put as input of the model. These descriptors will encode the symmetries of the system (often global translation/rotation) and only take into account the most important degrees of freedom: for example for the folding process of a protein in water, one may only use the interatomic distances between heavy atoms. This phase is highly system dependent, and the complexity of choosing descriptors grows with the system size. Once descriptors are chosen, a __representation is learned__ and the derivatives present in the loss function are computed with respect to atomic positions. After learning the representation, __the resolvent is regressed__ on the learned dictionary of functions. If the initial dataset does not have enough transitions, it might be that the eigenfunctions obtained are not good enough. In that case, one can use the learned eigenfunction as a collective variable for biasing the simulation to collect more transitions and further improve our model of the generator, and hence the final estimation of its eigenfunctions. Note, however, that the last point is based on the fact that the eigenfunctions of the generator are optimal CVs, see  [A]. While we did not need to implement it in our experiments, future work will study more complex systems where this iterative process can be beneficial.

## New experiment 
Finally, to showcase the power of our method, we implement it for the Chignolin miniprotein: we first performed a 1 microsecond biased simulation using as CV the deep-TDA one [B]. This allowed us to gather transitions. Then, this trajectory was used to train our method, which is, to the best of our knowledge, the largest scale work on generator learning. To validate our results, we trained our method on a very long unbiased trajectory of 106 microseconds from D.E. Shaw research [C]. The results are presented in Fig. 2 of the attached pdf. The two presented eigenfunctions are completely similar, but one is from only 1 microsecond long simulation, while the other is from 106 microseconds (so a 106 times more costly computations). This  shows that __our work on biased simulations paves the way towards scalable and reliable estimation of dynamical quantities in realistic problems of molecular dynamics__.  

We once more thank all the reviewers whose comments inspired the discussion above which, we believe, further strengthens our paper, and demonstrates its impact.    
The Authors 

### Ref.
[A] Zhang & Schütte, Understanding recent deep‐learning techniques for identifying collective variables of molecular dynamics. PAMM (2023)

[B]  Trizio & Parrinello, From enhanced sampling to reaction profiles. J. Phy. Ch. Letters (2021)

[C] Lindorff-Larsen et al. “How Fast-Folding Proteins Fold”. Science (2011)

Meta-review: This paper introduces a biased potential to improve the efficiency of sampling and inferences of simulated SDE. Reviewers agree that the paper is tackling an important and difficult problem, especially in modeling protein dynamics. They introduce a new method supported by theory and experiments. If accepted, please make sure to include new experimental results and incorporate improvements that stemmed from the rebuttal discussion.
Decision: Accept (poster)