# Analyzing multimodal probability measures with autoencoders

Tony Lelièvre,1,2 Thomas Pigeon,2,3 Gabriel Stoltz1,2, and Wei Zhang4,5
1: CERMICS, École des Ponts ParisTech, 6-8 Avenue Blaise Pascal, 77455,Marne-
la-Vallée, France, France
2: MATHERIALS project-team, Inria Paris, 2 Rue Simone Iff, 75012 Paris, France
3: IFP Energies Nouvelles, Rond-Point de l’Echangeur de Solaize, 69360
Solaize, France
4: Department of Mathematics and Computer Science, Freie Universität Berlin,
Arnimallee 14, 14195 Berlin, Germany
5: Zuse Institute Berlin, Takustraße 7, 14195 Berlin, Germany

###### Abstract

Finding collective variables to describe some important coarse-grained
information on physical systems, in particular metastable states, remains a
key issue in molecular dynamics. Recently, machine learning techniques have
been intensively used to complement and possibly bypass expert knowledge in
order to construct collective variables. Our focus here is on neural network
approaches based on autoencoders. We study some relevant mathematical
properties of the loss function considered for training autoencoders, and
provide physical interpretations based on conditional variances and minimum
energy paths. We also consider various extensions in order to better describe
physical systems, by incorporating more information on transition states at
saddle points, and/or allowing for multiple decoders in order to describe
several transition paths. Our results are illustrated on toy two dimensional
systems and on alanine dipeptide.

##  1 Introduction

Molecular simulation has proved to be an effective computational approach to
understand complex systems in biophysics, chemistry, material science, etc.
[47, 26]. However, sampling complex physical systems in molecular dynamics
remains a computational challenge due to the metastability of the system’s
dynamics – by which we mean that the system is often stuck for a long time in
some state (called “metastable state”), before switching to another state
where it again remains stuck for a long time. Many efforts have been, and
still are pursued to enhance the sampling of molecular systems. One successful
strategy to this end relies on free energy biasing where (a fraction of) the
free energy is used as an importance sampling function [31]. This paradigm
underpins adaptive methods to compute the free energy such as the celebrated
metadynamics and its variants [58, 11], and the adaptive biasing force method
[23], to quote only two options (see [42] for a more extensive review). A key
element in these free energy biasing methods is the choice of collective
variables (CVs), which summarize some coarse-grained information on the
system, and describe in a lower dimensional space the various metastable
states.

CVs were chosen until recently based on expert knowledge and chemical
intuition. Following the many successes of machine learning approaches in data
science and computer vision, the last decade witnessed many propositions to
improve, complement, or even replace expert knowledge in molecular dynamics.
There are by now various reviews on the use of machine learning techniques in
molecular dynamics, including also the construction of empirical force fields
(see [95] which summarizes various other review papers). Dedicated review
articles are devoted to finding CVs [30, 84, 74, 38, 39, 16]; see also [82]
for more mathematically oriented elements. These methods mostly correspond to
unsupervised dimensionality reduction and manifold learning techniques,
although some of them are based on (semi)supervised approaches.

In order to make the best use of machine learning methods, one needs to
specify criteria which make CVs appropriate. Desirable requirements include:

  * •

the transition path in the collective variable space from one metastable state
to another should correspond to some form of minimum energy or free energy
path [90];

  * •

the CVs provide a parametrization of the eigenfunctions associated with the
dominant eigenvalues of the transfer operators (or the generator of the
dynamics), as these eigenfunctions allow to characterize metastable states
[78];

  * •

free energy biasing based on these CVs efficiently suppresses or reduces
metastability. From a mathematical viewpoint, this can be quantified through
improved constants in functional inequalities (see for instance [61] for a
pedagogical introduction).

Some CVs meet only some of these requirements. For instance, eigenfunctions
have sharp transitions between metastable states, and may therefore not be
used as such as CVs since otherwise the biasing forces used in free energy
biasing dynamics, which involve the gradient of the collective variable, would
be large and lead to unstable dynamics. The same remark applies to the
committor function.

Most machine learning methods to find CVs are unsupervised (see however [9]
for a supervised method based on Fisher’s linear discriminant analysis, and
[6] for a semi-supervised approach). There are two main classes of
unsupervised methods: those seeking high variances CVs, which aim at
reproducing overall features of the Boltzmann–Gibbs distribution at hand; and
those seeking slowly evolving CVs (also see [99] for related discussions).
Both classes can be separated into linear and nonlinear methods. Our focus in
this work is on methods that find high variance CVs, for which no notion of
dynamics is a priori used. It is useful to further distinguish between (i)
linear methods, e.g. principal component analysis (PCA) or factor analysis,
and (ii) non-linear algorithms, e.g. kernel methods, autoencoders, decision
trees and random forests (see for instance [68, 71] for introductory
references on these classes of methods). In the first case, CVs are
interpretable, but often too simple to give precise results.

#### Aims and scope.

We focus in this work on autoencoders, which have been successfully used in
molecular dynamics (see [14] and Section 2.1 for reviews). Our aim is to
better characterize the learning problem associated with autoencoders, and
relate it to various frameworks relevant for statistical methods (principal
curves and manifolds) and/or molecular dynamics (use of conditional
expectations as for free energy computations). We also aim at better
understanding what is learned when minimizing the reconstruction error, and
how/whether the so-obtained information can be leveraged. This can be done for
example by (i) adding extra terms in the loss function, to describe
transitions from one metastable state to another, following up on [79]; (ii)
by allowing for multiple decoders in order to represent multiple transition
paths connecting two metastable states; or (iii) by requiring the encoder to
parametrize dominant eigenfunctions.

#### Outline.

We briefly present autoencoders and some of their important properties in
Section 2. A focus on their interpretation in terms of conditional
expectations is provided in Section 3. We then first consider a supervised
setting where the probability measure on the data has only two modes which are
known a priori. Section 4 discusses options to construct a single path
relating the two modes, and incorporating useful information on the saddle
points and transition states observed during the transition. The situation
when the modes can be related via several paths is then addressed in Section
5. We consider in Section 6 an application to alanine dipeptide. We finally
discuss in Section 7 some perspectives of this work.

Let us emphasize that our baseline assumption in all this work is that the
dataset at hand is sufficiently rich to correctly describe all the modes of
the probability measure, and sometimes maybe also the transition zones. The
latter situation is somewhat unrealistic for actual systems of interest, but
can be achieved by an iterative procedure where one cycles between exploration
phases using a dynamics biased by the free energy, and a learning phase to
update the collective variables (as done in various works; see references in
Section 2.1).

##  2 Autoencoders and their properties

We present more precisely autoencoders in Section 2.1, with a focus on
bottleneck autoencoders and applications in molecular dynamics. We then
discuss various interpretations of physical relevance of the loss function in
Section 2.2.

###  2.1 Presentation of autoencoders

Autoencoders have been considered early on in the neural network literature,
where they were also called auto-associative neural networks [57]. The models
considered in these early works correspond to what is currently known as
bottleneck autoencoders, and were rather shallow. Deep autoencoders were used
later on with the advent of modern computing architectures [45]. Bottleneck
autoencoders were first introduced to provide a nonlinear generalization of
PCA, as it was shown that the linear neural networks obtained by minimizing
the mean-square error were essentially equivalent to PCA [10, 5]. We refer for
instance to [8, Section 12.4.2], [40, Chapter 14] and [71, Section 20.3] for
textbook presentations of autoencoders, which include some historical
perspectives, and present the many variations and extensions that were
considered.

Autoencoders fall into the class of unsupervised machine learning methods. For
a given input data point $x\in\mathcal{X}\subset\mathbb{R}^{D}$, we denote by
$f_{\theta}(x)$ the prediction of the neural network. The parameters
$\theta\in\Theta$ are chosen to minimize the loss function

$\mathscr{L}(\theta)=\mathbb{E}[\ell(X,f_{\theta}(X))],$ (1)

where $\ell$ is a given elementary loss function, and the expectation is over
the realizations of the input data $X$ distributed according to some
probability measure denoted by $\mu$. By default, in the sequel, expectations
are always taken with respect to the distribution of the data. The typical
choice for the latter elementary loss function is the square loss
$\ell(x,y)=\|x-y\|^{2}$, although other choices, such as the mean absolute
loss $\ell(x,y)=\|x-y\|$ could also be considered in order to give less weight
to outliers. In practice, the population loss $\mathscr{L}$ is replaced by the
empirical loss over a training set of $N_{\mathrm{data}}$ given input data
points $\\{x^{1},\dots,x^{N_{\mathrm{data}}}\\}$:

$\widehat{\mathscr{L}}(\theta)=\frac{1}{N_{\mathrm{data}}}\sum_{n=1}^{N_{\mathrm{data}}}\ell(x^{n},f_{\theta}(x^{n})).$

#### Families of autoencoders.

There are various classes of autoencoders. It is useful to distinguish between
undercomplete and overcomplete models. Undercomplete models have a limited
capacity that prevents them from achieving zero training loss. The most
prominent example is provided by bottleneck autoencoders for which

$f_{\theta}=f_{{\rm dec},\theta_{2}}\circ f_{{\rm enc},\theta_{1}},$ (2)

where the parameters $\theta=(\theta_{1},\theta_{2})$ have been decomposed
into parameters used in the encoder and decoder parts, respectively (see
Figure 1 below), and $\circ$ is the composition operator, namely

$h_{2}\circ h_{1}(x)=h_{2}(h_{1}(x)).$

The limitation in the capacity of the autoencoder arises from the fact that
the encoding function $f_{{\rm enc},\theta_{1}}$ has values in a latent space
$\mathcal{Z}\subset\mathbb{R}^{d}$ of dimension $d$ strictly smaller than the
dimension $D$ of the input/output space $\mathcal{X}$, usually much smaller in
fact. Overcomplete models can on the other hand achieve zero training error.
Such models are of course useless as such since they would simply copy the
input without extracting the salient features explaining the data at hand.
This is why some regularization process should be considered to limit the
capacity of the neural network. Regularization is however also useful for
undercomplete models. Standard examples of regularization mechanisms include:

  * •

resorting to regularization strategies commonly used for neural networks in
general, in particular early stopping, dropout, and standard weight decay to
name a few options (see for instance [71, Section 13.5] and [40, Chapter 7]
for a more precise presentation of these options);

  * •

sparse autoencoders where a penalization term is added to the loss function to
prevent too many neurons to be active [73];

  * •

denoising autoencoders [91], where outputs should be predicted from inputs
corrupted by some noise, which forces networks to learn the structure of the
distribution of the data [1];

  * •

contractive autoencoders, where the Jacobian of the encoder is penalized in
order to ensure smoother variations of this function, and limit its
sensitivity to small variations in the input [81];

  * •

variational autoencoders (VAEs) [52, 53, 37] can also be interpreted as some
regularized version of the usual autoencoder framework.

Other variations/extensions of autoencoders were also considered for more
specific purposes, and are still studied, in particular for manifold learning,
where the preservation of neighborhood relationships and/or geometric
information in the dimensionality reduction process are important [49, 25,
59].

In this work, we will consider bottleneck autoencoders only, as these are the
most relevant for molecular dynamics.

#### Autoencoders architectures.

Autoencoders are often symmetric in their structures. In some cases, tied
weights are being used, _i.e._ the weights $\theta_{2}$ are the transpose of
the weights $\theta_{1}$ when writing the prediction function as (2). This
choice reproduces at the level of autoencoders the symmetric structure of PCA,
where the decoder matrix is the transpose of the encoder matrix (see (20)
below and the discussion around this equation). In this case, a regularization
on the encoder part only, as in contractive autoencoders, in fact also
regularizes the decoder part. However, there is no particular motivation to
use symmetric architectures, and various works, such as [86], studied the
impact of an asymmetric architecture on the performance of the model.

Another important consideration, more specific to bottleneck autoencoders, is
the choice of the bottleneck dimension. In PCA, the number of dimensions to
retain usually corresponds to some “knee” in the plot of eigenvalues of the
empirical covariance matrix. Similar plots can be performed for bottleneck
autoencoders, for which the reconstruction error can be reported as a function
of the bottleneck dimension (the remainder of the architecture being fixed).
The optimal dimension should ideally coincide with the intrinsic dimension of
the data set in manifold learning [22], which can be quantified using the
Frechet inception distance (see for instance the method described in [48] in
the context of VAEs).

#### Autoencoders in molecular dynamics.

In molecular dynamics, only bottleneck autoencoders are considered. In the
remainder of this work, we simply refer to them as ”autoencoders”. They have
been used for a few years now to find CVs, by using the encoder part as a CV.
In this situation, the decoder is not useful as such, but is still required to
learn a good encoder.

Autoencoders were first considered in a static framework, to extract
information from samples distributed according to the Boltzmann–Gibbs
distribution of the system. In fact, since the systems under consideration are
metastable, it is often difficult to directly obtain a good sampling of the
target distribution, and one should therefore turn to an iterative procedure
alternating between a free energy biased sampling to obtain new data points,
and a subsequent update of the collective variable given by the encoder. Such
strategies were considered for applications in biophysics and materials
science in unsupervised [17, 20, 7, 15, 4] and semi-supervised [6] settings,
possibly also with variational autoencoders [80, 93].

More recently, autoencoders were also considered to learn the dynamics of the
systems under consideration. Besides methods for generic dynamical systems
[13], dedicated approaches were introduced for the dynamics considered in
molecular simulation. Examples include time-lagged autoencoders [94, 18] and
time-lagged variational autoencoders [44]; VAMPnets and their subsequent
extensions and modifications [67, 96, 19, 85], which can be seen as versions
of variational approaches to conformational dynamics [75, 56] using a basis of
functions represented by neural networks (in fact, only their encoder parts),
instead of using linear functions as in tICA [69, 72]; simultaneous
reconstruction of the probability measure and learning of the committor
function [33]; approaches based on recurrent neural networks, in particular
long short term memory networks [46], such as [88], or mixture density
networks [92]. All these methods differ in the details of the architecture and
choice of loss function (e.g. time-lagged reconstruction loss for time-lagged
autoencoders vs. autocorrelation of slow modes for VAMPnets and their
extensions/modifications).

It is often the case that no explicit regularization is introduced as the
networks under consideration are rather shallow and narrow. Some works however
explicitly introduce such regularization terms: for instance, an elastic net
penalization is considered in [50, 51], and dropout plus weight decay in [67].
In addition, it is sometimes suggested to add physically motivated terms to
the loss function, for instance to favor the exploration of new metastable
modes [50, 51], or find minimum energy paths [79].

######  Remark 1.

In contrast to applications of ML for which the distribution of the data is
typically unknown (see for instance [1] where it is shown how to estimate the
gradient of the log-density using autoencoders), an important difference in
the use of autoencoders in MD is that the distribution of the training data
generally has a known density (typically, the Boltzmann–Gibbs distribution).
Although standard losses to train autoencoders do not make use of the
expression of the density of the data, it might be possible to leverage this
information in the MD context to find better encodings/decodings than in
standard ML contexts.

###  2.2 Interpretations of the loss function

We discuss in this section various reformulations and reinterpretations of the
loss function (1) for bottleneck autoencoders (2) when the elementary loss
function is the square loss, and discuss in particular the relationship with
principal curves/manifolds [41, 87].

#### Three viewpoints on the loss function.

We consider an ideal setting where we minimize upon all measurable functions
$f_{\mathrm{enc}}:\mathcal{X}\to\mathcal{Z}$ and
$f_{\mathrm{dec}}:\mathcal{Z}\to\mathcal{X}$. We denote by
$\mathcal{F}_{\mathrm{enc}}$ and $\mathcal{F}_{\mathrm{dec}}$ the sets of
measurable functions from $\mathcal{X}$ to $\mathcal{Z}$ and from
$\mathcal{Z}$ to $\mathcal{X}$, respectively; and by $\mathcal{F}$ the set of
measurable functions obtained by composing functions of
$\mathcal{F}_{\mathrm{enc}}$ with functions of $\mathcal{F}_{\mathrm{dec}}$:

$\mathcal{F}=\left\\{f=f_{\mathrm{dec}}\circ f_{\mathrm{enc}},\
f_{\mathrm{enc}}\in\mathcal{F}_{\mathrm{enc}},\
f_{\mathrm{dec}}\in\mathcal{F}_{\mathrm{dec}}\right\\}.$

Minimizing the reconstruction error over the set of functions in $\mathcal{F}$
can be then rewritten as

$\inf_{f\in\mathcal{F}}\mathbb{E}\left[\left\|X-f(X)\right\|^{2}\right].$ (3)

Note that we do not consider a regularization term here, so that overfitting
may occur in practice (for instance, even with $\mathcal{Z}$ of dimension 1,
$f_{\mathrm{dec}}$ can parametrize a space-filling curve).

As discussed in [34] (which complements [36] which was already hinting at
autoencoders), the unsupervised least-square problem (3) can be thought of in
various way. In particular, there is some duality in the way the minimization
over $f\in\mathcal{F}$ is performed, as one can decide to either

  1. (i)

simultaneously minimize over $f_{\mathrm{enc}}$ and $f_{\mathrm{dec}}$, which
is the standard way to proceed when training neural networks;

  2. (ii)

minimize first over the encoder part, which allows to reformulate the
minimization as the well-known problem of finding principal manifolds;

  3. (iii)

minimize first over the decoder part, which is natural when thinking of the
reconstruction error as some total variance to be decomposed using a
conditioning on the values of the encoder.

We discuss in the remainder of this section option (ii), leaving the
discussion of (iii) to Section 3.

