# From Biased to Unbiased Dynamics:
An Infinitesimal Generator Approach

Timothee Devergne
Atomistic Simulations
Computational Statistics and Machine Learning
Italian Institute of Technology Vladimir R. Kostic
Computational Statistics and Machine Learning, Italian Institute of Technology
Department of Mathematics and Informatics, University of Novi Sad Michele
Parrinello
Atomistic Simulations, Italian Institute of Technology Massimiliano Pontil
Computational Statistics and Machine Learning, Italian Institute of Technology
Department of Computer Science, University College London

###### 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.

##  1 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. 44, 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}=-\nabla
U(X_{t})dt+\sigma 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 $\sigma$ 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 [33]. 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, 42, 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, 28].
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, 40, 14, 9, 10, 4, 27]. A key idea is to use available
dynamical information to construct the CVs [39, 30, 7, 40]. 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 [45]. 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, 47, 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.

##  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, [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.
32], 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 equation

$dX_{t}=-\nabla U(X_{t})dt+\sqrt{2\beta^{-1}}dW_{t}\quad\text{and}\quad
X_{0}=x,$ (1)

describing dynamics in a (state) space $\mathcal{X}\subseteq\mathbb{R}^{d}$,
governed by the potential $V:\mathbb{R}^{d}\to\mathbb{R}$ at the temperature
$\beta^{-1}=k_{B}T$, where $W_{t}$ is a $\mathbb{R}^{d}$-dimensional standard
Brownian motion.

The SDE (1) admits a unique strong solution $X=(X_{t})_{\geqslant 0}$ that is
a Markov process to which we can associate the semigroup of Markov transfer
operators $(\mathcal{T}_{t})_{t\geq 0}$ defined, for every $t\geq 0$, as

$[\mathcal{T}_{t}f](x):=\mathbb{E}[f(X_{t})|X_{0}=x],\;\;x\in\mathcal{X},\,f\colon\mathcal{X}\to\mathbb{R}.$
(2)

For (1) the distribution of $X_{t}$ converges to the invariant measure $\pi$
on $\mathcal{X}$ called the Boltzmann distribution, given by $\pi(dx)\propto
e^{-\beta V(x)}dx$. In such cases, one can define the semigroup on
$L^{2}_{\pi}(\mathcal{X})$, and characterize the process by the infinitesimal
generator

$\mathcal{L}:=\lim_{t\to 0^{+}}\tfrac{\mathcal{T}_{t}-I}{t}$

defined on the Sobolev space $H^{1,2}_{\pi}(\mathcal{X})$ of functions in
$L^{2}_{\pi}(\mathcal{X})$ whose gradient are also in
$L^{2}_{\pi}(\mathcal{X})$. The transfer operator and the generator are linked
one to another by the formula $\mathcal{T}_{t}=\exp(t\mathcal{L})$. Moreover,
it can be shown (see Appendix A) that the generator $\mathcal{L}$ acts on
$f\colon\mathcal{X}\to\mathbb{R}$ as

$\mathcal{L}f=-{\langle\nabla U,\nabla f\rangle}+\beta^{-1}\Delta f,$ (3)

which, integrating by parts, gives
$\int\\!(\mathcal{L}f)g\,d\pi{=}{-}\beta^{-1}\\!\\!\int\\!\langle\nabla
f,\nabla g\rangle\,d\pi{=}\int\\!\\!f(\mathcal{L}g)d\pi$, showing that
$\mathcal{L}$ is self-adjoint. If $\mathcal{L}$ has only a discrete spectrum,
one can solve (1) by computing the spectral decomposition

$\mathcal{L}=\textstyle{\sum_{i\in\mathbb{N}}}\lambda_{i}f_{i}\otimes f_{i},$
(4)

Using (2) and the exponential relation between the transfer operator and the
generator, one can thus write

$[\mathcal{T}_{t}f](x):=\mathbb{E}[f(X_{t})|X_{0}=x]=\textstyle{\sum_{i\in\mathbb{N}}}{e^{t\lambda_{i}}f_{i}(x){\langle
f_{i},f\rangle}},\;\;x\in\mathcal{X},\,f\colon\mathcal{X}\to\mathbb{R}$ (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 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}:\mathcal{X}\rightarrow\mathbb{R}$,
$j\in[m]$,

$\mathcal{H}:=\Big{\\{}h_{u}=\textstyle{\sum_{j\in[m]}}u_{j}z_{j}\,\big{|}\,u\,{=}\,(u_{1},\dots,u_{m})\in\mathbb{R}^{m}\Big{\\}}.$
(6)

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, 29, 47]. The space $\mathcal{H}$ is naturally equipped with the
geometry induced by the norm
$\|h_{u}\|^{2}_{\mathcal{H}}\,{:=}\,{\sum_{j=1}^{m}u_{j}^{2}}.$ Moreover,
every operator $A\colon\mathcal{H}\to\mathcal{H}$ can be identified with
matrix $\textsc{A}\in\mathbb{R}^{m\times m}$ by
$Ah_{u}=z(\cdot)^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}\textsc{A}u$.
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 $\mathcal{T}_{t}$
can be simply seen as the vector-valued regression problem [22], in which the
action of $\mathcal{T}_{t}\colon L^{2}_{\pi}(\mathcal{X})\to
L^{2}_{\pi}(\mathcal{X})$ on the domain $\mathcal{H}\subseteq
L^{2}_{\pi}(\mathcal{X})$ is estimated by an operator
$\widehat{T}_{t}\colon\mathcal{H}\to\mathcal{H}$. This aims to minimize the
mean square error (MSE) w.r.t. the invariant distribution. Given a dataset
$\mathcal{D}_{n}:=(x_{i},y_{i}=x_{i+1})_{i=1}^{n}$ 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
$\widehat{\textsc{T}}_{\gamma}{:=}\widehat{\textsc{C}}_{\gamma}^{-1}\widehat{\textsc{C}}_{t}$,
where the empirical covariance matrices are
$\widehat{\textsc{C}}\,{=}\,\textstyle{\tfrac{1}{n}\sum_{i\in[n]}}\,z(x_{i})z(x_{i})^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}$
and
$\widehat{\textsc{C}}_{t}\,{=}\,\textstyle{\tfrac{1}{n}\sum_{i\in[n]}}\,z(x_{i})z(y_{i})^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}$.
We then estimate the eigenpairs $(\lambda_{i},f_{i})$ in (4) by the eigenpairs
$(\widehat{\mu}_{i},\widehat{u}_{i})$ of $\widehat{\textsc{T}}_{\gamma}$ as
$\widehat{\lambda}_{i}{:=}\ln(\widehat{\mu}_{i}/t)$ and
$\widehat{f}_{i}{:=}z(\cdot)^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}\widehat{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.

##  3 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^{\prime}(x):=U(x)+V(x),\;\;x\in\mathcal{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 it challenging to recover the unbiased dynamics from biased data.
Denoting the invariant measure of the perturbed process by $\pi^{\prime}$ and
its generator by $\mathcal{L}^{\prime}\colon
H^{1,2}_{\pi^{\prime}}(\mathcal{X})\to H^{1,2}_{\pi^{\prime}}(\mathcal{X})$,
our principal objective is thus to:

  * Gather data from simulations generated by $\mathcal{L}^{\prime}$ to learn the spectral decomposition
of the unperturbed generator $\mathcal{L}$.

To tackle this problem, we note that since the eigenfunctions of the generator
$\mathcal{L}$ are also eigenfunctions of every transfer operator
$\mathcal{T}_{t}=e^{t\mathcal{L}}$, 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
$\mathcal{L}^{\prime}=\mathcal{L}+{\langle\nabla V,\nabla(\cdot)\rangle},$
this knowledge is not easily transferred to the perturbation of the transfer
operator. Indeed, recalling that
$\mathcal{T}:=\mathcal{T}_{1}=e^{\mathcal{L}}$, and since the differential
operator ${\langle\nabla V,\nabla(\cdot)\rangle}$ in general does not share
the same eigenstructure of $\mathcal{L}$, one has that

$\mathcal{T}^{\prime}:=e^{\mathcal{L}^{\prime}}=e^{\mathcal{L}-{\langle\nabla
V,\nabla(\cdot)\rangle}}\neq\mathcal{T}e^{-{\langle\nabla
V,\nabla(\cdot)\rangle}}.$

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

$\pi(dx)=\frac{e^{-\beta U(x)}dx}{\int e^{-\beta
U(x)}dx},\;\pi^{\prime}(dx)=\frac{e^{-\beta U^{\prime}(x)}dx}{\int e^{-\beta
U^{\prime}(x)}dx}\;\text{ and }\;\frac{d\pi}{d\pi^{\prime}}(x)=\frac{e^{\beta
V(x)}}{\int e^{\beta V(x)}\pi^{\prime}(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
features

$\textsc{C}\,{=}\mathbb{E}_{X^{\prime}\sim\pi^{\prime}}\left[\tfrac{d\pi}{d\pi^{\prime}}(X^{\prime})z(X^{\prime})z(X^{\prime})^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}\right].$
(8)

However, since the transition kernel of the process $(X_{t}^{\prime})_{t\geq
0}$ generated by $\mathcal{L}^{\prime}$ is different from that of the original
process, the above reasoning does not hold for the cross-covariance matrix,
that is,

$\textsc{C}_{t}\,{:=}\mathbb{E}_{X_{0}\sim\pi^{\prime}}\left[\tfrac{d\pi}{d\pi^{\prime}}(X_{0})\,z(X_{0})z(X_{t})^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}\right]\neq\mathbb{E}_{X_{0}^{\prime}\sim\pi^{\prime}}\left[\tfrac{d\pi}{d\pi^{\prime}}(X_{0}^{\prime})\,z(X_{0}^{\prime})z(X_{t}^{\prime})^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}\right]=:\textsc{C}^{\prime}_{t},$

