# Learning inducing points and uncertainty on molecular data by scalable
variational Gaussian processes

Mikhail Tsitsvero Institute for Chemical Reaction Design and Discovery (WPI-
ICReDD), Hokkaido University, Sapporo 001-0021, Japan,
email<EMAIL_ADDRESS>Mingoo Jin Institute for Chemical Reaction Design
and Discovery (WPI-ICReDD), Hokkaido University, Sapporo 001-0021, Japan
Andrey Lyalin Department of Chemistry, Faculty of Science, Hokkaido
University, Sapporo 060-0810, Japan; Research Center for Energy and
Environmental Materials (GREEN), National Institute for Materials Science,
Namiki 1-1, Tsukuba 305-0044, Japan.

###### Abstract

Uncertainty control and scalability to large datasets are the two main issues
for the deployment of Gaussian process (GP) models within the autonomous
machine learning-based prediction pipelines in material science and chemistry.
One way to address both of these issues is by introducing the latent inducing
point variables and choosing the right approximation for the marginal log-
likelihood objective. Here, we empirically show that variational learning of
the inducing points in a molecular descriptor space improves the prediction of
energies and atomic forces on two molecular dynamics datasets. First, we show
that variational GPs can learn to represent the configurations of the
molecules of different types that were not present within the initialization
set of configurations. We provide a comparison of alternative log-likelihood
training objectives and variational distributions. Among several evaluated
approximate marginal log-likelihood objectives, we show that predictive log-
likelihood provides excellent uncertainty estimates at the slight expense of
predictive quality. Furthermore, we extend our study to a large molecular
crystal system, showing that variational GP models perform well for predicting
atomic forces by efficiently learning a sparse representation of the dataset.

##  1  Introduction

Machine learning methods can significantly speed up the computational
prediction of the physical and chemical properties of atomic systems by
substituting the time-consuming electronic structure calculations with a
surrogate predictive model. As predictive models for physical properties, such
as energies and force fields, tend to become the active-learning components of
the larger pipelines for material and chemical reaction design [38, 41, 46,
30, 8], it is essential for such models to autonomously reason about the
uncertainty of predictions. Although significant progress has been made in
designing accurate models for approximating the potential energy surfaces [5,
47, 12, 10, 11, 39], robust uncertainty control is still in an under-developed
state to be suitable for direct use in the systems with autonomous reasoning.
Uncertainty quantification is typically addressed either by Bayesian model
averaging, see e.g. [50, 23] or by directly fitting the uncertainty
quantification parameters of the predictive distribution, for example, within
the framework of Gaussian processes [33, 40, 3, 45]. In this work, we follow
the second route and address the problem of learning the uncertainty
parameters by scalable Gaussian processes. For an extensive review and
applications of Gaussian processes to materials and molecules, we refer the
reader to Ref. [14].

Gaussian processes (GP) are often deemed the gold standard in machine
learning, however, their application is hindered by poor scalability and over-
estimation of a true uncertainty [24]. Parameters of the GPs are typically
learned by maximizing the marginal log-likelihood objective function by the
backpropagation method [18] and stochastic gradient descent [22], [37].
Evaluation of exact marginal log-likelihood objective is computationally
expensive for large datasets, therefore, several scalable approximations of
the marginal log-likelihood were proposed recently [27, 24] that are based on
the Evidence Lower BOund (ELBO) inequalities. One way to address the
scalability issue of the exact GPs is to sparsify them by introducing latent
inducing point variables [43]. In the following, we refer to the GP model
augmented by the variational parameters that are learned by the stochastic
gradient descent optimization of some objective function as the Stochastic
Variational Gaussian Process (SVGP) model.

The main contribution of our work is the application of SVGP framework to fit
the potential energies and atomic forces of molecular systems with continuous
molecular descriptors [4]. We start by describing how a GP model can be
augmented with inducing point latent variables to achieve scalability and
improved uncertainty control. Then we present several options for
approximating the marginal log-likelihood objective as well as outline options
for modeling the variational distributions over the inducing variables. We
illustrate the performance of the variational GP models with two numerical
examples. First, we fit the potential energies of molecular configurations
from Molecular Dynamics (MD) trajectories of three structural isomers of a
$12$-atom molecule. We evaluate a collection of scalable GP models in terms of
both the prediction quality on the test sets and the ability to separate the
uncertainty for the training and test configurations. We demonstrate that
learning the inducing points significantly improves prediction quality and
show that inducing points initialized on the conformations of one kind of
molecule can learn to represent the conformations of different kinds of
molecules. Second, we illustrate that the SVGP framework also performs well
for fitting the atomic forces. We provide an example of fitting the atomic
forces of a large $528$-atom molecular crystal [25] by the SVGP framework and
show that predictions by the sparsified variational models are on par with the
predictions by the exact GP models.

Additionally, in the last section, we briefly discuss whether a GP model makes
predictions on test configurations in high-dimensional descriptor space by
interpolating within the training set of data points. By following the
arguments by [1] we show, similarly to [53], that most test samples indeed
fall outside of the convex hull of densely sampled training points in the
high-dimensional descriptor space, and therefore conclude that the GP model
almost always makes predictions in the extrapolation regime.

##  2 Fitting molecular data with scalable Gaussian processes

###  2.1 Invariant Descriptors

To construct the rotation and translation invariant representations of atomic
configurations in 3D space various methods were proposed where the data
representation is either learned jointly with model parameters or constructed
deterministically. Here we employ the widely used Smooth Overlap of Atomic
Positions (SOAP) continuous and differentiable descriptors [4, 13]. We briefly
descirbe the construction of the invariant vectors. The atomic structure
$\mathcal{A}$, composed of atoms with atomic numbers $Z$, is described by the
set of atomic positions $\\{\mathbf{R}_{i}\\}_{i\in\mathcal{A}}$ with respect
to some selected origin which is referred to as the center. To build a
representation of the entire atomic structure, we first select a collection of
centers, capture the local environments around each center and combine them
into global representation, see e.g. Ref. [9] for details. The local atomic
structure can be mapped to smooth atomic density by
$\rho^{Z}(\mathbf{r})=\sum_{i}^{\left|Z_{i}\right|}e^{-1/2\sigma^{2}\left|\mathbf{r}-\mathbf{R}_{i}\right|^{2}}f_{cut}(\mathbf{r})$,
where $f_{cut}$ is a cutoff function allowing to take into account only the
atoms within short range $[0,r_{cut}]$ from the center. To build the simplest
invariant representation vector, the atomic density is first projected onto
the set of radial and spherical basis functions by

(1)
$c_{nlm}^{Z}=\int_{\mathbb{R}^{3}}g_{n}(r)Y_{lm}(\theta,\phi)\rho^{Z}(\mathbf{r})\mathrm{~{}d}V.$