The chosen numerical approach has a natural impact on the topology of the
networks which are considered: in situation (i), encoders and decoders are
treated on an equal footing, and it is therefore natural to consider them to
be of a similar complexity; whereas options (ii) and (iii) suggest to consider
asymmetric autoencoders. For instance, in option (iii), the minimization over
the decoder part, which is performed first, could be done more carefully, with
more expressive networks in order to better approximate the optimal decoder
for a given encoder. We numerically investigate this point in Section 3.3.

#### Principal manifold reformulation.

We start by minimizing the reconstruction error (3) over the encoder part for
a given decoder:

$\displaystyle\inf_{f\in\mathcal{F}}\mathbb{E}\left[\left\|X-f(X)\right\|^{2}\right]$
$\displaystyle=\inf_{f_{\mathrm{dec}}\in\mathcal{F}_{\mathrm{dec}}}\left\\{\inf_{f_{\mathrm{enc}}\in\mathcal{F}_{\mathrm{enc}}}\mathbb{E}\left[\left\|X-f_{\mathrm{dec}}\circ
f_{\mathrm{enc}}(X)\right\|^{2}\right]\right\\}$
$\displaystyle=\inf_{f_{\mathrm{dec}}\in\mathcal{F}_{\mathrm{dec}}}\mathbb{E}\left[\left\|X-f_{\mathrm{dec}}\circ
h^{\star}_{f_{\mathrm{dec}}}(X)\right\|^{2}\right],$ (4)

where the optimal encoder
$h^{\star}_{f_{\mathrm{dec}}}:\mathcal{X}\to\mathcal{Z}$ for a given decoder
$f_{\mathrm{dec}}:\mathcal{Z}\to\mathcal{X}$ is defined pointwise as

$h^{\star}_{f_{\mathrm{dec}}}(x)\in\mathop{\mathrm{argmin}}_{z\in\mathcal{Z}}\|x-f_{\mathrm{dec}}(z)\|,$

provided that this minimization problem admits a solution. When
$f_{\mathrm{dec}}$ is smooth and has a non-vanishing Jacobian, the principal
manifold is then the set $f_{\mathrm{dec}}(\mathcal{Z})\subset\mathcal{X}$.
For any $x\in\mathcal{X}$, $h^{\star}_{f_{\mathrm{dec}}}(x)\in\mathcal{Z}$
gives the coordinates in the latent space of the projection of $x$ on the
principal manifold (the so-called projection index).

The reformulation of the minimization problem in terms of the decoder function
leads to a minimization problem similar to the one encountered when searching
for principal curves and manifolds. These mathematical concepts generalize in
some sense PCA to curves and surfaces rather than lines and hyperplanes, as
already discussed in the work introducing principal curves [41]. The
minimization problems associated with finding principal manifolds are in
general difficult to solve as these manifolds correspond to saddle-points of
the loss functional [24]. In practice, this means that it is possible to move
away from a critical point without increasing the test loss. This leads to
overfitting issues and prevents traditional cross-validation procedures to
tune regularization hyperparameters.

In the small temperature regime, principal curves, associated with decoders
starting from a one-dimensional space $\mathcal{Z}$, reduce to minimum energy
paths when the measure under consideration is the Boltzmann–Gibbs measure; see
the discussion in [90]. Another perspective on this statement is provided at
the end of Section 3.1, using the interpretation relying on conditional
expectations.

##  3 Reformulating autoencoders with conditional expectations

We first show in Section 3.1 how to reformulate the training of autoencoders
with conditional expectations, and provide alternative interpretations to the
reconstruction error. Such reformulations can in fact be considered for other
unsupervised machine learning methods such as principal component analysis and
clustering, as we discuss in Section 3.2. We finally numerically illustrate in
Section 3.3, on a simple two-dimensional example, the interpretation of the
reconstruction loss with conditional expectations.

###  3.1 Reformulation using conditional expectations

In contrast to (4), we minimize here the reconstruction error (3) by first
minimizing over the decoder part for a given encoding function, as already
considered in [35]. This approach is natural in molecular dynamics, as it is
reminiscent of free energy computations [21, 62] where average quantities are
computed for a fixed value of the collective variable $f_{\mathrm{enc}}$. From
a mathematical viewpoint, it corresponds to introducing conditional averages
associated with fixed values of the encoder. More generally, various
unsupervised learning methods can be understood in this manner, see Section
3.2.

#### Rewriting the reconstruction error with conditional expectations.

The loss function for unsupervised least-squares can be rewritten as

$\displaystyle\inf_{f\in\mathcal{F}}\mathbb{E}\left[\left\|X-f(X)\right\|^{2}\right]$
$\displaystyle=\inf_{f_{\mathrm{enc}}\in\mathcal{F}_{\mathrm{enc}}}\left\\{\inf_{f_{\mathrm{dec}}\in\mathcal{F}_{\mathrm{dec}}}\mathbb{E}\left[\left\|X-f_{\mathrm{dec}}\circ
f_{\mathrm{enc}}(X)\right\|^{2}\right]\right\\}$
$\displaystyle=\inf_{f_{\mathrm{enc}}\in\mathcal{F}_{\mathrm{enc}}}\mathbb{E}\left[\left\|X-g^{\star}_{f_{\mathrm{enc}}}\circ
f_{\mathrm{enc}}(X)\right\|^{2}\right],$ (5)

where the ideal decoder $g^{\star}_{f_{\mathrm{enc}}}$ for a given encoder
$f_{\mathrm{enc}}$ is the Bayes predictor associated with the least square
regression problem (see [3, Section 2.2.3]):

$g^{\star}_{f_{\mathrm{enc}}}(z)=\mathbb{E}[\,X\,|\,f_{\mathrm{enc}}(X)=z].$
(6)

Let us recall that, in all these expressions, expectations are taken with
respect to the probability distribution $\mu$ of the input data (which is not
necessarily the Boltzmann–Gibbs distribution). Equations (5)-(6) show that the
question of finding the best autoencoder can be reduced to finding the best
encoding function, provided that one is able to compute good approximations of
the conditional expectation.

#### Derivation of (6).

For completeness, let us recall here the derivation of (6), as this derivation
also allows to reformulate the minimization problem under consideration as
some variance maximization for conditional expectations. We start by
introducing the quantity $g^{\star}_{f_{\mathrm{enc}}}\circ
f_{\mathrm{enc}}(x)$ in the reconstruction error: for any
$f_{\mathrm{dec}}\in\mathcal{F}_{\mathrm{dec}}$ and
$f_{\mathrm{enc}}\in\mathcal{F}_{\mathrm{enc}}$,

$\displaystyle\mathbb{E}\left[\left\|X-f_{\mathrm{dec}}\circ
f_{\mathrm{enc}}(X)\right\|^{2}\right]=\mathbb{E}\left[\left\|\left[X-g^{\star}_{f_{\mathrm{enc}}}\circ
f_{\mathrm{enc}}(X)\right]+\left[g^{\star}_{f_{\mathrm{enc}}}\circ
f_{\mathrm{enc}}(X)-f_{\mathrm{dec}}\circ
f_{\mathrm{enc}}(X)\right]\right\|^{2}\right]$
$\displaystyle\qquad=\mathbb{E}\left[\left\|X-g^{\star}_{f_{\mathrm{enc}}}\circ
f_{\mathrm{enc}}(X)\right\|^{2}\right]+\mathbb{E}\left[\left\|g^{\star}_{f_{\mathrm{enc}}}\circ
f_{\mathrm{enc}}(X)-f_{\mathrm{dec}}\circ
f_{\mathrm{enc}}(X)\right\|^{2}\right],$ (7)

where we used the following identity, obtained by conditioning on the values
of the random variable $Z=f_{\mathrm{enc}}(X)$:

$\displaystyle\mathbb{E}\left[\left(X-g^{\star}_{f_{\mathrm{enc}}}\circ
f_{\mathrm{enc}}(X)\right)\cdot\left(g^{\star}_{f_{\mathrm{enc}}}\circ
f_{\mathrm{enc}}(X)-f_{\mathrm{dec}}\circ f_{\mathrm{enc}}(X)\right)\right]$
$\displaystyle\qquad=\mathbb{E}\left[\left(\mathbb{E}\left[X\middle|Z\right]-g^{\star}_{f_{\mathrm{enc}}}(Z)\right)\cdot\left(g^{\star}_{f_{\mathrm{enc}}}(Z)-f_{\mathrm{dec}}(Z)\right)\right]=0,$

in view of the definition (6) of $g^{\star}_{f_{\mathrm{enc}}}$. It is clear
from (7) that the decoding function which minimizes the reconstruction error
for a given encoder $f_{\mathrm{enc}}\in\mathcal{F}_{\mathrm{enc}}$ is indeed
$g^{\star}_{f_{\mathrm{enc}}}$, as defined by (6).

#### Alternative interpretation of the reconstruction error.

The reconstruction error (5) can be reinterpreted in terms of variances.
Indeed, on the one hand,

$\displaystyle\mathbb{E}\left[\left\|X-g^{\star}_{f_{\mathrm{enc}}}\circ
f_{\mathrm{enc}}(X)\right\|^{2}\right]$
$\displaystyle=\mathbb{E}\left[\left\|X-\mathbb{E}\left[X|f_{\mathrm{enc}}(X)\right]\right\|^{2}\right]=\mathbb{E}\left[\mathbb{E}\left(\left.\left\|X-\mathbb{E}\left[X|f_{\mathrm{enc}}(X)\right]\right\|^{2}\right|f_{\mathrm{enc}}(X)\right)\right]$
$\displaystyle=\mathbb{E}\left[\mathrm{Var}(X|f_{\mathrm{enc}}(X))\right].$
(8)

On the other hand,

$\displaystyle\mathbb{E}\left[\left\|X-g^{\star}_{f_{\mathrm{enc}}}\circ
f_{\mathrm{enc}}(X)\right\|^{2}\right]$
$\displaystyle=\mathbb{E}\left[\left\|X-\mathbb{E}\left[X|f_{\mathrm{enc}}(X)\right]\right\|^{2}\right]=\mathbb{E}\left(\|X\|^{2}\right)-\mathbb{E}\left(\mathbb{E}\left[X|f_{\mathrm{enc}}(X)\right]^{2}\right)$
$\displaystyle=\mathrm{Var}(X)-\mathrm{Var}\left[\mathbb{E}(X|f_{\mathrm{enc}}(X))\right].$
(9)

These two equalities yield the well-known formula for the total variance
decomposition, namely
$\mathrm{Var}(X)=\mathbb{E}\left[\mathrm{Var}(X|f_{\mathrm{enc}}(X))\right]+\mathrm{Var}\left[\mathbb{E}(X|f_{\mathrm{enc}}(X))\right]$.

A consequence of (9) is that the minimization problem (5) can be reformulated
as the following equivalent maximization problem:

$\sup_{f_{\mathrm{enc}}\in\mathcal{F}_{\mathrm{enc}}}\mathrm{Var}\left[\mathbb{E}(X|f_{\mathrm{enc}}(X))\right].$
(10)

In words, this reformulation translates the equivalence between (the ”classes”
referring here to the level sets of $f_{\mathrm{enc}}$)

  * •

minimizing the intraclass dispersion (5): the distribution of configurations
$x\in\mathcal{X}$ for a fixed value $z$ of $f_{\mathrm{enc}}$ should
concentrate around the mean value $g^{\star}_{f_{\mathrm{enc}}}(z)$ by having
a variance as small as possible;

  * •

maximizing the interclass dispersion (10): the values of the conditional
averages for a fixed value of $f_{\mathrm{enc}}$ should be as spread out as
possible over the range of $f_{\mathrm{enc}}$.

#### Formal characterization of the optimal encoding function.

A key equality to establish (9), namely

$\mathbb{E}\left[\left\|X-g^{\star}_{f_{\mathrm{enc}}}\circ
f_{\mathrm{enc}}(X)\right\|^{2}\right]=\mathbb{E}\left[\left\|X\right\|^{2}\right]-\mathbb{E}\left[\left\|g^{\star}_{f_{\mathrm{enc}}}\circ
f_{\mathrm{enc}}(X)\right\|^{2}\right],$

shows that the minimization of the reconstruction error is equivalent to the
maximization of the second moment of the conditional expectation:

$\sup_{f_{\mathrm{enc}}\in\mathcal{F}_{\mathrm{enc}}}\mathbb{E}\left[\left\|g^{\star}_{f_{\mathrm{enc}}}\circ
f_{\mathrm{enc}}(X)\right\|^{2}\right]\,.$ (11)

This alternative viewpoint allows to characterize the optimal encoding
function $f_{\mathrm{enc}}$ by some orthogonality condition similar to the
self-consistency condition of principal curves, see [36, Section 2]. In fact,
it can be formally shown that critical points of (11) satisfy

$\forall j\in\\{1,\dots,d\\},\quad\forall
x\in\mathrm{Supp}(\mu),\qquad\left[x-g^{\star}_{f_{\mathrm{enc}}}(f_{\mathrm{enc}}(x))\right]^{\top}\partial_{z_{j}}g^{\star}_{f_{\mathrm{enc}}}(f_{\mathrm{enc}}(x))=0,$
(12)

where $\mathrm{Supp}(\mu)$ is the support of the probability measure $\mu$.
The derivation of this condition, which can be read in Appendix A, can be seen
as a variation of derivations of optimality conditions for principal curves,
as written already in [41]; see also [36] where (12) is used to construct a
new objective function to minimize to find $f_{\mathrm{enc}}$. From a
technical viewpoint, an originality of our approach to obtain (12) is that we
rely on the co-area formula [29, 2] together with the use of weak derivatives;
see Appendix A.

An interesting implication of (12) is that the intersection of
$\mathrm{Supp}(\mu)$ and the submanifold

$\Sigma_{z}=f_{\mathrm{enc}}^{-1}\\{z\\}=\left\\{x\in\mathcal{X}\,\middle|\,f_{\mathrm{enc}}(x)=z\right\\}$
(13)

is in fact included in the $(D-d)$-dimensional hyperplane containing the point
$g^{\star}_{f_{\mathrm{enc}}}(z)$ and orthogonal to the vectors
$\partial_{z_{1}}g^{\star}_{f_{\mathrm{enc}}}(z),\dots,\partial_{z_{d}}g^{\star}_{f_{\mathrm{enc}}}(z)$
(recalling that $\mathcal{X}$ and $\mathcal{Z}$ have dimensions $D$ and $d$,
respectively). As such hyperplanes generally have a non-empty intersection,
this condition can be true only for distributions $\mu$ which have a support
sufficiently concentrated around the principal manifold. The condition (12)
can be generalized when several principal manifolds are considered, as in
Section 5 here. The issue of having hyperplanes intersecting can be seen as
the counterpart in the context we consider here of the concept of ambiguity
points for principal curves [41]. Let us mention that a condition similar to
(12) is derived in [90] to reinterpret the string method using principal
curves, followed by a discussion on hyperplanes corresponding to isosurfaces
of fixed values of the assignment function (the counterpart of the encoder
function in our context).

#### Minimum energy paths and the small temperature limit.

We finally study the limit $\beta\to+\infty$ in (12), when the probability
measure $\mu$ under consideration is the Boltzmann–Gibbs measure with a
density proportional to $\mathrm{e}^{-\beta V(x)}$, with $V$ the potential
energy function.

Transition paths are one-dimensional curves. It would be possible to work in a
latent space $\mathcal{Z}$ of dimension $d\geqslant 2$, and to extract some
one dimensional path relating the initial and end configurations in this
space. Here, for simplicity, we restrict ourselves to a one-dimensional latent
space $\mathcal{Z}$, _i.e._ $d=1$.

We consider two local minima $x_{A}$ and $x_{B}$ of the potential energy
function $V$, located respectively on $\Sigma_{z_{A}}$ and $\Sigma_{z_{B}}$
with $z_{A}=f_{\mathrm{enc}}(x_{A})$ and $z_{B}=f_{\mathrm{enc}}(x_{B})$
(assuming $z_{A}\leqslant z_{B}$ without loss of generality). Similarly to the
discussion in [90] for principal curves, it can formally be shown in the low
temperature limit that the path
$\\{g^{\star}_{f_{\mathrm{enc}}}(z)\\}_{z\in[z_{A},z_{B}]}$ converges to a
minimum energy path (MEP). A MEP $\\{\gamma(z)\\}_{z\in[z_{A},z_{B}]}$ between
$x_{A}$ and $x_{B}$ is characterized by the boundary conditions
$\gamma(z_{A})=x_{A}$ and $\gamma(z_{B})=x_{B}$, and the fact that
$\gamma^{\prime}(z)$ is collinear to $\nabla V(\gamma(z))$ for
$z\in(z_{A},z_{B})$ (except at critical points along the path); see for
instance [43] as well as [64] for a rigorous mathematical formulation.

In view of (6) and (33) in Appendix A, the conditional expectation
concentrates on a minimum of $V$ on $\Sigma_{z}$ the limit $\beta\to+\infty$.
More precisely, recalling the definition (13) of $\Sigma_{z}$,

$\forall z\in(z_{A},z_{B}),\qquad
g^{\star}_{f_{\mathrm{enc}}}(z)\xrightarrow[\beta\to+\infty]{}g^{\infty}_{f_{\mathrm{enc}}}(z)\in\mathop{\mathrm{argmin}}_{x\in\Sigma_{z}}V(x).$
(14)

We assume in the sequel that the minimum in (14) is well defined (_i.e._ it
exists and is unique). Then, (14) implies that