Consequently, the estimator $\widehat{T}_{t}$ obtained by minimizing the
reweighed risk functional
$\mathcal{R}^{\prime}(\widehat{T}_{t})\\!:=\\!\mathbb{E}_{X_{0}\sim\pi^{\prime}}\big{[}\,\tfrac{d\pi}{d\pi^{\prime}}(X_{0})\,\lVert
z(X_{t}^{\prime})\\!-\\!\widehat{\textsc{T}}_{t}^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}z(X_{0})\rVert^{2}_{2}\,\big{]}$
does not minimize the true risk since
$\mathcal{R}^{\prime}(\widehat{T}_{t})\\!\neq\\!\mathcal{R}(\widehat{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, 30, 46]. 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 $\textsc{C}_{t}$
by $\textsc{C}^{\prime}_{t}$ and minimize
$\mathcal{R}^{\prime}(T)\approx\mathcal{R}(T)$ over
$T\colon\mathcal{H}\to\mathcal{H}$.

##  4 Infinitesimal generator learning

In this section, we address generator learning. While there has been
significant progress on this topic [24, 15, 34, 47, 20], we follow the recent
approach in [21] for learning the generator $\mathcal{L}$ on an a priori fixed
hypothesis space $\mathcal{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 $(\eta I\\!-\\!\mathcal{L})^{\\!-\\!1}$, where $\eta\,{>}\,0$ is
a prescribed parameter. To this end, we observe that the action of the
resolvent in $\mathcal{H}$ can be expressed as $((\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}h_{u})(x)\,{=}\,\chi_{\eta}(x)^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}u$,
where $\chi_{\eta}$ is the embedding of the resolvent $(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}$ into $\mathcal{H}$, given by
$\chi_{\eta}(x)\,{=}\,\int_{0}^{\infty}\mathbb{E}[z(X_{t})e^{-\eta
t}\,|\,X_{0}{=}x]dt,$ $x\,{\in}\,\mathcal{X}$, see [21]. We then aim to
approximate $\chi_{\eta}(x)\approx\textsc{G}^{*}z(x)$ by a matrix
$\textsc{G}\in\mathbb{R}^{m\times 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
$\mathfrak{E}^{\eta}_{\pi}\colon
H^{1,2}_{\pi}(\mathcal{X})\,{\times}\,H^{1,2}_{\pi}(\mathcal{X})\,{\to}\,\mathbb{R}$,
given by
$\mathfrak{E}^{\eta}_{\pi}[f,g]\,{=}\,\mathbb{E}_{x\sim\pi}\left[\eta{f(x)g(x)}\,{-}\,f(x)[\mathcal{L}g](x)\right]$,
which using (3) becomes

$\mathfrak{E}^{\eta}_{\pi}[f,g]{=}\mathbb{E}_{x\sim\pi}\left[\eta
f(x)g(x){+}f(x)\nabla
U(x)^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}\nabla
g(x){-}\tfrac{1}{\beta}f(x)\Delta g(x)\right],$ (9)

and, due to the identity $\int f\mathcal{L}gd\pi=-\beta^{-1}\int(\nabla
f)^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}(\nabla g)d\pi$, also

$\mathfrak{E}^{\eta}_{\pi}[f,g]{=}\mathbb{E}_{x\sim\pi}\left[\eta
f(x)g(x){+}\tfrac{1}{\beta}\textstyle{\sum_{k\in[d]}}\partial_{k}f(x)\partial_{k}g(x)\right].$
(10)

Since $\mathcal{L}$ is negative semi-definite, the above kernel induces the
regularized squared energy norm $\mathfrak{E}^{\eta}_{\pi}\colon
H^{1,2}_{\pi}(\mathcal{X})\to[0,+\infty)$ by
$\mathfrak{E}^{\eta}_{\pi}[f]:=\mathfrak{E}^{\eta}_{\pi}[f,f]=\mathbb{E}_{x\sim\pi}\left[\eta{f^{2}(x)}-f(x)[\mathcal{L}f](x)\right]$.
It counteracts the resolvent and balances the transient dynamics (energy) of
the process with the invariant distribution $\pi$. In a nutshell, instead of
using the mean square error of
$f(x):=\lVert\chi_{\eta}(x)-\textsc{G}^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}z(x)\rVert_{2}$
to define the risk, we "fight fire with fire" and penalize the energy to
formulate the generator regression problem

$\min_{G\colon\mathcal{H}\to\mathcal{H}}\mathcal{R}_{\partial}(G)\equiv\mathcal{R}_{\partial}(\textsc{G}):=\mathfrak{E}^{\eta}_{\pi}\big{[}\lVert\chi_{\eta}(\cdot)\\!-\\!\widehat{\textsc{G}}^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}z(\cdot)\rVert_{2}\big{]}.$
(11)

Indeed, this risk overcomes the difficulty of not knowing $\chi_{\eta}$. To
show this, let us define the space
$\mathcal{H}^{\eta}_{\pi}(\mathcal{X}):=\\{f\in
H^{1,2}_{\pi}(\mathcal{X})\,|\,\mathfrak{E}^{\eta}_{\pi}[f]<\infty\\}$
associated to the energy norm $\lVert
f\rVert_{\mathcal{H}^{\eta}_{\pi}}:=\sqrt{\mathfrak{E}^{\eta}_{\pi}[f]}$, and
recalling that the operator $G\colon\mathcal{H}\to\mathcal{H}$ is identified
with a matrix $\textsc{G}\in\mathbb{R}^{m\times m}$ via
$Gh_{u}=z(\cdot)^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}(\textsc{G}u)$,
define the (injection) operator
$\mathcal{Z}\colon\mathbb{R}^{m}\to\mathcal{H}^{\eta}_{\pi}$ by
$\mathcal{Z}u=z(\cdot)^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}u$,
for every $u\in\mathbb{R}^{m}$. Then, since
${\rm{HS}}\left(\mathbb{R}^{m},\mathcal{H}^{\eta}_{\pi}\right)\equiv{\rm{HS}}\left(\mathcal{H},\mathcal{H}^{\eta}_{\pi}\right)$,
the norm is the sum of squared $\mathcal{H}^{\eta}_{\pi}$ norm over the
standard basis in $\mathbb{R}^{m}$, and one obtains

$\displaystyle\mathcal{R}_{\partial}(G)$ $\displaystyle=$
$\displaystyle\lVert(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}-G\rVert^{2}_{{\rm{HS}}\left(\mathcal{H},\mathcal{H}^{\eta}_{\pi}\right)}$
(12) $\displaystyle=$ $\displaystyle\underbrace{\lVert P_{\mathcal{H}}(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}-G\rVert^{2}_{{\rm{HS}}\left(\mathcal{H},\mathcal{H}^{\eta}_{\pi}\right)}}_{\text{projected
problem}}+\underbrace{\lVert(I-P_{\mathcal{H}})(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}\rVert^{2}_{{\rm{HS}}\left(\mathcal{H},\mathcal{H}^{\eta}_{\pi}\right)}}_{\text{
representation error }\rho(\mathcal{H})},$

where $P_{\mathcal{H}}$ is the orthogonal projector in
$\mathcal{H}^{\eta}_{\pi}(\mathcal{X})$ onto $\mathcal{H}$. In learning theory
$\rho(\mathcal{H})$ is known as the approximation error of the hypothesis
space $\mathcal{H}$ [see e.g. 41]. While this error may vanish for infinite-
dimensional spaces, when $\mathcal{H}$ is finite dimensional, controlling
$\rho(\mathcal{H})$ is crucial to achieving statistical consistency. This can
be accomplished by minimizing (11), which is equivalent to

$\min_{G\colon\mathcal{H}\to\mathcal{H}}\lVert P_{\mathcal{H}}(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}{-}G\rVert^{2}_{{\rm{HS}}\left(\mathcal{H},\mathcal{H}^{\eta}_{\pi}\right)}{=}\lVert\mathcal{Z}(\mathcal{Z}^{*}\mathcal{Z})^{\dagger}\mathcal{Z}^{*}(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}\mathcal{Z}-\mathcal{Z}\textsc{G}\rVert^{2}_{{\rm{HS}}\left(\mathbb{R}^{m},\mathcal{H}^{\eta}_{\pi}\right)}$
(13)

where $(\cdot)^{\dagger}$ is the Moore-Penrose’s pseudoinverse. Using the
covariance matrices

$\mathcal{Z}^{*}(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}\mathcal{Z}\,{=}\,\textsc{C}\,{=}\,\big{(}\mathbb{E}_{x\sim\pi}[z_{i}(x)z_{j}(x)]\big{)}_{i,j\in[m]},\quad\text{
and
}\quad\textsc{W}=\mathcal{Z}^{*}\mathcal{Z}\,{=}\,\big{(}\mathfrak{E}^{\eta}_{\pi}[z_{i},z_{j}]\big{)}_{i,j\in[m]},$
(14)

w.r.t. the invariant distribution and energy, respectively, gives the ridge
regularized (RR) solution
$\textsc{G}=(\textsc{W}+\gamma\textsc{I})^{-1}\textsc{C}$, $\gamma>0$. The
induced RR estimator of the resolvent,
$G_{\eta,\gamma}:\mathcal{H}\rightarrow\mathcal{H}$ is given, for every
$h_{u}\in\mathcal{H}$, by
$G_{\eta,\gamma}h_{u}:=\mathcal{Z}(\textsc{W}+\gamma\textsc{I})^{-1}\textsc{C}u=z(\cdot)^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}(\textsc{W}+\gamma\textsc{I})^{-1}\textsc{C}u$,
and it can be estimated given data from $\pi$ by replacing expectation and the
energy in (14) 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 eigenfunctions (corresponding to the slowest time scales) of the
infinitesimal generator from biased data. First, we address regressing the
generator on an a priori fixed hypothesis space $\mathcal{H}$. Then we
introduce our deep-learning method to either build a suitable space
$\mathcal{H}$, or even directly learn the eigenfunctions.