Then, the simplest symmetry-invariant representation for a chosen center,
referred to as SOAP power spectrum, is given by summation over the magnetic
quantum number indices $m$:

(2)
$p_{nn^{\prime}l}^{Z_{1}Z_{2}}=\pi\sqrt{\frac{8}{2l+1}}\sum_{m}\left(c_{nlm}^{Z_{1}}\right)^{*}\left(c_{n^{\prime}lm}^{Z_{2}}\right).$

As a radial basis set one can choose, for example, the Gaussian Type Orbital
(GTO) basis functions. Depending on the complexity of the molecular dataset
and size of the atomic structures, only the expansion coefficients with radial
and spherical indices of the orders up to $n_{max}$ and $l_{max}$ are
considered. A simple way to build a vector representing the entire structure
is just to average the local power spectra over all the centers indexed by $i$

(3)
$\mathbf{x}\propto\sum_{i}\sum_{m}\left(c_{nlm}^{i,Z_{1}}\right)^{*}\left(c_{n^{\prime}lm}^{i,Z_{2}}\right).$

Since all maps required to produce the feature vector $\mathbf{x}$ from atomic
coordinates are continuous and differentiable, the derivatives of the feature
vector components with respect to atomic positions can also be computed. The
feature vectors $\mathbf{x}$ are then used as inputs to a GP model to fit the
global molecular property such as the potential energy of a molecular
configuration. On the other hand, the derivatives of the power spectra with
respect to atomic positions are used as inputs to a GP model to fit local
properties such as atomic forces. Since recently, a variety of efficient
implementations for computing the invariants emerged, such as Librascal [51],
Quip [6], and Dscribe [21]. While we briefly described the simplest invariant
representation that captures the $3$-body correlations, a more expressive
representations capturing higher-order correlations sometimes are required, we
refer the reader to [51] for a detailed review.

###  2.2 Scalable Variational Gaussian processes

Gaussian processes offer rich non-parametric models for regression and
classification tasks [52]. A GP is fully specified by a mean function
$\mu:\mathbb{R}^{d}\rightarrow\mathbb{R}$ and kernel or covariance function
$k:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$. Choosing a
different kernel and a mean allows us to encode the prior information on the
process that generated data. A prototypical example of a kernel that we will
use across this paper is the scaled Radial Basis Function (RBF) kernel with
Automatic Relevance Determination (ARD) of features [52], given by

(4)
$k\left(\mathbf{x}-\mathbf{x}^{\prime}\mid\theta\right)=\theta_{s}\exp\left(-\frac{1}{2}\sum_{i=1}^{d}\frac{1}{\theta_{i}^{2}}\left(x_{i}-x_{i}^{\prime}\right)^{2}\right),$

where $\theta_{s}>0$ is a scale parameter and $\theta_{i}>0$ are the
lengthscale parameters; this type of kernel allows to assign short
lengthscales over the irrelevant dimensions in the descriptor space.

In the most common case of a Gaussian likelihood, the predictive distribution
of the target value $y^{*}$ evaluated at the test location $\mathbf{x}^{*}$
can be computed in closed form

(5)
$p\left(y^{*}\mid\mathbf{x}^{*}\right)=\mathcal{N}\left(y^{*}\mid\mu_{\mathbf{f}}\left(\mathbf{x}^{*}\right),\sigma_{\mathbf{f}}\left(\mathbf{x}^{*}\right)^{2}+\sigma_{\text{obs
}}^{2}\right),$

where $\sigma_{\mathbf{f}}(\mathbf{x}^{*})^{2}$ is the local latent function
variance parameter and $\sigma_{\text{obs}}^{2}$ is a parameter representing
the global observation noise. The variance of the prediction $y^{*}$ evaluated
at the test point $\mathbf{x}^{*}$ is given by the sum of the global
observation noise parameter $\sigma^{2}_{\mathrm{obs}}$ and the local latent
function variance $\sigma_{\mathbf{f}}(\mathbf{x}^{*})^{2}$:

(6)
$\sigma\left[y^{*}\mid\mathbf{x}^{*}\right]=\sigma_{\mathrm{obs}}^{2}+\sigma_{\mathbf{f}}\left(\mathbf{x}^{*}\right)^{2}.$

Typically, all the parameters of the GP are jointly learned with
backpropagation by maximizing the marginal likelihood

(7) $p(\mathbf{y}\mid\mathbf{X})=\int
p\left(\mathbf{y}\mid\mathbf{f},\sigma_{\mathrm{obs}}^{2}\right)p(\mathbf{f}\mid\mathbf{X})\mathrm{~{}d}\mathbf{f},$

where $\mathbf{y}$ is the vector with real-valued targets, $\mathbf{f}$ are
the latent function values and
$\mathbf{X}=\left\\{\mathbf{x}_{i}\right\\}_{i=1}^{N}$ are the $N$ input
locations with $\mathbf{x}_{i}\in\mathbb{R}^{d}$. Unfortunately, expression
(7) is too expensive to compute when number of training samples $N$ is
becoming large, since computation of the inverse of the kernel matrix
$\mathbf{K}_{NN}^{-1}$ scales as $\mathcal{O}\left(N^{3}\right)$, where
$\mathbf{K}_{NN}=k(\mathbf{X},\mathbf{X})$. This motivates the design of the
approximate methods allowing for scalable GP inference.

Inducing point methods [42, 43, 19] offer a scalable extension to the exact
GPs. The main idea behind the inducing point methods is to replace the large
kernel matrix $\mathbf{K}_{NN}$ with its low rank approximation
$\mathbf{K}_{MM}$ by introducing the latent variables
$\mathbf{u}\in\mathbb{R}^{M}$ which are evaluated over the set of $M$ inducing
points $\mathbf{Z}=\left\\{\mathbf{z}_{m}\right\\}_{m=1}^{M}$, with each
$\mathbf{z}_{m}\in\mathbb{R}^{d}$. Inducing point locations can be set to be
the learnable variational parameters of the model. In this scenario, the GP
prior in (7) is augmented by

(8) $p(\mathbf{f}\mid\mathbf{X})\rightarrow
p(\mathbf{f}\mid\mathbf{u},\mathbf{X},\mathbf{Z})p(\mathbf{u}\mid\mathbf{Z}).$

Using Jensen’s inequality we can lower bound the log joint density over the
targets and inducing variables

(9) $\displaystyle\log p(\mathbf{y},\mathbf{u}\mid\mathbf{X},\mathbf{Z})$
$\displaystyle=\log\int
p(\mathbf{y}\mid\mathbf{f})p(\mathbf{f}\mid\mathbf{u})p(\mathbf{u})\mathrm{~{}d}\mathbf{f}$
$\displaystyle\geq\mathbb{E}_{p(\mathbf{f}\mid\mathbf{u})}[\log
p(\mathbf{y}\mid\mathbf{f})+\log p(\mathbf{u})].$

