From Biased to Unbiased Dynamics: An Infinitesimal Generator Approach

Authors: timothée  devergne
vladimir r kostic
michele  parrinello
massimiliano  pontil
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 17:50:11

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.

Main paper: 
 introduction: 
 Dynamical systems and stochastic differential equations (SDEs) provide a general mathematical framework to study natural phenomena, with broad applications in science and engineering. Langevin SDEs, the main focus of this paper, are widely used to simulate physical processes such as protein folding or catalytic reactions [see e.g. 47, and references therein]. A main objective is to describe the dynamics of the process, forecast its evolution from a starting state, ultimately gaining insights on macroscopic properties of the system.In molecular dynamics, the motion of a molecule is sampled according to a potential energy U (x), where the state vector x represents the positions of all the atoms. Specifically, the Langevin equation dX t = -∇U (X t )dt + σdW t describes the stochastic behavior of the system at thermal equilibrium, where X t is the random position of the state at time t, the scalar σ is a multiple of the square root of the system's temperature, and W t is a 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 a simulation. For instance, for a protein the free energy barrier between the folded and unfolded states is larger than thermal agitation, making the transition between the two states a rare event. Consequently, long trajectories need to be simulated before such interesting events are observed. In fact, one needs to observe many events to get the relevant thermodynamics (free energy) and kinetics (transition rates) information [34] . Beyond molecular dynamics, the slow mixing behavior of many systems modeled by SDEs is a major bottleneck in 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, 11] . The main idea is to add to the potential energy a 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 is often chosen as a function of a few wisely selected variables called collective variables (CVs). For instance, if a chemical reaction involves a 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 [8, 42, 14, 9, 10, 4, 27] . A key idea is to use available dynamical information to construct the CVs [41, 31, 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 at 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.More recently, there has been growing interest in learning the infinitesimal generator of the process [15, 1, 50, 20] , which allows one to 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, 21] , we introduce a 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 to 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 state-of-the-art transfer operator and recent generator learning approaches in biased simulations. Remarkably, even with datasets containing only a few relevant transitions due to sub-optimal biasing, our method effectively recovers crucial information about the transition mechanism.Paper organization In Section 2, we introduce the learning problem. Section 3 explores limitations of transfer operator approaches. In Section 4, we review a recent generator learning approach [21] and adapt it to nonlinear regression with a finite dictionary of functions. Section 5 presents our method for learning from biased dynamics. Finally, in Section 6, we report our experimental findings.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 6 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 be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
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, [22] [23] [24] , contributing in particular to the existence of physical priors and feasibility of data acquisition for successful learning.Stochastic differential equations (SDEs) and evolution operators While our observations in the paper naturally extend to general forms of SDEs [see e.g. 33], to 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 equationEQUATIONdescribing dynamics in a (state) space X ⊆ R d , governed by the potential V : R d → R at the temperature β -1 = k B T , where W t is a R d -dimensional standard Brownian motion.The SDE (1) admits a unique strong solution X = (X t ) ⩾0 that is a Markov process to which we can associate the semigroup of Markov transfer operators (T t ) t≥0 defined, for every t ≥ 0, asEQUATIONFor (1) the distribution of X t converges to the invariant measure π on X called the Boltzmann distribution, given by π(dx) ∝ e -βV (x) dx. In such cases, one can define the semigroup on L 2 π (X ), and characterize the process by the infinitesimal generatorL := lim t→0 + (T t -I)/t defined on the Sobolev space H 1,2π (X ) of functions in L 2 π (X ) whose gradient are also in L 2 π (X ). The transfer operator and the generator are linked one to another by the formula T t = exp(tL). Moreover, it can be shown (see Appendix A) that the generator L acts on f : X → R asEQUATIONwhich, integrating by parts, gives (Lf )g dπ=-β -1 ⟨∇f, ∇g⟩ dπ= f (Lg)dπ, showing that L is self-adjoint. If L has only a discrete spectrum, one can solve (1) by computing the spectral decompositionEQUATIONUsing (2) and the exponential relation between the transfer operator and the generator, one can writeEQUATIONwhere 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 algorithms for the transfer operator and the generator lies in the nature of the data used. While for the transfer operator we can only observe a 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 an unbounded differential operator. Consequently, learning methods for the former align with vector-valued regression in function spaces [22] , whereas methods for the latter, as discussed in the following section, are more akin to physics-informed regression algorithms.In both settings, we learn operators defined on a function (hypothesis) space, formed by the linear combinations of a prescribed set of basis functions (dictionary) z j : X → R, j ∈ [m],EQUATION)The choice of the dictionary, instrumental in designing successful learning algorithms, may be based on prior knowledge on the process or learned from data [24, 30, 50] . The space H is naturally equipped with the geometry induced by the norm ∥h u ∥ 2 H := m j=1 u 2 j . Moreover, every operator A : H → H can be identified with matrix A ∈ R m×m by Ah u = z(•) T Au. In the following, we will refer to A and A as the same object, explicitly stating the difference when necessary.Transfer operator learning Learning the transfer operator T t can be simply seen as the vector-valued regression problem [22] , in which the action of T t : L 2 π (X ) → L 2 π (X ) on the domain H ⊆ L 2 π (X ) is estimated by an operator T t : H → H. This aims to minimize the mean square error (MSE) w.r.t. the invariant distribution. Given a dataset D n := (x i , y i = x i+1 ) n i=1 of time-lag t > 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 Tγ := C -1 γ C t , where the empirical covariance matrices are C = 1 n i∈[n] z(x i )z(x i ) T and C t = 1 n i∈[n] z(x i )z(y i ) T .We then estimate the eigenpairs (λ i , f i ) in (4) by the eigenpairs ( µ i , u i ) of T γ as λ i := ln( µ i /t) and f i :=z(•) T u i .We stress that transfer operator approaches crucially relies on the definition of the time-lag t from which dynamics is observed. Setting this value is a delicate task, depending on the events one wants to study. If t is chosen too small, the cross-covariance matrices will be too noisy for slowly mixing processes. On the other hand, if t is too large, because the relevant phenomena occur at large time scales, a 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.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 are not discussed in the paper. • The authors are encouraged to create a 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 on 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 a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not 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 to address problems of privacy and fairness. • 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 acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
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 ′ (x) := U (x) + V (x), x ∈ Xwhere 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 it challenging to recover the unbiased dynamics from biased data. Denoting the invariant measure of the perturbed process by π ′ and its generator byL ′ : H 1,2 π ′ (X ) → H 1,2 π ′ (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 L are also eigenfunctions of every transfer operator T t = e tL , we can address the related problem of learning the transfer operator from perturbed dynamics. Unfortunately, there is an inherent difficulty in doing so. While one typically knows the perturbation in the generator, that is L ′ = L + ⟨∇V, ∇(•)⟩, this knowledge is not easily transferred to the perturbation of the transfer operator. Indeed, recalling that T := T 1 = e L , and since the differential operator ⟨∇V, ∇(•)⟩ in general does not share the same eigenstructure of L, one has that T ′ := e L ′ = e L-⟨∇V,∇(•)⟩ ̸ = T e -⟨∇V,∇(•)⟩ .Simply put, the generator depends linearly on 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 is the Boltzmann distribution defined by the potential. Hence, we have that π(dx) = e -βU (x) dx e -βU (x) dx , π ′ (dx) = e -βU ′ (x) dx e -βU ′ (x) dx and dπ dπ ′ (x) = e βV (x) e βV (x) π ′ (dx) (7) where the last term is the Radon-Nikodym derivative, which exposes the data-distribution change.Consequently, we can express the covariance operators for the unperturbed process as weighted expectations of the perturbed data featuresEQUATIONHowever, since the transition kernel of the process (X ′ t ) t≥0 generated by L ′ is different from that of the original process, the above reasoning does not hold for the cross-covariance matrix, that is,C t :=E X0∼π ′ dπ dπ ′ (X 0 ) z(X 0 )z(X t ) T ̸ = E X ′ 0 ∼π ′ dπ dπ ′ (X ′ 0 ) z(X ′ 0 )z(X ′ t ) T =: C ′ t ,Consequently, the estimator T t obtained by minimizing the reweighed risk functionalR ′ ( T t ) := E X0∼π ′ dπ dπ ′ (X 0 ) ∥z(X ′ t ) -T T t z(X 0 )∥ 2 2does not minimize the true risk since R ′ ( T t ) ̸ = R( T t ). Despite this difference, whenever the perturbation is small or controlled and the time-lag t is small enough, estimating the true transfer operator of the process from the perturbed dynamics via reweighed covariance/cross-covariance operators has been systematically used as the state-of-the art approach in the field of atomistic simulations [7, 9, 31, 49] . The (limited) success of such approaches is based on a delicate balance of a small enough lag-time and biased potential, since for small t > 0 one can approximate C t by C ′ t and minimize R ′ (T ) ≈ R(T ) over T : H → H.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer Justification: please see Appendix E.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or 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). 
infinitesimal generator learning: 
 In this section, we address generator learning. While there has been significant progress on this topic [24, 15, 35, 50, 20] , we follow the recent approach in [21] for learning the generator L on an a priori fixed hypothesis space H through its resolvent. Leveraging on its strong statistically guarantees, we adapt it from kernel regression to nonlinear regression over a 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 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 (ηI-L) -1 , where η > 0 is a prescribed parameter. To this end, we observe that the action of the resolvent in H can be expressed as ((ηI -L) -1 h u )(x) = χ η (x) T u, where χ η is the embedding of the resolvent [21] . We then aim to approximate χ η (x) ≈ G * z(x) by a matrix G ∈ R m×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(ηI -L) -1 into H, given by χ η (x) = ∞ 0 E[z(X t )e -ηt | X 0 =x]dt, x ∈ X , seeEQUATIONand, due to the identity f Lgdπ = -β -1 (∇f ) T (∇g)dπ, alsoEQUATION)Since L is negative semi-definite, the above kernel induces the regularized squared energy normE η π : H 1,2 π (X ) → [0, +∞) by E η π [f ] := E η π [f, f ] = E x∼π ηf 2 (x) -f (x)[Lf ](x). It counteracts the resolvent and balances the transient dynamics (energy) of the process with the invariant distribution π. In a nutshell, instead of using the mean square error of f (x) := ∥χ η (x) -G T z(x)∥ 2 to define the risk, we "fight fire with fire" and penalize the energy to formulate the generator regression problem minEQUATIONIndeed, this risk overcomes the difficulty of not knowing χ η . To show this, let us define the spaceH η π (X ) := {f ∈ H 1,2 π (X ) | E η π [f ] < ∞} associated to the energy norm ∥f ∥ H η π := E η π [f ],and recalling that the operator G : H → H is identified with a matrix G ∈ R m×m via Gh u = z(•) T (Gu), define the (injection) operator Z : R m → H η π by Zu = z(•) T u, for every u ∈ R m . Then, since HS (R m , H η π ) ≡ HS (H, H η π ), the norm is the sum of squared H η π norm over the standard basis in R m , and one obtainsEQUATIONwhere P H is the orthogonal projector in H η π (X ) onto H. In learning theory ρ(H) is known as the approximation error of the hypothesis space H [see e.g. 43]. While this error may vanish for infinite-dimensional spaces, when H is finite dimensional, controlling ρ(H) is crucial to achieving statistical consistency. This can be accomplished by minimizing (11) , which is equivalent toEQUATIONwhere (•) † is the Moore-Penrose's pseudoinverse. Using the covariance matricesEQUATIONw.r.t. the invariant distribution and energy, respectively, gives the ridge regularized (RR) solutionG = (W + γI) -1 C, γ > 0. The induced RR estimator of the resolvent, G η,γ : H → H is given, for every h u ∈ H, by G η,γ h u := Z(W + γI) -1 Cu = z(•) T (W + γI) -1Cu, and it can be estimated given data from π by replacing expectation and the energy in (14) with their empirical counterparts. 