$\forall z\in(z_{A},z_{B}),\quad\forall t_{z}\in
T\Sigma_{z}\left(g^{\infty}_{f_{\mathrm{enc}}}(z)\right),\qquad
t_{z}^{\top}\nabla V\left(g^{\infty}_{f_{\mathrm{enc}}}(z)\right)=0,$ (15)

where we denote by $T\Sigma_{z}(x)$ the tangent plane to $\Sigma_{z}$ at
$x\in\Sigma_{z}$.

Moreover, for a perfectly converged autoencoder model, we obtain, by passing
in (12) to the limit $x\rightarrow g^{\star}_{f_{\mathrm{enc}}}(z)$ for
$x\in\Sigma_{z}$:

$\forall z\in\mathcal{Z},\quad\forall\,t_{z}\in
T\Sigma_{z}\left(g^{\star}_{f_{\mathrm{enc}}}(z)\right),\qquad
t_{z}^{\top}\left(g^{\star}_{f_{\mathrm{enc}}}\right)^{\prime}(z)=0\,,$ (16)

where $(g^{\star}_{f_{\mathrm{enc}}})^{\prime}(z)$ is the derivative of the
decoder. By passing to the limit $\beta\to+\infty$ in the latter condition
(assuming that the derivatives in $z$ and the limit $\beta\to+\infty$ can be
exchanged), it follows that

$\forall z\in(z_{A},z_{B}),\quad\forall t_{z}\in
T\Sigma_{z}\left(g^{\infty}_{f_{\mathrm{enc}}}(z)\right),\qquad
t_{z}^{\top}\left(g^{\infty}_{f_{\mathrm{enc}}}\right)^{\prime}(z)=0.$

The comparison between the latter condition and (15) shows that $\nabla
V(g^{\infty}_{f_{\mathrm{enc}}}(z))$ is collinear to
$(g^{\infty}_{f_{\mathrm{enc}}})^{\prime}(z)$ (provided $\nabla
V(g^{\infty}_{f_{\mathrm{enc}}}(z))\neq 0$), which shows that the curve
$\\{g^{\infty}_{f_{\mathrm{enc}}}(z)\\}_{s\in[z_{A},z_{B}]}$ is a MEP.

###  3.2 Conditional expectations for other unsupervised dimensionality
reduction methods

We have seen in Section 3.1 how the unsupervised least-square reconstruction
error can be reformulated in terms of conditional expectations as follows (see
(5) and (6)):

$\inf_{f_{\mathrm{enc}}\in\mathcal{F}_{\mathrm{enc}}}\mathbb{E}\left[\Big{\|}X-\mathbb{E}\left(X\,\middle|\,f_{\mathrm{enc}}(X)\right)\Big{\|}^{2}\right],$
(17)

where the random variable $X$ is distributed according to $\mu$. We discuss in
this section how two other unsupervised learning methods, $K$-means clustering
and PCA, can be understood in terms of a similar minimization problem
involving conditional expectations. Such interpretations have already been
discussed in the literature (see for instance [32, Section 4]), but are of
interest for discussions in later sections, in particular Section 5 where we
combine autoencoders with ideas from clustering.

####  $K$-means clustering.

For simplicity of exposition and to be closer to traditional expositions in
machine learning textbooks, we discuss here a reformulation of clustering
methods for sample data, and not at the level of probabililty distributions
(corresponding to so-called population losses, as we have considered until
now). We therefore introduce a sample data set
$\mathcal{D}_{N_{\mathrm{data}}}=\\{x^{1},\dots,x^{N_{\mathrm{data}}}\\}\subset\mathcal{X}$
of cardinality $N_{\mathrm{data}}$, and denote by
$\widehat{\mu}_{N_{\mathrm{data}}}$ the empirical probability measure

$\widehat{\mu}_{N_{\mathrm{data}}}=\frac{1}{N_{\mathrm{data}}}\sum_{n=1}^{N_{\mathrm{data}}}\delta_{x^{n}}.$

The data is partitioned into $K\leqslant N_{\mathrm{data}}$ classes indexed by
an integer $1\leqslant k\leqslant K$, using an assignment function
$f_{\mathrm{ass}}:\mathcal{D}_{N_{\mathrm{data}}}\to\\{1,\dots,K\\}$. For a
given assignment function, we denote by

$\mathscr{C}_{k}=\left\\{x\in\mathcal{D}_{N_{\mathrm{data}}}\,\middle|\,f_{\mathrm{ass}}(x)=k\right\\}=f_{\mathrm{ass}}^{-1}\\{k\\}$

the various clusters identified by the clustering method. The cardinality of
$\mathscr{C}_{k}$ is denoted by $|\mathscr{C}_{k}|$. With this notation, the
objective function of the $K$-means algorithm reads (see for instance [71,
Section 21.3] or [83, Section 22.2])

$\sum_{k=1}^{K}\sum_{i\in\mathscr{C}_{k}}\left\|x^{i}-\frac{1}{|\mathscr{C}_{k}|}\sum_{j\in\mathscr{C}_{k}}x^{j}\right\|^{2}.$
(18)

It is easy to see that (18) is indeed of the form (17) with expectations taken
with respect to $\widehat{\mu}_{N_{\mathrm{data}}}$ and $f_{\mathrm{enc}}$
replaced by $f_{\mathrm{ass}}$ since

$\forall
k\in\\{1,\dots,K\\},\qquad\mathbb{E}_{\widehat{\mu}_{N_{\mathrm{data}}}}(X\,|\,f_{\mathrm{ass}}(X)=k)=\frac{1}{|\mathscr{C}_{k}|}\sum_{j\in\mathscr{C}_{k}}x^{j}.$

The latter two equalities show that autoencoders whose encoders have values in
the discrete space $\\{1,\dots,K\\}$ aim at solving the same minimization
problem as $K$-means, the assignment function being the encoder function.

#### Principal component analysis.

We denote by $K\in\\{1,\dots,D-1\\}$ the number of dimensions kept in the
dimensionality reduction process performed by PCA. For a given a sample data
set $\mathcal{D}_{N_{\mathrm{data}}}$, we denote by

$\overline{x}_{N_{\mathrm{data}}}=\frac{1}{N_{\mathrm{data}}}\sum_{n=1}^{N_{\mathrm{data}}}x^{n}\in\mathbb{R}^{D}$

the empirical mean of the data, and by

$\widehat{\Sigma}_{N_{\mathrm{data}}}=\frac{1}{N_{\mathrm{data}}}\sum_{n=1}^{N_{\mathrm{data}}}\left(x^{n}-\overline{x}_{N_{\mathrm{data}}}\right)\left(x^{n}-\overline{x}_{N_{\mathrm{data}}}\right)^{\top}\in\mathbb{R}^{D\times
D}$

the empirical covariance matrix. The PCA algorithm performs a dimensionality
reduction by computing normalized eigenvectors
$\widehat{u}_{1},\cdots,\widehat{u}_{K}$ associated with the $K$ largest
eigenvalues of the empirical covariance matrix, and proposing the following
reconstruction of a test point $x\in\mathbb{R}^{D}$:

$\widetilde{x}=\overline{x}_{N_{\mathrm{data}}}+\sum_{k=1}^{K}\left[\widehat{u}_{k}^{\top}\left(x-\overline{x}_{N_{\mathrm{data}}}\right)\right]\widehat{u}_{k}=\overline{x}_{N_{\mathrm{data}}}+\widehat{U}\widehat{U}^{\top}\left(x-\overline{x}_{N_{\mathrm{data}}}\right),$
(19)

where the matrix $\widehat{U}\in\mathbb{R}^{D\times K}$ has columns
$\widehat{u}_{1},\cdots,\widehat{u}_{K}$. The maximal variance formulation of
PCA is equivalent to a minimum reconstruction error formulation. In fact, as
made precise in [83, Section 23.1] for instance, the above matrix
$\widehat{U}$ is a solution to the minimization problem

$\min_{\begin{subarray}{c}W\in\mathbb{R}^{D\times K}\\\
W^{\top}W=\mathrm{I}_{K}\end{subarray}}\mathbb{E}_{\widehat{\mu}_{N_{\mathrm{data}}}}\left(\left\|X-WW^{\top}X\right\|^{2}\right).$
(20)

This formulation is the one relevant to write PCA in the autoencoder framework
[10, 5].

PCA can be reformulated as the minimization task (17) in situations when the
affine transformation in (19) can be interpreted in terms of some
conditioning. The only situation we are aware of is when a population loss is
considered, with $X$ a Gaussian random variable. Without loss of generality,
to simplify the notation, the Gaussian distribution can be assumed to be
centered (upon subtracting the mean of the distribution to the random variable
under consideration). To state the result, we introduce the set
$\mathcal{F}_{\mathrm{LOP}}$ of linear orthogonal projections:

$\mathcal{F}_{\mathrm{LOP}}=\left\\{f_{\mathrm{LOP}}:\left\\{\begin{aligned}
\mathbb{R}^{D}&\rightarrow\mathbb{R}^{K}\\\ x&\mapsto
W^{\top}x\end{aligned}\right.\ \mbox{for some}~{}W\in\mathbb{R}^{D\times
K}\mbox{ such that}~{}W^{\top}W=\mathrm{I}_{K}\right\\}.$ (21)

######  Lemma 1.

Assume that $X$ is a Gaussian random variable with mean 0. Then, the PCA
minimization problem (20) is equivalent to

$\min_{f_{\mathrm{LOP}}\in\mathcal{F}_{\mathrm{LOP}}}\mathbb{E}\left[\Big{\|}X-\mathbb{E}\left(X\,\middle|\,f_{\mathrm{LOP}}(X)\right)\Big{\|}^{2}\right].$
(22)

###### Proof.

For a given a matrix $W\in\mathbb{R}^{D\times K}$ such that
$W^{\top}W=\mathrm{I}_{K}$, let us compute the conditional expectation
$\mathbb{E}(X\,|\,f_{\mathrm{LOP}}(X))$ for the linear map
$f_{\mathrm{LOP}}(x)=W^{\top}x$. Introduce to this end a matrix
$U\in\mathbb{R}^{D\times(D-K)}$ such that $U^{\top}W=0$ and
$U^{\top}U=\mathrm{I}_{D-K}$. Since the columns of $W$ and $U$ form an
orthonormal basis of $\mathbb{R}^{D}$, the identities
$(WW^{\top}+UU^{\top})W=W$ and $(WW^{\top}+UU^{\top})U=U$ imply that
$WW^{\top}+UU^{\top}=\mathrm{I}_{D}$.

Given $z\in\mathbb{R}^{K}$, we can then write (with $\mu$ the distribution of
the centered Gaussian measure under consideration)

$\displaystyle\mathbb{E}\left(X\,\middle|\,f_{\mathrm{LOP}}(X)=z\right)$
$\displaystyle=\mathbb{E}\left(X\,\middle|\,W^{\top}X=z\right)$
$\displaystyle=Z_{z}^{-1}\int_{\\{x\in\mathbb{R}^{D}:W^{\top}x=z\\}}x\,\mu(dx)$
$\displaystyle=Z_{z}^{-1}\int_{\\{x\in\mathbb{R}^{D}:W^{\top}x=z\\}}\left(WW^{\top}x+UU^{\top}x\right)\mu(dx)$
$\displaystyle=Wz+Z_{z}^{-1}U\int_{\\{x\in\mathbb{R}^{D}:W^{\top}x=z\\}}U^{\top}x\,\mu(dx),$
(23)

where

$Z_{z}=\int_{\\{x\in\mathbb{R}^{D}:W^{\top}x=z\\}}\mu(dx).$

Note that, since $X$ is a centered Gaussian random variable and $U^{\top}W=0$,
the two vectors $U^{\top}X$ and $W^{\top}X$ are independent centered Gaussian
random variables in $\mathbb{R}^{D-K}$ and $\mathbb{R}^{K}$, respectively.
Therefore, the second term in (23) vanishes, and

$\mathbb{E}\left(X\,\middle|\,f_{\mathrm{LOP}}(X)\right)=WW^{\top}X.$

This shows that the two minimization tasks (20) and (22) are equivalent. ∎

###  3.3 Numerical illustration

In this section, we illustrate the results obtained in Section 3.1 by
considering a concrete example in $\mathbb{R}^{2}$. In particular, the formula
(6) suggests to consider decoders which are sufficiently expressive, in order
to correctly approximate the conditional expectation. This motivates working
with asymmetric autoencoders, where decoders have a larger complexity than
encoders. In contrast, it is important to keep relatively simple encoders
(computationally not too expensive) in molecular dynamics since they are
intended to be used as CVs, for instance to perform free energy computations.

#### Architecture and training of autoencoders.

In the following, the architectures of autoencoders are characterized using
the formulation of Figure 1. We systematically use hyperbolic tangents as
activation functions for all layers, except the one leading to the bottleneck
and the output, for which linear activation functions are considered. As
discussed in Section 2.1, autoencoders with complex architectures can overfit
when the dataset is not large enough, which is why various regularization
strategies are considered. In our situation, as the models we use typically
contain a moderate number of parameters, it is sufficient to rely on early
stopping for regularization. Concretely, in the present case, the training was
stopped once the validation loss did not decrease after a certain number of
epochs $n_{\mathrm{wait}}$. The model giving the lowest validation loss was
kept. Obviously, increasing $n_{\mathrm{wait}}$ leads to potential
overfitting, and the appropriate value of this parameter, which depends on the
chosen architecture, has to be chosen with some form of cross–validation.

Figure 1: Schematic representation of a symmetric autoencoder. Blue neurons
correspond to hidden layers, while the input and output layers are
respectively in green and red. The encoder and decoder are respectively
referred to as (5, 4, 3, 2) and (2, 3, 4, 5). Activation functions are
hyperbolic tangents, except for the bottleneck and output layers, for which
linear activation functions are considered in order not to restrict the range
of values.

#### Measure of quality for conditional expectations.

For simplicity, we choose here the latent space $\mathcal{Z}$ to be
one–dimensional, _i.e._ $d=1$; but our analysis can be extended to more
general settings. Since $\nabla
f_{\mathrm{enc}}(g^{\star}_{f_{\mathrm{enc}}}(z))$ is orthogonal to
$\Sigma_{z}$ at $g^{\star}_{f_{\mathrm{enc}}}(z)\in\Sigma_{z}$ (as
$\Sigma_{z}$ is a level set of $f_{\mathrm{enc}}$), and recalling (16), we
conclude that the cosine of the angle between vectors
$(g^{\star}_{f_{\mathrm{enc}}})^{\prime}(z)$ and $\nabla
f_{\mathrm{enc}}(g^{\star}_{f_{\mathrm{enc}}}(z))$, namely

$\rho(z)=\frac{\nabla
f_{\mathrm{enc}}(g^{\star}_{f_{\mathrm{enc}}}(z))^{\top}(g^{\star}_{f_{\mathrm{enc}}})^{\prime}(z)}{\left\|\nabla
f_{\mathrm{enc}}(g^{\star}_{f_{\mathrm{enc}}}(z))\right\|\left\|(g^{\star}_{f_{\mathrm{enc}}})^{\prime}(z)\right\|}$
(24)

should ideally be $1$. We monitor this quantity in the numerical tests below.

#### Müller–Brown potential.

We consider the Müller-Brown potential [70] defined for
$x=(x_{1},x_{2})\in\mathbb{R}^{2}$ as

$V(x_{1},x_{2})=\sum_{i=1}^{4}A_{i}\exp\left(a_{i}\left(x_{1}-u_{i}\right)^{2}+b_{i}\left(x_{1}-u_{i}\right)\left(x_{2}-v_{i}\right)+c_{i}\left(x_{2}-v_{i}\right)^{2}\right),$
(25)

with the parameters $A=(-200,-100,-170,15)$, $a=(-1,-1,-6.5,0.7)$,
$b=(0,0,11,0.6)$, $c=(-10,-10,-6.5,0.7)$, $u=(1,0,-0.5,-1)$ and
$v=(0,0.5,1.5,1)$. The dataset is obtained by sampling the Boltzmann–Gibbs
measure $\mu$ associated with $V$ for $\beta=0.05$ (_i.e._ the density of
$\mu$ is proportional to $\mathrm{e}^{-\beta V}$), using overdamped Langevin
dynamics discretized with a Euler–Maruyama scheme and a time step $\Delta
t=10^{-4}$. A dataset of $2\times 10^{7}$ points covering all the relevant
parts of the potential landscape was generated by running two trajectories of
$10^{7}$ points starting from two minima of the potential, located
respectively at $(-0.56,1.44)$ and $(0.62,0.03)$.

The training and validation datasets were built by first drawing randomly
without replacement $N_{\mathrm{data}}=10^{4}$ points for each of them from
the full dataset. The hyperparameters for the training are chosen using a
$K$-fold cross–validation procedure. The test set, composed of the remaining
points in the initial dataset, was used to assess the quality of the results.
All the models were trained using the Adam algorithm [54] with a learning rate
of $\eta=0.005$. The remaining training hyperparameters to select are the
model architecture, the minibatch size and the early stopping criterion
$n_{\mathrm{wait}}$. We first use a model with a given architecture (namely
(2, 5, 5, 1, 20, 20, 2)) to select the batch size and $n_{\mathrm{wait}}$, and
then study the impact of the model architecture in a second stage. For all
these trainings, the same random seed was used to initialize the parameters of
the model. We give typical orders of magnitudes for the number of epochs in
the training in the last column of Table 1.

Table 1: Test reconstruction error (3) for the training of the AE model (2, 5, 5, 1, 20, 20, 2) with fixed learning rate $\eta=0.005$, for various values of $n_{\mathrm{wait}}$ and batch sizes. $n_{\mathrm{wait}}$ | batch size |  relative test reconstruction error |  number of epochs (order of magnitude)
---|---|---|---
25 | $1\times 10^{3}$ | $1.98\times 10^{-2}$ | 500
50 | $1\times 10^{3}$ | $1.93\times 10^{-2}$ | 700
100 | $1\times 10^{3}$ | $1.86\times 10^{-2}$ | 850
25 | $5\times 10^{2}$ | $1.90\times 10^{-2}$ | 300
50 | $5\times 10^{2}$ | $1.83\times 10^{-2}$ | 400
100 | $5\times 10^{2}$ | $1.82\times 10^{-2}$ | 700

Table 1 shows that, on the range of hyperparameter values under consideration,
the test reconstruction error (3) slightly decreases as the batch size
increases and $n_{\mathrm{wait}}$ increases. The variation is however not very
significant as the relative test reconstruction error (test reconstruction
error divided by the empirical variance of the test data) is quite small for
all the values of the parameters that are considered. The choice of batch size
and early stopping criterion are not very critical as the quality of the
results is rather robust to changes in these parameters.

In view of (9), the fraction of unexplained variance, which corresponds to the
reconstruction error divided by the total variance, is (for a perfectly
converged decoder)

$\frac{\mathbb{E}\left[\mathrm{Var}(X|f_{\mathrm{enc}}(X))\right]}{\mathrm{Var}(X)}=1-\frac{\mathrm{Var}\left[\mathbb{E}(X|f_{\mathrm{enc}}(X))\right]}{\mathrm{Var}(X)}.$

This quantity can be estimated by replacing expectations with empirical
averages on the test set. For the settings considered in Table 1, the
empirical variance of the test set is $0.86$, and the model fails to explain
less than 2% of the variance of the test set. A more careful determination of
the batch size or $n_{\mathrm{wait}}$ may even further decrease this
percentage, but such levels of error are sufficiently low to illustrate the
theoretical points discussed in the previous sections. In the following the
training of AE models are done with a batchsize of $500$ and
$n_{\mathrm{wait}}=100$.

We now illustrate the impact of the decoder’s architecture on its capacity to
approximate conditional expectations (see (6)). We trained two AEs, with
respective architectures (2, 5, 5, 1, 5, 5, 2) and (2, 5, 5, 1, 20, 20, 2),
using the same training and validation sets. To compute the conditional
averages, the range of values of the encoder in the one dimensional latent
space $\mathcal{Z}$ was split into $100$ non-overlapping intervals (bins) of
constant lengths. Conditional averages are then approximated as the mean
values of the data points whose encoded values are in the corresponding bins.
Conditional averages and values of the decoder are plotted on the potential
energy heatmap in Figures 2(a) and 2(b) for both architectures, together with
the MEP and isolevels of the encoder $f_{\mathrm{enc}}$. These pictures
qualitatively demonstrate that an increase in the complexity of the decoder
allows it to better approximate conditional averages. Note also that the
decoder path seems to be orthogonal to the isolevels of the encoder, which is
in agreement with the discussion above (24). The quality of the agreement can
be made more quantitative by computing the Euclidean distance between the
values of the decoder and the conditional expectations for the 100 bins
considered in $\mathcal{Z}$ space, and the alignment criterion (24) between
the gradient of the encoder and the derivative of the decoder; see Figures 3
and 4. These pictures show respectively that the Euclidean distance is small,
and that the alignment factor is close to 1. Note also that more expressive
decoders lead to better values of these criteria.

(a)

(b)

Figure 2: Minimum energy path, conditional averages and decoder on the
potential energy heatmap with isolevels of the encoder $f_{\mathrm{enc}}$ for
the AE models (a) (2, 5, 5, 1, 5, 5, 2) and (b) (2, 5, 5, 1, 20, 20, 2).
Conditional averages are computed using 100 uniformly spaced bins in encoder
space.

(a)

(b)

Figure 3: Euclidean distance between the conditional averages and the values
of the decoder, where conditional averages are computed using 100 uniformly
spaced bins in encoder space $\mathcal{Z}$ for the AE models (a) (2, 5, 5, 1,
5, 5, 2) and (b) (2, 5, 5, 1, 20, 20, 2).

(a)

(b)

Figure 4: Alignment criterion (24) evaluated at the values of the conditional
averages computed using 100 uniformly spaced bins in encoder space
$\mathcal{Z}$ for the AE models (a) (2, 5, 5, 1, 5, 5, 2) and (b) (2, 5, 5, 1,
20, 20, 2).

Note moreover from Figure 2 that the decoder path with the more complex
architecture is closer to the MEP, although there are some discrepancies. The
agreement would presumably be better for larger values of $\beta$. Although
more complex decoders allow to better approximate conditional expectations,
let us mention that the associated test losses may not be smaller than those
for less complex decoders. The value of the test loss alone may therefore be
misleading to assess the convergence of the decoder to the conditional
average, as a trained model with a slightly smaller test loss might lead to an
encoder/decoder pair less close to the conditional average (in the sense of
equation (24)). Here, the relative test reconstruction error for the trained
model with the smaller decoder was $1.71\times 10^{-2}$, which is comparable
to the relative test reconstruction error $1.83\times 10^{-2}$ obtained with
the larger decoder (the fact that this reconstruction error can be larger may
seem surprising, but can be attributed to some variability in the training
procedure and choice of data sets). On the other hand, the decoder path
qualitatively looks nicer for the more complex model, and this can be
quantified using our measure of alignment and/or the distance between the
decoder path and conditional averages, as illustrated in Figures 3 and 4.

##  4 Better describing the transition

As discussed at the end of Section 3.1, the decoder parametrizes the MEP in
the small temperature regime, and thus gives information on the transition
between local minima. We discuss in this section ways to improve the ability
of autoencoders to describe such transitions.

The loss function of autoencoders focuses the reconstruction effort on more
likely regions under the probability measure $\mu$ of the data points. For
physical systems described by the Boltzmann–Gibbs measure, this means that the
reconstruction should be accurate in metastable states around local minima of
the potential energy function. In contrast, transition states between
metastable regions are typically only scarcely sampled, and hence a large
reconstruction error on such configurations may not have a large impact on the
overall reconstruction error. However, from a physical viewpoint, transition
states are very important to understand how the system can undergo a
transition from one metastable state to another.

We explore in this section one idea to better take into account transition
states, by changing the reference measure from the Boltzmann–Gibbs measure to
a measure putting more mass on regions between metastable states, such as the
reactive trajectory measure [27, 66] (which is the distribution of
configurations sampled by portions of trajectories switching from one
metastable state to another). We explore in Section 6 another idea, namely the
use of extra physical information encoded via additional terms in the loss
functions (as considered in [79] for instance).

We illustrate the impact of changing the reference probability measure $\mu$
of the data distribution on the quality of the autoencoder training on the
two-dimensional Müller–Brown potential considered in Section 3.3. When
training an autoencoder on this model with any kind of architecture, there is
always a chance that the training converges to a “wrong” solution, depending
on the random initialization of the weights and the subsequent randomness in
the training procedure. Such a situation is illustrated in Figure 6(a).
Although this “wrong” solution allows to capture the variance, the decoder
path between the probability mode on the left and the one on the right is far
away from the MEP. Somehow, the mode on the left is described in the wrong
direction. We use such spurious solutions to study the impact of the change of
the reference probability measure by observing whether the wrong solution
subsists after training.

We consider a family of probability measures obtained by a convex combination
of the usual Boltzmann–Gibbs measure and the reactive trajectory distribution.
The latter one is sampled by running adaptive multilevel splitting (AMS) [12,
65, 77] to sample the transition paths for the overdamped Langevin dynamics
described in Section 3.3. The initial state $A$ and final state $B$ are
defined as small discs of radius $0.1$ centered on the two local minima of the
potential, and the one dimensional reaction coordinate used for AMS is the
finite element approximation of the committor function (as done in [28] for
instance). A dataset of $\widetilde{N}_{\mathrm{data}}=2\times 10^{5}$
configurations $\\{\widetilde{x}^{n}\\}_{1\leqslant
n\leqslant\widetilde{N}_{\mathrm{data}}}$ distributed according to the
reactive trajectory measure is obtained by drawing at random without
replacement configurations obtained by running two AMS simulations: one
forward and one backward, with 1000 replicas each. To compare visually this
distribution with the Boltzmann–Gibbs distribution, two sets of $10^{4}$
points distributed according to these two distributions are plotted in Figure
5. The probability distribution $\mu$ is obtained by considering a fraction
$\lambda\in[0,1]$ of the Boltzmann–Gibbs measure, and a fraction $1-\lambda$
of the reaction path measure. Denoting by $N_{\mathrm{data}}$ the number of
data points $\\{x^{n}\\}_{1\leqslant n\leqslant N_{\mathrm{data}}}$
distributed according to the Boltzmann–Gibbs measure (the same set as in
Section 3.3), the associated empirical training loss is then (see (2))

$\widehat{\mathscr{L}}_{\lambda}(\theta)=\frac{\lambda}{N_{\mathrm{data}}}\sum_{n=1}^{N_{\mathrm{data}}}\left\|x^{n}-f_{\theta}(x^{n})\right\|^{2}+\frac{1-\lambda}{\widetilde{N}_{\mathrm{data}}}\sum_{n=1}^{\widetilde{N}_{\mathrm{data}}}\left\|\widetilde{x}^{n}-f_{\theta}(\widetilde{x}^{n})\right\|^{2}.$

(a)

(b)

Figure 5: Minimum energy path (purple line) and configurations (orange dots)
distributed according to (a) the Boltzmann–Gibbs distribution, and (b) the
reactive trajectory distribution.

The differences between the results obtained with the trained models in
Figures 2(b) and 6(a) come from the values of the initial weights of the model
and the changes in the chosen minibatches during the training. Here, the
initial weights of the models are chosen at random, using the default
initialization of PyTorch (see [40, Section 8.4] for further background). This
initialization and the choices made in the minibatching procedures are
determined by the seed of the random number generators, which is set to
different values for Figures 2(b) and 6(a).

(a)

(b)

Figure 6: Minimum energy path, conditional averages and decoder on the
potential energy heatmap with isolevels of the encoder $f_{\mathrm{enc}}$ for
the AE models (2, 5, 5, 1, 20, 20, 2) trained using (a) the Boltzmann–Gibbs
distribution ($\lambda=1$) and (b) the reactive trajectory distribution
($\lambda=0)$ as reference measures. Conditional averages are computed using
(a) the Boltzmann-Gibbs distribution or (b) the reactive trajectory
distribution with 100 uniformly spaced bins in encoder space $\mathcal{Z}$.

(a)

(b)

Figure 7: Minimum energy path, conditional averages and decoder on the
potential energy heatmap with isolevels of the encoder $f_{\mathrm{enc}}$ for
the AE models (2, 5, 5, 1, 20, 20, 2) trained with (a) $\lambda=0.5$ and (b)
$\lambda=\frac{1}{11}$. Conditional averages are consistently computed using
$\lambda$-mixture of both measures with 100 uniformly spaced bins in encoder
space $\mathcal{Z}$.

Depending on the relative weight of the two reconstruction errors in the
definition of $\widehat{\mathscr{L}}_{\lambda}(\theta)$, the “wrong solution”
can be avoided since configurations distributed according to the reactive
measure emphasize reconstruction efforts in the transition region. In Figures
6 and 7, the training results starting from the same initial conditions are
shown for various values of $\lambda$. For $\lambda=0$ (training on the
Boltzmann–Gibbs distribution only), the “wrong solution” is obtained. For
$\lambda=1$ (only configurations from the reactive path measure), transition
regions are well described but the two important modes in the probability
measure are not well taken into account. Additional numerical tests (not
reported here) indicate that the best results are obtained for $\lambda$ close
to 1, as this is a situation where the addition of more points in the
transition regions allows to avoid converging to the “wrong solution”, while
still focusing reconstruction efforts to the important regions of the
configuration space. The value of the parameter $\lambda$ could be fine-tuned
(using a $K$-fold cross-validation procedure for instance) to obtain the most
satisfactory results in terms of the values of the alignment criterion (24)
and distance between decoder values and conditional averages.

##  5 Parametrizing multiple transition paths in two state systems

In a situation where multiple transition paths link two metastable states, the
autoencoder may fail to properly represent the system in the transition region
between local minima for a one-dimensional latent space $\mathcal{Z}$, as it
constructs only a single curve for the conditional expectations. An idea to
address this issue is to consider multiple decoders associated with a common
encoder, and to choose for a given configuration the decoder which best
reconstructs the state through some assignment function reminiscent of the one
considered for clustering (recall Section 3.2). We illustrate such a strategy
in this section. Related strategies were suggested for principal curves (see
[90, Section II.B] and [32, Section 3.5] for instance), and actually tested in
[55].

#### Loss function for multiple decoders.

In order to take into account the presence of multiple decoders, the
reconstruction loss in (1) for a given encoder function $f_{\mathrm{enc}}$ and
$K$ decoder functions $f_{{\rm dec},k}$ is modified as

$\mathbb{E}\left[\underset{k\in 1,..,K}{\min}\big{\|}x-f_{{\rm dec},k}\circ
f_{\mathrm{enc}}(x)\big{\|}^{2}\right].$

The associated training loss for a given dataset $\\{x^{n}\\}_{1\leqslant
n\leqslant N_{\mathrm{data}}}$ reads

$\widehat{\mathscr{L}}_{K}(\theta)=\frac{1}{N_{\mathrm{data}}}\sum_{n=1}^{N_{\mathrm{data}}}\underset{k\in
1,..,K}{\min}\left\|x^{n}-f_{{\rm dec},\theta_{k+1}}\circ
f_{\mathrm{enc},\theta_{1}}(x^{n})\right\|^{2}.$ (26)

######  Remark 2.

Although we did not explore this option here, it would be possible to replace
the hard assignment encoded by taking the minimum over $k$ by a probabilistic
assignment arising from a softmin function, similarly to what is considered in
the probabilistic approach to classification problems.

#### Numerical illustration.

To illustrate the approach described above, we first consider a double well
potential on a circle. More precisely, we consider the potential energy
function in $\mathbb{R}^{2}$

$\displaystyle
V(x_{1},x_{2})=2x_{2}^{2}+10\left(x_{1}^{2}+x_{2}^{2}-1\right)^{2}.$ (27)

A dataset of $2\times 10^{5}$ points distributed according to the
Boltzmann–Gibbs measure was generated using a Euler–Maruyama discretization of
the overdamped Langevin dynamics with $\beta=2$ and $\Delta t=0.01$. The early
stopping waiting time $n_{\mathrm{wait}}=50$, learning rate $0.005$ and
minibatch size $500$ were kept constant for all the subsequent training
experiments with (26). We compare in Figure 8 results obtained with $K=2$
decoders to results with a single decoder. When a single decoder is
considered, the decoder path describes accurately both modes, but only one
transition path is correctly indexed, while on the second path, in the zone of
lowest weight according to the Boltzmann–Gibbs measure, the description is
less precise and the encoder function varies very fast. When $K=2$ decoders
are considered, each decoder takes care of a mode of the probability measure,
with some interface more or less in the middle. The encoder has rather
horizontal isolevels in the modes, constructed so that the decoder path is
close to some principal curve. Of course, this does not allow to infer
anything about the transition mechanism for configurations switching from one
probability mode to another.

(a)

(b)

Figure 8: (a): Minimum energy path, conditional averages and decoders on the
potential energy heatmap with isolevels of the encoder $f_{\mathrm{enc}}$ for
models with a single decoder ($K=1$). (b): Same plots for a model with $K=2$
decoders on top of which the training dataset points are colored according to
the index of the decoder which leads to the minimal reconstruction error.
Conditional averages are computed using 100 bins uniformly spaced in the
encoded dimension.

Further physical insights need to be injected in the model in order to obtain
a reconstruction of the system which is also relevant to describe the
transition paths switching from one metastable region to another. One possible
approach to this end would be to train the model on configurations sampled
according to the reactive path measure, already considered in Section 4, as
each decoder would concentrate on one mode of this probability measure
(located around the transition points in the vicinity of $(0,-1)$ and
$(0,1)$).

An alternative approach, which we consider in the remainder of this section,
is to add extra terms to the loss (26), so that the encoder better captures
the transition paths. A first option is

$\displaystyle\widehat{\mathscr{L}}_{K,\lambda_{0},\lambda_{1},\lambda_{2}}(\theta)=\lambda_{0}\widehat{\mathscr{L}}_{K}(\theta)+\lambda_{1}\left(\widehat{\mathrm{Var}}_{N_{\mathrm{data}}}(f_{\rm{enc},\theta_{1}})-1\right)^{2}+\frac{\lambda_{2}}{N_{\mathrm{data}}}\sum_{n=1}^{N_{\mathrm{data}}}\left|\nabla
f_{\rm{enc},\theta_{1}}(x^{n})\right|^{2},$ (28)

where $\lambda_{0},\lambda_{1},\lambda_{2}\geqslant 0$ allow to tune the extra
loss terms on the encoder and

$\widehat{\mathrm{Var}}_{N_{\mathrm{data}}}(f)=\frac{1}{N_{\mathrm{data}}}\sum_{n=1}^{N_{\mathrm{data}}}|f(x^{n})|^{2}-\left|\frac{1}{N_{\mathrm{data}}}\sum_{n=1}^{N_{\mathrm{data}}}f(x^{n})\right|^{2}$
(29)

denotes the empirical variance of a function $f$. Note that the extra terms in
the loss function are insensitive to shifts in the values of the encoder
function.

In the particular case $\lambda_{0}=0$, one recovers in the limit
$N_{\mathrm{data}}\to+\infty$ and $\lambda_{1}\to+\infty$ the loss function
characterizing eigenfunctions associated with the first non-zero eigenvalue of
the generator of the overdamped Langevin dynamics. Indeed, the optimization
problem converges in the limit $N_{\mathrm{data}}\to+\infty$ to the
minimization of

$\lambda_{1}\left(\int_{\mathcal{X}}\left(f_{\mathrm{enc}}(x)-\int_{\mathcal{X}}f_{\mathrm{enc}}\,d\mu\right)^{2}\,\mu(dx)-1\right)^{2}+\lambda_{2}\int_{\mathcal{X}}\left|\nabla
f_{\mathrm{enc}}(x)\right|^{2}\mu(dx).$

When $\lambda_{1}\to+\infty$, the first term in the expression above enforces
a normalization constraint on $f_{\mathrm{enc}}$. The minimization of
$\widehat{\mathscr{L}}_{K,0,\lambda_{1},\lambda_{2}}$ in the limits
$N_{\mathrm{data}}\to+\infty$ and $\lambda_{1}\to+\infty$ can be then
reformulated as

$\inf\left\\{\int_{\mathcal{X}}\left|\nabla
f_{\mathrm{enc}}(x)\right|^{2}\mu(dx)\ \middle|\
\int_{\mathcal{X}}\left|f_{\mathrm{enc}}(x)-\int_{\mathcal{X}}f_{\mathrm{enc}}\,d\mu\right|^{2}\mu(dx)=1\right\\},$

which characterizes (up to constant shifts) eigenfunctions associated with the
first non-zero eigenvalue of the generator of the overdamped Langevin
dynamics; see [98] for a more detailed derivation and a numerical strategy
based on neural networks.

(a)

(b)

Figure 9: Minimum energy path, conditional averages and decoders on the
potential energy heatmap with isolevels of the encoder $f_{\mathrm{enc}}$ for
models trained to minimize the loss (28) (a) $\lambda_{0}=1$,
$\lambda_{1}=10^{2}$ and $\lambda_{2}=10^{-2}$ and (b) $\lambda_{0}=0$,
$\lambda_{1}=10^{2}$ and $\lambda_{2}=10^{-2}$. These results should be
compared to the one obtained in Figure 8(b), which corresponds to the choice
$\lambda_{0}=1$ and $\lambda_{1}=\lambda_{2}=0$.

Including the regularization terms to the multiple decoder loss as in (28)
allows to obtain a model for which each decoder correctly indexes a different
transition paths (see Figure 9(a)). In fact using only the added penalization
terms ($\lambda_{0}=0$) allows to obtain a collective variable able to index
both paths (see Figure 9(b)). The reconstruction error with the two decoders
allows to perform a clustering of the dataset corresponding to the two
transition paths.

In addition to the extra terms in (28), one can also fix the values of the
encoder at specific points, such as local mimima or centers of metastable
modes. Here for instance, one can add penalization terms
$f_{\mathrm{enc}}(-1,0)^{2}$ and $[f_{\mathrm{enc}}(1,0)-1]^{2}$ to force the
encoder to have values close to 0 in the left mode and close to 1 in the right
mode, respectively. This requires the encoder to describe the transition from
values 0 to 1. In essence, this can be seen as some form of semi-supervised
approach where specific points are identified as belonging to the left or
right mode, and the values of the encoder are softly enforced there. Of
course, this approach relies on some prior information on the metastable
states.

In view of the above discussion, the loss finally reads

$\displaystyle\widehat{\mathscr{L}}_{K,\lambda_{0},\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4}}(\theta)=\lambda_{0}\widehat{\mathscr{L}}_{K}(\theta)$
$\displaystyle+\lambda_{1}\left(\widehat{\mathrm{Var}}_{N_{\mathrm{data}}}(f_{\rm{enc},\theta_{1}})-1\right)^{2}$
$\displaystyle+\lambda_{2}\frac{1}{N_{\mathrm{data}}}\sum_{n=1}^{N_{\mathrm{data}}}\left|\nabla
f_{\rm{enc},\theta_{1}}(x^{n})\right|^{2}$
$\displaystyle+\lambda_{3}\frac{1}{N_{\mathrm{pen}}}\sum_{n=1}^{N_{\mathrm{pen}}}\sum_{k=1}^{K}\big{(}\widetilde{x}^{n}-f_{{\rm
dec},\theta_{k+1}}\circ
f_{\mathrm{enc},\theta_{1}}(\widetilde{x}^{n})\big{)}^{2}$
$\displaystyle+\lambda_{4}\frac{1}{N_{\mathrm{pen}}}\sum_{n=1}^{N_{\mathrm{pen}}}\left(\widetilde{y}^{n}-f_{\rm{enc},\theta_{1}}(\widetilde{x}^{n})\right)^{2},$
(30)