Finally, by introducing the variational distribution $q(\mathbf{u})$ over the
inducing variables, we come up with the following lower bounds for the
evidence $p(\mathbf{y})$ [20, 19]

(10) $\displaystyle\log p(\mathbf{y})$
$\displaystyle\geq\mathbb{E}_{q(\mathbf{u})}[\log
p(\mathbf{y},\mathbf{u}\mid\mathbf{X},\mathbf{Z})]+H[q(\mathbf{u})]$
$\displaystyle\geq\mathbb{E}_{q(\mathbf{u})}\left[\mathbb{E}_{p(\mathbf{f}\mid\mathbf{u})}[\log
p(\mathbf{y}\mid\mathbf{f})]\right]-\mathrm{KL}[q(\mathbf{u})\|p(\mathbf{u})],$

where $H[q(\mathbf{u})]$ is the entropy term and
$\mathrm{KL}[q(\mathbf{u})\|p(\mathbf{u})]$ is the Kullback-Leibler divergence
between the distributions $q(\mathbf{u})$ and $p(\mathbf{u})$.

In the most common scenario, the variational distribution over the inducing
variables $\mathbf{u}$ is given by the multivariate normal
$q(\mathbf{u})=\mathcal{N}(\mathbf{m},\mathbf{S})$ parametrized by mean vector
$\mathbf{m}\in\mathrm{R}^{M}$ and covariance matrix
$\mathbf{S}\in\mathrm{R}^{M\times M}$. Depending on the desired expressivity
of the GP model and complexity of the data, we can choose to learn either the
full positive semidefinite symmetric covariance matrix with non-zero non-
diagonal elements, use the mean-field approximation by restricting
$\mathbf{S}$ to be diagonal, or even collapse the multivariate normal
$q(\mathbf{u})$ to a Dirac delta distribution by taking the formal limit
$\mathbf{S}\rightarrow 0$.

Next, we briefly present several ELBO-based approximations to the high-cost
marginal log-likelihood objective that enable the scalable GP inference. Note
how all three objectives will be naturally written as the sum of the data-fit
term and the regularization term, thus explicitly targeting the model close to
the desired prior $p(\mathbf{u})$ that fits the observed data. It is crucially
important that all three objectives allow for data subsampling, therefore,
enabling us to use the stochastic gradient-based optimization methods [22].
The first approximation directly extends (10) by adding a regularization
parameter $\beta$, we refer to it as variational ELBO objective:

(11)
$\displaystyle\mathcal{L}_{\mathrm{ELBO}}=\sum_{i=1}^{N}\mathbb{E}_{q\left(f_{i}\right)}\left[\log
p\left(y_{i}\mid
f_{i}\right)\right]-\beta\mathrm{KL}[q(\mathbf{u})\|p(\mathbf{u})],$

where $q\left(f_{i}\right):=\int
p\left(f_{i}\mid\mathbf{u},\mathbf{x}_{i}\right)q(\mathbf{u})\mathrm{~{}d}\mathbf{u}$
and $\beta$ is a regularization parameter.

Another objective based on ELBO, referred to as $\gamma$-robust [27], aims at
providing robustness to the training of GPs with respect to model
misspecification and parameters learning, where the log-likelihood term in
ELBO is replaced with $\gamma$-divergence:

(12)
$\displaystyle\mathcal{L}_{\gamma}=\sum_{i=1}^{N}\mathbb{E}_{q(f_{i})}\left[-\frac{\gamma}{\gamma-1}\frac{p\left(y_{i}\mid
f_{i}\right)^{\gamma-1}}{\int p\left(y\mid
f_{i}\right)^{\gamma}\mathrm{~{}d}y}\right]-\beta\mathrm{KL}[q(\mathbf{u})\|p(\mathbf{u})],$

where $\gamma$ is a hyperparameter and is set to be equal to $1.03$ for the
rest of the paper.

The third objective, proposed by [24], aims at improving the uncertainty
estimates of the latent function and is referred to as predictive log-
likelihood objective:

(13)
$\mathcal{L}_{\mathrm{PLL}}=\sum_{i=1}^{N}\log\mathbb{E}_{q(f_{i})}\left[p\left(y_{i}\mid
f_{i}\right)\right]-\beta\mathrm{KL}[q(\mathbf{u})\|p(\mathbf{u})].$

For the presented objectives that approximate the evidence, we can calculate
the derivatives efficiently and perform the parameter learning in a black-box
manner [36].

In all of the above approximations of ELBO, the first term, which ensures
fitness to the training data, is tightly balanced by the second regularization
term, which additionally biases the search towards variational distributions
with higher entropy, as can be seen from the first bound in (10). The KL
divergence term may also be replaced by the more general Rényi divergence
[28], thus allowing for wider classes of variational distributions to be
explored.

##  3 Experiments

In this section, we provide numerical experiments aiming to show the benefits
of using variational GPs for fitting potential energy and atomic forces for
two molecular systems. The first system is given by the artificial example of
isomers of a $12$-atom molecule, illustrated in Fig. 1. For this system we
evaluate the performance of the variational GP models with respect to the
choice of the training objective, approximations of the variational
distributions, and number of inducing points. The second system is given by
the large molecular crystal [25], for which we show that variational GP models
can efficiently learn a sparse representation of the derivatives of the SOAP
feature vectors and provide predictions of atomic forces on par with exact GP
models.

###  3.1 Fitting energies of molecular isomers

Here we show that using the SVGP model augmented by the variationally learned
inducing points can improve the prediction quality, especially on the
configurations that have not been represented by the initialization selection
of data points from the training set. We proceed with the example of fitting
the potential energies of the conformations from the MD trajectories of three
structural isomer molecules, as illustrated in Fig. 1, with a collection of
variational GPs. The target values $\mathbf{y}$ that we aim to fit in this
example are the normalized potential energies of molecules from Fig. 1. For
each conformer, we generated the short $20$ ps MD trajectory, which was then
divided into the contiguous train and test parts. Each molecular conformation
was mapped to the $720$-dimensional feature vector in molecular descriptor
space, i.e. the input dimension for the model is $d=720$. The details of the
preparation of the training dataset and details of the training are given in
Appendix A.1. Here we present our main findings.

To compare different variational GP models and training strategies described
in the previous section, we built a collection of $3\times 3\times 3\times
2=54$ models corresponding to: three objective functions given by the ELBO
approximations by Eqs. (11), (12) and (13); three types of normal variational
distributions: one with learning the full covariance matrix, one with diagonal
covariance (mean-field approximation), and one with delta covariance, where
only the mean vector of the normal distribution is learned; three numbers of
inducing points of $100$, $500$ and $1000$ points; two regimes defining
whether the inducing point vector components are learnable parameters or not.
Each of the $54$ models was trained $3$ times to ensure the stability of
convergence. Empirically we found little dependence on the learning quality
and uncertainty estimates of varying the regularization parameter $\beta$
within the approximations of ELBO and we used $\beta=0.1$ in all the
experiments.