Unbiasing generator regression  Whenever $\pi$ is absolutely continuous w.r.t.
$\pi^{\prime}$, the regularized energy kernel (9) satisfies the simple
identity

$\mathfrak{E}^{\eta}_{\pi}[f,g]=\mathfrak{E}^{\eta}_{\pi^{\prime}}\big{[}f\sqrt{\tfrac{d\pi}{d\pi^{\prime}}},g\sqrt{\tfrac{d\pi}{d\pi^{\prime}}}\big{]},\quad
f,g\in H^{1,2}_{\pi}(\mathcal{X}),$ (15)

which, recalling the rightmost equation in (7), implies that when the bias $V$
and the diffusion coefficient $\beta$ are known, the energy kernel can be
empirically estimated through samples from $\pi^{\prime}$ via (10). Moreover,
when the potential $U$ is known too, we can use (9). Now, leveraging on (15)
we directly obtain that

$\mathcal{R}_{\partial}(G)\equiv\mathcal{R}_{\partial}(\textsc{G}){=}\mathfrak{E}^{\eta}_{x\sim\pi^{\prime}}\big{[}\lVert\chi_{\eta}(x)\\!-\\!\widehat{\textsc{G}}^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}z(x)\rVert_{2}\sqrt{\tfrac{d\pi}{d\pi^{\prime}}(x)}\big{]}\leq\kappa_{V}\mathfrak{E}^{\eta}_{x\sim\pi^{\prime}}\lVert\chi_{\eta}(x)\\!-\\!\widehat{\textsc{G}}^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}z(x)\rVert_{2}$
(16)

where $\kappa_{V}=\operatorname*{\text{\rm
ess\,sup}}_{x\sim\pi^{\prime}}\frac{d\pi}{d\pi^{\prime}}(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
$\chi_{\eta}(x)$ can be estimated, one can derive principled estimators of the
true generator $\mathcal{L}$’s dominant eigenpairs from the biased dynamics
generated by $\mathcal{L}^{\prime}$. This is established by the following
proposition, the proof of which is presented in Appendix B.

######  Theorem 1.

Let $\mathcal{D}_{n}=(x_{i}^{\prime})_{i\in[n]}$ be the biased dataset
generated from $\pi^{\prime}$. Let $w(x)=e^{\beta V(x)}$ and define the
empirical covariances w.r.t. the empirical distribution
$\widehat{\pi}^{\prime}\,{=}\,n^{{-}1}\sum_{i{\in}[n]}\delta_{x_{i}^{\prime}}$
by

$\widehat{\textsc{C}}=\big{(}\mathbb{E}_{x^{\prime}\sim\widehat{\pi}^{\prime}}[w(x^{\prime})z_{i}(x^{\prime})z_{j}(x^{\prime})]\big{)}_{i,j\in[m]}\quad\text{
and
}\quad\widehat{\textsc{W}}=\big{(}\mathfrak{E}^{\eta}_{\widehat{\pi}^{\prime}}[\sqrt{w}z_{i},\sqrt{w}z_{j}]\big{)}_{i,j\in[m]}.$
(17)

Compute the eigenpairs $(\nu_{i},v_{i})_{i\in[m]}$ of the RR estimator
$\widehat{\textsc{G}}_{\eta,\gamma}\,{=}\,(\widehat{\textsc{W}}\,{+}\,\eta\gamma\textsc{I})^{-1}\widehat{\textsc{C}}$,
and estimate the eigenpairs in (4) as
$(\widehat{\lambda}_{i},\widehat{f}_{i})\,{=}\,(\eta\,{-}\,1/\nu_{i},z(\cdot)^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}v_{i})$.
If the elements of $\mathcal{H}$ and their gradients are essentially bounded,
and $\displaystyle{\lim_{m\to\infty}}\rho(\mathcal{H})\,{=}\,0$, then for
every $\varepsilon>0$, there exist
$(m,n,\gamma)\,{\in}\,\mathbb{N}\,{\times}\,\mathbb{N}\,{\times}\,\mathbb{R}_{+}$,
such that, for every $i\in[m]$,
$\lvert\lambda_{i}\,{-}\,\widehat{\lambda}_{i}\rvert\leq\varepsilon$ and
$\sin_{L^{2}_{\pi}}(\sphericalangle(f_{i},\widehat{f}_{i}))\,{\leq}\,\varepsilon$,
with high probability.