unbiased learning of the infinitesimal generator from biased simulations: 
 In this section, we present the main contributions of this work: approximating the leading eigenfunctions (corresponding to the slowest time scales) of the infinitesimal generator from biased data.While the general pipeline for the method can be found in figure 1 , in the following, we first address regressing the generator on an a priori fixed hypothesis space H. Then we introduce our deep-learning method to either build a suitable space H, or even directly learn the eigenfunctions.Figure 1 : Pipeline of our method: from biased simulations to timescales and metastable states.Unbiasing generator regression Whenever π is absolutely continuous w.r.t. π ′ , the regularized energy kernel ( 9) satisfies the simple identityEQUATIONwhich, recalling the rightmost equation in (7) , implies that when the bias V and the diffusion coefficient β are known, the energy kernel can be empirically estimated through samples from π ′ via (10) . Moreover, when the potential U is known too, we can use (9) . Now, leveraging on (15) we directly obtain thatR ∂ (G) ≡ R ∂ (G)=E η x∼π ′ ∥χ η (x)-G T z(x)∥ 2 dπ dπ ′ (x) ≤ κ V E η x∼π ′ ∥χ η (x)-G T z(x)∥ 2 (16)where κ V = ess sup x∼π ′ dπ dπ ′ (x), which recalling ( 7) is finite whenever the bias V is essentially bounded. Therefore, in sharp contrast to transfer operator learning, whenever the true embedding χ η (x) 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 n = (x ′ i ) i∈[n]be the biased dataset generated from π ′ . Let w(x) = e βV (x) and define the empirical covariances w.r.t. the empirical distribution πEQUATIONCompute the eigenpairs(ν i , v i ) i∈[m] of the RR estimator G η,γ = ( W + ηγI) -1 C, and estimate the eigenpairs in (4) as( λ i , f i ) = (η -1/ν i , z(•) T v i ).If the elements of H and their gradients are essentially bounded, andlim m→∞ ρ(H) = 0, then for every ε > 0, there exist (m, n, γ) ∈ N × N × R + , such that, for every i ∈ [m], |λ i -λ i | ≤ ε and sin L 2 π (∢(f i , f i )) ≤ ε, with high probability.Note that, due to the form of the estimator, the normalizing constant w(x)dx does not need be computed. Moreover, relying on the upper bound in (16) we can alternatively compute C and W without the weights w and still ensure that the above result holds true.Neural network based learning Theorem 1 guarantees successful estimation of the eigendecomposition of the generator in (4) whenever the energy-based representation error ρ(H) in ( 12) is controlled. It is therefore natural to minimize ρ(H) by choosing an appropriate basis function z i 's.Inspired by the recent work [24] , we parameterize them by a neural network, and optimize them to span the leading invariant subspace of the generator.Let z θ = (z θ i ) i∈[m] : X → R m be a neural network (NN) embedding parameterized by θ ∈ Θ weights with continuously differentiable activation functions, and let λ θ i , i ∈ [m], be real non-positive (trainable) weights. We propose to optimize the NN to find the slowest time-scales λ θ i that solve the eigenvalue equationLz θ i = λ θ i z θ i , i ∈ [m]. Letting Z θ : R m → H η π (X ) be the (pa- rameterized) injection operator, given, for every u ∈ R m by Z θ u = i∈[m] z θ i u i , and denot- ing Λ η θ = (ηI -diag(λ θ 1 , . . . , λ θ m )) -1, the eigenvalue equations for the resolvent then become(ηI -L) -1 Z θ = Z θ Λ η θ .In other words, we aim to find the best rank-m decomposition of resolvent(ηI -L) -1 ≈ Z θ Λ η θ Z * θ .Therefore, for some hyperparameter α ≥ 0 we introduce the lossE α (θ) := ∥(ηI -L) -1 -Z θ Λ θ η Z * θ ∥ 2 HS(H η π ) -∥(ηI -L) -1 ∥ 2 HS(H η π ) + α i,j∈[m] (⟨z θ i , z θ j ⟩ L 2 π -δ i,j ) 2 .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 (ηI -L) -1 ∈ HS (H η π (X )) (see Appendix C for a discussion). The third term is optional; specifically, if the goal is not only to identify the proper invariant subspace of the generator (α = 0), but also to optimize the neural network to extract eigenfunctions as features, then this last term (α > 0) encourages the orthonormality of features in L 2 π (X ), an idea successfully exploited in machine learning and computational chemistry [see e.g. 24, and references therein].Recalling (14) and denoting by C θ and W θ the covariance matrices associated to the parameterized features, after some algebra, we obtain thatEQUATION)In turn, this can be estimated from biased data by two independent samples π ′ 1 and π ′ 2 asE π ′ 1 , π ′where C k θ and W k θ are the empirical covariances given by ( 17) for distribution π ′ k , while [2] . Importantly, the computational complexity of the loss ( 19) is of the order O(nm 2 d), where d is the state dimension and n the sample size, however it can be reduced to O(nmd) (see Appendix C) allowing its application to learn large dictionaries for high-dimensional problems with big amounts of (biased) data.w k = E x ′ ∼ π ′ k x w(x ′ ), k∈The following result, linked to controlling of the representation error as detailed in Theorem 1, provides theoretical guarantees for our approach. The proof and discussion are provided in Appendix B.EQUATIONwhere w = E x∼π ′ [w(x)]. Moreover, if α > 0 and λ m+2 < λ m+1 , then the equality holds if and only if This theorem provides a justification for minimizing the loss in (19) , which can be achieved by stochastic optimization algorithms, to obtain an approximation of either the leading invariant subspace of the resolvent (ηI -L) -1 (without orthonormality loss, i.e. α = 0), on which the estimator in Theorem 1 can be computed, or even the individual eigenpairs (α > 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 β are known, recalling (10) , 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 1 to compute the approximate eigendecomposition of L.(λ θ i , z θ i ) = (λ i , f i ) π-a. 