First, in Table 1, we show that SVGP framework augmented by variational
inducing points $\mathbf{Z}$ can successfully learn to represent different
types of molecules. We initialized the inducing points by random selection of
the training data points corresponding only to the first molecule and then
trained the SVGP model using all the training data points including also
second and third molecules. In Table 1 we compare two types of training
strategies: the first one is given by the fixed inducing points (FIP) where
the inducing points are initialized by the configurations of the first
molecule and are not learned during the training, and the second strategy is
given by the variational inducing points (VIP) where the inducing points are
learned during the training. In Table 1 we provide the best root mean square
error (RMSE) score among all the training runs for each number of inducing
points. Clearly, the variational learning of the inducing points significantly
improved the predictions for all three molecules. The gain in the prediction
performance by learning of the inducing points is greater for the low number
of inducing points, which could be simply because too few initial inducing
configurations failed to represent a large part of the training dataset; while
introducing an additional $100\times 720=72,000$ parameters within the model
makes it more expressive. Although the inducing points were initialized by the
configurations of the first molecule, the learned high-dimensional inducing
point vectors could represent well configurations of the second and third
molecules.

Figure 1: Structural isomers of C3H8O. | FIP/VIP
---|---
IP | $RMSE_{1}$ | $RMSE_{2}$ | $RMSE_{3}$
100 | 0.049/0.027 | 0.057/0.028 | 0.100/0.044
500 | 0.032/0.018 | 0.046/0.018 | 0.079/0.038
1000 | 0.025/0.018 | 0.039/0.018 | 0.081/0.037
Table 1:  Best root mean square error (RMSE) score among training runs for
test sets from different conformers (lower is better). IP stands for the
number of Inducing Points, FIP and VIP stand for the Fixed and Variational
Inducing Points, respectively. Inducing points were initialized by random
sampling from the configurations of the first molecule. Units are the
normalized potential energies.

(a) Figure 2:

Figures (a) and (c) show the share of the training runs with best uncertainty
estimates ($R_{\Sigma}<\epsilon_{R}$) for all three conformers, while figures
(b) and (d) show the share of the training runs with best RMSE scores
($RMSE_{\Sigma}<\epsilon_{RMSE}$) for all three conformers.

(b)

(c)

(d)

To check how well GP models can distinguish uncertainties on the train and
test datasets, we introduce a simple metric measuring the ratio of the mean
variance of the training subset of configurations $T$ to the mean variance of
the test subset of configurations $S$:

(14) $R(T,S)=\frac{|S|\sum_{i\in
T}\sigma\left(\mathbf{x}^{*}_{i}\right)^{2}}{|T|\sum_{j\in
S}\sigma\left(\mathbf{x}^{*}_{j}\right)^{2}}.$

By denoting the $T_{1}$, $T_{2}$, $T_{3}$ and $S_{1}$, $S_{2}$, and $S_{3}$
the training and test subsets of configurations corresponding to the first,
the second, and the third molecule, we computed the ratios
$R_{i}=R(T_{i},S_{i})$. We have chosen the test sets to contain the
continuations of the trajectories because, on the one hand, test
configurations are somewhat similar to the subset of training configurations,
therefore, making predictions on them makes sense, on the other hand, test
configurations should still be different from the training configurations
because they lie outside of the intervals of short MD trajectories, and,
therefore, measuring uncertainty is also reasonable for them.

In Figs. 2(a) (a) and (c), we provide the share of training runs with
corresponding training objectives and variational distributions that satisfied
the criterion $R_{\Sigma}=\sum_{i=1}^{3}R(T_{i},S_{i})<\epsilon_{R}$, with
$\epsilon_{R}=0.87\cdot 3=2.61$, corresponding approximately to the top $30\%$
of the score among all runs. In simple terms, we filtered training runs with
the best uncertainty ratios for all three conformers. In Figs. 2(a) (b) and
(d), we provide the share of training runs with best RMSE scores, i.e. runs
for which $RMSE_{\Sigma}=\sum_{i=1}^{3}RMSE_{i}<\epsilon_{RMSE}$, with
$\epsilon_{RMSE}=0.028\cdot 3=0.084$, also corresponding approximately to the
top $30\%$. The histograms of the $RMSE_{\Sigma}$ and $R_{\Sigma}$ scores for
variational objectives and variational distributions for all the training runs
are provided in Appendix, Fig. 6. While the predictive log-likelihood
objective provides the best uncertainty separation, it comes at the slight
expense of the predictive quality. Numerical experiments additionally showed
that mean-field variational distribution provides the best uncertainty
separation, while the full covariance variational distribution provides the
best predictive quality.

The improved uncertainty estimates by the predictive log-likelihood objective
can be justified theoretically by the following argument, as given in [24].
For all three objectives given by eq. (11), (12) and (13) the variance
parameter $\sigma$ is given by (6). The variational ELBO (11) and
$\gamma$-robust (12) objectives tend to underestimate the latent function
variance by attributing most of the uncertainty to the global observation
noise $\sigma^{2}_{\mathrm{obs}}$. This is because the data fit term in (11)
contains the term proportional to
$\frac{1}{\sigma_{\mathrm{obs}}^{2}}\left|y_{i}-\mu_{\mathbf{f}}\left(\mathbf{x}_{i}\right)\right|^{2}$,
while in the predictive log-likelihood objective the data fit term contains
the term proportional to
$\frac{1}{\sigma_{\mathrm{obs}}^{2}+\sigma_{\mathbf{f}}\left(\mathbf{x}_{i}\right)^{2}}\left|y_{i}-\mu_{\mathbf{f}}\left(\mathbf{x}_{i}\right)\right|^{2}$,
thus equilibrating the scales of the learned global and local contributions to
the variance $\sigma$.

###  3.2 Fitting atomic forces of molecular crystal

Figure 3: a) Machine learning scheme for modeling atomic forces in a molecular
crystal (b). A collection of 14 GP models was trained. Each GP model
correspond to atomic group with similar local atomic environments within a
molecular crystal. (b) Molecular crystal system with dynamic gearing parts.

(a) Absolute error for atomic forces on all atoms predicted by the set of SVGP
models, each parametrized with 500 inducing points and trained on 1000 points,
and for two exact GP models trained with 500 and 1000 points. (b) Distribution
of the normalized variance parameter learned by the SVGP and exact GP models.

(a)

(b)

Figure 4:

In this section, we provide numerical results aiming to illustrate the
performance of the SVGP framework for fitting the atomic forces and compare
the variational modeling framework with the exact GPs. Based on the
performance results from the example of molecular conformers, here we chose to
showcase the performance of the predictive log-likelihood objective for
fitting the atomic forces of a much larger molecular system. Next we provide
the numerical results for fitting the atomic forces of a large crystalline
system with dynamic gearing parts [25].

Based on the symmetries of local chemical environments within the molecular
crystal, we split the 528 atoms of the molecular crystal within a unit cell
into 14 groups and trained 14 SVGP and exact GP models as illustrated in Fig.
3. For each SVGP model, we initialized 500 inducing points by randomly
sampling a total of 1000 training points for each GP model. To compare the
performance of SVGP models trained with predictive log-likelihood objective
(13) and exact GP trained with objective (7), we additionally trained two sets
of exact GP models: with 500 and 1000 force samples. The feature vectors that
describe the local chemical environment of each atom and their derivatives
with respect to atomic positions were represented by a $1200$-dimensional SOAP
vectors [4], i.e. the input dimension for the model is $d=1200$. We used the
derivatives of the SOAP vectors of an atom with respect to the position of the
same atom as input feature vectors for the model. In this example, we used the
full covariance matrix to represent the variational distribution. The details
on the generation of the training dataset and the training procedure are given
in Appendix A.2.

The kernel function used in this example was the scaled Matérn kernel with
Automatic Relevance Determination (ARD) of features, defined by [52]

(15) $k_{\text{Matern
}}\left(\mathbf{x}_{1},\mathbf{x}_{2}\right)=\theta_{s}\frac{2^{1-\nu}}{\Gamma(\nu)}(\sqrt{2\nu}d^{\prime})^{\nu}K_{\nu}(\sqrt{2\nu}d^{\prime})$

where $\Gamma$ is a Gamma function, $K_{\nu}$ modified Bessel function of the
second kind,
$d^{\prime}=\frac{1}{2}\sum_{i=1}^{d}\frac{1}{\theta_{i}^{2}}\left(x_{i}-x_{i}^{\prime}\right)^{2}$
is a weighted distance between points, $\theta_{s}>0$ is a scale parameter,
$\theta_{i}>0$ are the lengthscale parameters, and in this example we used
$\nu=2.5$.

To train the SVGP and exact GP models, we generated the training dataset by
performing ab initio molecular dynamics simulation (simulation and training
details are given in Appendix A.2). To evaluate the performance of the SVGP
and exact GP models, as a test trajectory we took the trajectory and
corresponding atomic forces of correlated barrier crossing by all the 12 gears
within a unit cell as illustrated in Fig. 3 (b). This trajectory was used in
the paper [25] for computational estimates of the rotational barriers of
gears.

Model |  RMSE | Outliers |
---|---|---|---
GP: 500 points | 0.07822 | 4858 |
GP: 1000 points | 0.05045 | 1977 |
SVGP: 500 ind. points | 0.05008 | 1518 |
Table 2:  Performance of GP models for atomic force predictions of a molecular
crystal trained with variational and exact GPs. Outliers refer to the number
of test predictions with an absolute error greater than $0.8$
$\mathrm{e}\mathrm{V}$/$\mathrm{\SIUnitSymbolAngstrom}$.

In Fig. 4 we provide numerical results for three modeling settings: (1)
collection of 14 SVGP models, each with 500 inducing points and trained on
1000 points; (2) collection of 14 exact GP models trained on 1000 randomly
sampled points; (3) collection of 14 exact GP models trained on the 500
points. The inducing points of SVGP models were initialized by randomly
sampling 500 points out of 1000 training data points. The SVGP and exact GP
models were trained similarly with backpropagation by optimizing the objective
functions (7) and (13) respectively.

In Fig. 4 (a) we provide the absolute errors of the predicted atomic forces
for three modeling setups. Even though SVGP used 500 variationally learned
inducing points, the absolute error for the test dataset for the SVGP models
is on par with the exact GPs trained on 1000 points. This could be related to
the fact that each SVGP model has additional $500\times 1200=600,000$
variationally learned parameters allowing the SVGP models to find a more
optimal representation of the dataset. In Fig. 4 (b) we provide the histogram
of the learned GP variance parameter (normalized to be in $[0,1]$ interval)
for all atoms within a unit cell. Distributions of the variance parameters
differ significantly for the SVGP and exact GP models, with exact GP models
trained on 1000 points tending to have relatively higher normalized variance
than SVGP models. Additionally, in Table 2 we provide the RMSE scores for the
atomic force predictions and the number of outliers for the SVGP and exact GP
models. The SVGP models with 500 inducing points and trained on 1000 points
gave a similar RMSE score to the exact GP models trained on 1000 points,
however, with fewer outlier predictions.

###  3.3 Initialization of inducing points and stochasticity of the training

Although the training of the variational GP model has been cast in a black-box
manner, we still need to comment on possible pitfalls. One possible problem
while training SVGPs is the initialization of the inducing points
$\mathbf{Z}$. For example, we observed that optimization of the components of
high-dimensional inducing point vectors, initialized by random draws from the
normal distribution $\mathcal{N}(0,1)$, is typically stuck in some local
minimum and the trained model performs poorly. Therefore proper initialization
of inducing point vectors is needed. In this work, we initialized the inducing
points by randomly sampling from the training data points. A better
initialization strategy could be to sample the farthest points from the
training set [32] or use the $k$-means clustering algorithm. Even better
methods for choosing samples could also take into account the intrinsic
topology of the dataset [15] since descriptors of conformers of molecular
dynamics trajectory in principle follow a continuous path in a high-
dimensional descriptor space. Additionally, since training datasets are
frequently obtained by the MD simulations, the configurations are often
sampled non-uniformly. To mitigate this issue, alternative approaches have
been proposed for generating the training dataset [48, 49].

Another issue is whether inducing points can be learned by gradient descent,
and related to what kind of stochastic gradient descent algorithm is used for
iterations. In this paper, for both examples, we used full-batch training with
Adam optimizer [26] and multi-step gradient optimization with $2$ gradient
step sizes. The necessity of mini-batch training for generalization is
currently an open question: while some works [17] claim that generalization is
achievable in a full-batch regime, others [34] claim that mini-batch training
is necessary for generalization.

##  4 Interpolation vs. extrapolation