Note that, due to the form of the estimator, the normalizing constant $\int
w(x)dx$ does not need be computed. Moreover, relying on the upper bound in
(16) we can alternatively compute $\widehat{\textsc{C}}$ and
$\widehat{\textsc{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 $\rho(\mathcal{H})$ in (12) is controlled. It is
therefore natural to minimize $\rho(\mathcal{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^{\theta}\,{=}\,(z^{\theta}_{i})_{i\in[m]}\colon\mathcal{X}\,{\to}\,\mathbb{R}^{m}$
be a neural network (NN) embedding parameterized by $\theta\,{\in}\,\Theta$
weights with continuously differentiable activation functions, and let
$\lambda^{\theta}_{i}$, $i\in[m]$, be real non-positive (trainable) weights.
We propose to optimize the NN to find the slowest time-scales
$\lambda^{\theta}_{i}$ that solve the eigenvalue equation
$\mathcal{L}z^{\theta}_{i}\,{=}\,\lambda^{\theta}_{i}z^{\theta}_{i}$,
$i\,{\in}\,[m]$. Letting
$\mathcal{Z}_{\theta}\colon\mathbb{R}^{m}\to\mathcal{H}^{\eta}_{\pi}(\mathcal{X})$
be the (parameterized) injection operator, given, for every
$u\in\mathbb{R}^{m}$ by
$\mathcal{Z}_{\theta}u\,{=}\,\sum_{i\in[m]}z^{\theta}_{i}u_{i}$, and denoting
$\Lambda^{\eta}_{\theta}\,{=}\,(\eta I\,{-}\,\operatorname{\text{\rm
diag}}(\lambda^{\theta}_{1},\ldots,\lambda^{\theta}_{m}))^{-1}$, the
eigenvalue equations for the resolvent then become $(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}\mathcal{Z}_{\theta}\,{=}\,\mathcal{Z}_{\theta}\Lambda^{\eta}_{\theta}$.
In other words, we aim to find the best rank-$m$ decomposition of resolvent
$(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}\,{\approx}\,\mathcal{Z}_{\theta}\Lambda^{\eta}_{\theta}\mathcal{Z}_{\theta}^{*}$.
Therefore, for some hyperparameter $\alpha\,{\geq}\,0$ we introduce the loss

$\mathcal{E}_{\alpha}(\theta):=\lVert(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}-\mathcal{Z}_{\theta}\Lambda^{\theta}_{\eta}\mathcal{Z}_{\theta}^{*}\rVert^{2}_{{\rm{HS}}\left(\mathcal{H}^{\eta}_{\pi}\right)}-\lVert(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}\rVert^{2}_{{\rm{HS}}\left(\mathcal{H}^{\eta}_{\pi}\right)}+\alpha\\!\sum_{i,j\in[m]}({\langle
z^{\theta}_{i},z^{\theta}_{j}\rangle}_{L^{2}_{\pi}}{-}\delta_{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 $(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}\in{\rm{HS}}\left(\mathcal{H}^{\eta}_{\pi}(\mathcal{X})\right)$
(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 ($\alpha=0$), but also to optimize the neural network to extract
eigenfunctions as features, then this last term ($\alpha>0$) encourages the
orthonormality of features in $L^{2}_{\pi}(\mathcal{X})$, an idea successfully
exploited in machine learning and computational chemistry [see e.g. 24, and
references therein].

Recalling (14) and denoting by $\textsc{C}_{\theta}$ and $\textsc{W}_{\theta}$
the covariance matrices associated to the parameterized features, after some
algebra, we obtain that

$\mathcal{E}_{\alpha}(\theta)=\operatorname{\text{\rm
tr}}\left[\textsc{W}_{\theta}\Lambda^{\eta}_{\theta}\textsc{W}_{\theta}\Lambda^{\eta}_{\theta}-2\textsc{C}_{\theta}\Lambda^{\eta}_{\theta}+\alpha(\textsc{C}_{\theta}-\textsc{I})^{2}\right].$
(18)

In turn, this can be estimated from biased data by two independent samples
$\widehat{\pi}^{\prime}_{1}$ and $\widehat{\pi}^{\prime}_{2}$ as

$\mathcal{E}_{\alpha}^{\widehat{\pi}^{\prime}_{1},\widehat{\pi}^{\prime}_{2}}(\theta)=\operatorname{\text{\rm
tr}}\left[\widehat{\textsc{W}}_{\theta}^{1}\Lambda^{\eta}_{\theta}\widehat{\textsc{W}}^{2}_{\theta}\Lambda^{\eta}_{\theta}{-}\widehat{w}^{1}\widehat{\textsc{C}}_{\theta}^{2}\Lambda^{\eta}_{\theta}-\widehat{w}_{2}\widehat{\textsc{C}}_{\theta}^{1}\Lambda^{\eta}_{\theta}+\alpha(\widehat{\textsc{C}}_{\theta}^{1}-\widehat{w}_{1}\textsc{I})(\widehat{\textsc{C}}_{\theta}^{2}-\widehat{w}_{2}\textsc{I})\right],$
(19)

where $\widehat{\textsc{C}}_{\theta}^{k}$ and
$\widehat{\textsc{W}}_{\theta}^{k}$ are the empirical covariances given by
(17) for distribution $\widehat{\pi}^{\prime}_{k}$, while
$\widehat{w}^{k}=\mathbb{E}_{x^{\prime}\sim\widehat{\pi}^{\prime}_{k}x}w(x^{\prime})$,
$k{\in}[2]$. Importantly, the computational complexity of the loss (19) is of
the order $\mathcal{O}(n\hskip 0.85355ptm^{2}d)$, where $d$ is the state
dimension and $n$ the sample size, however it can be reduced to
$\mathcal{O}(n\hskip 0.85355pt{m}\hskip 0.85355ptd)$ (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 1, provides theoretical guarantees for our approach. The
proof and discussion are provided in Appendix B.

######  Theorem 2.

Given a compact operator $(\eta I\\!-\\!\mathcal{L})^{\\!-\\!1}$, $\eta>0$, if
$(z^{\theta})_{i\in[m]}\subseteq\mathcal{H}^{\eta}_{\pi}(\mathcal{X})$ for all
$\theta\\!\in\\!\Theta$, then

$\mathbb{E}\big{[}\mathcal{E}_{\alpha}^{\widehat{\pi}^{\prime}_{1},\widehat{\pi}^{\prime}_{2}}(\theta)\big{]}=\overline{w}^{2}\,\mathcal{E}_{\alpha}(\theta)\geq-{\textstyle{\sum_{i\in[m]}}\tfrac{\overline{w}^{2}}{(\eta-\lambda_{i})^{2}}},\quad\text{for
all }\theta\in\Theta,$ (20)

where $\overline{w}=\mathbb{E}_{x\sim\pi^{\prime}}[w(x)]$. Moreover, if
$\alpha\,{>}\,0$ and $\lambda_{m+2}\,{<}\,\lambda_{m+1}$, then the equality
holds if and only if
$(\lambda^{\theta}_{i},z^{\theta}_{i})\,{=}\,(\lambda_{i},f_{i})$ $\pi$-a.e.,
up to the ordering of indices and choice of eigenfunction signs for
$i\,{\in}\,[m]$.

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 $(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}$ (without orthonormality loss, i.e.
$\alpha=0$), on which the estimator in Theorem 1 can be computed, or even the
individual eigenpairs ($\alpha>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
$\beta$ 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 $\mathcal{L}$.

Algorithm 1 From biased to unbiased dynamics via infinitesimal generator

1:Parameters $\eta>0$ shift of the generator, $m$ number of wanted
eigenfunctions, $K$ number of optimization steps, $\gamma>0$ and $\alpha>0$
regression and NN hyperparameters

2:Inputs Dataset $\mathcal{D}_{n}=(x_{\ell})_{\ell\in[n]}$ gathered from a
simulation with bias potential $V$

3:Compute weights $w(x_{\ell})=\exp(\beta V(x_{\ell}))$, $\ell\in[n]$, to be
used in line 7

4:if the dictionary of function $z$ does not already exist then

5: Initialization: randomly initialize neural networks weights of
$\Lambda^{\theta}$ and $(z^{\theta}_{i})_{i\in[m]}$, set $k=0$

6: while $k<K$ do

7: Compute $\widehat{\textsc{C}}_{\theta}^{j}$ and
$\widehat{\textsc{W}}_{\theta}^{j}$, $j=1,2$, using (17) for two independent
batches $\widehat{\pi}^{\prime}_{1}$ and $\widehat{\pi}^{\prime}_{2}$

8: Compute loss
$\widehat{\mathcal{E}}^{\widehat{\pi}^{\prime}_{1},\widehat{\pi}^{\prime}_{2}}(\theta)$
using (19) and backpropagate

9: end while

10:end if

11:Compute $\widehat{\textsc{C}}$ and $\widehat{\textsc{W}}$ using (17) the
datatset $\mathcal{D}_{n}$

12:Compute the eigenpairs $(\nu_{i},v_{i})_{i\in[m]}$ of
$\widehat{\textsc{G}}_{\eta,\gamma}{=}(\widehat{\textsc{W}}{+}\eta\gamma\textsc{I})^{-1}\widehat{\textsc{C}}$

13:Output Estimated eigenpairs of $\mathcal{L}$ are
$(\widehat{\lambda}_{i},\widehat{f}_{i}){=}(\eta{-}1/\nu_{i},z^{\theta}(\cdot)^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}v_{i})$,
$i\in[m]$

##  6 Experiments

In this section, we test the method described above on well-established [14,
9, 31, 35] 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
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

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 5 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, 36, 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 1. We compare
the results with deepTICA and a state of the art generator learning approach
in [47]. 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 [47] 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.

Figure 1: 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. [47] (fourth column). x and y 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 $\lambda$ are also reported.

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 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, 47, 40]: it
presents a conformational change which is a rare event described by the angles
$\phi$ and $\psi$. In the studies made on this system, the angle $\phi$ 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 $\psi$ 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 $\psi$ 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 $\phi$ 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 [47, 4, 10] to transform
the positions of the atoms. The results are presented in Figure 2. 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 $\phi$ and
$\theta$. 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].

Figure 2: 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 $\phi$ and $\psi$ dihedral angles. c) first eigenfunction represented on
dataset 2, in the plane of the $\phi$ and $\theta$ 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]

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 2, when
comparing to the committor model obtained in [19] indicating the good
performance of our method.

##  7 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 sub-optimal 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.

## References

  * A Cabannes and Bach, [2024] A Cabannes, V. and Bach, F. (2024).  The Galerkin method beats graph-based approaches for spectral algorithms.  In Proceedings of The 27th International Conference on Artificial Intelligence and Statistics, volume 238, pages 451–459. PMLR.
  * Abraham et al., [2015] 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.
  * Bakry et al., [2014] Bakry, D., Gentil, I., and Ledoux, M. (2014).  Analysis and Geometry of Markov Diffusion Operators.  Springer.
  * Belkacemi et al., [2022] 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–78.
  * Berezhkovskii and Szabo, [2004] 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):9186–9187.
  * Bolhuis et al., [2000] 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.
  * Bonati et al., [2021] 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):e2113533118.
  * Bonati et al., [2023] Bonati, L., Trizio, E., Rizzi, A., and Parrinello, M. (2023).  A unified framework for machine learning collective variables for enhanced sampling simulations: mlcolvar.  The Journal of Chemical Physics, 159(1):014801.
  * Chen et al., [2023] 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.
  * Chen et al., [2018] 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.
  * Comer et al., [2015] 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 to ask.  The Journal of Physical Chemistry B, 119(3):1129–1151.
  * Dama et al., [2014] Dama, J. F., Parrinello, M., and Voth, G. A. (2014).  Well-tempered metadynamics converges asymptotically.  Phys. Rev. Lett., 112:240602.
  * Davis and Kahan, [1970] Davis, C. and Kahan, W. M. (1970).  The rotation of eigenvectors by a perturbation. iii.  SIAM Journal on Numerical Analysis, 7(1):1–46.
  * France-Lanord et al., [2024] France-Lanord, A., Vroylandt, H., Salanne, M., Rotenberg, B., Saitta, A. M., and Pietrucci, F. (2024).  Data-driven path collective variables.  Journal of Chemical Theory and Computation, 20(8):3069–3084.
  * Hou et al., [2023] Hou, B., Sanjari, S., Dahlin, N., Bose, S., and Vaidya, U. (2023).  Sparse learning of dynamical systems in RKHS: An operator-theoretic approach.  In Proceedings of the 40th International Conference on Machine Learning, volume 202, pages 13325–13352. PMLR.
  * Invernizzi and Parrinello, [2020] Invernizzi, M. and Parrinello, M. (2020).  Rethinking metadynamics: From bias potentials to probability distributions.  The Journal of Physical Chemistry Letters, 11(7):2731–2736.
  * Invernizzi and Parrinello, [2022] Invernizzi, M. and Parrinello, M. (2022).  Exploration vs convergence speed in adaptive-bias enhanced sampling.  Journal of Chemical Theory and Computation, 18(6):3988–3996.
  * Kabsch, [1976] Kabsch, W. (1976).  A solution for the best rotation to relate two sets of vectors.  Acta Crystallographica Section A, 32(5):922–923.
  * Kang et al., [2024] Kang, P., Trizio, E., and Parrinello, M. (2024).  Computing the committor with the committor: an anatomy of the transition state ensemble.  Preprint, arXiv:2401.05279.
  * Klus et al., [2020] 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, 406:132416.
  * [21] Kostic, V., Lounici, K., Halconruy, H., Devergne, T., and Pontil, M. (2024a).  Learning the infinitesimal generator of stochastic diffusion processes.  Preprint, arXiv:2405.12940.
  * Kostic et al., [2022] 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.
  * Kostic et al., [2023] Kostic, V. R., Lounici, K., Novelli, P., and Pontil, M. (2023).  Sharp spectral rates for koopman operator learning.  In 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-homogeneous dynamical systems.  In The Twelfth International Conference on Learning Representations.
  * Laio and Parrinello, [2002] Laio, A. and Parrinello, M. (2002).  Escaping free-energy minima.  Proceedings of the National Academy of Sciences, 99(20):12562–12566.
  * Leitold et al., [2020] 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).
  * Lelièvre et al., [2024] Lelièvre, T., Pigeon, T., Stoltz, G., and Zhang, W. (2024).  Analyzing multimodal probability measures with autoencoders.  The Journal of Physical Chemistry B, 128(11):2607–2631.
  * Magrino et al., [2022] 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):8887–8900.
  * Mardt et al., [2018] Mardt, A., Pasquali, L., Wu, H., and Noé, F. (2018).  VAMPnets for deep learning of molecular kinetics.  Nature Communications, 9(1).
  * McCarty and Parrinello, [2017] 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.
  * Mori et al., [2020] 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, 153(5):054115.
  * Oksendal, [2013] Oksendal, B. (2013).  Stochastic Differential Equations: an Introduction with Applications.  Springer Science & Business Media.
  * Pietrucci, [2017] Pietrucci, F. (2017).  Strategies for the exploration of free energy landscapes: Unity in diversity and challenges ahead.  Reviews in Physics, 2:32–45.
  * Pillaud-Vivien and Bach, [2023] Pillaud-Vivien, L. and Bach, F. (2023).  Kernelized diffusion maps.  In The Thirty Sixth Annual Conference on Learning Theory, pages 5236–5259. PMLR.
  * Plainer et al., [2023] Plainer, M., Stark, H., Bunne, C., and Günnemann, S. (2023).  Transition path sampling with Boltzmann generator-based MCMC moves.  In NeurIPS 2023 Generative AI and Biology (GenBio) Workshop.
  * Ray et al., [2023] Ray, D., Trizio, E., and Parrinello, M. (2023).  Deep learning collective variables from transition path ensemble.  The Journal of Chemical Physics, 158(20):204102.
  * Reed and Simon, [1980] Reed, M. and Simon, B. (1980).  Functional Analysis (Modern Mathematical Physics, Volume 1).  Academic Press.
  * Salomon-Ferrer et al., [2013] Salomon-Ferrer, R., Case, D. A., and Walker, R. C. (2013).  An overview of the amber biomolecular simulation package.  WIREs Computational Molecular Science, 3(2):198–210.
  * Schwantes and Pande, [2015] Schwantes, C. R. and Pande, V. S. (2015).  Modeling molecular kinetics with tica and the kernel trick.  Journal of Chemical Theory and Computation, 11(2):600–608.
  * Shmilovich and Ferguson, [2023] 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, 127(15):3497–3517.
  * Steinwart and Christmann, [2008] Steinwart, I. and Christmann, A. (2008).  Support Vector Machines.  Springer New York.
  * Torrie and Valleau, [1977] 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.
  * Tribello et al., [2014] Tribello, G. A., Bonomi, M., Branduardi, D., Camilloni, C., and Bussi, G. (2014).  Plumed 2: New feathers for an old bird.  Computer Physics Communications, 185(2):604–613.
  * Tuckerman, [2023] Tuckerman, M. E. (2023).  Statistical Mechanics: Theory and Molecular Simulation.  Oxford university press.
  * Vanden-Eijnden, [2006] Vanden-Eijnden, E. (2006).  Transition Path Theory, pages 453–493.  Springer Berlin Heidelberg, Berlin, Heidelberg.
  * Yang and Parrinello, [2018] Yang, Y. I. and Parrinello, M. (2018).  Refining collective coordinates and improving free energy representation in variational enhanced sampling.  Journal of Chemical Theory and Computation, 14(6):2889–2894.
  * Zhang et al., [2022] Zhang, W., Li, T., and Schütte, C. (2022).  Solving eigenvalue pdes of metastable diffusion processes using artificial neural networks.  Journal of Computational Physics, 465:111377.