experiments: 
 In this section, we test the method described above on well-established [14, 9, 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 Muller-Brown potential which is a 2D potential, where this time, sampling is Compute loss E π ′ 1 , π ′ 2 (θ) using ( 19) and backpropagate 9:end while 10: end if 11: Compute C and W using (17) the datatset D n 12: Compute the eigenpairs(ν i , v i ) i∈[m] of G η,γ =( W+ηγI) -1 C 13: Output Estimated eigenpairs of L are ( λ i , f i )=(η-1/ν i , z θ (•) T v i ), i ∈ [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-like approach and showcase the efficiency of our approach, even with a small number of transitions in the training set. The codes used to train the models can be found in the following repository: https://github.com/DevergneTimothee/GenLearn One dimensional double well potential We first showcase the efficiency of our method on a 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 are shown Figure 7 in the appendix, where our method clearly outperforms transfer operator approaches and recovers the true underlying dynamics.Muller Brown potential with metadynamics biasing Muller Brown is a 2 dimensional potential presenting metastable states often studied in the context of enhanced sampling [19, 37, 14, 9] . It 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 [25] . The results of the training procedure are presented in Figure 2 . We compare the results with deepTICA 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 a bias potential. As expected, it is clearly 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.We next treat a more concrete molecular dynamics example with the study of the conformational change of alanine dipeptide in gas phase. It is a 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 a 30 dimensional system. This molecule has widely been used to test methods in enhanced sampling [9, 50, 42] : it presents a conformational change which is a rare event described by the angles ϕ and ψ. In the studies made on this system, the angle ϕ has been shown to be a 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 a non optimal CV. We generated biased dataset using a variation of metadynamics called on the fly probability enhanced sampling (OPES) [16] , which allows a more extensive and faster exploration of the state space than metadynamics [25] :Dataset 1: 800ns simulation, biasing on the ψ 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 two 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 1 mimics situations where one has only a basic prior knowledge of the system: only a "suboptimal" CV is used yielding to only a few transitions between the metastable states within the affordable simulation time. In order to ensure translational and rotational invariant vectors, we use Kabsch [18] algorithm, which has been used in previous studies [50, 4, 10] to transform the positions of the atoms. The results are presented in Figure 3 . Panels a) and b) display the first and second eigenfunctions learned by our method respectively. Notice that, even though only 2 transitions are 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 as above, but this time visualized on points from dataset 2 and in the plane of dihedral angles ϕ and θ. Interestingly, we discover that a linear relationship is present in the transition region, in agreement with recent findings in the molecular dynamics literature [6, 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 it on Dataset 2. One key quantity in molecular dynamics is the committor function for metastable states A and B, which is defined as 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, a 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.In this section, we report the results of our method obtained on a larger scale experiment: the folding/unfolding mechanism of the chignolin miniprotein. This system has c ) 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 text for more discussion). On all three panels, points are colored according to the value of the eigenfunction. d) Comparison of our method with the committor (q) of [19] extensively been studied [46, 19, 39, 7] . We first performed a 1 µs biased simulation using the deep-TDA collective variable [46, 37] 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 106µs unbiased trajectory provided by D.E. Shaw research [28] . The results are presented in figure 4 , showing a very good agreement between the training on an unbiased trajectory and on a biased one. 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 (Gly7O) and the distance between ASP3N and the oxygen atom of residue 8: THR (THR8) which allow visualizing the folded and unfolded states. 
conclusions: 
 We presented a method to 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 a recent generator learning approach based on Rayleigh quotients. In experiments, we observed that our approach is effective even when trained from suboptimal biased simulations. In the future our method could be applied to larger-scale simulations to discover rare events such as protein-ligand binding or catalytic processes. A main limitation of our method is that, in its current form, it is formulated for time-homogeneous bias potentials. However, the proposed framework could be naturally extended to time-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. 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? 
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 is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or 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.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a 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 is 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 to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] his paper is a 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 a 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 a 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 of use in the 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 a dataset. 
 
 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. In this work we consider stochastic differential equations (SDE) of the formEQUATIONThe special case of Langevin equation considered in the main body of the paper corresponds to a(x) = -∇U (x) and b(x) = 2 β I. Equation ( 21) describes the dynamics of the random vector X t in the state space X ⊆ R d , governed by the drift a : R d → R d and the diffusion b : R d → R d×p coefficients, where W t is a R p -dimensional standard Brownian motion. Under the usual conditions [see e.g. 33] that a and b are globally Lispchitz and sub-linear, the SDE (21) admits an unique strong solution X = (X t ) ⩾0 that is a Markov process to which we can associate the semigroup of Markov transfer operators (T t ) t≥0 defined, for every t ≥ 0, asEQUATION)For stable processes, the distribution of X t converges to an invariant measure π on X , such that X 0 ∼ π implies that X t ∼ π for all t ≥ 0. In such cases, one can define the semigroup on L 2 π (X ), and characterize the process by the infinitesimal generator of the semi-group (T t ) t⩾0 ,EQUATIONdefined on the Sobolev space H 1,2 π (X ) of functions in L 2 π (X ) whose gradient are in L 2 π (X ), too, i.e. L : L 2 π (X ) → L 2 π (X ) and dom(L) = H 1,2 π (X ). The transfer operator and the generator are linked to each other by the formula T t = exp(tL).After defining the infinitesimal generator for Markov processes by (23) , we provide its explicit form for solution processes of equations like (21) . Given a smooth function f ∈ C 2 (X , R), Itô's formula [see for instance 3, p. 495] provides for t ∈ R + ,f (X t ) -f (X 0 ) = t 0 d i=1 ∂ i f (X s )dX i s + 1 2 t 0 d i,j=1 ∂ 2 ij f (X s )d⟨X i , X j ⟩ s = t 0 ∇f (X s ) T dX s + 1 2 t 0 Tr X T s (∇ 2 f )(X s )X s ds.Recalling ( 21), we getEQUATIONProvided f and b are smooth enough, the expectation of the stochastic integral vanishes so that we getE[f (X t )|X 0 = x] = f (x) + t 0 E ∇f (X s ) T a(X s ) + 1 2 Tr b(X s ) T (∇ 2 f (X s ))b(X s ) X 0 = x dsRecalling that L = lim t→0 + (T t f -f )/t, we get for every x ∈ X ,EQUATIONwhich provides the closed formula for the IG associated with the solution process of (21) . In particular, for Langevin dynamics this reduces to (3).Next, recalling that for a bounded linear operator A on some Hilbert space H the resolvent set of the operator A is defined as ρ(A) = {λ ∈ C | A -λI is bijective}, and its spectrum Sp(A) = C\{ρ(A)}, let λ ⊆ Sp(A) be the isolated part of the spectra, i.e. both λ and µ = Sp(A) \ λ are closed in Sp(A). Then, the Riesz spectral projector P λ : H → H is defined byEQUATION)where Γ is any contour in the resolvent set Res(A) with λ in its interior and separating λ from µ. Indeed, we have that P 2 λ = P λ and H = Im(P λ ) ⊕ Ker(P λ ) where Im(P λ ) and Ker(P λ ) are both invariant under A, and we have Sp(A | Im(P λ ) ) = λ and Sp(A | Ker(P λ ) ) = µ. Moreover, P λ + P µ = I and P λ P µ = P µ P λ = 0.Finally if A is a compact operator, then the Riesz-Schauder theorem [see e.g. 38] assures that Sp(T ) is a discrete set having no limit points except possibly λ = 0. Moreover, for any nonzero λ ∈ Sp(T ), then λ is an eigenvalue (i.e. it belongs to the point spectrum) of finite multiplicity, and, hence, we can deduce the spectral decomposition in the formEQUATIONwhere the geometric multiplicity of λ, r λ = rank(P λ ), is bounded by the algebraic multiplicity of λ. If additionally A is a normal operator, i.e. AA * = A * A, then P λ = P * λ is an orthogonal projector for each λ ∈ Sp(A) and P λ = r λ i=1 ψ i ⊗ ψ i , where ψ i are normalized eigenfunctions of A corresponding to λ and r λ is both algebraic and geometric multiplicity of λ.We conclude this section by stating the well-known Davis-Kahan perturbation bound for eigenfunctions of self-adjoint compact operators.). Let A be compact self-adjoint operator on a separable Hilbert space H. Given a pair ( µ, f ) ∈ C × H such that ∥ f ∥ = 1, let λ be the eigenvalue of A that is closest to µ and let f be its normalized eigenfunction.If g = min{| µ -λ| | λ ∈ Sp(A) \ {λ}} > 0, then sin(∢( f , f )) ≤ ∥A f -µ f ∥/ g.In this section, we prove Theorem 1 relying on recently developed statistical theory of generator learning [21] . To that end, let ϕ(x) := z(•) T z(x) ∈ H be a feature map of the RKHS space H of dimension dim(H) = m. Let Wu j = σ 2 j u j be the eigenvalue decomposition of W and let v j := u j /σ j . This induces the SVD of the injection operator, Zu j = σ j zj for zj := z(•) T v j .Since H ⊆ H 1,∞ π (X ), we have thatc τ = ess sup x∼π j∈N [η|z j (x)| 2 -z j (x)[Lz j ](x)] < +∞,γ := W + ηγI tr[W -1 γ W γ ] ≤ m.In addition denote the empirical version of W γ as W γ := W + ηγI. Now, we can apply the following propositions from [21] to our setting, recalling the notation for normalizing constant w := E x∼π ′ [w(x)] for which w(•)/w = dπ/dπ ′ . Proposition 2. Given δ > 0, with probability in the i.i.d. draw of (x i ) n i=1 from π, it holds thatP{∥ W -W∥ ≤ ε n (δ)} ≥ 1 -δ,whereEQUATION)Proposition 3. Given δ > 0, with probability in the i.i.d. draw of (x i ) n i=1 from π, it holds thatEQUATIONwhereEQUATIONandL 1 (γ, δ) = ln 4 δ + ln tr(W -1 γ W) ∥W -1 γ W∥ .Moreover,EQUATION)Proposition 4. With probability in the i.i.d. draw of (x i ) n i=1 from π, it holdsP ∥W -1/2 γ ( C -C)∥ F ≤ ε 2 n (γ, δ) ≥ 1 -δ,whereEQUATION)We are now ready to prove Theorem 1, which we restate here for convenience.Theorem 1. Let D n = (x ′ i ) i∈[n]be the biased dataset generated from π ′ . Let w(x) = e βV (x) and define the empirical covariances w.r.t. the empirical distribution πEQUATION)Compute the eigenpairs (ν i , v i ) i∈[m] of the RR estimator G η,γ = ( W + ηγI) -1 C, and estimate the eigenpairs in (4) as( λ i , f i ) = (η -1/ν i , z(•) T v i ).If the elements of H and their gradients are essentially bounded, and lim m→∞ ρ(H) = 0, then for every ε > 0, there exist (m, n, γ)∈ N × N × R + , such that, for every i ∈ [m], |λ i -λ i | ≤ ε and sin L 2 π (∢(f i , f i )) ≤ ε, with high probability.Proof. We first show that R ∂ ( G η,γ ) < ε for big enough m, n ∈ N and small enough γ > 0.Observe that W = E[ W]/w and C = E[ C]/w (33) and, henceG η,γ := W -1 γ C = (E[ W γ ]) -1 (E[ C]), due to cancellation of w. Given ε > 0, let m ∈ N be such that ρ(H) = ∥(I -P H )(ηI -L) -1 ∥ 2 HS(H,H η π ) < ε/3. Next, since P H (ηI -L) -1 Z-Z G η,γ =(I-Z W -1 γ Z * )(ηI -L) -1 Z = Z(W † C -W -1 γ C), we have that ∥P H (ηI-L) -1 Z-Z G η,γ ∥ HS(H,H η π ) = ∥W 1/2 (W † C-W -1 γ C)∥ F = ∥W 1/2 (W x † -W -1 γ )C∥ F → 0, as γ → 0. Hence, let γ > 0 be such that ∥P H (ηI -L) -1 Z-Z G η,γ ∥ H→H η π < ε/3. Finally, using the decomposition of the risk R ∂ ( G η,γ ) ≤ ρ(H)+∥P H (ηI -L) -1 Z-Z G η,γ ∥ HS(R m ,H η π ) +∥Z( G η,γ -G η,γ )∥ HS(R m ,H η π )it remains to show that for large enough n we have∥Z(G η,γ -G η,γ )∥ HS(R m ,H η π ) ≤ ε/3. To that end observe that W 1/2 γ ( G η,γ -G η,γ ) = W 1/2 γ W -1 γ ( C -W γ W -1 γ C ± C) = W 1/2 γ W -1 γ W 1/2 γ W -1/2 γ ( C -C) -W -1/2 γ ( W -W)W -1/2 γ (W -1/2 γ C) .Thus, by multiplying the above expression by w and applying Propositions 3 and 4, we obtain that there exists n ∈ N such that ∥Z(G η,γ -G η,γ )∥ HS(R m ,H η π ) ≤ ε/3. Next, assuming that ∥ W -W∥ is small, for the normalization of the estimated eigenfunctions we have that∥v j ∥ 2 2 ∥ f j ∥ 2 H η π = v T j v j v T j Wv j ≤ v T j v j v T j Wv j -v T j (W -W)v j ≤ 1 λ + min ( W) -∥ W -W∥ ≤ 1 λ m (W) -2∥ W -W∥ .where we have that λ m (W) > 0 due to fact that (z j ) are linearly independent.Therefore, to conclude the proof, we apply [21, Proposition 2] which directly relying on Proposition 1 yields the result.At last we remark, based on the observation that W = E[ W]/w and C = E[ C]/w, one can readily obtain stronger version of Theorem 1 in the general RKHS setting of [21] .In this section, we provide details on our DNN method and prove Theorem 2. To that end, let us denote the terms in the loss asEQUATIONand recall that the injection operator is Z θ =(z θ ) T (•) : R m → H η π (X ). It is easy to show from the definition of the adjoint that, for every f ∈ H η π (X ), we haveEQUATIONThus,EQUATIONwherew = E x∼π ′ [w(x)]. Moreover, if α > 0 and λ m+2 < λ m+1 , then the equality holds if and only if Proof. Let P k : H η π (X ) → H η π (X ) be spectral projector of L corresponding to the k largest eigenvalues. Now, consider(λ θ i , z θ i ) = (λ i , f i ) π-a.E k (θ)=∥P k (ηI -L) -1 -Z θ Λ θ η Z * θ ∥ 2 HS(H η π ) -∥P k (ηI -L) -1 ∥ 2 HS(H η π ). Due to Eckhart-Young theorem, we have that for every θ ∈ Θ, the best rank-m approximation of (ηI -L) -1 is (ηI -L) -1 P m = (ηI -L) -1 P k P m , for k > m, and it holdsE k (θ)≥ k j=m+1 (η -λ i ) -1 - k j=1 (η -λ i ) -1 = - m j=1 (η -λ i ) -1 .As before, expanding in E k the HS norm via the trace, we obtainE k (θ)=∥Z θ Λ θ η Z * θ ∥ 2 HS(H η π ) -2 tr[Z θ Λ θ η Z * θ (ηI -L) -1 P k ]=E(θ)+2 tr[Z θ Λ θ η Z * θ (ηI -L) -1 (I-P k )],and, hence, by Cauchy-Schwartz inequality, we have that|E(θ) -E k (θ)| ≤ ∥Z θ ∥ 2 HS(H,H η π ) ∥Λ θ η ∥∥(ηI -L) -1 (I-P k )∥ = 1 η i∈[m] E η π [z θ i ](η-λ k+1 ) -1 .Observing that 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 Pz θ i ∈ H η π (X ), i.e. E η π [z θ i ] < ∞, we conclude that, for every θ ∈ Θ, lim k→∞ E k (θ) = E(θ). Therefore, noting that E α (θ) ≥ E(θ), inequality in (m (ηI -L) -1 , i.e. (z θ i ) i∈[m] span the leading invariant subspace (f i ) i∈[m] of the generator. So, if E α (θ) = E(θ) = m j=1 (η -λ i ) -1and α > 0, we have that E on (θ) = 0, implying that (z i ) i∈[m] is an orthonormal basis, and, henceP m (ηI -L) -1 , i.e. (z θ i ) i∈[m] = Z θ Λ η θ Z * θ . The result follows. To show that E γ (θ) = E E π ′ 1 , π ′ 2 γ (θ) /w 2 ,we rewrite (18) to encounter the distribution change, noting that the empirical covariances are reweighted but not normaliezed by w. So, we have thatE γ (θ) = tr w -2 E[ C θ ]Λ η θ E[ W θ ]Λ η θ -2w -1 E[C θ ]Λ η θ + w -2 α(E[ C θ ] -wI) 2 .But, sinceE x∼π [f (x)g(x)]=E x ′ ∼π ′ dπ dπ ′ (x ′ ) f (x ′ ) g(x ′ ) =E x ′ 1 ∼π ′ dπ dπ ′ (x ′ 1 )f (x ′ 1 ) E x ′ 2 ∼π ′ dπ dπ ′ (x ′ 2 )g(x ′ 2 ) ,where x ′ 1 and x ′ 2 are two independent r.v. with a law π ′ , the proof is completed. It is interesting to note that during the training process, the loss reaches progressively lower plateaus. This is due to the fact that the NN has found a novel eigenfunction orthogonal to the previously found ones, starting with the constant one. Then during the plateau phase, the subspace is being explored 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 was observed in all the training sessions. This nice property is a handy tool in properly optimizing the loss and understanding the proper stopping time.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 1. In Figure 6 , we illustrate such situation on a simulation of the Muller Brown 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. In all of our experiments we used the hyperbolic tangent activation function. This choice was made because it is a widely used, bounded function with continuous derivative. It thus satisfies all the criteria needed for this method. Finally, when looking for m eigenpairs, instead of having a neural network with m outputs, we choose to have m neural networks with one output.Besides common hyperparameters such as learning rate, neural network architecture and activation function, our method requires only two hyperparameters: η and α. Other methods such as [50] do not require η, but on the other hand requires one weight per searched eigenfunction.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 a workstation with a 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 a training and a validation dataset. The proportion were set to 80% for training and 20% for validation. The training of deepTICA models was performed using the mlcolvar package [8] .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 U tg (x) = 4(-1.5 exp(-80x 2 ) + x 8 ), which has a form of two wells separated by a high barrier which can hardly be crossed during a 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: U sim = 4(-0.5 exp(-80x 2 ) + x 8 ), which thus makes a bias potential: V bias (x) = U sim (x) -U tg (x) = -4(exp(-80x 2 )). 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 deepTICA [7] which is a state-of-the-art method for molecular dynamics simulations. For this experiment we have used a Gaussian kernel of lengthscale 0.1, η = 0.1 and a regularization parameter of 10 -5For this experiment, the dataset was generated using an in-house code implementing the Euler-Maruyama scheme to discretize the overdamped Langevin equation. The simulation was performed at The bias potential is built according to the following equationEQUATIONwhere the centers x i are built on the fly: every 500 timestep one more center is added to this list. This kind of bias potential is called metadynamics [25] and allows reducing the height of the barrier for a better exploration of the space. In order to have a static potential, no more centers are added after 300000 timesteps. We use a learning rate of 5.10 -3 , the architecture of the neural network used is a multilayer perceptron with layers of size 2 (inputs), 20, 20 and 1. The parameter η was chosen to be 0.05.Simulation details All the simulations are run 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 a timestep of 2fs with a damping coefficent γ i = m i /0.05 at a target temperture of 300K. For both dataset, in order to make proper comparison, we used the explore version of OPES [17] , with a barrier parameter of 20 kJ/mol and a pace of deposition of 500 timesteps.Neural networks training For our neural networks training, we assume overdamped Langevin dynamics with the same value of the friction coefficient as in the simulation, as done in other works [50] . We use a learning rate of 10 -3 , and the architecture of the neural network used is a multilayer perceptron with layers of size 30 (inputs), 20, 20 and 1. The parameter η was chosen to be 0.1.Simulation details The simulations are run 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, we 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 deepTDA model used for biasing the simulation is the one obtained in ref [46] Neural networks training For our neural networks training, we assume overdamped Langevin dynamics. We use a learning rate of 5.10 -4 , and the architecture of the neural network used is a multilayer perceptron with layers of size 210 (inputs), 50, 50 and 1. The parameter η was chosen to be 0.2.Justification: We enjoy research and we respect other people work.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS 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., if there is a special consideration due to laws or regulations in their jurisdiction).38th Conference on Neural Information Processing Systems (NeurIPS 2024). 


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: 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.
questions: 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}$? 
I 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.