where $\lambda_{0}$, $\lambda_{1}$, $\lambda_{2}$, $\lambda_{3}$ and
$\lambda_{4}$ are nonnegative numbers which allow to tune the strength of the
constraints, and $\\{\widetilde{x}^{n}\\}_{1\leqslant n\leqslant N_{\rm pen}}$
is the set of points at which the value of the encoder is softly constrained.
These points are typically not points of the training data set, but extra
points where the value of the encoder is a priori known. The term in factor of
$\lambda_{3}$ ensures that the reconstruction is correct on the set
$\\{\widetilde{x}^{n}\\}_{1\leqslant n\leqslant N_{\rm pen}}$ for all $K$
decoders, while the term in factor of $\lambda_{4}$ forces the encoder values
to be close to $\widetilde{y}^{n}$ for the configuration $\widetilde{x}^{n}$.
The chosen values of $\widetilde{y}^{n}$ depend on the probability mode to
which the configuration belongs. When there are two modes,
$\widetilde{y}^{n}\in\\{0,1\\}$. In practice, for the following illustration,
only two points ($N_{\mathrm{pen}}=2$) corresponding to the two local minima
on the considered potential will be penalized in such a manner with
$\widetilde{y}^{1}=0$ and $\widetilde{y}^{2}=1$.

Let us emphasize that the variance term in factor of $\lambda_{1}$ in the
above loss function is incompatible with the reconstruction term in factor of
$\lambda_{4}$, as both terms determine in some sense the range of values taken
by the encoder, but not in the same manner. Therefore, one should choose
either $\lambda_{1}>0$ or $\lambda_{4}>0$, but both parameters cannot be
simultaneously positive. Overall, the learning problem is semi-supervised as
it is the concatenation of unsupervised learning on $\\{x^{n}\\}_{1\leqslant
n\leqslant N_{\mathrm{data}}}$ and supervised learning on
$\\{\widetilde{x}^{n}\\}_{1\leqslant n\leqslant N_{\rm pen}}$.

We illustrate the behavior of models obtained by minimizing the loss function
(30) for the entropic switch potential [76]:

$\displaystyle V(x_{1},x_{2})$
$\displaystyle=3\,\mathrm{e}^{-x_{1}^{2}}\left(\mathrm{e}^{-\left(x_{2}-\frac{1}{3}\right)^{2}}-\mathrm{e}^{-\left(x_{2}-\frac{5}{3}\right)^{2}}\right)-5\,\mathrm{e}^{-x_{2}^{2}}\left(\mathrm{e}^{-\left(x_{1}-1\right)^{2}}-\mathrm{e}^{-\left(x_{1}+1\right)^{2}}\right)$
(31) $\displaystyle\;+0.2\,x_{1}^{4}+0.2\left(x_{2}-\frac{1}{3}\right)^{4}.$

The dataset of points distributed according to the associated Boltzmann–Gibbs
measure was generated following the same procedure as for the Müller–Brown
potential (see Section 3.3) with a time step $\Delta t=0.01$ and $\beta=2$.
There are $N_{\mathrm{pen}}=2$ points at which the values of the encoder are
fixed, respectively to 0 for $\widetilde{x}^{1}=(-1.05,0.04)$ and to 1 for
$\widetilde{x}^{2}=(1.05,-0.04)$. The training results of an AE model (2, 5,
5, 1, 20, 20, 2) with two decoders using the loss (30) are presented in Figure
10. Using the semi–supervised approach to constrain both decoders in each well
(Figure 10(b)) allows to differentiate correctly the top and bottom paths
linking the right and left well. Depending on the initial parameters in the
model, when the variance term is used instead (Figure 10(a)), these paths are
not as nicely separated. Such a behavior cannot be detected by observing the
convergence of the decoder to the conditional averages, or the angle between
the gradient of the encoder and the derivative of the decoder. One needs to
rely on the semi-supervised nature of the problem, and analyze decoder paths
to see whether they go from one probability mode to the other.

(a)

(b)

Figure 10: Minimum energy path, conditional averages and decoders on the
potential energy heatmap with isolevels of the encoder for models trained to
minimize the loss (28) with (a) $\lambda_{0}=1$, $\lambda_{1}=1$,
$\lambda_{2}=10^{-4}$ and $\lambda_{3}=\lambda_{4}=0$; (b) $\lambda_{1}=0$,
$\lambda_{2}=10^{-4}$ and $\lambda_{3}=\lambda_{4}=0.1$. The dataset points
are colored according to the decoders leading the the minimal reconstruction
error.

Finally, one may wonder how a model with more decoders than transition paths
would behave. We tested this on a system for which there is only one
transition path, namely the Müller–Brown potential (with the same dataset as
the one used in Section 3.3), for which we are looking for $K=2$ paths
corresponding to different decoders. In such a situation, the two paths build
for each decoders are almost parallel except in the high energy region (see
Figure 11). An analysis of the transition paths therefore reveals that there
is in fact a single transition mechanism here. We expect that such an analysis
can be similarly carried out in more complicated physical systems in order to
detect redundant decoders.

Figure 11:  Minimum energy path, conditional averages and decoders on the
potential energy heatmap with isolevels of the encoder $f\mathrm{enc}$. The
training dataset points are colored according to which decoder leads to the
minimal reconstruction error for models trained to minimize the loss (28) with
$\lambda_{0}=1$, $\lambda_{1}=1$ and $\lambda_{2}=10^{-4}$.

##  6 Alanine dipeptide in water

We study in this section the use of autoencoders on the molecular system
alanine dipeptide. We compare the quality of unregularized and regularized
autoencoders. For the latter ones, a regularization term is added to the
reconstruction loss in order to require autoencoders to parametrize the
leading eigenfunctions (corresponding to the largest non-trivial eigenvalues)
of the transfer operator of the underlying dynamics; see (32) below. Here, we
consider transfer operators rather than generators since the loss function of
the transfer operator framework does not require knowing the equation of the
underlying dynamics. This is advantageous when applied to MD applications
where data comes from sampling of an underlying dynamics (via an MD package)
whose potential function often involves many parameters and is too complicated
to be written down explicitly. Alternatively, one can consider similar
regularization terms involving eigenfunctions associated to the generator when
the underlying dynamics is governed by the overdamped Langevin equation (as in
Section 5).

#### Description of the system.

Alanine dipeptide is composed of $22$ atoms, $10$ of which are non-hydrogen
atoms. It is well known that the conformations of the system, in particular
the metastable states, can be well described by two backbone dihedral angles
$\phi$ and $\psi$. A training data consisting of $1.5\times 10^{5}$
configurations sampled along a trajectory was obtained from MD simulations of
the molecule solvated in water, using the GROMACS package [89]; see Appendix B
for further precisions. As Figure 12(a) indicates, the metastable
conformations of the system are well sampled by the training data. Figure
12(b) shows the conditional variances of the Cartesian coordinates of the $10$
non-hydrogen atoms, given the values of the two dihedral angles, computed
using trajectory data binned on a uniform grid of size $100\times 100$
according to the values of the dihedral angles. One can observe that the
conditional variance is small in the metastable region where the $\phi$ angle
is negative, but it is larger when the configuration is near the boundary of
the metastable regions or in the metastable region where $\phi$ is positive.
In essence, the conditional variance is smaller in regions of the phase-space
which have a larger weight under the target Boltzmann–Gibbs measure. The total
variance of the trajectory data is $3.988$ nm2, whereas the average of the
conditional variances (represented in Figure 12(b)) with respect to the
empirical distribution of the two dihedral angles (shown in Figure 12(a)) is
$0.167$ nm2. The latter is an approximation over a grid in latent space of the
quantity in (8) (and therefore an approximation of the reconstruction error)
in the case where the encoder $f_{\mathrm{enc}}$ is given by the two dihedral
angles $\phi,\psi$ and the decoder is the conditional mean.

(a)

(b)

Figure 12: Alanine dipeptide. (a) Histogram of dihedral angles along the
$1.5~{}\mu$s-long trajectory. (b) Conditional variances of the Cartesian
coordinates of the $10$ non-hydrogen atoms given the value of two dihedral
angles.

#### Conditional variances for unregularized autoencoders.

With the data prepared above, we started to train an autoencoder with the
standard reconstruction loss. The dimension of the state space was chosen to
be $D=30$, since we are only interested in the configurations of the $10$ non-
hydrogen atoms in alanine dipeptide and the coordinate of each atom is in
$\mathbb{R}^{3}$. The encoded dimension was set to $d=2$. For the encoder, we
employed a neural network that consists of a transformation layer (a map from
$\mathbb{R}^{30}$ to $\mathbb{R}^{30}$ without training parameters; see [98]
for implementation details) in order to guarantee both translation and
rotation invariance, and a feedforward neural network of size
$(30,20,15,10,2)$, _i.e._ an input layer of size $30$, an output layer of size
$2$, and three hidden layers of sizes $20,15,10$, respectively. For the
decoder, we employed a feedforward neural network of size $(2,10,15,20,30)$.
We used tanh as activation functions in both encoder and decoder. To train the
neural network, the dataset of size $1.5\times 10^{5}$ was split into a
training set and a test set with ratio of sizes $4:1$. We employed the Adam
optimizer [54] in PyTorch with a (large) batch size $20,000$ and learning rate
$2\times 10^{-3}$. The random seed was fixed in the training experiments
below.

(a)

(b)

(c)

Figure 13: Autoencoder learned by training the standard autoencoder on alanine
dipeptide with reconstruction loss. (a) Histogram of the encoded values along
the $1.5~{}\mu$s-long trajectory. (b) Conditional variances of coordinates of
the $10$ non-hydrogen atoms given the encoded values. (c) Mean squared
differences between the decoded coordinates (by the autoencoder) and the
conditional mean coordinates of the $10$ non-hydrogen atoms, given the encoded
values.