## 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.

Table 1: Summary of used notations.

notation meaning notation meaning $[\,\cdot\,]$ set $\\{1,2\ldots,\cdot\\}$
$[\,\cdot\,]_{+}$ nonnegative part of a number $\mathcal{X}$ state space of
the Markov process $(X_{t})_{t\geq 0}$ time-homogeneous Markov process
$dW_{t}$ Brownian motion $\beta$ inverse temperature of the system $U$
potential energy of the system $V$ bias potential $\pi$ Boltzmann distribution
of potential $U$ $\pi^{\prime}$ Boltzmann distribution of potential $U{+}V$
$\widehat{\pi}$ empirical version of $\pi$ $\widehat{\pi}^{\prime}$ empirical
version of $\pi^{\prime}$ $d\pi/d\pi^{\prime}$ density of $\pi$ w.r.t.
$\pi^{\prime}$ $w$ exponential weights of $d\pi/d\pi^{\prime}$
$L^{2}_{\pi}(\mathcal{X})$ $L^{2}$ space on $\mathcal{X}$ w.r.t. the measure
$\pi$ $H^{1,2}_{\pi}(\mathcal{X})$ Sobolev space w.r.t. $\pi$ on $\mathcal{X}$
$\mathcal{T}_{t}$ transfer operator with lag time $t$ $\hat{T}_{t}$ empirical
estimator of $\mathcal{T}_{t}$ $\mathcal{L}$ generator of the true process
$\mathcal{L}^{\prime}$ generator of the biased process $\eta$ shift parameter
$(\eta I\\!-\\!\mathcal{L})^{\\!-\\!1}$ resolvent at $\eta$ $\mathcal{H}$
hypothetical Hilbert space $z$ basis functions
$\mathfrak{E}^{\eta}[\cdot,\cdot]$ regularized energy kernel
$\mathfrak{E}^{\eta}[\cdot]$ regularized energy norm
$\mathcal{H}^{\eta}_{\pi}(\mathcal{X})$ space associated to the energy norm
$\mathcal{Z}$ injection operator $\mathcal{R}_{\partial}$ risk functional
$\widehat{\mathcal{R}}_{\partial}$ empirical risk functional C covariance
$\widehat{\textsc{C}}$ empirical covariance W covariance in energy space
$\widehat{\textsc{W}}$ empirical covariance in energy space $\lambda_{i}$
generator eigenvalue $f_{i}$ generator eigenfunction $\widehat{\lambda}_{i}$
empirical estimator of $\lambda_{i}$ $\widehat{f}_{i}$ empirical estimator of
$f_{i}$ $z^{\theta}$ neural network embedding $\theta$ neural network weights
$\mathcal{Z}_{\theta}$ neural network injection operator
$\Lambda_{\theta}^{\eta}$ neural network weights $\mathcal{E}_{\alpha}$
regularized loss function $\widehat{\mathcal{E}}_{\alpha}$ regularized
empirical loss function

##  Appendix A Background

In this work we consider stochastic differential equations (SDE) of the form

$dX_{t}=a(X_{t})dt+b(X_{t})dW_{t}\quad\text{and}\quad X_{0}=x.$ (21)

The special case of Langevin equation considered in the main body of the paper
corresponds to $a(x)=-\nabla U(x)$ and $b(x)=\sqrt{\frac{2}{\beta}}I$.
Equation (21) describes the dynamics of the random vector $X_{t}$ in the state
space $\mathcal{X}\subseteq\mathbb{R}^{d}$, governed by the drift
$a:\mathbb{R}^{d}\to\mathbb{R}^{d}$ and the diffusion
$b:\mathbb{R}^{d}\to\mathbb{R}^{d\times p}$ coefficients, where $W_{t}$ is a
$\mathbb{R}^{p}$-dimensional standard Brownian motion. Under the usual
conditions [see e.g. 32] that $a$ and $b$ are globally Lispchitz and sub-
linear, the SDE (21) admits an unique strong solution $X=(X_{t})_{\geqslant
0}$ that is a Markov process to which we can associate the semigroup of Markov
transfer operators $(\mathcal{T}_{t})_{t\geq 0}$ defined, for every $t\geq 0$,
as

$[\mathcal{T}_{t}f](x):=\mathbb{E}[f(X_{t})|X_{0}=x],\;\;x\in\mathcal{X},\,f\colon\mathcal{X}\to\mathbb{R}.$
(22)

For stable processes, the distribution of $X_{t}$ converges to an invariant
measure $\pi$ on $\mathcal{X}$, such that $X_{0}\sim\pi$ implies that
$X_{t}\sim\pi$ for all $t\geq 0$. In such cases, one can define the semigroup
on $L^{2}_{\pi}(\mathcal{X})$, and characterize the process by the
infinitesimal generator of the semi-group $(\mathcal{T}_{t})_{t\geqslant 0}$,

$\mathcal{L}:=\lim_{t\to 0^{+}}\frac{\mathcal{T}_{t}-I}{t}$ (23)

defined on the Sobolev space $H^{1,2}_{\pi}(\mathcal{X})$ of functions in
$L^{2}_{\pi}(\mathcal{X})$ whose gradient are in $L^{2}_{\pi}(\mathcal{X})$,
too, i.e. $\mathcal{L}\colon L^{2}_{\pi}(\mathcal{X})\to
L^{2}_{\pi}(\mathcal{X})$ and $\operatorname*{\text{\rm
dom}}(\mathcal{L})=H^{1,2}_{\pi}(\mathcal{X})$. The transfer operator and the
generator are linked to each other by the formula
$\mathcal{T}_{t}=\exp(t\mathcal{L})$.

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\in\mathcal{C}^{2}(\mathcal{X},\mathbb{R})$, Itô’s formula
[see for instance 3, p. 495] provides for $t\in\mathbb{R}_{+}$,

$\displaystyle f(X_{t})-f(X_{0})$
$\displaystyle=\int_{0}^{t}\sum_{i=1}^{d}\partial_{i}f(X_{s})dX_{s}^{i}+\tfrac{1}{2}\int_{0}^{t}\sum_{i,j=1}^{d}\partial_{ij}^{2}f(X_{s})d\langle
X^{i},X^{j}\rangle_{s}$ $\displaystyle=\int_{0}^{t}\nabla
f(X_{s})^{{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}}dX_{s}+\tfrac{1}{2}\int_{0}^{t}\mathrm{Tr}\big{[}X_{s}^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}(\nabla^{2}f)(X_{s})X_{s}\big{]}ds.$