- 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?

- Potential typos: Page 2, line 67: "principle" -> "principled".
Page 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)$.
Page 8, line 282, "approqch"

- 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: Limitations are discussed.
paper_summary: paper_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.
strength_weakness: strength_weakness: I find it difficult to read the paper and follow it, perhaps because I am too far from the field
questions: 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: limitations: Authors discuss correctly the limitations
paper_summary: paper_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.

The 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.
strength_weakness: strength_weakness: My points contain many understanding questions and some "criticizing" questions. It would be great if you could also answer the understanding questions.
1. 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. 
	1. 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?
	2. **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?
	3. Why do you omit the deep learning baselines from the ALDP experiments?
2. 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?
3. 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?
Minor:
1. 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? 
   Once you trained your neural network, how do you obtain the eigenfunctions from it.
2. 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?
questions: questions: I would appreciate any time you can take to answer some of my questions.
1.  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"?
2. 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: 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.

Comments: title: Acknowledgement of rebuttal by authors
Acknowledgement to the reviewer
Acknowledgement of reviewer's comments
Further discussion
comment: Thank you for the detailed response. They are very clear.
We are happy that our replies were helpful. We thank reviewer for all their comments and discussions, which I will improve our paper. We commit to incorporate them all in the revision.
Thank you very much for confirming some points, and for addressing the reference. I would like to take the freedom to raise my score from 6 to 8 and to increase my confidence.
We thank the reviewer for appreciating our rebuttal and suggesting our paper for acceptance. In the following we address the reviewer's additional questions.