Figure 13 presents the results of the trained autoencoder after $2,000$
training epochs. The histogram in Figure 13(a) shows that, similarly to the
dihedral angles, the learned encoder is able to map the system’s metastable
conformations into different clusters in latent space. Figure 13(b) shows the
conditional variances of the $10$ non-hydrogen atoms given the encoded values,
computed using trajectory data binned on a uniform grid of size $100\times
100$ according to the encoded values. A comparison between Figure 12(b) and
Figure 13(b) reveals that the learned encoder indeed outperforms dihedral
angles in reducing the conditional variances. In fact, as shown in Table 2,
the average conditional variance (with respect to the empirical distribution
of the encoder shown in Figure 13(a) is $0.126$ nm2, which is smaller than the
average conditional variance in the case of dihedral angles. Figure 13(c)
shows the mean square difference between the decoded coordinates and the
conditional average coordinates of the $10$ non-hydrogen atoms given the
encoded values. It can be seen that the decoded coordinates are close to the
conditional averages in the metastable regions, while the difference is larger
in certain low-density regions.

Figure 14: Illustration of an autoencoder architecture with a regularization
component to represent the leading eigenfunction.

#### Regularized autoencoders.

We next numerically study the training of autoencoders with regularization.
The idea is to learn an encoder that is optimized not only for configuration
reconstruction but also for parametrizing the leading eigenfunctions
(corresponding to the largest non-trivial eigenvalues) of the transfer
operator of the underlying dynamics. This is achieved by adding to the
standard autoencoder architecture another neural network, which takes the
output of the encoder as input and provides the values of eigenfunctions as
output; see Figure 14 for an illustration of the architecture where the
regularization part involves the first eigenfunction only. The extension to
regularization involving $K\geqslant 2$ eigenfunctions is straightforward.

Recall the notation $f_{\mathrm{enc},\theta_{1}}$ and $f_{\rm dec,\theta_{2}}$
from (2) for the encoder and the decoder with parameters
$\theta=(\theta_{1},\theta_{2})$, and denote by
$\\{\widetilde{\varphi}_{\mathrm{reg},\widetilde{\theta}_{i}}:\mathbb{R}^{d}\rightarrow\mathbb{R}\\}_{1\leqslant
i\leqslant K}$ the functions corresponding to the $K$ neural networks with
parameters
$\widetilde{\theta}=(\widetilde{\theta}_{1},\dots,\widetilde{\theta}_{K})$
introduced to the regularization part of the architecture. The leading
eigenfunctions
$\\{\varphi_{i}:\mathbb{R}^{D}\rightarrow\mathbb{R}\\}_{1\leqslant i\leqslant
K}$ are represented as compositions of the encoder and regularizers, _i.e._
$\varphi_{i}$ is approximated by
$\widetilde{\varphi}_{\mathrm{reg},\widetilde{\theta}_{i}}\circ
f_{\mathrm{enc},\theta_{1}}$ for $i=1,\dots,K$ (recall Figure 14). We consider
the training loss

$\displaystyle\widehat{\mathscr{L}}_{\lambda_{0},\lambda_{1},\lambda_{2},\\{\omega_{i}\\}_{1\leqslant
i\leqslant K},\tau}\left(\theta,\widetilde{\theta}\right)$ (32)
$\displaystyle\qquad=\frac{\lambda_{0}}{N_{\mathrm{data}}}\sum_{n=1}^{N_{\mathrm{data}}}\big{\|}x^{n}-f_{\rm
dec,\theta_{2}}\circ
f_{\mathrm{enc},\theta_{1}}(x^{n})\big{\|}^{2}+\frac{\lambda_{1}}{\tau}\sum_{i=1}^{K}\omega_{i}\frac{\widehat{\mathscr{E}}_{N_{\mathrm{data}},\tau}(\widetilde{\varphi}_{\mathrm{reg},\widetilde{\theta}_{i}}\circ
f_{\mathrm{enc},\theta_{1}})}{\widehat{\mathrm{Var}}_{N_{\mathrm{data}}}(\widetilde{\varphi}_{\mathrm{reg},\widetilde{\theta}_{i}}\circ
f_{\mathrm{enc},\theta_{1}})}$ $\displaystyle\qquad\ \
+\lambda_{2}\sum_{1\leqslant i\leqslant j\leqslant
K}\left(\widehat{\mathrm{Cov}}_{N_{\mathrm{data}}}\left(\widetilde{\varphi}_{\mathrm{reg},\widetilde{\theta}_{i}}\circ
f_{\mathrm{enc},\theta_{1}},\widetilde{\varphi}_{\mathrm{reg},\widetilde{\theta}_{j}}\circ
f_{\mathrm{enc},\theta_{1}}\right)-\delta_{ij}\right)^{2},$

where $\tau=\ell\Delta t$ is a lag-time for some integer $\ell\geqslant 1$,
$\\{\omega_{i}\\}_{1\leqslant i\leqslant K}$ are fixed weights,
$\widehat{\mathrm{Var}}_{N_{\mathrm{data}}}(\cdot)$ is the empirical variance
defined in (29) and, for functions $f,g:\mathcal{X}\rightarrow\mathbb{R}$, we
have used the shorthand notation

$\displaystyle\widehat{\mathscr{E}}_{N_{\mathrm{data}},\tau}(g)$
$\displaystyle=\frac{1}{2(N_{\mathrm{data}}-\ell)}\sum_{n=1}^{N_{\mathrm{data}}-\ell}|g(x^{n+\ell})-g(x^{n})|^{2},$
$\displaystyle\widehat{\mathrm{Cov}}_{N_{\mathrm{data}}}(g,h)$
$\displaystyle=\frac{1}{N_{\mathrm{data}}}\sum_{n=1}^{N_{\mathrm{data}}}g(x^{n})h(x^{n})-\left(\frac{1}{N_{\mathrm{data}}}\sum_{n=1}^{N_{\mathrm{data}}}h(x^{n})\right)\left(\frac{1}{N_{\mathrm{data}}}\sum_{n=1}^{N_{\mathrm{data}}}h(x^{n})\right)\,.$

Note that (32) is a sum of the reconstruction loss and the loss in [99] for
learning eigenfunctions of transfer operators (also see [67, 19]). The term in
factor of $\lambda_{2}$ in the loss function encourages the approximate
eigenfunctions
$\\{\widetilde{\varphi}_{\mathrm{reg},\widetilde{\theta}_{i}}\circ
f_{\mathrm{enc},\theta_{1}}\\}_{1\leqslant i\leqslant K}$ to form an
orthonormal basis since

$\displaystyle\widehat{\mathrm{Cov}}_{N_{\mathrm{data}}}(\widetilde{\varphi}_{\mathrm{reg},\widetilde{\theta}_{i}}\circ
f_{\mathrm{enc},\theta_{1}},\widetilde{\varphi}_{\mathrm{reg},\widetilde{\theta}_{j}}\circ
f_{\mathrm{enc},\theta_{1}})$
$\displaystyle\qquad\xrightarrow[N_{\mathrm{data}}\to+\infty]{}\int_{\mathcal{X}}\Pi\left(\widetilde{\varphi}_{\mathrm{reg},\widetilde{\theta}_{i}}\circ
f_{\mathrm{enc},\theta_{1}}\right)\Pi\left(\widetilde{\varphi}_{\mathrm{reg}\widetilde{\theta}_{j}}\circ
f_{\mathrm{enc},\theta_{1}}\right)d\mu,$

where we introduced the centering operator $\Pi$ acting as

$(\Pi g)(x)=g(x)-\int_{\mathcal{X}}g\,d\mu.$

The various terms in factor of $\lambda_{1}$ encode the property that
$\\{\widetilde{\varphi}_{\mathrm{reg},\widetilde{\theta}_{i}}\circ
f_{\mathrm{enc},\theta_{1}}\\}_{1\leqslant i\leqslant K}$ are approximations
of the eigenfunctions of the transfer operator $P_{\tau}$ since, assuming that
$P_{\tau}$ is self-adjoint (see [99, Lemma 1 and Appendix A]),

$\frac{\widehat{\mathscr{E}}_{N_{\mathrm{data}},\tau}(g)}{\widehat{\mathrm{Var}}_{N_{\mathrm{data}}}(g)}\xrightarrow[N_{\mathrm{data}}\to+\infty]{}\frac{\displaystyle\int_{\mathcal{X}}g\left(g-P_{\tau}g\right)d\mu}{\displaystyle\int_{\mathcal{X}}\left(\Pi
g\right)^{2}d\mu},$

where $P_{\tau}g(x)=\mathbb{E}(g(x^{\ell})|x^{0}=x)$ encodes the average
evolution of the system over the time $\tau$. The quantity on the right-hand
side of the previous equality is equal to $1-\mathrm{e}^{-\rho\tau}$ when $f$
is an eigenfunction associated with the eigenvalue $\mathrm{e}^{-\rho\tau}$ of
the transfer operator $P_{\tau}=\mathrm{e}^{\tau L}$ corresponding to the
generator $L$ of the dynamics (which admits $\rho$ as dominant non-zero
eigenvalue). This quantity is minimized when $\rho$ is minimized, _i.e._ the
largest eigenvalues of the transfer operator are recovered. This expression
also motivates the factor $1/\tau$ multiplying $\lambda_{1}$, so that the
second term of the loss is of order 1 even for $\tau$ small.

In our numerical experiments, we chose $K=3$,
$(\omega_{1},\omega_{2},\omega_{3})=(1.0,0.4.0.2)$,
$(\lambda_{0},\lambda_{1},\lambda_{2})=(1.0,0.2,10.0)$, and a lag-time
$\tau=10\,\mathrm{ps}$. We adopted the same encoder, decoder and training
parameters as for the unregularized networks considered later on in this
section. For the regularizers, we employed three neural networks with
identical sizes $(2,20,20,20,1)$ in order to represent the three leading
eigenfunctions.

Figure 15 shows the jointly learned encoder $f_{\mathrm{enc},\theta_{1}}$ and
the three leading eigenfunctions $\varphi_{1},\varphi_{2},\varphi_{3}$. As one
can see in Figures 15(a)-15(b) and Table 2, similarly to the encoder trained
with the standard reconstruction loss, the encoder trained with the loss (32)
is able to achieve dimensionality reduction with a smaller mean conditional
variance (_i.e._ $0.134$ nm2) then the one obtained with dihedral angles.
Figure 15(c) shows the mean square difference between the decoded coordinates
and the conditional mean coordinates of the $10$ non-hydrogen atoms given the
encoded values. Similarly to the case of the standard reconstruction loss
(Figure 13(c)), the decoded coordinates are close to the conditional averages
in the metastable regions, while the difference is visible in certain
low–density regions. Figures 15(d) to 15(i) show the three leading
eigenfunctions jointly learned in the training. Each of these eigenfunctions
reveals information about certain transition events of the system on a large
time-scale (related to the eigenvalues in Table 2, where larger eigenvalues
correspond to transition events on longer time-scales; see [78]). The
eigenfunctions are close to being constant inside metastable regions, and
exhibit a sign change in the transition regions between metastable states.
From these results, we see that training with the loss (32) indeed yields an
encoder that is capable of reconstructing conformations and at the same time
parametrizing the first three eigenfunctions. We have also performed
experiments with $K=1$ and $K=2$ in the loss (32), for which we have obtained
similar results in terms of the encoder’s quality for both conformation
reconstruction and eigenfunction parametrization.

To summarize, training with the loss (32) yields for this example results
seemingly similar to the ones obtained with the standard reconstruction loss.
This may be due to the simplicity of the alanine dipeptide system, for which
the autoencoder trained to minimize the reconstruction error is already
sufficient to characterize the dynamics, _i.e._ the leading eigenfunctions.
Nevertheless, the numerical experiments validate the capability of the
training algorithm with the loss (32). We expect that this can provide more
prominent regularizing effects in training autoencoders on MD systems with
more complex kinetics.

CV map | cond. var. | $\nu_{1}$ | $\nu_{2}$ | $\nu_{3}$
---|---|---|---|---
dihedrals | $0.167$ | $0.99$ | $0.88$ | $0.74$
standard AE | $0.126$ | $0.99$ | $0.88$ | $0.74$
reg. AE | $0.134$ | $0.99$ | $0.88$ | $0.76$
Table 2: Comparison of different CV maps for alanine dipeptide. The CV maps
are built using either the two dihedral angles (“dihedrals”), the encoder
learned using the standard reconstruction loss (“standard AE”), and the
encoder learned using the loss (32) (“reg. AE”). Column “cond. var.” records
the average conditional variances (unit: $\mathrm{nm}^{2}$) computed using
different CV maps. Each row of the columns “$\nu_{1}$”, “$\nu_{2}$” and
“$\nu_{3}$” records the estimations of the first three non-trivial eigenvalues
of the transfer operator with $\tau=10$ ps. For the CV map “reg. AE”, these
eigenvalues and the eigenfunctions are estimated jointly during the training
with the loss (32), whereas for the CV maps “dihedrals” and “standard AE” they
are estimated by performing a subsequent training with the loss (32) where
$\lambda_{0}=0$ and the corresponding CV map is fixed.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 15: Autoencoder and first three eigenfunctions learned by training the
autoencoder (see Figure 14) with the loss (32). (a) Histogram of the encoded
values along the $1.5~{}\mu$s-long trajectory. (b) Conditional variances of
coordinates of the $10$ non-hydrogen atoms given the encoded values. (c) Mean
squared differences between the decoded coordinates and the conditional mean
coordinates of the 10 non-hydrogen atoms, given the encoded values. (d)-(f):
Scatter plots of the first three eigenfunctions $\varphi_{1}$, $\varphi_{2}$,
and $\varphi_{3}$ (as compositions of the learned encoder and regularizer)
evaluated on the trajectory data in the space of dihedral angles. (g)-(i):
Scatter plots of the first three eigenfunctions evaluated on the trajectory
data in the space of the learned encoder.

##  7 Conclusions and perspectives

We have shown in this work that the usual setting of autoencoder training in
terms of the reconstruction error can be reformulated in various ways. This
allows to relate the problem to the well-studied class of techniques for
computing principal curves and manifolds. Somewhat more interestingly in our
view, it also allows to interpret the decoder as some approximation of the
Bayes predictor associated with a given encoder, which corresponds to a
conditional expectation. The quality of the decoding can be assessed by
comparing actual conditional expectations and the output of the decoder for a
given value of the latent variable corresponding to the value of the encoder.
We also suggested various extensions beyond the usual training of autoencoders
to capture some information on transitions from one metastable state to
another, by (i) emphasizing transition points through some modified
distribution of training points, (ii) allowing for multiple decoders to
parametrize multiple transition paths, and (iii) adding extra terms to the
loss function to take into account dynamical properties.

There are various natural follow-ups of this work. A first one is to combine
numerical methods to sample reactive trajectories from one metastable state to
another and the training of autoencoders. One instance of such an idea is to
start from some (partial) data set of configurations, train an autoencoder
with a one-dimensional bottleneck, and then use the corresponding encoder as a
one dimensional reaction coordinate in the adaptive multilevel splitting
algorithm [12, 65, 77] to sample reactive trajectories, in order to complement
the initial data set, as done in Section 4. A second follow-up is to extend
our analysis to slow CVs, _i.e._ genuinely consider dynamical aspects from the
start. This can be done by working with time-lagged autoencoders [94, 18],
relying on some data set where couples of configurations separated by a fixed
time lag are considered; see [99] for elements in this direction. An
alternative option is to search for CVs satisfying some optimality conditions
in terms of their effective dynamics [60, 97]. Finally, let us recall Remark
1, which highlights that the data distributions at hand in molecular dynamics
are not unknown, in contrast to the usual machine learning setting. This extra
information can probably be used to improve off-the-shelf machine learning
methods.

#### Acknowledgements.

The work of T.L and G.S. was funded in part by the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation
programme (grant agreement No 810367). The work of W.Z. was supported by the
DFG through grant CRC 1114 “Scaling Cascades in Complex Systems” (project ID:
235221301).

##  Appendix A Derivation of (12)

We formally derive in this appendix the necessary condition (12) satisfied by
(local) maxima of (11), for $\mathcal{Z}\subset\mathbb{R}^{d}$ and
$\mathcal{X}\subset\mathbb{R}^{D}$. We assume that $\mu$ has a smooth density
with respect to the Lebesgue measure, still denoted by $\mu$ with some abuse
of notation. We denote by $\mathrm{Supp}(\mu)$ the support of this measure,
namely the closure of the set of points at which the density is positive. A
key assumption in our derivation is that the encoder is a smooth function,
more precisely $f_{\mathrm{enc}}\in C^{2}(\mathcal{X},\mathcal{Z})$, and
$\nabla f_{\mathrm{enc}}(x)\in\mathbb{R}^{D\times d}$ has full rank $d$ for
all $x\in\mathrm{Supp}(\mu)$ (we use here the convention that the columns of
$\nabla f_{\mathrm{enc}}$ are the gradients of the components of the function
$f_{\mathrm{enc}}$).

#### Analytical reformulation of the maximization problem.

Our aim here is to rewrite the problem (11) in an analytic form, where
expectations are replaced by integrals. We start by writing it as the
minimization of

$\mathscr{S}(f_{\mathrm{enc}})=\int_{\mathcal{X}}\left\|g^{\star}_{f_{\mathrm{enc}}}(f_{\mathrm{enc}}(x))\right\|^{2}\mu(x)\,dx,$

where we recall that
$g^{\star}_{f_{\mathrm{enc}}}(z)=\mathbb{E}[X\,|\,f_{\mathrm{enc}}(X)=z]$ is
defined in (6). From an analytical viewpoint, the conditional expectation can
be rewritten as

$g^{\star}_{f_{\mathrm{enc}}}(z)=\frac{\displaystyle\int_{\Sigma_{z}}x\,\mu(x)\,\delta_{f_{\mathrm{enc}}(x)-z}(dx)}{\displaystyle\int_{\Sigma_{z}}\mu(x)\,\delta_{f_{\mathrm{enc}}(x)-z}(dx)},$
(33)

where $\Sigma_{z}\subset\mathcal{X}$ is the submanifold defined in (13) (the
definition makes sense in view of the assumptions on $f_{\mathrm{enc}}$), and
the delta measure $\delta_{f_{\mathrm{enc}}(x)-z}(dx)$ on $\Sigma_{z}$ is
defined by the decomposition formula

$dx=\delta_{f_{\mathrm{enc}}(x)-z}(dx)\,dz.$

The latter equality should be understood as follows: for any bounded
measurable functions $\varphi:\mathcal{X}\to\mathbb{R}$ and
$u:\mathcal{Z}\to\mathbb{R}$ with compact supports,

$\int_{\mathcal{X}}\varphi(x)u(f_{\mathrm{enc}}(x))\,dx=\int_{\mathcal{Z}}u(z)\left(\int_{\Sigma_{z}}\varphi(x)\,\delta_{f_{\mathrm{enc}}(x)-z}(dx)\right)dz.$
(34)

As discussed in [62, Section 3.2.1.1], the delta measure
$\delta_{f_{\mathrm{enc}}(x)-z}(dx)$ can be related to the surface measure
$\sigma_{\Sigma_{z}}(dx)$ induced by the Lebesgue measure in the ambient
Euclidean space $\mathbb{R}^{D}$ through the co-area formula [29, 2] as

$\delta_{f_{\mathrm{enc}}(x)-z}(dx)=\left(\det
G(x)\right)^{-1/2}\sigma_{\Sigma_{z}}(dx),$

where we introduced the Gram matrix $G(x)=\nabla
f_{\mathrm{enc}}(x)^{\top}\nabla f_{\mathrm{enc}}(x)\in\mathbb{R}^{d\times d}$
in order to alleviate the notation.

In order to make sense of derivatives of the conditional expectation, one
possible strategy, considered for instance in [36, Section 2], is to consider
some mollified version of the conditional expectation, namely

$g^{\varepsilon}_{f_{\mathrm{enc}}}(z)=\frac{\displaystyle\int_{\mathcal{X}}x\chi_{\varepsilon}(f_{\mathrm{enc}}(x)-z)\,\mu(x)\,dx}{\displaystyle\int_{\mathcal{X}}\chi_{\varepsilon}(f_{\mathrm{enc}}(x)-z)\,\mu(x)\,dx},$

where

$\chi_{\varepsilon}(z)=\frac{1}{\varepsilon^{d}}\chi\left(\frac{z}{\varepsilon}\right)$

is an approximation of the Dirac mass at $0\in\mathbb{R}^{d}$ for
$\varepsilon\to 0$ (the function $\chi$ being a smooth nonnegative function of
integral 1). We follow here an alternative route, where we do not consider
regularization by convolutions. More precisely, we provide results ensuring
that the functions under consideration are differentiable, by relying on a
weak formulation to compute derivatives of integrals with respect to the delta
measure; which then allows us to establish (12). For completeness, we also
give an expression for the derivatives of these functions, although this will
not be used in the proof of (12).

#### Differentiability results.

A first differentiability result makes precise derivatives with respect to
$z$.

######  Lemma 2.

Fix a function $\varphi\in L^{1}(\mathcal{X})$ such that
$\mathrm{div}_{x}\left(\varphi\nabla f_{\mathrm{enc}}G^{-1}\right)\in
L^{1}(\mathcal{X})$. Then, for almost every $z\in\mathcal{Z}$,

$\nabla_{z}\left(\int_{\Sigma_{z}}\varphi(x)\,\delta_{f_{\mathrm{enc}}(x)-z}(dx)\right)=\int_{\Sigma_{z}}\mathrm{div}_{x}\left(\varphi(x)\nabla
f_{\mathrm{enc}}(x)G(x)^{-1}\right)\,\delta_{f_{\mathrm{enc}}(x)-z}(dx),$

where the divergence of the $D\times d$ matrix $\varphi(x)\nabla
f_{\mathrm{enc}}(x)G(x)^{-1}$ on the right-hand side of the previous equality
is applied column by column.

The above equality should be understood as follows: for any
$\omega\in\mathbb{R}^{d}$,

$\omega^{\top}\nabla_{z}\left(\int_{\Sigma_{z}}\varphi(x)\,\delta_{f_{\mathrm{enc}}(x)-z}(dx)\right)=\int_{\Sigma_{z}}\mathrm{div}_{x}\left(\varphi(x)\nabla
f_{\mathrm{enc}}(x)G(x)^{-1}\omega\right)\,\delta_{f_{\mathrm{enc}}(x)-z}(dx).$

The formula established in Lemma 2 allows to compute the derivative of (33)
for a given function $f_{\mathrm{enc}}$, by considering
$\varphi(x)=x_{i}\,\mu(x)$ with $1\leqslant i\leqslant D$ for the integral in
the numerator (where $x_{i}$ is the $i$-th component of $x\in\mathbb{R}^{D}$),
and $\varphi(x)=\mu(x)$ for the one in the denominator.

###### Proof.

We establish the result for a smooth bounded function $\varphi$ with compact
support, the general case following by density. We start from (34), written
for a smooth bounded function $u:\mathcal{Z}\to\mathbb{R}$ with compact
support:

$\int_{\mathcal{Z}}\nabla_{z}u(z)\left(\int_{\Sigma_{z}}\varphi(x)\,\delta_{f_{\mathrm{enc}}(x)-z}(dx)\right)dz=\int_{\mathcal{X}}\varphi(x)(\nabla_{z}u)(f_{\mathrm{enc}}(x))\,dx.$

Since $\nabla_{x}(u\circ
f_{\mathrm{enc}})=(\nabla_{x}f_{\mathrm{enc}})[(\nabla_{z}u)\circ
f_{\mathrm{enc}}]$ with our convention for $\nabla_{x}f_{\mathrm{enc}}$, we
can rewrite the derivative in $z$ as a derivative in $x$ as

$(\nabla_{z}u)(f_{\mathrm{enc}}(x))=G(x)^{-1}\nabla
f_{\mathrm{enc}}(x)^{\top}\nabla_{x}\left[u(f_{\mathrm{enc}}(x))\right].$ (35)

Therefore,

$\displaystyle\int_{\mathcal{Z}}\omega^{\top}\nabla
u(z)\left(\int_{\Sigma_{z}}\varphi(x)\,\delta_{f_{\mathrm{enc}}(x)-z}(dx)\right)dz$
$\displaystyle\qquad\qquad=\int_{\mathcal{X}}\varphi(x)\,\omega^{\top}G(x)^{-1}\nabla
f_{\mathrm{enc}}(x)^{\top}\nabla_{x}\left[u(f_{\mathrm{enc}}(x))\right]\,dx$
$\displaystyle\qquad\qquad=-\int_{\mathcal{X}}\mathrm{div}_{x}\left(\varphi(x)\nabla
f_{\mathrm{enc}}(x)G(x)^{-1}\omega\right)u(f_{\mathrm{enc}}(x))\,dx$
$\displaystyle\qquad\qquad=-\int_{\mathcal{Z}}u(z)\left(\int_{\Sigma_{z}}\mathrm{div}_{x}\left(\varphi(x)\nabla
f_{\mathrm{enc}}(x)G(x)^{-1}\omega\right)\,\delta_{f_{\mathrm{enc}}(x)-z}(dx)\right)dz.$

This allows to identify the desired derivative. ∎

In order to understand variations with respect to $f_{\mathrm{enc}}$, we
consider a smooth and bounded perturbation $\widetilde{f}_{\mathrm{enc}}$ of
$f_{\mathrm{enc}}$. We assume that there exists $\eta_{\star}>0$ (which
depends on $f_{\mathrm{enc}}$ and $\widetilde{f}_{\mathrm{enc}}$) such that,
for any $\eta\in(-\eta_{\star},\eta_{\star})$, the smooth function
$f_{\mathrm{enc}}+\eta\widetilde{f}_{\mathrm{enc}}$ has a gradient of full
rank $d$ at every $x\in\mathcal{X}$. This can be ensured for instance when
$\widetilde{f}_{\mathrm{enc}}$ has compact support. We can then introduce the
functional
$I_{\varphi}:\mathcal{Z}\times(-\eta_{\star},\eta_{\star})\to\mathbb{R}$
defined as

$I_{\varphi}(z,\eta)=\int_{\widetilde{\Sigma}_{z,\eta}}\varphi(x)\,\delta_{(f_{\mathrm{enc}}+\eta\widetilde{f}_{\mathrm{enc}})(x)-z}(dx),\qquad\widetilde{\Sigma}_{z,\eta}=\left(f_{\mathrm{enc}}+\eta\widetilde{f}_{\mathrm{enc}}\right)^{-1}\\{z\\}.$

Note that

$\mathscr{S}\left(f_{\mathrm{enc}}+\eta\widetilde{f}_{\mathrm{enc}}\right)=\int_{\mathcal{X}}\sum_{i=1}^{D}\left|\frac{\displaystyle
I_{\phi_{i}}\left((f_{\mathrm{enc}}+\eta\widetilde{f}_{\mathrm{enc}})(x),\eta\right)}{I_{\mu}\left((f_{\mathrm{enc}}+\eta\widetilde{f}_{\mathrm{enc}})(x),\eta\right)}\right|^{2}\mu(x)\,dx,\qquad\phi_{i}(x)=x_{i}\,\mu(x).$
(36)

The next result establishes the differentiability of $I_{\varphi}$ with
respect to $\eta$.

######  Lemma 3.

Fix a function $\varphi\in L^{1}(\mathcal{X})$ such that
$\mathrm{div}_{x}\left(\varphi\nabla
f_{\mathrm{enc}}G^{-1}\widetilde{f}_{\mathrm{enc}}\right)\in
L^{1}(\mathcal{X})$. Then, for almost every $z\in\mathcal{Z}$,

$\partial_{\eta}I_{\varphi}(z,0)=-\int_{\Sigma_{z}}\mathrm{div}_{x}\left(\varphi\nabla
f_{\mathrm{enc}}G^{-1}\widetilde{f}_{\mathrm{enc}}\right)(x)\,\delta_{f_{\mathrm{enc}}(x)-z}(dx).$

###### Proof.

As for the proof of Lemma 2, we establish the result for a smooth bounded
function $\varphi$ with compact support, the general case following by
density. We introduce a smooth function $u:\mathcal{Z}\to\mathbb{R}$ with
compact support and write

$\int_{\mathcal{Z}}I_{\varphi}(z,\eta)u(z)\,dz=\int_{\mathcal{X}}\varphi(x)u\left((f_{\mathrm{enc}}+\eta\widetilde{f}_{\mathrm{enc}})(x)\right)dx.$

The right-hand side of the previous equality is a smooth function of $\eta$,
which can be differentiated. Using again (35) to rewrite $(\nabla_{z}u)\circ
f_{\mathrm{enc}}$ in terms of the gradient of $u\circ f_{\mathrm{enc}}$,

$\displaystyle\left.\frac{d}{d\eta}\left(\int_{\mathcal{Z}}I_{\varphi}(z,\eta)u(z)\,dz\right)\right|_{\eta=0}$
$\displaystyle=\int_{\mathcal{X}}\varphi(x)\widetilde{f}_{\mathrm{enc}}(x)^{\top}\nabla_{z}u(f_{\mathrm{enc}}(x))\,dx$
(37)
$\displaystyle=\int_{\mathcal{X}}\varphi(x)\widetilde{f}_{\mathrm{enc}}(x)^{\top}\left(G(x)^{-1}\nabla
f_{\mathrm{enc}}(x)^{\top}\nabla_{x}[u\circ f_{\mathrm{enc}}](x)\right)dx$
$\displaystyle=-\int_{\mathcal{X}}\mathrm{div}_{x}\left(\varphi\nabla
f_{\mathrm{enc}}G^{-1}\widetilde{f}_{\mathrm{enc}}\right)(u\circ
f_{\mathrm{enc}})\,dx$
$\displaystyle=-\int_{\mathcal{Z}}u(z)\left(\int_{\Sigma_{z}}\mathrm{div}_{x}\left(\varphi\nabla
f_{\mathrm{enc}}G^{-1}\widetilde{f}_{\mathrm{enc}}\right)(x)\,\delta_{f_{\mathrm{enc}}(x)-z}(dx)\right)dz,$

which leads to the claimed result. ∎

#### Necessary condition for $f_{\mathrm{enc}}$ to be a maximizer.

We finally establish (12). We rely on the differentiablity results established
in Lemmas 2 and 3, and assume to this end that the functions $\phi_{i}$ (for
$1\leqslant i\leqslant d$) all belong to $L^{1}(\mathcal{X})$, and also that
$\mathrm{div}_{x}\left(\mu(\nabla f_{\mathrm{enc}})G^{-1}\right)$,
$\mathrm{div}_{x}\left(\mu\nabla
f_{\mathrm{enc}}G^{-1}\widetilde{f}_{\mathrm{enc}}\right)$,
$\mathrm{div}_{x}\left(\phi_{i}\nabla f_{\mathrm{enc}}G^{-1}\right)$ and
$\mathrm{div}_{x}\left(\phi_{i}\nabla
f_{\mathrm{enc}}G^{-1}\widetilde{f}_{\mathrm{enc}}\right)$ all belong to
$L^{1}(\mathcal{X})$ (recalling (36)). Denoting by
$\mathscr{I}_{\varphi}(z)=I_{\varphi}(z,0)$, it holds

$\displaystyle\left.\frac{d}{d\eta}\left[\frac{1}{2}\mathscr{S}\left(f_{\mathrm{enc}}+\eta\widetilde{f}_{\mathrm{enc}}\right)\right]\right|_{\eta=0}$
$\displaystyle\qquad\qquad=\int_{\mathcal{X}}\sum_{i=1}^{D}\left(\widetilde{f}_{\mathrm{enc}}(x)^{\top}\nabla_{z}\left(\frac{\mathscr{I}_{\phi_{i}}}{\mathscr{I}_{\mu}}\right)(f_{\mathrm{enc}}(x))+\mathscr{D}_{\phi_{i},\widetilde{f}_{\mathrm{enc}}}(f_{\mathrm{enc}}(x))\right)\frac{\displaystyle\mathscr{I}_{\phi_{i}}(f_{\mathrm{enc}}(x))}{\mathscr{I}_{\mu}(f_{\mathrm{enc}}(x))}\mu(x)\,dx,$

where

$\mathscr{D}_{\phi_{i},\widetilde{f}_{\mathrm{enc}}}(z)=\frac{\partial_{\eta}I_{\phi_{i}}(z,0)}{\mathscr{I}_{\mu}(z)}-\frac{\mathscr{I}_{\phi_{i}}(z)\partial_{\eta}I_{\mu}(z,0)}{\mathscr{I}_{\mu}(z)^{2}}.$

Note also that

$\frac{\displaystyle\mathscr{I}_{\phi_{i}}(z)}{\mathscr{I}_{\mu}(z)}=g^{\star}_{f_{\mathrm{enc}},i}(z).$

The aim is to factor out the dependence on $\widetilde{f}_{\mathrm{enc}}$ in
the integral involving $\mathscr{D}_{\phi_{i},\widetilde{f}_{\mathrm{enc}}}$.
We consider the first term on the right hand side of the above equality, and
use (37) to write

$\displaystyle\int_{\mathcal{X}}\partial_{\eta}I_{\phi_{i}}(f_{\mathrm{enc}}(x),0)\frac{\displaystyle\mathscr{I}_{\phi_{i}}(f_{\mathrm{enc}}(x))}{\mathscr{I}_{\mu}(f_{\mathrm{enc}}(x))^{2}}\mu(x)\,dx$
$\displaystyle=\int_{\mathcal{Z}}\partial_{\eta}I_{\phi_{i}}(z,0)\frac{\displaystyle\mathscr{I}_{\phi_{i}}(z)}{\mathscr{I}_{\mu}(z)^{2}}\left(\int_{\Sigma_{z}}\mu(x)\,\delta_{f_{\mathrm{enc}}(x)-z}(dx)\right)dz$
$\displaystyle=\int_{\mathcal{Z}}\partial_{\eta}I_{\phi_{i}}(z,0)\frac{\displaystyle\mathscr{I}_{\phi_{i}}(z)}{\mathscr{I}_{\mu}(z)}\,dz$
$\displaystyle=\int_{\mathcal{Z}}\partial_{\eta}I_{\phi_{i}}(z,0)g^{\star}_{f_{\mathrm{enc}},i}(z)\,dz$
$\displaystyle=\int_{\mathcal{X}}\phi_{i}(x)\widetilde{f}_{\mathrm{enc}}(x)^{\top}\nabla_{z}g^{\star}_{f_{\mathrm{enc}},i}(f_{\mathrm{enc}}(x))\,dx.$

Similarly,

$\int_{\mathcal{X}}\partial_{\eta}I_{\mu}(f_{\mathrm{enc}}(x),0)\frac{\displaystyle\mathscr{I}_{\phi_{i}}(f_{\mathrm{enc}}(x))^{2}}{\mathscr{I}_{\mu}(f_{\mathrm{enc}}(x))^{3}}\mu(x)\,dx=\int_{\mathcal{X}}\mu(x)\widetilde{f}_{\mathrm{enc}}(x)^{\top}\nabla_{z}\left[\left(g^{\star}_{f_{\mathrm{enc}},i}\right)^{2}\right](f_{\mathrm{enc}}(x))\,dx.$

A necessary condition for $f_{\mathrm{enc}}$ to be a critical point is
therefore

$\displaystyle\int_{\mathcal{X}}\widetilde{f}_{\mathrm{enc}}(x)^{\top}\left(\sum_{i=1}^{D}\nabla_{z}g^{\star}_{f_{\mathrm{enc}},i}(f_{\mathrm{enc}}(x))g^{\star}_{f_{\mathrm{enc}},i}(f_{\mathrm{enc}}(x))\right)\mu(x)\,dx$
$\displaystyle\qquad\qquad+\sum_{i=1}^{D}\int_{\mathcal{X}}\widetilde{f}_{\mathrm{enc}}(x)^{\top}\left(\frac{\phi_{i}(x)}{\mu(x)}\nabla_{z}g^{\star}_{f_{\mathrm{enc}},i}(f_{\mathrm{enc}}(x))-\nabla_{z}\left[\left(g^{\star}_{f_{\mathrm{enc}},i}\right)^{2}\right](f_{\mathrm{enc}}(x))\right)\mu(x)\,dx=0.$

Since $\widetilde{f}_{\mathrm{enc}}$ is arbitrary, this shows that the
following equality must hold for all $x\in\mathcal{X}$ and $1\leqslant
j\leqslant d$:

$\sum_{i=1}^{D}\left(\frac{\phi_{i}(x)}{\mu(x)}\partial_{z_{j}}g^{\star}_{f_{\mathrm{enc}},i}(f_{\mathrm{enc}}(x))-\partial_{z_{j}}g^{\star}_{f_{\mathrm{enc}},i}(f_{\mathrm{enc}}(x))g^{\star}_{f_{\mathrm{enc}},i}(f_{\mathrm{enc}}(x))\right)=0.$

This indeed gives (12) since $\phi_{i}(x)=x_{i}\mu(x)$ by (36).

##  Appendix B Details of the molecular dynamics simulations of Section 6

Alanine dipeptide was put into a cubic simulation box with periodic boundary
condition, with box sizes set to the diameter of the system plus $2.4$ nm.
Water molecules were added into the box, resulting in $3,469$ atoms in total.
The force field AMBER99SB-ILDN [63] and TIP3P water models were adopted for
alanine dipeptide and water, respectively. Energy minimization was performed
using steepest descent minimization, with initial step-size set to $0.01$ nm
in the GROMACS input file and maximum number of minimization steps set to
$500$. The convergence was reached after $289$ steps when the maximum of the
absolute values of the force components were below $500$ kJ mol-1 nm-1.

The MD simulation to obtain the training data was conducted using a leap-frog
stochastic dynamics integrator with step-size $1$ fs. The linear constraint
solver (LINCS) of order $4$ was adopted to impose holonomic constraints
(fixing bond lengths involving H-atoms). Verlet lists were used to compute
interactions between neighboring atoms, the frequency to update the neighbor
list being initially set to $20$ steps and the cut-off radius to $1.0$ nm.
Both Coulomb and van der Waals potentials were shifted to zero at the cut-off
distance $1.0$ nm. Electrostatic interactions were computed using the fourth
order fast smooth particle-mesh Ewald (SPME) method with grid of size $0.16$
nm. A modified Berendsen thermostat was used with time constant $1$ ps and
reference temperature $300$ K. After equilibrating the system by a NVT
simulation for $100$ ps followed by a NPT simulation for $200$ ns, we
simulated the system for $1.5~{}\mu$s and, by recording the state every $10$
ps, we obtained a dataset consisting of $1.5\times 10^{5}$ states along the
system’s trajectory.

## References

  * [1] G. Alain and Y. Bengio.  What regularized auto-encoders learn from the data-generating distribution.  J. Mach. Learn. Res, 15(110):3743–3773, 2014.
  * [2] L. Ambrosio, N. Fusco, and D. Pallara.  Functions of Bounded Variation and Free Discontinuity Problems.  Oxford Science Publications, 2000.
  * [3] F. Bach.  Learning Theory from First Principles.  MIT Press, 2023.
  * [4] J. Baima, A. M. Goryaeva, T. D. Swinburne, J.-B. Maillet, M. Nastar, and M.-C. Marinica.  Capabilities and limits of autoencoders for extracting collective variables in atomistic materials science.  Phys. Chem. Chem. Phys., 24:23152–23163, 2022.
  * [5] P. Baldi and K. Hornik.  Neural networks and principal component analysis: Learning from examples without local minima.  Neural Networks, 2:53–58, 1989.
  * [6] Z. Belkacemi, M. Bianciotto, H. Minoux, T. Lelièvre, G. Stoltz, and P. Gkeka.  Autoencoders for dimensionality reduction in molecular dynamics: Collective variable dimension, biasing, and transition states.  J. Chem. Phys., 159(2):024122, 2023.
  * [7] Z. Belkacemi, P. Gkeka, T. Lelièvre, and G. Stoltz.  Chasing collective variables using autoencoders and biased trajectories.  J. Chem. Theory and Comput., 18(1):59–78, 2022.
  * [8] C. M. Bishop.  Pattern Recognition and Machine Learning.  Information Science and Statistics. Springer-Verlag, 2006.
  * [9] L. Bonati, V. Rizzi, and M. Parrinello.  Data-driven collective variables for enhanced sampling.  J. Phys. Chem. Lett., 11(8):2998–3004, 2020.
  * [10] H. Bourlard and Y. Kamp.  Auto-association by multilayer perceptrons and singular value decomposition.  Biol. Cybern., 59(4-5):291–294, 1988.
  * [11] G. Bussi and A. Laio.  Using metadynamics to explore complex free-energy landscapes.  Nat. Rev. Phys., 2:200–212, 2020.
  * [12] F. Cérou and A. Guyader.  Adaptive multilevel splitting for rare event analysis.  Stoch. Anal. Appl., 25(2):417–443, 2007.
  * [13] K. Champion, B. Lusch, J. N. Kutz, and S. L. Brunton.  Data-driven discovery of coordinates and governing equations.  Proc. Natl. Acad. Sci. U. S. A., 116(45):22445–22451, 2019.
  * [14] H. Chen and C. Chipot.  Chasing collective variables using temporal data-driven strategies.  QRB Discovery, 4:e2, 2023.
  * [15] H. Chen, H. Liu, H. Feng, H. Fu, W. Cai, X. Shao, and C. Chipot.  MLCV: Bridging machine-learning-based dimensionality reduction and free-energy calculation.  J. Chem. Inf. Model., 62(1):1–8, 2022.
  * [16] M. Chen.  Collective variable-based enhanced sampling and machine learning.  Eur. Phys. J. B, 94:211, 2021.
  * [17] W. Chen and A.L. Ferguson.  Molecular enhanced sampling with autoencoders: On-the-fly collective variable discovery and accelerated free energy landscape exploration.  J. Comput. Chem., 39(25):2079–2102, 2018.
  * [18] W. Chen, H. Sidky, and A. L. Ferguson.  Capabilities and limitations of time-lagged autoencoders for slow mode discovery in dynamical systems.  J. Chem. Phys., 151:064123, 2019.
  * [19] W. Chen, H. Sidky, and A. L. Ferguson.  Nonlinear discovery of slow molecular modes using state-free reversible VAMPnets.  J. Chem. Phys., 150(21):214114, 2019.
  * [20] W. Chen, A.R. Tan, and A.L. Ferguson.  Collective variable discovery and enhanced sampling using autoencoders: Innovations in network architecture and error function design.  J. Chem. Phys., 149(7):072–312, 2018.
  * [21] C. Chipot and A. Pohorille, editors.  Free Energy Calculations, volume 86 of Springer Series in Chemical Physics.  Springer, 2007.
  * [22] B. Dai and D. Wipf.  Diagnosing and enhancing VAE models.  In International Conference on Learning Representations, 2019.
  * [23] E. Darve, D. Rodríguez-Gómez, and A. Pohorille.  Adaptive biasing force method for scalar and vector free energy calculations.  J. Chem. Phys., 128(14):144120, 2008.
  * [24] T. Duchamp and W. Stuetzle.  Extremal properties of principal curves in the plane.  Ann. Statist., 24(4):1511–1520, 1996.
  * [25] A. F. Duque, S. Morin, G. Wolf, and K. Moon.  Extendable and invertible manifold learning with geometry regularized autoencoders.  In 2020 IEEE International Conference on Big Data (Big Data), Atlanta, GA, USA, pages 5027–5036, 2020.
  * [26] J. D. Durrant and J. A. McCammon.  Molecular dynamics simulations and drug discovery.  BMC Biol., 9:71, 2011.
  * [27] W. E and E. Vanden-Eijnden.  Towards a theory of transition paths.  J. Stat. Phys., 123(3):503–523, 2006.
  * [28] L. Evans, M. K. Cameron, and P. Tiwary.  Computing committors via Mahalanobis diffusion maps with enhanced sampling data.  J. Chem. Phys., 157(21):214107, 2022.
  * [29] L. C. Evans and R. F. Gariepy.  Measure Theory and Fine Properties of Functions.  Studies in Advanced Mathematics. CRC Press, 1992.
  * [30] A. L. Ferguson.  Machine learning and data science in soft materials engineering.  J. Phys. Condens. Matter, 30(4):043002, 2017.
  * [31] G. Fiorin, M. L. Klein, and J. Hénin.  Using collective variables to drive molecular dynamics simulations.  Mol. Phys., 111(22-23):3345–3362, 2013.
  * [32] A. Fischer.  Deux méthodes d’apprentissage non supervisé : synthèse sur la méthode des centres mobiles et présentation des courbes principales.  J. SFdS, 155(2):2–35, 2014.
  * [33] M. Frassek, A. Arjun, and P. G. Bolhuis.  An extended autoencoder model for reaction coordinate discovery in rare event molecular dynamics datasets.  J. Chem. Phys., 155(6):064103, 2021.
  * [34] S. Gerber.  Saddlepoints in unsupervised least squares.  arXiv preprint, 2104.05000, 2021.
  * [35] S. Gerber, T. Tasdizen, and R. Whitaker.  Dimensionality reduction and principal surfaces via kernel map manifolds.  In 2009 IEEE 12th International Conference on Computer Vision, pages 529–536, 2009.
  * [36] S. Gerber and R. Whitaker.  Regularization-free principal curve estimation.  J. Mach. Learn. Res., 14:1285–1302, 2013.
  * [37] P. Ghosh, M. S. M. Sajjadi, A. Vergari, M. Black, and B. Scholkopf.  From variational to deterministic autoencoders.  In International Conference on Learning Representations, 2020.
  * [38] P. Gkeka, G. Stoltz, A. Barati Farimani, Z. Belkacemi, M. Ceriotti, J. D. Chodera, A. R. Dinner, A. L. Ferguson, J.-B. Maillet, H. Minoux, C. Peter, F. Pietrucci, A. Silveira, A. Tkatchenko, Z. Trstanova, R. Wiewiora, and T. Lelièvre.  Machine learning force fields and coarse-grained variables in molecular dynamics: Application to materials and biological systems.  J. Chem. Theory and Comput., 16(8):4757–4775, 2020.
  * [39] A. Glielmo, B. E. Husic, A. Rodriguez, C. Clementi, F. Noé, and A. Laio.  Unsupervised learning methods for molecular simulation data.  Chem. Rev., 121(16):9722–9758, 2021.
  * [40] I. Goodfellow, Y. Bengio, and A. Courville.  Deep Learning.  MIT Press, 2016.
  * [41] T. Hastie and W. Stuetzle.  Principal curves.  J. Amer. Statist. Assoc., 84(406):502–516, 1989.
  * [42] J. Hénin, T. Lelièvre, M. R. Shirts, O. Valsson, and L. Delemotte.  Enhanced sampling methods for molecular dynamics simulations [article v1.0].  Living J. Comp. Mol. Sci., 4(1):1583, 2022.
  * [43] G. Henkelman, G. Jóhannesson, and H. Jónsson.  Methods for finding saddle points and minimum energy paths.  In S. D. Schwartz, editor, Theoretical Methods in Condensed Phase Chemistry, pages 269–302. Springer Netherlands, Dordrecht, 2002.
  * [44] C. X. Hernández, H. K. Wayment-Steele, M. M. Sultan, B. E. Husic, and V. S. Pande.  Variational encoding of complex dynamics.  Phys. Rev. E, 97:062412, 2018.
  * [45] G. E. Hinton and R. R. Salakhutdinov.  Reducing the dimensionality of data with neural networks.  Science, 313(5786):504–507, 2006.
  * [46] S. Hochreiter and J. Schmidhuber.  Long Short-Term Memory.  Neural Comput., 9(8):1735–1780, 1997.
  * [47] S. A. Hollingsworth and R. O. Dror.  Molecular dynamics simulation for all.  Neuron, 99(6):1129–1143, 2018.
  * [48] K. Iyer, R. Bhalodia, and S. Elhabian.  RENs: Relevance encoding networks.  arXiv preprint, 2205.13061, 2022.
  * [49] K. Jia, L. Sun, S. Gao, Z. Song, and B. E. Shi.  Laplacian auto-encoders: An explicit learning of nonlinear data manifold.  Neurocomputing, 160:250–260, 2015.
  * [50] R. Ketkaew, F. Creazzo, and S. Luber.  Machine learning-assisted discovery of hidden states in expanded free energy space.  J. Phys. Chem. Lett., 13(7):1797–1805, 2022.
  * [51] R. Ketkaew and S. Luber.  DeepCV: A deep learning framework for blind search of collective variables in expanded configurational space.  J. Chem. Inf. Model., 62(24):6352–6364, 2022.
  * [52] D. P. Kingma and M. Welling.  Auto-encoding variational Bayes.  In Yoshua Bengio and Yann LeCun, editors, 2nd International Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April 14-16, 2014, Conference Track Proceedings, 2014.
  * [53] D. P. Kingma and M. Welling.  An introduction to variational autoencoders.  Foundations and Trends® in Machine Learning, 12(4):307–392, 2019\.
  * [54] D.P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  3rd International Conference for Learning Representations, 2014\.
  * [55] S. Kirov and D. Slepcev.  Multiple penalized principal curves: Analysis and computation.  J. Math. Imaging Vis., 59:234–256, 2017.
  * [56] S. Klus, F. Nüske, P. Koltai, H. Wu, I. Kevrekidis, C. Schütte, and F. Noé.  Data-driven model reduction and transfer operator approximation.  J. Nonlinear Sci., 28(3):985–1010, 2018.
  * [57] M. A. Kramer.  Nonlinear principal component analysis using autoassociative neural networks.  AIChE Journal, 37(2):233–243, 1991.
  * [58] A. Laio and M. Parrinello.  Escaping free-energy minima.  Proc. Natl. Acad. Sci. USA, 99(20):12562–12566, 2002.
  * [59] Y. Lee, S. Yoon, M. Son, and F. C. Park.  Regularized autoencoders for isometric representation learning.  In International Conference on Learning Representations, 2022.
  * [60] F. Legoll and T. Lelièvre.  Effective dynamics using conditional expectations.  Nonlinearity, 23(9):2131–2163, 2010.
  * [61] T. Lelièvre.  Two mathematical tools to analyze metastable stochastic processes.  In A. Cangiani, R. L. Davidchack, E. Georgoulis, A. N. Gorban, J. Levesley, and M. V. Tretyakov, editors, Proceedings of the ENUMATH 2011, the 9th European Conference held at the University of Leicester, Leicester, September 5–9, 2011, pages 791–810. Springer, 2013\.
  * [62] T. Lelièvre, M. Rousset, and G. Stoltz.  Free-energy Computations: A Mathematical Perspective.  Imperial College Press, 2010.
  * [63] K. Lindorff-Larsen, S. Piana, K. Palmo, P. Maragakis, J. L. Klepeis, R. O. Dror, and D. E. Shaw.  Improved side-chain torsion potentials for the amber ff99sb protein force field.  Proteins: Struct. Funct., 78(8):1950–1958, 2010.
  * [64] X. Liu, H. Chen, and C. Ortner.  Stability of the minimum energy path.  arXiv preprint, 2204.00984, 2022.
  * [65] L. J. S. Lopes and T. Lelièvre.  Analysis of the adaptive multilevel splitting method with the alanine di-peptide’s isomerization.  J. Comput. Chem., 40(11):1198–1208, 2019.
  * [66] J. Lu and J. Nolen.  Reactive trajectories and the transition path process.  Probab. Theory Related Fields, 161(1-2):195–244, 2015.
  * [67] A. Mardt, L. Pasquali, H. Wu, and F. Noé.  VAMPnets for deep learning of molecular kinetics.  Nat. Commun., 9:5, 2018.
  * [68] P. Mehta, M. Bukov, C.-H. Wang, A. G. R. Day, C. Richardson, C. K. Fisher, and D. J. Schwab.  A high-bias, low-variance introduction to machine learning for physicists.  Physics Reports, 810:1–124, 2019.
  * [69] L. Molgedey and H. G. Schuster.  Separation of a mixture of independent signals using time delayed correlations.  Phys. Rev. Lett., 72:3634–3637, 1994.
  * [70] K. Müller and L. D. Brown.  Location of saddle points and minimum energy paths by a constrained simplex optimization procedure.  Theoretica Chimica Acta, 53:75–93, 1979.
  * [71] K. P. Murphy.  Probabilistic Machine Learning: An introduction.  MIT Press, 2022.
  * [72] Yusuke Naritomi and Sotaro Fuchigami.  Slow dynamics of a protein backbone in molecular dynamics simulation revealed by time-structure based independent component analysis.  J. Chem. Phys., 139(21):215102, 2013.
  * [73] A. Ng.  Sparse autoencoder.  In CS294A Lecture notes, volume 72, pages 1–19, 2011.
  * [74] F. Noé, A. Tkatchenko, K.-R. Müller, and C. Clementi.  Machine learning for molecular simulation.  Annual Review of Physical Chemistry, 71(1):361–390, 2020.
  * [75] F. Nüske, B. G. Keller, G. Pérez-Hernández, A. S. J. S. Mey, and F. Noé.  Variational approach to molecular kinetics.  J. Chem. Theory and Comput., 10(4):1739–1752, 2014.
  * [76] S. Park, M. K. Sener, D. Lu, and K. Schulten.  Reaction paths based on mean first-passage times.  J. Chem. Phys., 119(3):1313–1319, 2003.
  * [77] T. Pigeon, G. Stoltz, M. Corral-Valero, A. Anciaux-Sedrakian, M. Moreaud, T. Lelièvre, and P. Raybaud.  Computing surface reaction rates by adaptive multilevel splitting combined with machine learning and ab initio molecular dynamics.  J. Chem. Theory Comput., 19(12):3538–3550, 2023.
  * [78] J.-H. Prinz, H. Wu, M. Sarich, B. Keller, M. Senne, M. Held, J. D. Chodera, Ch. Schütte, and F. Noé.  Markov models of molecular kinetics: Generation and validation.  J. Chem. Phys., 134(17):174105, 2011.
  * [79] M. Ramil, C. Boudier, A. M. Goryaeva, M.-C. Marinica, and J.-B. Maillet.  On sampling minimum energy path.  J. Chem. Theory and Comput., 18(10):5864–5875, 2022.
  * [80] J. M. L. Ribeiro, P. Bravo, Y. Wang, and P. Tiwary.  Reweighted autoencoded variational Bayes for enhanced sampling (RAVE).  J. Chem. Phys., 149(7):072301, 2018.
  * [81] S. Rifai, P. Vincent, X. Muller, X. Glorot, and Y. Bengio.  Contractive auto-encoders: Explicit invariance during feature extraction.  In L. Getoor and T. Scheffer, editors, ICML, pages 833–840. Omnipress, 2011.
  * [82] C. Schütte, S. Klus, and C. Hartmann.  Overcoming the timescale barrier in molecular dynamics: Transfer operators, variational principles and machine learning.  Acta Numerica, 32:517–673, 2023.
  * [83] S. Shalev-Shwartz and S. Ben-David.  Understanding Machine Learning: From Theory to Algorithms.  Cambridge University Press, 2014.
  * [84] H. Sidky, W. Chen, and A. L. Ferguson.  Machine learning for collective variable discovery and enhanced sampling in biomolecular simulation.  Mol. Phys., 118(5):e1737742, 2020.