Recalling (21), we get

$\displaystyle f(X_{t})$ $\displaystyle=f(X_{0})+\int_{0}^{t}\bigg{[}\nabla
f(X_{s})^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}a(X_{s})+\tfrac{1}{2}\mathrm{Tr}\big{[}b(X_{s})^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}(\nabla^{2}f(X_{s}))b(X_{s})\big{]}\bigg{]}ds$
$\displaystyle\quad+\int_{0}^{t}\nabla
f(X_{s})^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}b(X_{s})dW_{s}.$
(24)

Provided $f$ and $b$ are smooth enough, the expectation of the last stochastic
integral vanishes so that we get

$\mathbb{E}[f(X_{t})|X_{0}=x]=f(x)+\int_{0}^{t}\mathbb{E}\Big{[}\nabla
f(X_{s})^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}a(X_{s})+\tfrac{1}{2}\mathrm{Tr}\big{[}b(X_{s})^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}(\nabla^{2}f(X_{s}))b(X_{s})\big{]}\Big{|}X_{0}=x\Big{]}ds$

Recalling that $\mathcal{L}=\underset{t\rightarrow
0^{+}}{\lim}(\mathcal{T}_{t}f-f)/t$, we get for every $x\in\mathcal{X}$,

$\displaystyle\mathcal{L}f(x)$ $\displaystyle=\lim_{t\rightarrow
0}\frac{\mathbb{E}[f(X_{t})|X_{0}=x]-f(x)}{t}$
$\displaystyle=\lim_{t\rightarrow
0}\frac{1}{t}\bigg{[}\int_{0}^{t}\mathbb{E}\Big{[}\nabla
f(X_{s})^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}a(X_{s})+\tfrac{1}{2}\mathrm{Tr}\big{[}(X_{s})^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}(\nabla^{2}f(X_{s}))b(X_{s})\big{]}\Big{]}ds\Big{|}X_{0}=x\bigg{]}$
$\displaystyle=\nabla
f(x)^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}a(x)+\tfrac{1}{2}\mathrm{Tr}\big{[}b(x)^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}(\nabla^{2}f(x))b(x)\big{]},$
(25)

which 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
$\mathcal{H}$ the resolvent set of the operator $A$ is defined as
$\rho(A)=\left\\{\lambda\in\mathbb{C}\,|\,A-\lambda I\text{ is
bijective}\right\\}$, and its spectrum $\operatorname*{\text{\rm
Sp}}(A)=\mathbb{C}\setminus\\{\rho(A)\\}$, let
$\lambda\subseteq\operatorname*{\text{\rm Sp}}(A)$ be the isolated part of the
spectra, i.e. both $\lambda$ and $\mu=\operatorname*{\text{\rm
Sp}}(A)\setminus\lambda$ are closed in $\operatorname*{\text{\rm Sp}}(A)$.
Then, the Riesz spectral projector
$P_{\lambda}\colon\mathcal{H}\to\mathcal{H}$ is defined by

$P_{\lambda}=\frac{1}{2\pi}\int_{\Gamma}(zI-A)^{-1}dz,$ (26)

where $\Gamma$ is any contour in the resolvent set $\operatorname*{\text{\rm
Res}}(A)$ with $\lambda$ in its interior and separating $\lambda$ from $\mu$.
Indeed, we have that $P_{\lambda}^{2}=P_{\lambda}$ and
$\mathcal{H}=\operatorname*{\text{\rm
Im}}(P_{\lambda})\oplus\operatorname{\text{\rm Ker}}(P_{\lambda})$ where
$\operatorname*{\text{\rm Im}}(P_{\lambda})$ and $\operatorname{\text{\rm
Ker}}(P_{\lambda})$ are both invariant under $A$, and we have
$\operatorname*{\text{\rm Sp}}(A_{|_{\operatorname*{\text{\rm
Im}}(P_{\lambda})}})=\lambda$ and $\operatorname*{\text{\rm
Sp}}(A_{|_{\operatorname{\text{\rm Ker}}(P_{\lambda})}})=\mu$. Moreover,
$P_{\lambda}+P_{\mu}=I$ and $P_{\lambda}P_{\mu}=P_{\mu}P_{\lambda}=0$.

Finally if $A$ is a compact operator, then the Riesz-Schauder theorem [see
e.g. 37] assures that $\operatorname*{\text{\rm Sp}}(T)$ is a discrete set
having no limit points except possibly $\lambda=0$. Moreover, for any nonzero
$\lambda\in\operatorname*{\text{\rm Sp}}(T)$, then $\lambda$ is an eigenvalue
(i.e. it belongs to the point spectrum) of finite multiplicity, and, hence, we
can deduce the spectral decomposition in the form

$A=\sum_{\lambda\in\operatorname*{\text{\rm Sp}}(A)}\lambda\,P_{\lambda},$
(27)

where the geometric multiplicity of $\lambda$,
$r_{\lambda}=\operatorname*{\text{\rm rank}}(P_{\lambda})$, is bounded by the
algebraic multiplicity of $\lambda$. If additionally $A$ is a normal operator,
i.e. $AA^{*}=A^{*}A$, then $P_{\lambda}=P_{\lambda}^{*}$ is an orthogonal
projector for each $\lambda\in\operatorname*{\text{\rm Sp}}(A)$ and
$P_{\lambda}=\sum_{i=1}^{r_{\lambda}}\psi_{i}\otimes\psi_{i}$, where
$\psi_{i}$ are normalized eigenfunctions of $A$ corresponding to $\lambda$ and
$r_{\lambda}$ is both algebraic and geometric multiplicity of $\lambda$.

We conclude this section by stating the well-known Davis-Kahan perturbation
bound for eigenfunctions of self-adjoint compact operators.

######  Proposition 1 ([13]).

Let $A$ be compact self-adjoint operator on a separable Hilbert space
$\mathcal{H}$. Given a pair
$(\widehat{\mu},\widehat{f})\in\mathbb{C}\times\mathcal{H}$ such that
$\lVert\widehat{f}\rVert=1$, let $\lambda$ be the eigenvalue of $A$ that is
closest to $\widehat{\mu}$ and let $f$ be its normalized eigenfunction. If
$\widehat{g}=\min\\{\lvert\widehat{\mu}-\lambda\rvert\,|\,\lambda\in\operatorname*{\text{\rm
Sp}}(A)\setminus\\{\lambda\\}\\}>0$, then
$\sin(\sphericalangle(\widehat{f},f))\leq\lVert
A\widehat{f}-\widehat{\mu}\widehat{f}\rVert/\widehat{g}$.

##  Appendix B Unbiased generator regression

In this section, we prove Theorem 1 relying on recently developed statistical
theory of generator learning [21]. To that end, let
$\phi(x):=z(\cdot)^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}z(x)\in\mathcal{H}$
be a feature map of the RKHS space $\mathcal{H}$ of dimension
$\dim(\mathcal{H})=m$. Let $\textsc{W}u_{j}=\sigma_{j}^{2}u_{j}$ be the
eigenvalue decomposition of W and let $v_{j}:=u_{j}/\sigma_{j}$. This induces
the SVD of the injection operator, $\mathcal{Z}u_{j}=\sigma_{j}\tilde{z}_{j}$
for
$\tilde{z}_{j}:=z(\cdot)^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}v_{j}$.

Since $\mathcal{H}\subseteq H^{1,\infty}_{\pi}(\mathcal{X})$, we have that

$c_{\tau}=\operatorname*{\text{\rm
ess\,sup}}_{x\sim\pi}\textstyle{\sum_{j\in\mathbb{N}}}[\eta|\tilde{z}_{j}(x)|^{2}-z_{j}(x)[\mathcal{L}\tilde{z}_{j}](x)]<+\infty,$

and, denoting $\textsc{W}_{\gamma}:=\textsc{W}+\eta\gamma\textsc{I}$

$\operatorname{\text{\rm tr}}[\textsc{W}_{\gamma}^{-1}\textsc{W}_{\gamma}]\leq
m.$

In addition denote the empirical version of $\textsc{W}_{\gamma}$ as
$\widehat{\textsc{W}}_{\gamma}:=\widehat{\textsc{W}}+\eta\gamma\textsc{I}$.

Now, we can apply the following propositions from [21] to our setting,
recalling the notation for normalizing constant
$\overline{w}:=\mathbb{E}_{x\sim\pi^{\prime}}[w(x)]$ for which
$w(\cdot)/\overline{w}=d\pi/d\pi^{\prime}$.

######  Proposition 2.

Given $\delta>0$, with probability in the i.i.d. draw of $(x_{i})_{i=1}^{n}$
from $\pi$, it holds that

$\mathbb{P}\\{\lVert\widehat{\textsc{W}}-\textsc{W}\rVert\leq\varepsilon_{n}(\delta)\\}\geq
1-\delta,$

where

$\varepsilon_{n}(\delta)=\frac{2\lVert\textsc{W}\rVert}{3n}\mathcal{L}(\delta)+\sqrt{\frac{2\lVert\textsc{W}\rVert}{n}\mathcal{L}(\delta)}\quad\text{
and }\quad\mathcal{L}(\delta)=\log\frac{4\operatorname{\text{\rm
tr}}(\textsc{W})}{\delta\,\lVert\textsc{W}\rVert}.$ (28)

######  Proposition 3.

Given $\delta>0$, with probability in the i.i.d. draw of $(x_{i})_{i=1}^{n}$
from $\pi$, it holds that

$\mathbb{P}\left\\{\lVert\textsc{W}_{\gamma}^{-1/2}(\widehat{\textsc{W}}-\textsc{W})\textsc{W}_{\gamma}^{-1/2}\rVert\leq\varepsilon^{1}_{n}(\gamma,\delta)\right\\}\geq
1-\delta,$ (29)