- __Concerning the reference.__ Thanks for bringing this paper to our attention, we will include it in the revision. In this work the authors build their model from D.E. Shaw dataset formed by a long trajectory containing already all the necessary information for training, without the need for biasing. They learn a transition kernel (i.e. a conditional probability to go from $X_t$ to $X_{t+\Delta t}$), which, as discussed in Sec. 3, has an inherent difficulty to be adapted to biased simulations. 

- __Scaling of generator algorithms.__ First, note that transfer operator (TO) based methods apply only to equally spaced data and that the sampling frequency  $1/\Delta t$ must be high enough to distinguish all relevant time-scales in the dynamics. Otherwise, since TO eigenvalues are $e^{\lambda_i \Delta t}$, small spectral gaps complicate learning (see Thm. 3 [23]). Conversely, our IG method, which uses gradient information, is time-scale independent, handles irregularly spaced measurements, and does not rely on time discretizations. This important aspect allows one to learn from biased simulations without notorious time-lag bottlenecks inherent to TO methods.  Alas, as there is no free lunch,  this incurs higher computational complexity. However, his complexity is to an extent mitigated through our representation learning, by exploiting automatic differentiation tools in deep learning.  

- __IG’s eigenfunctions.__ Indeed, the learned eigenfunctions of the IG can be used as CVs, which in fact are optimal CVs, as discussed in the general reply, see  [A]. Moreover, the leading eigenfunctions of IG encode the true dynamics. Hence, once they are properly learned, no more biasing is needed. For example, they can be used to infer the transition mechanism, see for instance figure 2, c) in the alanine dipeptide experiment, where we show that even with few transitions, we manage to recover a linear relationship between $\theta$ and $\phi$. Another important aspect is that one can build good approximations of all transfer operators $\cal{T}_{\Delta t}$, $\Delta t>0$, from the leading eigenpairs of IG, enabling forecasting of system’s observables see e.g. [21]. 