It is often tacitly assumed that simple non-variational GPs make predictions
on the test data points by interpolating within the training set of points,
which, for the low-dimensional datasets, is mostly an obvious statement. At
the same time, looking at the simulations of the MD of small molecules, like
the ones from Fig. 1, makes us think that conformational diversity during
dynamics is quite restricted – if some test snapshot is taken in between two
closely sampled in time neighboring snapshots, intuitively it “interpolates”
between its neighbors. Hence, if we want to predict some physical property,
like potential energy, for a collection of test configurations, we might be
lured into thinking that the GP generally will interpolate in between the
training configurations also in the descriptor space. By “in between,” it is
typically assumed that a test point will be within the convex hull of the
training points. Here we follow the arguments by [1], where it was shown that
for high-dimensional datasets, like the image data, learning always amounts to
extrapolation. Recently, a similar kind of analysis was performed on datasets
of bulk structures and crystals [53], where authors additionally propose an
alternative protocol to quantify the similarity of a test point to the
training set of points based on the local estimated density of training
points.

Next, we demonstrate, similarly to [1] and [53], that for the MD datasets with
high-dimensional descriptor representations like SOAP, the situation is quite
similar, even if we form the dataset by very similar molecular configurations,
the randomly sampled test configuration data point will most likely fall
outside the convex hull of the training set of points. This phenomenon is due
to the curse of dimensionality, it clearly shows that the meaning of
interpolation and extrapolation needs to be properly adjusted to the problem
at hand when working with high-dimensional data and alternative similarity
measures need to be constructed.

(a) Figure 5:

(a) Fraction of the total variance explained by the first $i$ singular
vectors. (b) Fraction of points from the test set within the convex hull of
the training points for random selections of $N_{dim}$ dimensions; the
standard deviations across the runs are shown with vertical bars.

(b)

We illustrate this phenomenon on the dataset of MD trajectories of the
molecules from Fig. 1. As before, for these molecules, we generated a large
set of conformations simply by simulating a $100$ $\mathrm{p}\mathrm{s}$ MD
trajectories at $300$ K temperature using the same parameters as in the
previous section. To get an estimate of the variability of the invariant
vectors within the dataset, we applied the Principal Component Analysis (PCA)
to the dataset containing invariant descriptor vectors for $10\,000$
consecutive snapshots. In Fig. 5(a) (a), we illustrate what fraction of the
total variance is explained by the first $i$ largest singular values
$\sigma_{j}$, i.e., we plot
$r_{i}=\sum_{j=0}^{i-1}\sigma_{j}^{2}/\sum_{j=0}^{d-1}\sigma_{j}^{2}$. For
example, the first $15$ principal components explain $96.5\%$, $94.7\%$ and
$97.1\%$ of total variance for the first, the second, and the third molecule.

Next, we show what fraction of test points falls within the convex hull of
training points for randomly selected subspaces in the descriptor space.
First, we formed a randomized collection of test vectors by randomly selecting
$200$ configurations from the pool of $10\,000$ configurations corresponding
to the $100$ $\mathrm{p}\mathrm{s}$ trajectory for each molecule. Then for a
collection of train and test vectors, we randomly subsampled $N_{dim}$
dimensions.

Constructing the convex hull for high-dimensional datapoints is
computationally expensive, however, checking hull membership for a test point
can be cast as a linear programming problem with positivity constraints. By
successfully solving the linear program we can check whether a test point can
be represented by a linear combination of the training points with positive
weights summing up to one.

In Fig. 5(a) (b), we illustrate what fraction of test points falls within the
convex hull of the training points for randomly sampled $N_{dim}$ dimensions
in the descriptor space. For each number of dimensions $N_{dim}$ we made $10$
runs with different random selections of test configurations and dimensions
for random projections. To conclude, we see that due to the curse of
dimensionality, in the high-dimensional molecular descriptor space, only a
small fraction of predictions will be performed in the interpolation regime.

##  5 Conclusion

In this work, we performed an empirical study of applying scalable variational
GPs for fitting the potential energies and atomic forces for two molecular
systems. First, we demonstrated that during the training, the SVGP models can
learn to represent different types of molecular configurations that were not
present within the initialization set of inducing points. Then we provided a
comparison of various alternatives for the SVGPs training, including different
variational objectives and variational distributions. Our experiments showed
that introducing and learning the latent inducing point variables in a high-
dimensional molecular descriptor space with a predictive log-likelihood
objective provided improved uncertainty estimates, however at a slight expense
of accuracy. Finally, we illustrated on the example of a large molecular
crystal system that a sparse SVGP model trained with a predictive log-
likelihood objective can provide predictions of atomic forces on par with the
exact GP models.

Although augmenting a GP model with a large number of variational parameters
clearly showed its advantages, it also made the training process more complex.
Moreover, even though the learned inducing points could represent different
types of configurations, their proper initialization was necessary, and the
extent to which inducing points can adapt to represent different molecular
configurations is still unexplored. In future work, we aim to address the
above issues as well as apply scalable variational GP models to large-scale
problems.