where

$\varepsilon_{n}^{1}(\gamma,\delta)=\frac{2c_{\tau}}{3n}\mathcal{L}^{1}(\gamma,\delta)+\sqrt{\frac{2\,c_{\tau}}{n}\mathcal{L}^{1}(\gamma,\delta)},$
(30)

and

$\mathcal{L}^{1}(\gamma,\delta)=\ln\frac{4}{\delta}+\ln\frac{\operatorname{\text{\rm
tr}}(\textsc{W}_{\gamma}^{-1}\textsc{W})}{\lVert\textsc{W}_{\gamma}^{-1}\textsc{W}\rVert}.$

Moreover,

$\mathbb{P}\left\\{\lVert\textsc{W}_{\gamma}^{1/2}\widehat{\textsc{W}}_{\gamma}^{-1}\textsc{W}_{\gamma}^{1/2}\rVert\leq\frac{1}{1-\varepsilon_{n}^{1}(\gamma,\delta)}\right\\}\geq
1-\delta.$ (31)

######  Proposition 4.

With probability in the i.i.d. draw of $(x_{i})_{i=1}^{n}$ from $\pi$, it
holds

$\mathbb{P}\left\\{\lVert\textsc{W}_{\gamma}^{-1/2}(\widehat{\textsc{C}}-\textsc{C})\rVert_{F}\leq\varepsilon^{2}_{n}(\gamma,\delta)\right\\}\geq
1-\delta,$

where

$\varepsilon_{n}^{2}(\gamma,\delta)=\frac{4\,\sqrt{2\,m\lVert\textsc{W}\rVert}}{\eta}\,\ln\frac{2}{\delta}\,\sqrt{\frac{c_{\beta}}{n}+\frac{c_{\tau}}{n^{2}}}.$
(32)

We are now ready to prove Theorem 1, which we restate here for convenience.

See 1

###### Proof.

We first show that
$\mathcal{R}_{\partial}(\widehat{\textsc{G}}_{\eta,\gamma})<\varepsilon$ for
big enough $m,n\in\mathbb{N}$ and small enough $\gamma>0$.

Observe that

$\textsc{W}=\mathbb{E}[\widehat{\textsc{W}}]/\overline{w}\quad\text{ and
}\quad\textsc{C}=\mathbb{E}[\widehat{\textsc{C}}]/\overline{w}$ (33)

and, hence

$\textsc{G}_{\eta,\gamma}:=\textsc{W}_{\gamma}^{-1}\textsc{C}=(\mathbb{E}[\widehat{\textsc{W}}_{\gamma}])^{-1}(\mathbb{E}[\widehat{\textsc{C}}]),$

due to cancellation of $\overline{w}$.

Given $\varepsilon>0$, let $m\in\mathbb{N}$ be such that
$\rho(\mathcal{H})=\lVert(I-P_{\mathcal{H}})(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}\rVert^{2}_{{\rm{HS}}\left(\mathcal{H},\mathcal{H}^{\eta}_{\pi}\right)}<\varepsilon/3$.
Next, since

$P_{\mathcal{H}}(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}\mathcal{Z}{-}\mathcal{Z}\textsc{G}_{\eta,\gamma}{=}(I{-}\mathcal{Z}\textsc{W}_{\gamma}^{-1}\mathcal{Z}^{*})(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}\mathcal{Z}=\mathcal{Z}(\textsc{W}^{\dagger}\textsc{C}-\textsc{W}_{\gamma}^{-1}\textsc{C}),$

we have that

$\lVert P_{\mathcal{H}}(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}\mathcal{Z}{-}\mathcal{Z}\textsc{G}_{\eta,\gamma}\rVert_{{\rm{HS}}\left(\mathcal{H},\mathcal{H}^{\eta}_{\pi}\right)}=\lVert\textsc{W}^{1/2}(\textsc{W}^{\dagger}\textsc{C}-\textsc{W}_{\gamma}^{-1}\textsc{C})\rVert_{F}=\lVert\textsc{W}^{1/2}(Wx^{\dagger}-\textsc{W}_{\gamma}^{-1})\textsc{C}\rVert_{F}\to
0,$

as $\gamma\to 0$. Hence, let $\gamma>0$ be such that $\lVert
P_{\mathcal{H}}(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}\mathcal{Z}{-}\mathcal{Z}\textsc{G}_{\eta,\gamma}\rVert_{\mathcal{H}\to\mathcal{H}^{\eta}_{\pi}}<\varepsilon/3$.

Finally, using the decomposition of the risk

$\mathcal{R}_{\partial}(\widehat{\textsc{G}}_{\eta,\gamma})\leq\rho(\mathcal{H}){+}\lVert
P_{\mathcal{H}}(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}\mathcal{Z}{-}\mathcal{Z}\textsc{G}_{\eta,\gamma}\rVert_{{\rm{HS}}\left(\mathbb{R}^{m},\mathcal{H}^{\eta}_{\pi}\right)}{+}\lVert\mathcal{Z}(\widehat{\textsc{G}}_{\eta,\gamma}-\widehat{\textsc{G}}_{\eta,\gamma})\rVert_{{\rm{HS}}\left(\mathbb{R}^{m},\mathcal{H}^{\eta}_{\pi}\right)}$

it remains to show that for large enough $n$ we have
$\lVert\mathcal{Z}(\textsc{G}_{\eta,\gamma}-\widehat{\textsc{G}}_{\eta,\gamma})\rVert_{{\rm{HS}}\left(\mathbb{R}^{m},\mathcal{H}^{\eta}_{\pi}\right)}\leq\varepsilon/3$.

To that end observe that

$\displaystyle\textsc{W}_{\gamma}^{1/2}(\widehat{\textsc{G}}_{\eta,\gamma}-\textsc{G}_{\eta,\gamma})$
$\displaystyle\,{=}\,\textsc{W}_{\gamma}^{1/2}\widehat{\textsc{W}}_{\gamma}^{-1}(\widehat{\textsc{C}}-\widehat{\textsc{W}}_{\gamma}\textsc{W}_{\gamma}^{-1}\textsc{C}\pm\textsc{C})$
$\displaystyle\,{=}\,\textsc{W}_{\gamma}^{1/2}\widehat{\textsc{W}}_{\gamma}^{-1}\textsc{W}_{\gamma}^{1/2}\left(\textsc{W}_{\gamma}^{-1/2}(\widehat{\textsc{C}}-\textsc{C})-\textsc{W}_{\gamma}^{-1/2}(\widehat{\textsc{W}}-\textsc{W})\textsc{W}_{\gamma}^{-1/2}(\textsc{W}_{\gamma}^{-1/2}\textsc{C})\right).$

Thus, by multiplying the above expression by $\overline{w}$ and applying
Propositions 3 and 4, we obtain that there exists $n\in\mathbb{N}$ such that
$\lVert\mathcal{Z}(\textsc{G}_{\eta,\gamma}-\widehat{\textsc{G}}_{\eta,\gamma})\rVert_{{\rm{HS}}\left(\mathbb{R}^{m},\mathcal{H}^{\eta}_{\pi}\right)}\leq\varepsilon/3$.

Next, assuming that $\lVert\widehat{\textsc{W}}-\textsc{W}\rVert$ is small,
for the normalization of the estimated eigenfunctions we have that

$\frac{\lVert
v_{j}\rVert^{2}_{2}}{\lVert\widehat{f}_{j}\rVert^{2}_{\mathcal{H}^{\eta}_{\pi}}}=\frac{v_{j}^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}v_{j}}{v_{j}^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}\textsc{W}v_{j}}\leq\frac{v_{j}^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}v_{j}}{v_{j}^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}\widehat{\textsc{W}}v_{j}-v_{j}^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}(\textsc{W}-\widehat{\textsc{W}})v_{j}}\leq\frac{1}{\lambda_{\min}^{+}(\widehat{\textsc{W}})-\lVert\widehat{\textsc{W}}-\textsc{W}\rVert}\leq\frac{1}{\lambda_{m}(\textsc{W})-2\lVert\widehat{\textsc{W}}-\textsc{W}\rVert}.$