If there are any remaining concerns and/or questions, we are happy to discuss more.
Thank you for these nice explanations! I fully agree that studying and recovering dynamics from biased simulations is an important problem and appreciate any ML work done for it. I think the paper should be accepted.


I appreciate the chignolin experiment.
1. Minor: Could you put your work a bit into context with "Implicit Transfer Operator Learning: Multiple Time-Resolution Surrogates for Molecular Dynamics" https://arxiv.org/abs/2305.18046 . They only learn the transfer operator. However, they also provide experiments on Chignolin. Are there inherent reasons why learning the generator might scale worse than learning the operator?

My main concern (maybe it was naive) was in thinking that CVs (which you are learning by learning the generator and using its eigenfunctions) main or only value is in using them for further biased simulations. Thus, the question of why one biased simulation with your CVs would be better than the first biased simulation. This also seems valuable, and you suggest it in the general rebuttal, but do not try it. \ 
Important:\
However, as far as I understand, you suggest that I was wrong in assuming that learning the generator is only valuable for obtaining CVs and using them for biasing simulations. Learning them is also valuable to e.g. recover the free energy surface from a biased simulation or to recover dynamical quantities about transitions - please let me know if this understanding is correct. 

(sorry for the delayed reply - I will be quicker to respond in the remaining days)
We thank the reviewer for the prompt reply.
- __Q1__ We agree that Galerkin projection interpretation seems more convenient from an expository standpoint. So, based on the reviewer’s feedback, we propose to modify the content in lines 185-191, by briefly explaining why Galerkin projection should be done in energy space $\cal{H}^\eta_\pi$ instead of $\cal{L}^2_\pi$ space. Then, since the statistical learning risk formulation is the key to understand how to generalizes w.r.t. (unobserved) true distribution $\pi$ from biased simulations, we briefly discuss its equivalence with the Galerkin viewpoint. We believe that this change can help the reader to better grasp the main idea behind the energy, and connect it to the material that follows. Does the reviewer find this proposal adequate?