## References

  * [1] R. Balestriero, J. Pesenti, and Y. LeCun, Learning in high dimension always amounts to extrapolation, arXiv preprint arXiv:2110.09485, (2021).
  * [2] C. Bannwarth, S. Ehlert, and S. Grimme, GFN2-xTB—an accurate and broadly parametrized self-consistent tight-binding quantum chemical method with multipole electrostatics and density-dependent dispersion contributions, Journal of chemical theory and computation, 15 (2019), pp. 1652–1671.
  * [3] A. P. Bartók, J. Kermode, et al., Improved uncertainty quantification for gaussian process regression based interatomic potentials, arXiv preprint arXiv:2206.08744, (2022).
  * [4] A. P. Bartók, R. Kondor, and G. Csányi, On representing chemical environments, Physical Review B, 87 (2013), p. 184115.
  * [5] A. P. Bartók, M. C. Payne, R. Kondor, and G. Csányi, Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons, Physical review letters, 104 (2010), p. 136403.
  * [6] A. P. Bartók, M. C. Payne, R. Kondor, and G. Csányi, Gaussian approximation potentials: the accuracy of quantum mechanics, without the electrons, Phys. Rev. Lett., 104 (2010), p. 136403.
  * [7] P. E. Blöchl, Projector augmented-wave method, Physical review B, 50 (1994), p. 17953.
  * [8] M. Ceriotti, Beyond potentials: integrated machine-learning models for materials, arXiv preprint arXiv:2208.06139, (2022).
  * [9] M. Ceriotti, M. J. Willatt, and G. Csányi, Machine learning of atomic-scale properties based on physical principles, Handbook of Materials Modeling: Methods: Theory and Modeling, (2020), pp. 1911–1937.
  * [10] S. Chmiela, H. E. Sauceda, K.-R. Müller, and A. Tkatchenko, Towards exact molecular dynamics simulations with machine-learned force fields, Nature communications, 9 (2018), pp. 1–10.
  * [11] S. Chmiela, H. E. Sauceda, I. Poltavsky, K.-R. Müller, and A. Tkatchenko, sGDML: Constructing accurate and data efficient molecular force fields using machine learning, Computer Physics Communications, 240 (2019), pp. 38–45.
  * [12] S. Chmiela, A. Tkatchenko, H. E. Sauceda, I. Poltavsky, K. T. Schütt, and K.-R. Müller, Machine learning of accurate energy-conserving molecular force fields, Science advances, 3 (2017), p. e1603015.
  * [13] S. De, A. P. Bartók, G. Csányi, and M. Ceriotti, Comparing molecules and solids across structural and alchemical space, Physical Chemistry Chemical Physics, 18 (2016), pp. 13754–13769.
  * [14] V. L. Deringer, A. P. Bartók, N. Bernstein, D. M. Wilkins, M. Ceriotti, and G. Csányi, Gaussian process regression for materials and molecules, Chemical Reviews, 121 (2021), pp. 10073–10141.
  * [15] E. Dufresne, P. Edwards, H. Harrington, and J. Hauenstein, Sampling real algebraic varieties for topological data analysis, in 2019 18th IEEE International Conference On Machine Learning And Applications (ICMLA), IEEE, 2019, pp. 1531–1536.
  * [16] J. R. Gardner, G. Pleiss, D. Bindel, K. Q. Weinberger, and A. G. Wilson, GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration, in Advances in Neural Information Processing Systems, 2018.
  * [17] J. Geiping, M. Goldblum, P. E. Pope, M. Moeller, and T. Goldstein, Stochastic training is not necessary for generalization, arXiv preprint arXiv:2109.14119, (2021).
  * [18] A. Griewank and A. Walther, Evaluating derivatives: principles and techniques of algorithmic differentiation, SIAM, 2008.
  * [19] J. Hensman, N. Fusi, and N. D. Lawrence, Gaussian processes for big data, arXiv preprint arXiv:1309.6835, (2013).
  * [20] J. Hensman, A. Matthews, and Z. Ghahramani, Scalable variational gaussian process classification, in Artificial Intelligence and Statistics, PMLR, 2015, pp. 351–360.
  * [21] L. Himanen, M. O. J. Jäger, E. V. Morooka, F. Federici Canova, Y. S. Ranawat, D. Z. Gao, P. Rinke, and A. S. Foster, DScribe: Library of descriptors for machine learning in materials science, Computer Physics Communications, 247 (2020), p. 106949.
  * [22] M. D. Hoffman, D. M. Blei, C. Wang, and J. Paisley, Stochastic variational inference, Journal of Machine Learning Research, (2013).
  * [23] G. Imbalzano, Y. Zhuang, V. Kapil, K. Rossi, E. A. Engel, F. Grasselli, and M. Ceriotti, Uncertainty estimation for molecular dynamics and sampling, The Journal of Chemical Physics, 154 (2021), p. 074102.
  * [24] M. Jankowiak, G. Pleiss, and J. Gardner, Parametric gaussian process regressors, in International Conference on Machine Learning, PMLR, 2020, pp. 4702–4712.
  * [25] M. Jin, R. Kitsu, N. Hammyo, A. Sato-Tomita, M. Mizuno, A. S. Mikherdov, M. Tsitsvero, A. Lyalin, T. Taketsugu, and H. Ito, A steric-repulsion-driven clutch stack of triaryltriazines: Correlated molecular rotations and a thermoresponsive gearshift in the crystalline solid, Journal of the American Chemical Society, 145 (2023), pp. 27512–27520.
  * [26] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980, (2014).
  * [27] J. Knoblauch, Robust deep gaussian processes, arXiv preprint arXiv:1904.02303, (2019).
  * [28] J. Knoblauch, J. Jewson, and T. Damoulas, Generalized variational inference: Three arguments for deriving new posteriors, arXiv preprint arXiv:1904.02063, (2019).
  * [29] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Physical review B, 54 (1996), p. 11169.
  * [30] S. Maeda, Y. Harabuchi, M. Takagi, T. Taketsugu, and K. Morokuma, Artificial force induced reaction (afir) method for exploring quantum chemical potential energy surfaces, The Chemical Record, 16 (2016), pp. 2232–2248.
  * [31] A. G. d. G. Matthews, Scalable Gaussian process inference using variational methods, PhD thesis, University of Cambridge, 2017.
  * [32] C. Moenning and N. A. Dodgson, Fast marching farthest point sampling, tech. rep., University of Cambridge, Computer Laboratory, 2003.
  * [33] F. Musil, M. J. Willatt, M. A. Langovoy, and M. Ceriotti, Fast and accurate uncertainty estimation in chemical machine learning, Journal of chemical theory and computation, 15 (2019), pp. 906–915.
  * [34] S. W. Ober, C. E. Rasmussen, and M. van der Wilk, The promises and pitfalls of deep kernel learning, in Uncertainty in Artificial Intelligence, PMLR, 2021, pp. 1206–1216.
  * [35] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Physical review letters, 77 (1996), p. 3865.
  * [36] R. Ranganath, S. Gerrish, and D. Blei, Black box variational inference, in Artificial intelligence and statistics, PMLR, 2014, pp. 814–822.
  * [37] H. Robbins and S. Monro, A stochastic approximation method, The annals of mathematical statistics, (1951), pp. 400–407.
  * [38] L. Sbailò, Á. Fekete, L. M. Ghiringhelli, and M. Scheffler, The nomad artificial-intelligence toolkit: Turning materials-science data into knowledge and understanding, arXiv preprint arXiv:2205.15686, (2022).
  * [39] K. T. Schütt, S. Chmiela, O. A. von Lilienfeld, A. Tkatchenko, K. Tsuda, and K.-R. Müller, Machine learning meets quantum physics, Lecture Notes in Physics, (2020).
  * [40] G. N. Simm and M. Reiher, Error-controlled exploration of chemical reaction networks with gaussian processes, Journal of chemical theory and computation, 14 (2018), pp. 5238–5248.
  * [41] G. N. Simm, A. C. Vaucher, and M. Reiher, Exploration of reaction pathways and chemical transformation networks, The Journal of Physical Chemistry A, 123 (2018), pp. 385–399.
  * [42] E. Snelson and Z. Ghahramani, Sparse gaussian processes using pseudo-inputs, Advances in neural information processing systems, 18 (2005).
  * [43] M. Titsias, Variational learning of inducing variables in sparse gaussian processes, in Artificial intelligence and statistics, PMLR, 2009, pp. 567–574.
  * [44] M. Tsitsvero, Molecular dynamics trajectories of C3 H8 O molecule and its structural isomers, May 2022.
  * [45] M. Tsitsvero, J. Pirillo, Y. Hijikata, and T. Komatsuzaki, NMR spectrum prediction for dynamic molecules by machine learning: A case study of trefoil knot molecule, The Journal of Chemical Physics, 158 (2023).
  * [46] J. P. Unsleber, S. A. Grimmel, and M. Reiher, Chemoton 2.0: Autonomous exploration of chemical reaction networks, arXiv preprint arXiv:2202.13011, (2022).
  * [47] J. Vandermause, Y. Xie, J. S. Lim, C. J. Owen, and B. Kozinsky, Active learning of reactive bayesian force fields: Application to heterogeneous hydrogen-platinum catalysis dynamics, arXiv preprint arXiv:2106.01949, (2021).
  * [48] M. J. Waters and J. M. Rondinelli, Energy contour exploration with potentiostatic kinematics, Journal of Physics: Condensed Matter, 33 (2021), p. 445901.
  * [49] M. J. Waters and J. M. Rondinelli, Benchmarking structural evolution methods for training of machine learned interatomic potentials, Journal of Physics: Condensed Matter, 34 (2022), p. 385901.
  * [50] M. Wen and E. B. Tadmor, Uncertainty quantification in molecular simulations with dropout neural network potentials, npj Computational Materials, 6 (2020), pp. 1–10.
  * [51] M. J. Willatt, F. Musil, and M. Ceriotti, Atom-density representations for machine learning, The Journal of chemical physics, 150 (2019), p. 154110.
  * [52] C. K. Williams and C. E. Rasmussen, Gaussian processes for machine learning, vol. 2, MIT press Cambridge, MA, 2006.
  * [53] C. Zeni, A. Anelli, A. Glielmo, and K. Rossi, Exploring the robust extrapolation of high-dimensional machine learning potentials, Physical Review B, 105 (2022), p. 165141.