where we have that $\lambda_{m}(\textsc{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
$\textsc{W}=\mathbb{E}[\widehat{\textsc{W}}]/\overline{w}$ and
$\textsc{C}=\mathbb{E}[\widehat{\textsc{C}}]/\overline{w}$, one can readily
obtain stronger version of Theorem 1 in the general RKHS setting of [21].

##  Appendix C Unbiased deep learning of spectral features

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 as

$\mathcal{E}_{\gamma}(\theta){:=}\underbrace{\lVert(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}{-}\mathcal{Z}_{\theta}\Lambda^{\theta}_{\eta}\mathcal{Z}_{\theta}^{*}\rVert^{2}_{{\rm{HS}}\left(\mathcal{H}^{\eta}_{\pi}\right)}\\!{-}\lVert(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}\rVert^{2}_{{\rm{HS}}\left(\mathcal{H}^{\eta}_{\pi}\right)}}_{\text{expected
loss
}\mathcal{E}}\\!{+}\gamma\\!\underbrace{\textstyle{\sum_{i,j\in[m]}\\!({\langle
z^{\theta}_{i},z^{\theta}_{j}\rangle}_{L^{2}_{\pi}}\\!\\!\\!{-}\delta_{i,j})^{2}}}_{\text{orthonormality
loss }{\mathcal{E}}_{\rm on}},$ (34)

and recall that the injection operator is
$\mathcal{Z}_{\theta}{=}(z^{\theta})^{\raisebox{-0.75346pt}{\scalebox{0.6}{$\mathsf{T}$}}}(\cdot)\colon\mathbb{R}^{m}\to\mathcal{H}^{\eta}_{\pi}(\mathcal{X})$.
It is easy to show from the definition of the adjoint that, for every
$f\in\mathcal{H}^{\eta}_{\pi}(\mathcal{X})$, we have

$\mathcal{Z}_{\theta}^{*}f=\mathbb{E}_{x\sim\pi}[z^{\theta}(x)((\eta
I{-}\mathcal{L})f)(x)].$ (35)

Thus,
$\mathcal{Z}_{\theta}^{*}\mathcal{Z}_{\theta}=\mathfrak{E}^{\eta}[z^{\theta}_{i},z^{\theta}_{j}]_{i,j\in[m]}$
which we denote by $\textsc{W}_{\theta}$, while $\mathcal{Z}_{\theta}^{*}(\eta
I-\mathcal{L})\mathcal{Z}_{\theta}=\mathbb{E}_{x\sim\pi}[z^{\theta}_{i}(x)z^{\theta}_{j}(x)]_{i,j\in[m]}$
is denoted by $\textsc{C}_{\theta}$.

See 2

###### Proof.

Let
$P_{k}\colon\mathcal{H}^{\eta}_{\pi}(\mathcal{X})\to\mathcal{H}^{\eta}_{\pi}(\mathcal{X})$
be spectral projector of $\mathcal{L}$ corresponding to the $k$ largest
eigenvalues. Now, consider

$\mathcal{E}^{k}(\theta){=}\lVert P_{k}(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}{-}\mathcal{Z}_{\theta}\Lambda^{\theta}_{\eta}\mathcal{Z}_{\theta}^{*}\rVert^{2}_{{\rm{HS}}\left(\mathcal{H}^{\eta}_{\pi}\right)}{-}\lVert
P_{k}(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}\rVert^{2}_{{\rm{HS}}\left(\mathcal{H}^{\eta}_{\pi}\right)}.$

Due to Eckhart-Young theorem, we have that for every $\theta\in\Theta$, the
best rank-$m$ approximation of $(\eta I\\!-\\!\mathcal{L})^{\\!-\\!1}$ is
$(\eta I\\!-\\!\mathcal{L})^{\\!-\\!1}P_{m}=(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}P_{k}P_{m}$, for $k>m$, and it holds

$\mathcal{E}^{k}(\theta){\geq}\sum_{j=m+1}^{k}(\eta-\lambda_{i})^{-1}-\sum_{j=1}^{k}(\eta-\lambda_{i})^{-1}=-\sum_{j=1}^{m}(\eta-\lambda_{i})^{-1}.$

As before, expanding in $\mathcal{E}^{k}$ the HS norm via the trace, we obtain

$\mathcal{E}^{k}(\theta){=}\lVert\mathcal{Z}_{\theta}\Lambda^{\theta}_{\eta}\mathcal{Z}_{\theta}^{*}\rVert^{2}_{{\rm{HS}}\left(\mathcal{H}^{\eta}_{\pi}\right)}{-}2\operatorname{\text{\rm
tr}}[\mathcal{Z}_{\theta}\Lambda^{\theta}_{\eta}\mathcal{Z}_{\theta}^{*}(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}P_{k}]{=}\mathcal{E}(\theta){+}2\operatorname{\text{\rm
tr}}[\mathcal{Z}_{\theta}\Lambda^{\theta}_{\eta}\mathcal{Z}_{\theta}^{*}(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}(I{-}P_{k})],$

and, hence, by Cauchy-Schwartz inequality, we have that

$\lvert\mathcal{E}(\theta)-\mathcal{E}^{k}(\theta)\rvert\leq\lVert\mathcal{Z}_{\theta}\rVert^{2}_{{\rm{HS}}\left(\mathcal{H},\mathcal{H}^{\eta}_{\pi}\right)}\lVert\Lambda^{\theta}_{\eta}\rVert\lVert(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}(I{-}P_{k})\rVert=\tfrac{1}{\eta}\sum_{i\in[m]}\mathfrak{E}^{\eta}_{\pi}[z^{\theta}_{i}](\eta{-}\lambda_{k+1})^{-1}.$

Observing that $z^{\theta}_{i}\in\mathcal{H}^{\eta}_{\pi}(\mathcal{X})$, i.e.
$\mathfrak{E}^{\eta}_{\pi}[z^{\theta}_{i}]<\infty$, we conclude that, for
every $\theta\in\Theta$,
$\displaystyle{\lim_{k\to\infty}}\mathcal{E}^{k}(\theta)=\mathcal{E}(\theta)$.
Therefore, noting that $\mathcal{E}_{\alpha}(\theta)\geq\mathcal{E}(\theta)$,
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
$P_{m}(\eta I\\!-\\!\mathcal{L})^{\\!-\\!1}$, i.e.
$(z^{\theta}_{i})_{i\in[m]}$ span the leading invariant subspace
$(f_{i})_{i\in[m]}$ of the generator. So, if

$\mathcal{E}_{\alpha}(\theta)=\mathcal{E}(\theta)=\sum_{j=1}^{m}(\eta-\lambda_{i})^{-1}$

and $\alpha>0$, we have that ${\mathcal{E}}_{\rm on}(\theta)=0$, implying that
$(z_{i})_{i\in[m]}$ is an orthonormal basis, and, hence $P_{m}(\eta
I\\!-\\!\mathcal{L})^{\\!-\\!1}$, i.e.
$(z^{\theta}_{i})_{i\in[m]}=\mathcal{Z}_{\theta}\Lambda^{\eta}_{\theta}\mathcal{Z}_{\theta}^{*}$.
The result follows.

To show that

$\mathcal{E}_{\gamma}(\theta)=\mathbb{E}\big{[}\mathcal{E}_{\gamma}^{\widehat{\pi}^{\prime}_{1},\widehat{\pi}^{\prime}_{2}}(\theta)\big{]}/\overline{w}^{2},$

we rewrite (18) to encounter the distribution change, noting that the
empirical covariances are reweighted but not normaliezed by $\overline{w}$.
So, we have that

$\mathcal{E}_{\gamma}(\theta)=\operatorname{\text{\rm
tr}}\left[\overline{w}^{-2}\mathbb{E}[\widehat{\textsc{W}}_{\theta}]\Lambda^{\eta}_{\theta}\mathbb{E}[\widehat{\textsc{W}}_{\theta}]\Lambda^{\eta}_{\theta}-2\overline{w}^{-1}\mathbb{E}[\textsc{C}_{\theta}]\Lambda^{\eta}_{\theta}+\overline{w}^{-2}\alpha(\mathbb{E}[\widehat{\textsc{C}}_{\theta}]-\overline{w}\textsc{I})^{2}\right].$

But, since

$\mathbb{E}_{x\sim\pi}[f(x)g(x)]{=}\mathbb{E}_{x^{\prime}\sim\pi^{\prime}}\\!\\!\left[\tfrac{d\pi}{d\pi^{\prime}}(x^{\prime})\,f(x^{\prime})\,g(x^{\prime})\right]{=}\mathbb{E}_{x_{1}^{\prime}\sim\pi^{\prime}}\\!\big{[}\sqrt{\tfrac{d\pi}{d\pi^{\prime}}(x_{1}^{\prime})}f(x_{1}^{\prime})\big{]}\,\mathbb{E}_{x_{2}^{\prime}\sim\pi^{\prime}}\\!\big{[}\sqrt{\tfrac{d\pi}{d\pi^{\prime}}(x_{2}^{\prime})}g(x_{2}^{\prime})\big{]},$

where $x_{1}^{\prime}$ and $x_{2}^{\prime}$ are two independent r.v. with a
law $\pi^{\prime}$, the proof is completed. ∎

##  Appendix D Training of a neural network model

###  D.1 Evolution of the loss with eigenfunctions

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
3. 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.

Figure 3: Typical behavior of the loss function during a training.

###  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 1. In Figure
4, 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.

Figure 4: Comparison on Muller Brown potential with ground truth, learned
features and fitted eigenfunctions

###  D.3 Activation functions and structure of the neural network

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.

###  D.4 Hyperparameters

Besides common hyperparameters such as learning rate, neural network
architecture and activation function, our method requires only two
hyperparameters: $\eta$ and $\alpha$. Other methods such as [47] do not
require $\eta$, but on the other hand requires one weight per searched
eigenfunction.

##  Appendix 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 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].

###  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 $U_{\rm
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_{\rm sim}=4(-0.5\exp(-80x^{2})+x^{8})$, which thus makes a bias
potential: $V_{\rm bias}(x)=U_{\rm sim}(x)-U_{\rm tg}(x)=-4(\exp(-80x^{2}))$.
In Figure 5, 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, $\eta=0.1$ and a regularization parameter of
$10^{-5}$

Figure 5: 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)

###  E.2 Muller Brown potential

For 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 a temperature of 1 (arbitrary unit)
and a timestep of $10^{-3}$. The simulation was run for $10^{7}$ timesteps.
The bias potential is built according to the following equation

$U(x,t)=h\sum_{i=1}^{N_{t}}e^{-\frac{\|x-x_{i}\|^{2}}{2\sigma^{2}}}$ (36)

where 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 $\eta$ was chosen to be
$0.05$.

###  E.3 Alanine dipeptide

Figure 6: Time evolution of the $\phi$ angle in both datasets. Points are
colored according to the value of the first nontrivial eigenfunction

Simulation details  All the simulations are run with GROMACS 2022.3 [2] and
patched with plumed 2.10 [43] in order to perform enhanced sampling
simulations. We used the Amber99-SB [38] force field. The Langevin dynamics
was sampled with a timestep of 2fs with a damping coefficent
$\gamma_{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 [47]. 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 $\eta$
was chosen to be $0.1$.