- __Q2: about the loss:__ To answer this question, let us compare (19) to the two losses used in the most related works, [24] and [46]. While the former one (only tested on unbiased dynamics of a toy system) depends on one hyperparameter, the latter loss has $1+m$ hyperparameters, where $m$ is the latent dimension. Importantly, __both losses suffer from statistically biased estimation__ of gradients when optimized over the sample distributions, which may negatively impact the optimization.  Moreover, as reported in [46] (see also lines 287-289) additional $m$ hyperparameters in loss of [46] are delicate to tune. __In sharp contrast__, our empirical loss, as stated in Theorem 2, __is an unbiased estimate__ of the true loss in Equation (18), and we didn’t experience any particular difficulty in tuning our only hyperparameter $\alpha>0$. Note also that, in principle $\alpha>0$ is optional to speed up the training, and we may as well use $\alpha=0$. Concerning the computation, we respectfully disagree that the loss is hard to compute. Indeed, for two size $n$ batches of data, it just relies on computation of covariances of features $C$ and their gradients $W$, which is of the order $\cal{O}(n m^2 d)$. Since in the molecular dynamics setting $m$ is typically small, computing the loss is very efficient, similarly to losses in self-supervised learning based on canonical correlations, see e.g. [A]. The only computational bottleneck lies in using the gradients for $W$. However, we believe that this is a necessary price to pay to be able to generalize properly from biased simulations, as motivated in Sec. 3. Finally, we remark that optimization of our loss follows interpretable dynamics, see Figure 3 of the Appendix where plateaus reveal discovery of new relevant features, helping practitioners to decide when to stop the training.   