##  Appendix A Experiment details

###  A.1 Fitting energies of molecular isomers

For each molecule, we generated a short $20$ $\mathrm{p}\mathrm{s}$ MD
trajectory at $300$ $\mathrm{K}$ temperature with Berendsen thermostat and
GFN2-xTB tight-binding Hamiltonian [2]. Before simulation molecular geometries
were optimized with the same GFN2-xTB tight-binding Hamiltonian. The
equilibration period of $10$ $\mathrm{p}\mathrm{s}$ was used. All simulations
were performed in vacuum. Trajectories were calculated with $0.4$
$\mathrm{f}\mathrm{s}$ time step and then subsampled to $10$
$\mathrm{f}\mathrm{s}$ time step [44]. After subsampling each trajectory
consisted of $2000$ snapshots. We then selected the first $1600$
configurations within each trajectory to be included in a training dataset,
while the remaining $400$ configurations within each trajectory, were selected
for test sets. In total, the training set consists of $4800$ configurations
and there are three test sets with $400$ configurations each. For each
configuration of test/train sets, we computed the $720$-dimensional invariant
SOAP vector representation (3) with [21]. The invariant vectors were computed
with the following parameters: $n_{max}=5$, $l_{max}=5$, $r_{cut}=5.0$
$\mathrm{\SIUnitSymbolAngstrom}$ and $\sigma=0.5$
$\mathrm{\SIUnitSymbolAngstrom}$. Four atoms, namely, the three carbons and
the oxygen, were selected as centers for calculating the invariant vectors by
Eq. (3). For each configuration of MD trajectory, we extracted the potential
energy given by the GFN2-xTB tight-binding Hamiltonian [2]. The difference
between the maximum and minimum energies of configurations within the training
set (including all three molecules) was $0.8036$ $\mathrm{e}\mathrm{V}$, thus
we normalized the energies to fall into the interval $[0,1]$ scaling by this
number. In the figures and tables within the main text, we provide the error
metrics in the normalized units. Our goal was not to reach the ultimate
accuracy but to illustrate the concept, for the ultimate accuracy, one may
need to calculate the invariant descriptors with a greater number of radial
and spherical basis functions, design the optimal shape of the cutoff
function, choose optimal values of $\sigma$, as well as design more efficient
ways to combine the local invariants, other than summation (3) – a separate
set of problems. In Fig. 6 we provide the histograms of the RMSE and the
variance ratios for all training runs. Tables with scores for individual
isomers and used hyperparameters are provided in the supplementary material.

The estimated training time for all $54$ models was $54$ hours. All
variational models were trained with $20\,000$ gradient steps. The starting
learning rate was set to $0.01$ for the first $1000$ steps, then reduced to
$0.001$ until reaching the $10\,000$ steps, and finally reduced to $0.0001$
until reaching the $20\,000$ steps.

(a)

(b)

(c)

(d)

Figure 6:

Figures (a) and (c) show the distribution of the ratio score $R_{\Sigma}$ for
all training runs, while figures (b) and (d) provide the $RMSE_{\Sigma}$ score
for all training runs.

###  A.2 Fitting forces of molecular crystal

To cover the conformational diversity of the molecular crystal, we calculated
the ab initio molecular dynamics trajectory with $5000$ steps at temperature
$T=2000$ $\mathrm{K}$ with Langevin thermostat and a time-step of $1$
$\mathrm{f}\mathrm{s}$. The equilibration period of $2$ $\mathrm{p}\mathrm{s}$
was used. We sampled $1000$ configurations from the trajectory and computed
the SOAP vectors and their derivatives using the Librascal library [51]. The
parameters used to compute the SOAP vectors: $n_{max}=4$, $l_{max}=4$,
$r_{cut}=4.0$ $\mathrm{\SIUnitSymbolAngstrom}$ and $\sigma=0.3$
$\mathrm{\SIUnitSymbolAngstrom}$. The ab initio Density-functional-theory
(DFT) calculations were carried out using the Perdue-Burke-Ernzerhof (PBE)
exchange-correlation functional and the projector-augmented wave (PAW) method,
as implemented in the Vienna ab initio Simulation Package (VASP) [29, 7, 35].
A plane wave basis set with an energy cutoff of $500$ $\mathrm{e}\mathrm{V}$
was used. Gaussian smearing with a width of 0.05 eV was applied to treat
partial occupancies of the electronic orbitals, and an accuracy condition for
energy conversion of $10^{-6}$ eV was used. Only $\Gamma$-point calculations
were performed for the Brillouin zone sampling due to the large size of the
supercell.

The estimated training time for all 14 models was about 2 hours. Each model
was trained by $1000$ iterations with $0.1$ learning rate step, followed by
tuning $1000$ iterations with a reduced learning rate of $0.01$.

##  Appendix B Computational details

All the experiments were performed at NVIDIA RTX 8000 GPU. We used GPyTorch
library for modeling scalable Gaussian processes [16] and the DScribe [21] and
Librascal [51] libraries for computing SOAP invariant descriptors. The
“whitening” [31] transformation was used in both examples to accelerate the
training of variational parameters, as implemented within GyTorch library. The
models from both examples were trained with the Adam optimizer [26]. The
lengthscale and scale hyperparameters were initialized by standard
initialization procedures within GPyTorch library.