- __Q2: about the contribution:__ As the reviewer rightfully recognizes, introducing loss (19) for representation learning is an important contribution of this work. However, our main focus is on __biased simulations__, or, in other words,  __learning dynamics that ML algorithm did not observe__. With this in mind, our main contributions are summarized in lines 67-74. For the reviewer's convenience we rephrase and slightly expand them here. To the best of our knowledge, for the first time 
- __(1)__ we __formalize the problem of learning__ the spectral decomposition of the infinitesimal generator __from biased simulations__ as a regression, 
- __(2)__ we __derive empirical estimators with guaranteed statistical consistency__ (Theorem 1), and 
- __(3)__ we propose an end-to-end __representation learning + regression__ pipeline that efficiently scales to large problems. 
- __(4)__ we __empirically validate our methodology__ on molecular dynamics datasets of increasing complexity.  That said, __we will follow reviewer’s suggestion and emphasize better the important role of representation learning__ for the generator regression.  

Once more we thank the reviewer, and if any other doubt remains or clarification is needed, we are happy to discuss in more detail.

Ref.

[A] Balestriero et al., A Cookbook of Self-Supervised Learning. Arxiv:2304.12210 (2023).
Thank you for the detailed response. I have some follow-up points.

- **Q1**. I now agree that the estimator is a regularized Galerkin projection of the resolvent onto the space of basis functions, using the energy inner product. This interpretation seems more convenient from an expository standpoint. It would be helpful if the authors can incorporate this perspective into the paper and clarify the writing. In my experience, the current mathematical exposition, with its various operators, norms and formulas, did not provide a very pleasant reading experience.

- **Q2**. Yes, it is Figure 1. Thank you. The explanation and results are interesting. Do the authors have any insights regarding the superiority of the loss (19) for representation learning? It seems to me the loss is not that simple and easy to optimize or compute. There is a penalization term to make the $z_i$ normalized and one has the hyper parameter $\alpha$ to tune. Furthermore, based on the authors' response, do I understand correctly that the primary improvement stems from the loss function (19) for representation learning, which is more significant than the regression formulation proposed in the paper? If so I think these insights should be added to the paper.
I really appreciated the authors' explanations, including the new figures (pipeline and the experiment with the proteins). I will change my rating to accept.