# Efficient interpolation of molecular properties across chemical compound
space with low-dimensional descriptors

Yun-Wen Mao  Roman V. Krems  Department of Chemistry, University of British
Columbia, Vancouver, B.C. V6T 1Z1, Canada
Stewart Blusson Quantum Matter Institute, Vancouver, B.C. V6T 1Z4, Canada

###### Abstract

We demonstrate accurate data-starved models of molecular properties for
interpolation in chemical compound spaces with low-dimensional descriptors.
Our starting point is based on three-dimensional, universal, physical
descriptors derived from the properties of the distributions of the
eigenvalues of Coulomb matrices. To account for the shape and composition of
molecules, we combine these descriptors with six-dimensional features informed
by the Gershgorin circle theorem. We use the nine-dimensional descriptors thus
obtained for Gaussian process regression based on kernels with variable
functional form, leading to extremely efficient, low-dimensional interpolation
models. The resulting models trained with 100 molecules are able to predict
the product of entropy and temperature ($S\times T$) and zero point
vibrational energy (ZPVE) with the absolute error under 1 kcal mol-1 for $>78$
% and under 1.3 kcal mol-1 for $>92$ % of molecules in the test data. The test
data comprises 20,000 molecules with complexity varying from three atoms to 29
atoms and the ranges of $S\times T$ and ZPVE covering 36 kcal mol-1 and 161
kcal mol-1, respectively. We also illustrate that the descriptors based on the
Gershgorin circle theorem yield more accurate models of molecular entropy than
those based on graph neural networks that explicitly account for the atomic
connectivity of molecules.

Suggested keywords

##  I  Introduction

The number of molecules with properties suitable for medicinal applications is
estimated to be between $10^{23}$ and $10^{60}$ polishchuk2013estimation .
This large number underscores the importance of the development of efficient
methods for predicting molecular properties for any given combination of
chemically bonded atoms. As demonstrated by recent work Weinreich2021 ;
hansen2015 , prediction of molecular properties can be achieved by
interpolation in chemical compound space, defined as an abstract space with
individual molecules as elements of the space, by machine learning (ML)
models. A chemical space is often, although not always, restricted to a
particular number and combination of atoms (for example, a space of molecules
with chemical formula C7H10O2). Interpolation of molecular properties across
chemical space meets with four main challenges: (i) the vast number of
molecules; (ii) the requirement to represent molecules by numerical
descriptors that are universal, efficient and accurate; (iii) the sensitivity
of molecular properties to small structural changes in atomic arrangements;
(iv) lack of structure in chemical compound space. While interpolation of
molecular properties in chemical spaces is often considered as a big data
problem, requiring information from large sets of molecules, the focus of many
recent studies has been on data-efficient models Browning2017 ; Gubaev2018 ;
FCHL19 ; Bogojeski2020 . The efficiency of ML predictions of molecular
properties can be enhanced by building underlying symmetries into ML models
Gastegger2018 ; Rostami2018 ; Qiao2020 , by improving the molecular
descriptors musil2021physics ; Gallegos2021 ; Willatt_2018 ; Langer2022 ;
FCHL19 ; Faber2017 ; Faber2018  and by improving the ML models through kernel
optimization for kernel models duvenaud2013structure ; duvenaud2011additive ;
Deng2020  or optimization of neural network (NN) architecture and algorithm
hyperparameters for NN-based models Gilmer2017 ; schnet_2018 ; Smith2019 .

Key to any ML predictions of molecular properties are molecular descriptors
that serve as inputs into supervised learning models. The molecular
descriptors must be numerical, universal (i.e. describe any molecule), unique,
and embody all of the relevant features, such as the shape of molecules, the
arrangement of chemical bonds, bond strengths and bond order, vibrational
frequencies and conformational degrees of freedom. Many variations of
molecular descriptors have recently been proposed, including Coulomb matrices
PRL_CM_2011 , bag of bonds (BoB) features BoB_2015 , smooth overlap of atomic
positions (SOAP) SOAP_2012 , Faber-Christensen-Huang-Lilienfeld (FCHL19)
representation FCHL19 . Another possibility is to represent molecules by a
Graph Neural Network (GNN) kipf2016semi ; Wu2019  that explicitly accounts for
the spatial arrangement and connectivity of atoms in a molecule. GNN can be
combined with other features embodying molecular properties, yielding accurate
molecular descriptors Zhu2022 ; Yang2019 ; 9338432 ; Qiao2020 . However, for
complex molecules, these descriptors are necessarily high-dimensional. While
high-dimensional descriptors can be used as inputs into isotropic ML models
(for example ML models with kernels that do not discriminate between different
input dimensions) williams2006gaussian , such models do not account for
inherent differences in scales of different features. Willatt et. al.
Willatt_2018  showed that reducing the dimensionality of descriptors may
improve the accuracy of ML predictions. In addition, Musil et. al.
musil2021physics  argued that descriptors reflecting some fundamental
principles yield more robust, transferable and data-efficient ML models. In
addition to efficiency, low-dimensional descriptors enable the possibility of
Bayesian optimization of molecular properties in chemical compound space. The
present work demonstrates that it is possible to design low-dimensional
descriptors of large polyatomic molecules for extremely efficient
interpolation of entropy and zero point vibrational energy in chemical
compound space.

Different molecular properties exhibit vastly different variation across
chemical compound spaces. Therefore, it can be argued that the most data-
efficient predictions must be based on models and descriptors tailored for the
specific molecular property. In the present work, we test this conjecture by
first considering interpolation of molecular entropy across chemical compound
spaces and developing low-dimensional descriptors tailored for the prediction
of entropy. Our starting point is based on three-dimensional, universal,
physical descriptors derived from the properties of the distributions of the
eigenvalues of the Coulomb matrix. We show that the largest eigenvalue and the
two leading moments of the eigenvalue distribution provide the physical
descriptors for molecular entropy, as elucidated using the Gershgorin circle
theorem gershgorin1931uber ; saad2011numerical_gershgorin . We illustrate that
the accuracy of these descriptors can be enhanced by combining them with
either: (1) the leading principle components of the output of GNN tailored for
the prediction of entropy; or (2) with descriptors based on the Gershgorin
circle theorem. We show that the models of entropy using descriptors based on
the Gershgorin circle theorem perform better than the models using descriptors
based on GNN.

By using the resulting nine-dimensional descriptors in Gaussian process
regression with kernels optimized for pattern recognition
duvenaud2013structure ; duvenaud2011additive , we build ML models capable of
predicting entropy of 133,000 organic molecules with up to 29 atoms (HCNOF)
using as few as 100 molecules for training. The predictions have chemical
accuracy (with predicted entropy $\times$ temperature under 1 kcal mol-1),
illustrating, to the best of our knowledge, the most data-efficient models of
entropy to date. We then illustrate that these descriptors, even though
tailored for models of entropy, can be used to interpolate zero point
vibrational energy (ZPVE). The variation of ZPVE in the chemical space
considered here is much greater than the scale of values covered by entropy.
However, we show that the nine-dimensional descriptors can be used to produce
extremely data-starved yet accurate models of ZPVE, provided the kernels of
the underlying models are properly optimized. We illustrated that models
trained by only 100 molecules produce ZPVE with mean absolute error and
absolute error for $>92$ % of predictions under 1.3 kcal mol-1, in a chemical
compound space including molecules with up to 29 atoms and covering the ZPVE
range of 161 kcal mol-1

##  II Dataset summary

In this work, we focus on molecules from the QM9 dataset qm9_dataset , where a
chemical subspace consists of 133,885 stable organic molecules made up of a
combination of the following atoms: CHONF. Within the QM9 dataset, there are
526 subsets of constitutional isomers, where constitutional isomers are
molecules with the same chemical formula but different atom connectivity. Each
subset includes at least 2 molecules, while the largest constitutional isomer
subset includes 6059 molecules.

We reserve 20 $\%$ of these molecules for the hold-out test set, while samples
from the remaining distributions of molecules are used to train the ML models.
Since molecules in this chemical subspace range from 3 to 29 atoms in size, we
bin all molecules according to their size before randomly selecting the
distributions for the training and test sets. This ensures the ensuing ML
models sample the entire chemical subspace, both at the training and test
stage. Furthermore, 20 $\%$ of the training dataset is used as a model
validation set. The validation set is used for optimizing the hyperparameters
of the ML models.

The prediction accuracy of the ML models is quantified, for both the test and
validation sets, by the root-mean-square error (RMSE)

${\rm RMSE}=\sqrt{\frac{\Sigma_{i=1}^{N}(y_{i}-\hat{y_{i}})^{2}}{N}}$ (1)

or mean absolute error (MAE)

${\rm MAE}=\frac{1}{N}\sum_{i=1}^{N}|y_{i}-\hat{y_{i}}|$ (2)

evaluated over the number of samples $N$ in the test or validation set, as
appropriate. Here, $y_{i}$ and $\hat{y_{i}}$ are the observed entropy and ML
predicted entropy of the $i^{th}$ molecule in the test/validation set,
comprising $N$ molecules. Lower RMSE and MAE values indicate better accuracy
of ML predictions.

##  III Method

The main goal of the present work is to develop low-dimensional, general
descriptors of molecules that can be applied to both small and large molecules
and used for efficient interpolation of molecular properties across chemical
compound spaces with a great variety of molecules. We build these descriptors
incrementally, by starting with the properties of distributions describing the
Coulomb matrices PRL_CM_2011 , considering the enhancement of the descriptors
by Graph Neural Networks, and using physical insights offered by the
Gershgorin circle theorem Bellman+2021 ; saad2011numerical_gershgorin . Once
identified, the best low-dimensional descriptors are used as features of
inputs into Gaussian process models with kernels optimized for a specific
molecular property. The efficiency of the final models presented here is thus
determined by both the efficiency of the molecular descriptors and the
alignment of model kernels with the prediction targets.

###  III.1 Descriptor I: Distributions of eigenvalues of Coulomb matrices

Coulomb matrix is a molecular descriptor developed by Rupp et al. PRL_CM_2011
. Each molecule is represented by a $n_{\rm max}\times n_{\rm max}$ symmetric
matrix, where $n_{\rm max}$ is the number of atoms in the largest molecule in
the chemical subspace. Each matrix element $M_{ij}$ in matrix $M$ is defined
as

$M_{ij}=\begin{cases}0.5Z_{i}^{2.4}&\text{for $i=j$}\\\
\frac{Z_{i}Z_{j}}{R_{ij}}&\text{for $i\neq j$}\end{cases}$ (3)

where $Z_{i}$ and $Z_{j}$ are the atomic numbers of atoms $i$ and $j$,
respectively, and $R_{ij}$ is the separation between atoms $i$ and $j$. The
diagonal elements of the Coulomb matrix represent the polynomial fit relating
the atomic number to the total energies of the free atoms PRL_CM_2011 . The
off-diagonal elements describe the electrostatic interaction between atoms $i$
and $j$. If the size of a molecule is smaller than $n_{\rm max}$, the matrix
elements with $i$ or $j>n_{\rm max}$ are set to zero. The Coulomb matrix is
invariant to molecular translations and rotations, but not invariant to
permutations of the atoms. This problem can be overcome by diagonalizing the
Coulomb matrix PRL_CM_2011 ; schrier2020can .

In the present work, we begin by implementing the basic three-dimensional (3D)
descriptors of molecules using the largest eigenvalue $\epsilon_{\rm max}$,
the mean $\mu(\epsilon>1)$ and the standard deviation $\sigma(\epsilon>1)$ of
the distribution of the eigenvalues of the Coulomb matrix ($\epsilon$) with
values greater than 1. These 3D descriptors are used as inputs into the most
basic ML models of entropy in this work. To enhance the expressiveness of the
molecular descriptors and the accuracy of the ensuing models, we supplement
these 3D vectors with features based on either the outputs of Graph Neural
Networks, optimized to predict entropy as described in Descriptor II, or
features based on the Gershgorin circle theorem as described in Descriptor
III.

###  III.2 Descriptor II: Graph Convolution Networks

To enhance the accuracy of the resulting ML models as well as to explore the
role of explicit graph descriptors for interpolation of molecular properties,
we supplement the 3D descriptors derived from the Coulomb matrix by
descriptors based on the outputs of graph convolution networks (GCN). A
molecule can be mapped onto a graph $G=(\nu$, $e$$)$, where the atoms and
bonds are treated as nodes ($\nu$) and edges ($e$), respectively. Each graph
is represented by an adjacency matrix $A$ with matrix elements defined as

$A_{ij}=\begin{cases}1&\text{if $e_{ij}\in\varepsilon$}\\\
0&\text{otherwise,}\end{cases}$ (4)

where $e_{ij}$ is the edge between the $i^{th}$ and $j^{th}$ node. With the
adjacency matrix thus defined, the $(\lambda+1)^{\rm th}$ layer of a multi-
layer GCN is implemented using

$H^{\lambda+1}=\phi(\hat{D}^{-\frac{1}{2}}\hat{A}\hat{D}^{-\frac{1}{2}}H^{\lambda}W^{\lambda}),$
(5)

where $\hat{D}$ is the diagonal matrix with elements

$\hat{D}_{ii}=\sum_{j-1}^{m_{\rm max}}\hat{A}_{ij},$ (6)

$\hat{A}=A+I_{m_{\rm max}}$, $I_{m_{\rm max}}$ is $m_{\rm max}$-dimensional
identity matrix, $m_{\rm max}$ is the maximum number of heavy atoms (CONF)
that a molecule has for the dataset, $\phi$ is the activation function, and
$W^{\lambda}$ is the matrix of weights for layer $\lambda\geq 1$, producing
molecular descriptors illustrated in Fig. 1. The weight matrices are optimized
by training the GCN using the algorithm illustrated in Fig. 2.

The GCN are trained using the training and validation sets of molecules. For
$\lambda=1$, $H^{1}=X$ where $X$ are the node features in the one hot encoded
data format. The feature of each node is a vector including atomic number,
hybridization, degree, and chirality of each atom. Each node of the graph
aggregates information from its neighbouring node(s) as the graph passes
through layers of the GCN by Eq. (5). The weight matrices are initialized
randomly with a truncated normal distribution. The weight matrix of the last
GCN layer is set as a one-dimensional vector, so that the resulting final
layer $H^{\lambda+1}_{i}$ is a vector instead of a matrix. In this work, we
build GCN with $3+1$ layers, where $W_{1}$ is a $20\times 9$ matrix, $W_{2}$
is a $9\times 9$ matrix, and $W_{3}$ is a $9\times 1$ matrix.

Figure 1:  Schematic diagram illustrating the use of GCN to create molecular
descriptor ($H^{4}_{\rm{C_{2}H_{4}O}}$) with C2H4O as an example. The GCN used
in this work, as demonstrated in this figure, has four layers ($\lambda+1=4$).
The resulting molecular descriptor is nine dimensional, because the largest
molecule in QM9 dataset considered here has 9 heavy atoms.

The resulting $H^{\lambda+1}_{i}$ vector is treated as input into a Gaussian
process regression (GPR) model, with entropy of the corresponding molecule as
output. The RMSE over the validation set is used as the loss function for the
optimization of the GCN weight matrices. The outputs of GCN are converted into
a 9D vector by a densely connected NN layer, and the 9D vector is further
reduced to $k$ dimensions by principal component analysis (PCA), yielding GCN-
PC($k$) descriptors.

In the present work, the GCN model is trained with the cluster of
constitutional isomers with the chemical formula C7H10O2 from Ref. qm9_dataset
. The trained weight matrices are then used to compute the adjacency matrices
for other molecules in the QM9 database of Ref. qm9_dataset  to generate the
GCN-based descriptors. For smaller molecules with less than nine heavy atoms,
the adjacency matrix size is still $9\times 9$.

Figure 2:  Schematic diagram of the algorithm used in the present work for
training GCN by optimizing weight matrices ($W_{1}$, $W_{2}$, $W_{3}$). Each
molecule in A is processed as described in Fig. 1 to yield a unique GCN vector
($H^{4}_{\rm{molecule}}$), used as the molecular descriptor. A subset of
molecules in A, yielding a set of vectors depicted in C, is used to train a
Gaussian process (GP) model of entropy. The RMSE of the GP predicted values is
used in step E as the loss function to optimize the weight matrices ($W_{1}$,
$W_{2}$, $W_{3}$) of GCN by gradient descent.

###  III.3 Descriptor III: Gershgorin circle theorem

The Gershgorin circle theorem identifies a region in the complex plane that
contains all the eigenvalues of a complex square matrix Bellman+2021 ;
saad2011numerical_gershgorin . The theorem states that for an $m\times m$
matrix with the entries in complex plane $C$, the eigenvalues of matrix $M$
are in the range of $D_{1}$ $\cup$ $D_{2}$ $\cup$ …$\cup$ $D_{m}$. Each disc
$D_{i}$ is defined as

$D_{i}=\\{z\in C:|{z-M_{ii}}|<\sum_{j\neq i}|M_{ij}|\\}.$ (7)

Schrier schrier2020can  analyzed the eigenvalues of Coulomb matrices using the
Gershgorin circle theorem and showed that a highly substituted carbon in an
alkane has larger off-diagonal values, leading to higher upper bound for
$D_{i}$ and consequently higher eigenvalues. We use the results of the
Gershgorin circle theorem to define new descriptors of molecules, as follows.

Each Gershgorin circle disc ($D_{i}$) of the coulomb matrix ($M$) is used to
define a normal distribution function $f_{i}$

$f_{i}(x|\mu,\sigma^{2})=\frac{1}{\sigma}\phi\left(\frac{x-\mu}{\sigma}\right),$
(8)

where $\mu$ and $\sigma$ are the mean and the standard deviation of the normal
distribution, and $\phi(\cdot)$ is the probability density function (PDF)
defined as

$\phi(z)=\frac{e^{-z^{2}/2}}{\sqrt{2\pi}}.$ (9)

The mean of the normal distribution is defined with the value of the center of
each Gershgorin circle disc, $M_{ii}$. The standard deviation of the
distribution is determined by the radius of each Gershgorin circle disc,
$(\sum_{j\neq i}|M_{ij}|)^{\tau}$. Each molecule is then described by $f_{\rm
molecule}$, defined as the sum of all $f_{i}$

$f_{\rm molecule}=\sum_{i}^{n_{\rm max}}d_{i}\times
f_{i}(x|\mu=M_{ii},\sigma^{2}=(\sum_{j\neq i}|M_{ij}|)^{\tau}).$ (10)

where $d_{i}$ and $\tau$ are hyperparameters that allow adjustment on the
weight on and standard deviation of $f_{i}$, respectively. In this work, we
fix $d_{i}$ and $\tau$ to 1.

We further define an atomic reference PDF ($f_{\rm atom}$),

$f_{\rm atom}=d_{\rm atom}\times f_{i}(x|\mu=M_{\rm atom},\sigma^{2}_{\rm
atom}),$ (11)

where $M_{\rm atom}=0.5Z_{\rm atom}^{2.4}$, and $d_{\rm atom}$ is a weight
constant that is set to 10. By calculating the inner products between the
atomic reference PDF and the individual molecular PDF, we reduce the high
dimensional continuous PDF curve into a five dimensional descriptor, $\langle
f_{\rm HCNOF},f_{\rm molecule}\rangle$. We also use the area under each
$f_{\rm molecule}$ denoted hereafter by ${\rm AUC}(f_{\rm molecule})$ as an
additional molecular descriptor. The relative role of these descriptors is
illustrated in the following section.

Figure 3 shows $f_{\rm{molecule}}$ of two $\rm C_{7}H_{10}O_{2}$
constitutional isomers, where the blue curve is an isomer with higher entropy
and the orange curve is that with lower entropy. The variances ($\sigma^{2}$)
of $f_{i}$ are $\sum_{j\neq i}|M_{ij}|$. The atomic reference curves are shown
by black dotted lines.

Figure 3:  Left: Molecular probability density function for two isomers of
${\rm C_{7}H_{10}O_{2}}$ based on the Gershgorin circle theorem. The black
doted curves show the atomic reference curves (HCNOF). Right: Molecular
probability density function for an isomer of ${\rm C_{7}H_{10}O_{2}}$ and the
reference molecule based on the Gershgorin circle theorem. For all $f_{\rm
C_{7}H_{10}O_{2}}$ and $f_{\rm ref}$, $d_{i}=M_{ii}$, and $\tau=1$.

###  III.4 ML method: Gaussian Process Regression with variable kernels

Gaussian process regression (GPR) murphy2018machine  is a supervised ML method
that models an unknown target function, $f$, given input data $X$ and output
data $y$. The training dataset comprises
$\mathcal{D}=\\{(x_{i},f_{i}),i\in[1,N_{\rm train}]\\}$, where
$f_{i}=f(x_{i})$ is the noise-free observation of the function evaluated at
$x_{i}$. GPR infers a distribution over functions with $p(f|\mathcal{D})$. In
the context of property (e.g. entropy) interpolation across chemical space, a
probability distribution of all possible entropy values for a molecule is
inferred from the training dataset $\mathcal{D}$. The inputs $x_{i}$ are the
molecular descriptors and $f_{i}$ is the corresponding value of molecular
entropy. GPR defines a joint distribution to predict $f_{*}$ on the unseen
inputs $(X_{*})$ as

$\begin{pmatrix}f\\\ f_{*}\end{pmatrix}\sim N\left(\begin{pmatrix}\mu\\\
\mu_{*}\end{pmatrix},\begin{pmatrix}K&K_{*}\\\
K_{*}^{T}&K_{**}\end{pmatrix}\right)$ (12)

where $\mu=\left[m(x_{1}),...,m(x_{N_{\rm train}})\right]$, with $m(\cdot)$
representing the mean of a Gaussian distribution, and $K=\kappa(X,X)$,
$K_{*}=\kappa(X,X_{*})$, $K_{**}=\kappa(X_{*},X_{*})$ determine the covariance
matrix, defined by the kernel function $\kappa$. The efficiency of the GP
model is determined by the functional form of the kernel function.

In the present work, we follow Refs.duvenaud2011additive ;
duvenaud2013structure  to construct kernel functions with variable complexity
by forming linear combinations and products of the following basis kernel
functions:

$\kappa_{\rm{DP}}(x_{i},x_{j})=\sigma_{0}^{2}+\langle x_{i},x_{j}\rangle$ (13)

where $\sigma_{0}$ is a parameter which controls the inhomogenity of the
kernel,

$\kappa_{\rm{RBF}}(x_{i},x_{j})=\exp\left(-\frac{d(x_{i},x_{j})^{2}}{2l^{2}}\right)$
(14)

where $l$ is the length scale and $d(\cdot,\cdot)$ is the Euclidean distance
between two data points,

$\kappa_{\rm{RQ}}(x_{i},x_{j})=\left(1+\frac{d(x_{i},x_{j})^{2}}{2\alpha
l^{2}}\right)^{2}$ (15)

where $\alpha$ is the scale mixture parameter and $l$ is the length scale of
the kernel,

$\kappa_{\rm{MAT}}(x_{i},x_{j})=\frac{1}{\Gamma(\nu)2^{\nu-1}}\left(\frac{\sqrt{2\nu}}{l}d(x_{i},x_{j})\right)^{\nu}K_{\nu}\left(\frac{\sqrt{2\nu}}{l}d(x_{i},x_{j})\right)$
(16)

where $\nu$ and $l$ are positive parameters, $K_{\nu}(\cdot)$ is a modified
Bessel function and $\Gamma(\cdot)$ is the gamma function, as well as

$\kappa_{\rm{Periodic}}(x_{i},x_{j})=\exp\left(-\frac{2\sin^{2}(\pi
d(x_{i},x_{j}))/p}{l^{2}}\right)$ (17)

where $l$ is the length scale of the kernel, and $p$ is the periodicity of the
kernel.

We optimize the functional form of the kernel function using a greedy search
algorithm that combines simple functions into linear combinations and
products, as was described in Refs. duvenaud2011additive ;
duvenaud2013structure  and used for various applications in Refs.
torabian2023compositional ; vargas2020physical ; asnaashari2021gradient ;
dai2023neural ; dai2020interpolation ; vargas2018extrapolating . The kernel
construction algorithm is illustrated in Fig. 4.

Figure 4:  Schematic diagram of the kernel construction algorithm used in this
work to increase the complexity of the kernel function for optimal kernel
models.

All possible linear combinations and products of base kernels defined in Eqs.
(13) - (17) are created to yield more complex kernels. The model is selected
by the Bayesian information criterion (BIC) calculated as

${\rm BIC}(\mathcal{\kappa}_{i})=\log
P(y|{\kappa}_{i})-\frac{1}{2}|{\kappa}_{i}|\log(N_{\rm train})$ (18)

where $|\kappa_{i}|$ is the number of parameters of kernel $\kappa_{i}$, and
$\log P(y|{\kappa}_{i})$ is the maximum of log-likelihood marginalized over
the parameters of the functions in the Gaussian process, but not the kernel
hyperparameters. This is the same function that is maximized to train a
Gaussian process.

##  IV Results

###  IV.1 Analysis of different low-dimensional descriptors

Figure 5:  Test MAE of GP regression of entropy trained with different
combinations of molecular descriptors: circles – 3D CM($\epsilon_{\rm max}$,
$\mu(\epsilon>1)$, $\sigma(\epsilon>1)$) with $\rm{RQ+RQ}$ kernel; filled-
diamonds – CM + GCN-PC($k=5$) with $\rm{RQ+RQ}$ kernel; empty diamonds – CM +
GCN-PC($k=9$) with $\rm{RQ+RQ}$ kernel; up-triangles – the 4D
$\rm{CM}+\rm{AUC}[f_{\rm{molecule}}]$ with $\rm{RQ+RQ+MAT}$ kernel; down-
triangles – the 9D $\rm{CM}+\rm{AUC}[f_{\rm{molecule}}]+\langle f_{\rm
HCNOF},f_{\rm molecule}\rangle$ descriptor trained with $\rm{RQ+MAT\times DP}$
kernel. The kernel of each model is optimized as described in text with BIC as
the model selection metric. The lines represent the average of ten sets of
training using different distributions of molecules in the training set, and
the shaded area is the standard deviation of the results.

In order to analyze the relative performance of the different molecular
descriptors, we begin by building GP models of entropy with various
combinations of molecular descriptors described in the previous section.
Figure 5 compares the learning curves of GP models thus obtained. As shown in
Fig. 5, the best performing GP model is obtained with nine dimensional CM-GC
descriptor, which includes three descriptors stemming from the Coulomb matrix
denoted as CM($\epsilon_{\rm max}$, $\mu(\epsilon>1)$, $\sigma(\epsilon>1)$),
the area under the curve ${\rm AUC}(f_{\rm molecule})$, and five inner
products $\langle f_{\rm HCNOF},f_{\rm molecule}\rangle$. The figure shows
that models trained with CM-GC descriptor yield the most accurate
interpolation, reducing the MAE to below 1 $\rm{kcal~{}mol^{-1}}$ with less
than 100 molecules. This figure illustrates that physical descriptors derived
from combinations of CM-based quantities and properties based on the
Gershgorin circle theorem improve the accuracy of ensuing models much more
significantly than the descriptors based on graph neural networks.

To understand the role of the descriptors contributing to the best performing
combination in Figure 5, we examine correlations of molecular entropy with the
largest eigenvalue of the Coulomb matrix and the standard deviation $\sigma$
of the distribution of the eigenvalues $\epsilon>1$, both shown in Fig. (6).
While neither $\epsilon_{\rm max}$ nor $\mu$ are unique features for entropy,
the results in Figs. 6 illustrate that there exist strong correlations between
entropy and both of $\epsilon_{\rm max}$ and $\sigma$. These correlation can
be rationalized using the Gershgorin circle theorem.

Figure 6:  Entropy as a function of the largest eigenvalue ($\epsilon_{\rm
max}$) of the Coulomb matrix (left panel) and the standard deviation
($\sigma$) of the Coulomb matrix eigenvalues greater than 1 (right panel) for
constitutional isomers with the formula C7H10O2.

The Gershgorin circle theorem identifies a region in the complex plane that
contains all the eigenvalues of a complex square matrix Bellman+2021 ;
saad2011numerical_gershgorin . The theorem states that for an $m\times m$
matrix with the entries in complex plane $C$, the eigenvalues of matrix $M$
are in the range of $D_{1}$ $\cup$ $D_{2}$ $\cup$ …$\cup$ $D_{m}$. Each disc
$D_{i}$ is defined in Eq. 7. Schrier schrier2020can  analyzed the eigenvalues
of Coulomb matrices using the Gershgorin circle theorem and showed that a
highly substituted carbon in an alkane has larger off-diagonal values, leading
to higher upper bound for $D_{i}$ and consequently higher eigenvalues.

Figure 7:  Visual demonstration of the Gershgorin circle theorem on the
Coulomb matrices of two C7H10O2 constitutional isomers
1-hydroxyhept-5-yn-3-one (A, top) and 2,4-dioxatricyclo[$4.3.0.0^{3,8}$]nonane
(B, bottom), where the molecule A has the higher entropy. The entropy values
are in the unit of kcal mol-1.The Gershgorin discs of each isomer’s Coulomb
matrix are drawn as orange circles in the figures. The blue triangles are all
the eigenvalues of the Coulomb matrices. The red arrows emphasis the CM-EV1 of
each isomer’s Coulomb matrix.

To elucidate the connection with entropy, consider two constitutional isomers
with molecular formula C7H10O2 that have drastically different entropy. The
molecules with the highest and lowest entropy in the C7H10O2 isomeric cluster
are 1-hydroxyhept-5-yn-3-one (Fig. 7.A) and
2,4-dioxatricyclo[$4.3.0.0^{3,8}$]nonane (Fig. 7.B), respectively. For
molecule A: the entropy is $1.130\times 10^{-1}$ kcal mol-1K-1, $\epsilon_{\rm
max}=175.437$, $\sigma(\epsilon>1)=51.800$ and $\mu(\epsilon>1)=45.987$. For
molecule B: the entropy is $7.668\times 10^{-2}$ kcal mol-1K-1 (lower),
$\epsilon_{\rm max}=212.320$ (higher), $\sigma(\epsilon>1)=61.221$ (higher),
while $\mu(\epsilon>1)=45.980$ is similar to that of 1-hydroxyhept-5-yn-3-one.
Molecule B is more rigid, which leads to fewer microstates and lower entropy.

Atom 1 |  Atom 2 |  $Z_{1}Z_{2}$ |  $R$ (Å) |  $M_{ij}$
---|---|---|---|---
H |  H |  1 |  0.74 |  1.351
C |  H |  6 |  1.09 |  5.505
O |  H |  8 |  0.971 |  8.239
C |  C |  36 |  1.535 (single) |  23.453
|  |  |  1.339 (double) |  26.886
|  |  |  1.203 (triple) |  29.925
C |  O |  48 |  1.431 (single, ethanol) |  33.543
|  |  |  1.219 (double, ketone) |  39.377
Table 1: Off-diagonal entries ($M_{ij}$) of the Coulomb matrix for several
bonded atoms.

The more rigid molecule, 2,4-dioxatricyclo [$4.3.0.0^{3,8}$]nonane (Fig. 7B.),
has a lower entropy. This is the manifestation of the fact that entropy is
determined by the shape and rigidity of the molecule huyskens1989influence .
Rigid molecules tend to be restricted in motion, which leads to lower numbers
of available microstates and lower entropy. Thus, a good representation of
molecular rigidity should be a good descriptor for predicting entropy. Given
this premise, it is useful to analyze how the CM($\epsilon_{\rm max}$,
$\mu(\epsilon>1)$, $\sigma(\epsilon>1)$) descriptor reflects the rigidity of
molecules. The diagonal entries of CM (Table 1) for molecules considered here
are $M_{ii,{\rm H}}=0.5$, $M_{ii,{\rm C}}=36.858$, or $M_{ii,{\rm O}}=73.517$,
corresponding to hydrogen, carbon and oxygen. The off-diagonal elements
$M_{ij}$ are directly proportional to the atomic numbers $Z_{i}$ and $Z_{j}$,
and inversely proportional to the interatomic separation $R_{ij}$. This means
that the $M_{ij}$ value is larger for bonded pairs of heavy atoms. Table 1
lists some examples of $M_{ij}$ values, where the C$=$O double bond pair gives
the highest $M_{ij}$ value. For rigid molecules, the atoms are spatially
crowded, which leads to higher $M_{ij}$ values on average. We observe that the
Gershgorin circle theorem yields higher values of $|A_{ij}|$ in Eq. (7) for
rigid molecules on average, leading to larger values of $\sum_{j\neq
i}|A_{ij}|$. Thus, the range of eigenvalues defined by the union of $D_{i}$ in
the circle theorem is higher for rigid molecules. The opposite is true for
linear and less rigid molecules.

Fig. 7 visualizes the results of the Gershgorin circle theorem. The $y$\- and
$x$-axes are the imaginary and real axes of the complex plane containing the
eigenvalues. The orange discs are $D_{i}$ defined by Eq. (7). Fig. 7 shows
that the union of discs for the more rigid molecule on the right is
$D_{1}\cup...\cup D_{19}=\\{-152.77\leq z\leq 292.16\\}$, which covers a wider
area compared to that for the molecule on the left, $D_{1}\cup...\cup
D_{19}=\\{-128.78\leq z\leq 248.76\\}$. Since the range of eigenvalues is
wider for the more rigid molecule, it has a higher probability of having a
larger $\epsilon_{\rm max}$ value and bigger $\sigma(\epsilon>1)$. Since the
molecules from the same constitutional isomer cluster have the same atomic
composition, the distribution of the center of the discs, ${D_{i},i\in N}$,
are identical. This explains the similarity of $\mu(\epsilon>1)$ for molecules
from the same cluster observed in our analysis (not shown) of the dataset.

To illustrate the physical significance of the descriptors based on PDF
derived from the Gershgorin circle theorem as defined by Eq. (11), we examine
in Figure 8 the correlation of entropy with each of the individual overlaps
$\langle f_{i},f_{\rm molecule}\rangle$ and with ${\rm AUC}(f_{\rm
molecule}(x|\mu=M_{ii},\sigma^{2}=(\sum_{j\neq i}|M_{ij}|)^{2})$. Figure 8
shows that entropy is particularly well correlated with $\langle f_{\rm
H},f_{\rm molecule}\rangle$ and $\langle f_{\rm C},f_{\rm molecule}\rangle$
and strongly correlated with ${\rm AUC}(f_{\rm
molecule}(x|\mu=M_{ii},\sigma^{2}=(\sum_{j\neq i}|M_{ij}|)^{2})$.

Figure 8:  Correlation of entropy with AUC (upper left panel) and $\langle
f_{\rm atom},f_{\rm molecule}\rangle$ for constitutional isomers with the
chemical formula $\rm C_{7}H_{10}O_{2}$.

###  IV.2 Data-efficient models of entropy

Figure 9:  Left: BIC values for models with different kernels. Right: Learning
curves of GP models trained with the same 9D descriptor
($\rm{CM}+\rm{AUC}[f_{\rm{molecule}}]+\langle f_{\rm HCNOF},f_{\rm
molecule}\rangle$), but with the different kernels illustrated in the left
panel. The shaded area is the standard deviation computed with ten models
trained with randomly sampled data points. For the lines not showing the
shaded area, the standard deviation is too small to be visible on this scale.

All results in the present work show entropy multiplied by temperature at
$T=298.15$ K in the units of kcal mol-1. Our goal is to build models with
prediction accuracy below 1 kcal mol-1 using as few molecules on input as
possible. We aim to achieve this by exploiting the most performant descriptors
identified in the previous section in combination with kernel models based on
optimized kernels. When the dimensionality of descriptors is low, it is
possible to build anisotropic kernels, optimized for each dimension of the
input space separately. In the present work, we have found that including the
anisotropy into the kernels for the best low-dimensional descriptors of
molecules does not have any significant effect on the accuracy of the ensuing
models.

We use GP regression to model entropy. The kernel functions of the GP models
are chosen and optimized using greedy search in the space of kernel functions
guided by BIC as discussed in the ML method section. This aligns kernels with
available data as much as possible, given the constraints of the greedy
search, the descriptors and the uncertainties stemming from the extremely
limited size of the datasets considered here. To reduce bias in the
calculation of BIC used for the greedy search of the optimal kernel function,
1000 training data points are chosen randomly from the pool of 85,000 data
points reserved for training. Each of the randomly chosen training datasets is
passed through multiple layers of the BIC kernel selection algorithm. At each
level of kernel selection, we calculate the average value of BIC scores,
$\mu$(BIC), for the same kernel trained with ten different sets of randomly
selected data points. The base or complex kernel with the highest average BIC
score is retained for the next level of kernel function selection.

The effect of the kernel function on the accuracy of the resulting models of
entropy is shown in Fig. 9. Figure 9 compares the $\mu$(BIC) values of base
and complex kernels. Although the $\rm{RQ}+\rm{Matern}\times\rm{DP}$ kernel
has the highest BIC score, this BIC score is only slightly higher than that of
the $\rm{RQ}+\rm{Matern}$ kernel. The $\rm{RQ}+\rm{Matern}$ kernel is
therefore chosen as the optimal kernel for the best model of entropy in this
work due to its high BIC score and lower function complexity. The BIC
algorithm stops at level four because the BIC scores plateau and decrease with
further increase of the kernel complexity.

###  IV.3 Zero-point Vibrational Energy

The preceding sections illustrate that very efficient and accurate models of
entropy can be built with 9D descriptors
($\rm{CM}+\rm{AUC}[f_{\rm{molecule}}]+\langle f_{\rm HCNOF},f_{\rm
molecule}\rangle$). It is instructive to explore the performance of these
descriptors in models of other molecular properties. Here, we consider zero
point vibrational energy (ZPVE) of molecules. Figure 10 illustrates that ZPVE
is strongly correlated with ${\rm AUC}(f_{\rm molecule})$ and $\langle f_{\rm
ref},f_{\rm molecule}\rangle$.

Figure 10:  Zero point vibrational energy as a function of ${\rm AUC}(f_{\rm
molecule}(x|\mu=M_{ii},\sigma^{2}=(\sum_{j\neq i}|M_{ij}|)^{2})$ and $\langle
f_{\rm ref},f_{\rm molecule}\rangle$ for all molecules in the QM9 dataset
(left panels); and the molecules from ten largest isomer groups in the QM9
dataset (right panels). The reference molecule is ${\rm H_{3}C_{4}N_{2}OF}$
illustrated in Fig. 3.

Figure 11:  Comparison of the present results (squares) of ZPVE with the
results of Lilienfeld et al. (2023) khan2023quantum . The present results are
obtained with 9D molecular descriptors ${\rm CM}+{\rm AUC}+f_{\rm
molecule}+\langle f_{\rm HCNOF},f{\rm molecule}\rangle$ and the GP model with
the optimized ${\rm RQ+MAT\times DP}$ kernel. The results of Ref.
khan2023quantum  are based on kernel ridge regression with the following
descriptors: full CM including 435 features (orange circles), BOB with 1128
features (green triangles), SLATM 11960 features (red squares), FCHL with
10440 features (purple down triangles), and MBDF with 145 features (brown
diamonds).

Figure 11 shows the interpolation MAE for ZPVE illustrating the comparison of
the performance of the present 9D models with models developed in Ref.
khan2023quantum  . The results in Ref. khan2023quantum  use high-dimensional
descriptors: full CM including 435 features (shown in Fig. 11 orange circles),
BOB with 1128 features (shown by green triangles), SLATM 11960 features (shown
by red squares), FCHL with 10440 features (shown by purple down triangles),
and MBDF with 145 features (shown by brown diamonds). These high-dimensional
descriptors were used in Ref. khan2023quantum  to compute simple isotropic
kernels for kernel ridge regression models of ZPVE.

Figure 12:  Distribution of errors for predictions of entropy (upper) and ZPVE
(lower) with nine-dimensional descriptors of molecules identified as the best
descriptors in Fig. 5. Full bars – predictions of models trained with 100
molecules; open bars – models trained with 1000 molecules.

Our results shown in Fig. 11 by blue diamonds use nine-dimensional descriptors
with optimized, complex kernels given by the ${\rm RQ+MAT\times DP}$
combination of base kernels. These kernels are used for GP regression models
of ZPVE. The results illustrate that the 9D descriptors proposed here produce
ZPVE interpolation models with the prediction MAE under 1 kcal mol-1, when
trained by as few as 100 molecules.

While it is common to quantify the model performance by average errors, such
metrics do not illustrate fully the quality of the models. A combination of
average errors and largest errors is a better metric, which, however, can be
significantly affected by a small number of outliers. To quantify non-
ambiguously the performance of the bets models of entropy and ZPVE developed
in the present work, we show in Fig. 12 the distribution of errors of
predictions over a hold-out test set, containing 20,000 molecules of different
size. The shaded bars illustrate the performance of models trained with 100
molecules and the open bars – molecules trained with 1000 molecules. We note
that the range of $S\times T$ for the molecules considered here covers 26 kcal
mol-1, while the range of ZPVE covers 160 kcal mol-1. In both cases, we
observe that our models trained with 100 molecules are able to predict
$S\times T$ and ZPVE with absolute error under 1 kcal mol-1 for $>78$ % of the
test data and with error under 1.3 kcal mol-1 for $>92.4$ % of the test data.

##  V Summary

The goal of the present work is to build efficient, data-starved models of
molecular properties for interpolation in chemical compound spaces comprising
molecules of variable complexity with low-dimensional molecular descriptors.
The dataset considered here comprises 133,885 molecules with complexity
ranging from 3 atoms to 29 atoms. We consider two molecular properties:
entropy and zero point vibrational energy. We begin by designing low-
dimensional descriptors tailored for the prediction of entropy. Our results
show that entropy can be efficiently interpolated by nine-dimensional models
with molecular descriptors comprising three physical parameters describing the
distributions of eigenvalues of Coulomb matrices and six parameters accounting
for the composition and shape of molecules through quantities based on the
Gershgorin circle theorem.

In order to build efficient interpolation models of entropy, we combine the
nine-dimensional descriptors thus obtained with Gaussian process regression
based on kernels with variable functional form. The kernel functions are
optimized using a greedy search algorithm guided by the Bayesian information
criterion (BIC). In order to obtain meaningful, unbiased models based on
extremely restricted datasets, we use BIC averaged over multiple randomly
sampled instances of training data as a model selection metric. Nine-
dimensional Gaussian process models thus optimized are capable of predicting
molecular entropy across the entire dataset with the average error under 1
kcal mol-1, using as few as 100 molecules on input. We present a careful
analysis of model prediction errors, demonstrating that entropy models trained
with 100 molecules predict $S\times T$ with chemical accuracy under 1 kcal
mol-1 for $>80$ % of molecules in the test set, comprising 20,000 molecular
species.

The nine-dimensional descriptors based on the distributions of the Coulomb
matrix eigenvalues and the Gershgorin circle theorem are then shown to perform
well for the prediction of a completely different molecular property - ZPVE.
The error analysis illustrates that the Gaussian process models of ZPVE
trained with 100 molecules and using optimized kernels predict ZPVE with
absolute error under 1 kcal mol-1 for $>80$ % of the test set and with error
under 1.3 kcal mol-1 for $>90$ % of molecules in the test set. This is
remarkable, given that ZPVE of molecules in the test set covers the range from
10 to 171 kcal mol-1.

There are three main implications of the results presented here. First,
because the training sets of the demonstrated models are very small, accurate
predictions of molecular properties for a wide range of molecules can be based
on either experimental data or extremely accurate quantum chemistry
calculations. Second, the possibility of using low-dimensional descriptors for
accurate predictions with restricted datasets opens up (1) the possibility of
building models with anisotropic kernels that account for the difference in
scales of different features; and (2) the possibility of Bayesian optimization
in chemical compound space, which can be used to predict and design molecules
with desired properties. We will explore both of these possibilities in a
future work.

## Acknowldgments

This work was supported by NSERC of Canada.

## References

  * (1) Polishchuk P G, Madzhidov T I and Varnek A 2013 J. Comput. Aided Mol. Des. 27 675
  * (2) Weinreich J, Browning N J and von Lilienfeld O A 2021 J. Chem. Phys. 154 134113
  * (3) Hansen K, Biegler F, Ramakrishnan R, Pronobis W, von Lilienfeld O A, Müller K R and Tkatchenko A 2015 J. Phys. Chem. Lett. 6 2326
  * (4) Browning N J, Ramakrishnan R, von Lilienfeld O A and Roethlisberger U 2017 J. Phys. Chem. Lett. 8 1351
  * (5) Gubaev K, Podryabinkin E V and Shapeev A V 2018 J. Chem. Phys. 148 241727
  * (6) Christensen A S, Bratholm L A, Faber F A and von Lilienfeld O A 2020 J. Chem. Phys. 152 044107
  * (7) Bogojeski M, Vogt-Maranto L, Tuckerman M E, Müller K R and Burke K 2020 Nat. Commun. 11 5223
  * (8) Gastegger M, Schwiedrzik L, Bittermann M, Berzsenyi F and Marquetand P 2018 J. Chem. Phys. 148 241709
  * (9) Rostami S, Amsler M and Ghasemi S A 2018 J. Chem. Phys. 149 124106
  * (10) Qiao Z, Welborn M, Anandkumar A, Manby F R and Miller Thomas F I 2020 J. Chem. Phys. 153 124111
  * (11) Musil F, Grisafi A, Bartók A P, Ortner C, Csányi G and Ceriotti M 2021 Chem. Rev. 121 9759
  * (12) Gallegos L C, Luchini G, St John P C, Kim S and Paton R S 2021 Acc. Chem. Res. 54 827
  * (13) Willatt M J, Musil F and Ceriotti M 2018 Phys. Chem. Chem. Phys. 20 29661
  * (14) Langer M F, Goeßmann A and Rupp M 2022 npj Comput. Mater. 8 41
  * (15) Faber F A, Hutchison L, Huang B, Gilmer J, Schoenholz S S, Dahl G E, Vinyals O, Kearnes S, Riley P F and Von Lilienfeld O A 2017 J. Chem. Theory Comput. 13 5255
  * (16) Faber F A, Christensen A S, Huang B and von Lilienfeld O A 2018 J. Chem. Phys. 148 241717
  * (17) Duvenaud D, Lloyd J, Grosse R, Tenenbaum J and Zoubin G 2013 Proceedings of the 30th Int. Conf. on Machine Learning
  * (18) Duvenaud D K, Nickisch H and Rasmussen C 2011 Adv. Neural Inf. Process. Syst. 24
  * (19) Deng Z, Tutunnikov I, Averbukh I S, Thachuk M and Krems R V 2020 J. Chem. Phys. 153 164111
  * (20) Gilmer J, Schoenholz S S, Riley P F, Vinyals O and Dahl G E 2017 Proceedings of the 34th Int. Conf. on Machine Learning
  * (21) Schütt K T, Sauceda H E, Kindermans P J, Tkatchenko A and Müller K R 2018 J. Chem. Phys. 148 241722
  * (22) Smith J S, Nebgen B T, Zubatyuk R, Lubbers N, Devereux C, Barros K, Tretiak S, Isayev O and Roitberg A E 2019 Nat. Commun. 10 2903
  * (23) Rupp M, Tkatchenko A, Müller K R and von Lilienfeld O A 2012 Phys. Rev. Lett. 108 058301
  * (24) Hansen K, Biegler F, Ramakrishnan R, Pronobis W, von Lilienfeld O A, Müller K R and Tkatchenko A 2015 J. Phys. Chem. Lett. 6 2326
  * (25) Bartók A P, Kondor R and Csányi G 2013 Phys. Rev. B 87 184115
  * (26) Kipf T N and Welling M 2016 arXiv preprint arXiv:1609.02907
  * (27) Wu Z, Pan S, Chen F, Long G, Zhang C and Yu P S 2021 IEEE Trans. Neural Netw. Learn. Syst. 32 4
  * (28) Zhu W, Zhang Y, Zhao D, Xu J and Wang L 2022 J. Chem. Inf. Model. 63 43
  * (29) Yang K, Swanson K, Jin W, Coley C, Eiden P, Gao H, Guzman-Perez A, Hopper T, Kelley B, Mathea M, Palmer A, Settels V, Jaakkola T, Jensen K and Barzilay R 2019 J. Chem. Inf. Model. 59 3370
  * (30) Shui Z and Karypis G 2020 IEEE Int. Conf. on Data Mining
  * (31) Williams C K and Rasmussen C E 2006 Gaussian processes for machine learning 3 (Cambridge, MA: MIT press)
  * (32) Gershgorin S A 1931 News of the Russian Academy of Sciences. Mathematical series 749
  * (33) Saad Y 2011 Numerical Methods for Large Eigenvalue Problems: Revised Edition Classics in Applied Mathematics (Philadelphia, PA: Society for Industrial and Applied Mathematics)
  * (34) Ramakrishnan R, Dral P O, Rupp M and von Lilienfeld O A 2014 Sci. Data. 1 140022
  * (35) Bellman R 1966 Science 153 34
  * (36) Schrier J 2020 J. Chem. Inf. Model. 60 3804
  * (37) Murphy K P 2018 Machine learning: A probabilistic perspective (adaptive computation and machine learning series) (London, UK: The MIT Press)
  * (38) Torabian E and Krems R V 2023 Phys. Rev. Res. 5 013211
  * (39) Vargas-Hernández R A and Krems R V 2020 Mach. Learn. Meets Quantum Phys. 968 171
  * (40) Asnaashari K and Krems R V 2021 Mach. Learn.: Sci. Technol. 3 015005
  * (41) Dai J and Krems R V 2023 Mach. Learn.: Sci. Technol. 4 045027
  * (42) Dai J and Krems R V 2020 J. Chem. Theory Comput. 16 1386
  * (43) Vargas-Hernández R A, Sous J, Berciu M and Krems R V 2018 Phys. Rev. Lett. 121 255702
  * (44) Huyskens P, Vandevijvere P and Siegel G 1989 J. Mol. Struct. 200 555
  * (45) Khan D, Heinen S and von Lilienfeld O A 2023 J. Chem. Phys. 159 034106


11institutetext: Universite catholique de Louvain, 10 place Cardinal Mercier,
Louvain-la-Neuve, 1348, Belgium 11email<EMAIL_ADDRESS>
# Using CollGram to Compare Formulaic Language in Human and Neural Machine
Translation††thanks: The author is a Research Associate of the Fonds de la
Recherche Scientifique - FNRS (Federation Wallonie Bruxelles de Belgique). He
would like to warmly thank S. Granger and M. Dupont for access to the PLECI
corpus.

Yves Bestgen 11 0000-0001-7407-7797

###### Abstract

A comparison of formulaic sequences in human and neural machine translation of
quality newspaper articles shows that neural machine translations contain less
lower-frequency, but strongly-associated formulaic sequences, and more high-
frequency formulaic sequences. These differences were statistically
significant and the effect sizes were almost always medium or large. These
observations can be related to the differences between second language
learners of various levels and between translated and untranslated texts. The
comparison between the neural machine translation systems indicates that some
systems produce more formulaic sequences of both types than other systems.

###### Keywords:

Neural machine translation Multiword unit Lexical association indice DeepL
Google Translate Microsoft Translator.

##  1 Introduction

Neural machine translation (NMT) systems are currently considered to bridge
the gap between human and machine translation [22, 26]. However, little
research has been done to determine whether NMT systems are also very
effective in processing multiword units [20, 27], whereas the importance of
preformed units in language use is now well established, including in foreign
language learning and translation [1, 21, 24]. The present study addresses
this issue by comparing formulaic language in human and neural machine
translation. It focuses on a specific category of multiword units, the
”habitually occurring lexical combinations” [17], such as dramatic increase,
depend on, out of, which are not necessarily semantically non-compositional,
but are considered statistically typical of the language because they occur
”with markedly high frequency, relative to the component words or alternative
phrasings of the same expression” [2]. These formulaic sequences (FSs) are
analyzed by means of a technique proposed by [4], improved by [13] and
automated by [10] under the name CollGram111As one reviewer pointed out, this
automation is not an ”easily available plug-and-play implementation”. However,
there is a freely available system that implements it
(http://collgram.pja.edu.pl) [18, 25]. Some of the indices used can also be
easily obtained with the TAALES software [16] which allows the automatic
analysis of many other lexical indices. TAALES presents however an important
limitation because it only takes into account bigrams that occur at least 51
times in the reference corpus [9], a value much too high for the MI at the
heart of the CollGram approach..

CollGram rests on two lexical association indices that measure the strength of
attraction between the words that compose a bigram, mutual information (MI)
and t-score [14], calculated on the basis of the frequencies of occurrences in
a reference corpus [10, 13]. These two indices are complementary, MI favoring
lower-frequency, but strongly-associated, FSs such as self-fulfilling
prophecy, sparsely populated or sunnier climes while the t-score favors high-
frequency bigrams such as you know, out of or more than. A series of studies
have shown that, compared to native speakers, English as a foreign language
learners tend to underuse collocations with high MI scores, while overusing
those with high t-scores and that exactly the same differences are observed
between advanced learners and intermediate learners [10, 13]. These
observations are in agreement with usage-based models of language learning
which ”hold that a major determining force in the acquisition of formulas is
the frequency of occurrence and co-occurrence of linguistic forms in the
input” [13]. It is worth noting that the same differences were observed
between translated and untranslated texts, but the proposed explanation relies
on a tendency towards normalization in translation [5, 8]. Since neural models
also seem to be affected by frequency of use [15, 19], the hypothesis tested
in the present study is that the same effects could be observed when comparing
human translations (HTs) and NMTs, namely that NMTs will underuse high MI FSs
and overuse high t-score FSs.

##  2 Method

###  2.1 Translation Corpus

The texts used are taken from the journalistic section of the PLECI corpus
(uclouvain.be/en/research-institutes/ilc/cecl/pleci.html). It is a sentence-
aligned translation corpus of quality newspaper articles written in French and
published in Le Monde diplomatique and in English in one of the international
editions of this same newspaper. Two hundred and seventy-nine texts, published
between 2005 and 2012, were used for a total of 570,000 words in the original
version and of 500,000 words in the translation.

All original texts were translated into English by three well-known NMT
systems: DeepL (deepl.com/translator), Google Translate (translate.google.com)
and Microsoft Translator (microsoft.com/translator). Online translators were
used for the first two, while the version available in Office 365 was used for
the third. All these translations were performed between March 24 and April 6,
2021.

###  2.2 Procedure

Each translated text was tokenized and POS-tagged by CLAWS7 [23] and all
bigrams were extracted. Punctuation marks and any character sequences that did
not correspond to a word interrupted the bigram extraction. Each bigram, which
did not include any proper name or number according to CLAWS, was then
searched for in the 100 million word British National Corpus (BNC222[6] showed
that CollGram produces the same results if another reference corpus, such as
COCA (corpus.byu.edu/coca) or WaCKy [3], is used., www.natcorp.ox.ac.uk). When
it is present, the corresponding MI and t-score were used to decide whether it
is highly collocational or not. Based on [8] and [13], bigrams with a score
greater than or equal to 5 for the MI and 6 for the t-score were considered
highly collocational. The last step consisted in calculating, for each text
and for each association index, the percentage that the bigrams considered as
highly collocational represent compared to the total number of bigrams present
in the text.

##  3 Analyses and Results

Table 1 shows the average percentages of highly collocational bigrams for the
MI and t-score in the four type of translations. The four means for each
measure of association were analyzed using the Student’s test for repeated
measures since the same texts, which are the unit of analysis, were translated
by the four translators. All these comparisons were statistically significant
($p<0.0001$).

Table 1: Average percentages of highly collocational bigrams for the two indices in the four translation types. Measure | Human | DeepL | Google | Microsoft
---|---|---|---|---
High MI | 11.21 | 10.48 | 10.07 | 10.27
High t-score | 58.76 | 60.60 | 59.49 | 59.89

Table 2 presents the differences between the means as well as two effect
sizes. The first is Cohen’s $d$, which expresses the size of the difference
between the two means as a function of the score variability. According to
[11], a $d$ of 0.50 indicates a medium effect and that a $d$ of 0.80 a large
effect. The second effect size indicates the percentage of texts for which the
difference between the two translations has the same sign as the mean
difference. A value of 100 means that all texts produced by a given translator
have a higher score than those translated by the other one and a value of 50
means that there is no difference.

Figure 1: Average percentages of highly collocational bigrams for the two
indices.

As shown in these tables and in Figure 1, both hypotheses are verified.
Compared to HTs, texts translated by the three neural systems contain a
significantly smaller percentage of highly collocational bigrams for the MI
and a larger percentage of highly collocational bigrams for the t-score.
Cohen’s $d$s are almost always medium or large and the percentages of texts
for which differences are observed are greater than 70% except in one case.

Table 2: Differences (row translator minus column translator) and effect sizes for the two indices in the four translation types. | Human | DeepL | Google
---|---|---|---
| D | Es | % | D | Es | % | D | Es | %
| High MI
DeepL | -0.72 | 0.59 | 73.48 |  |  |  |  |  |
Google | -1.14 | 1.00 | 84.33 | -0.41 | 0.65 | 74.91 |  |  |
Microsoft | -0.94 | 0.77 | 81.00 | -0.21 | 0.35 | 62.72 | 0.20 | 0.36 | 64.52
| High t-score
DeepL | 1.83 | 0.84 | 80.65 |  |  |  |  |  |
Google | 0.72 | 0.32 | 62.37 | -1.11 | 0.98 | 84.59 |  |  |
Microsoft | 1.13 | 0.51 | 71.33 | -0.70 | 0.62 | 70.97 | -0.41 | 0.42 | 69.18

An analysis of the passages in which the differences between HT and NMT are
the largest suggests that the origin lies at least partially in the less
literal nature of human translations (see Table 3 for an example).

Table 3: Example of the four translation types and percentages of highly collocational bigrams for the two indices. Type | Phrase | %High MI | %High t-score
---|---|---|---
Original | A raison de huit heures par jour |  |
Human | In an eight-hour day | 67 | 33
DeepL | At eight hours a day | 25 | 100
Google & Microsoft | At the rate of eight hours a day | 14 | 100

The differences between the three NMT systems are smaller, but still
statistically significant. However, they require a different interpretation.
When NMTs are compared to HTs (see Figure 1), the patterns of differences are
reversed according to the MI or the t-score, as expected. For the NMT systems
(see Figure 2), these patterns are identical for both types of collocation.
The average percentages of highly collocational bigrams (see Table 1) are
always higher in texts translated by DeepL than in those translated by
Microsoft and also higher in the latter than in those translated by Google.
Only a detailed qualitative analysis could determine whether these results
indicate a difference in effectiveness.

Figure 2: Average percentages of highly collocational bigrams for the three
NMT systems.

##  4 Discussion and Conclusion

The reported analyses confirm the hypotheses and thus suggest that, compared
to HTs, NMTs more closely resemble texts written by intermediate learners than
by advanced learners of English as a foreign language, a result that could be
interpreted in the context of a usage-based model of language learning [13].
The NMTs also resemble translated texts more than untranslated texts, but it
is not clear that this can be explained by a normalization process.
Statistically significant differences, but smaller in terms of effect size,
were also observed between the three NMT systems.

It is important to keep in mind that the present study only considers global
quantitative properties of MWUs. At no point is the appropriateness in context
of the MWUs assessed. It is therefore a very partial approach. However, it has
the advantage of not requiring a human qualitative evaluation that is often
complicated and cumbersome to set up. Moreover, it is likely that the
appropriateness of a MWU is much more important for non-compositional
expressions than for the habitually occurring lexical combinations studied
here [12].

Another important feature of the approach is that it relies on a native
reference corpus to identify highly collocational bigrams for both indices. As
already mentioned, research on foreign language learning, but also on the
comparison of translated and untranslated texts, has shown that the use of
other large reference corpora such as COCA or WaCKy [3] did not change the
results [5, 6]. One can also wonder whether the use of a comparable reference
corpus, rather than a generic one, would have returned different results. In
the case of the comparison of translated and untranslated texts, [5] observed
that the use of a journalistic corpus, the Corpus Est Republicain (115 million
words) made available by the Centre National de Ressources Textuelles et
Lexicales, produced differences similar to those obtained with the WaCKy
corpus.

Before considering taking advantage of these observations to try to improve
NMT systems, a series of complementary analyses must be conducted. Indeed,
this study has many limitations, such as focusing only on a subcategory of
MWUs [20], on a single language pair, and on a single genre of texts.
Moreover, a thorough qualitative analysis is essential to better understand
the results and evaluate the proposed explanations. As it has been shown in
foreign language learning [7], it would also be interesting to verify that the
observed effects are not explained by differences in single-word lexical
richness. Finally, the differences between the three NMT systems also require
further analysis.

## References

  * [1] Baker, M.: Patterns of idiomaticity in translated vs. non-translated text. Belgian Journal of Linguistics 21, 11–21 (2007)
  * [2] Baldwin, T., Kim, S.N.: Multiword expressions. In: Indurkhya, N., Damerau, F.J. (eds.) Handbook of Natural Language Processing, pp. 267–292. CRC Press (2010)
  * [3] Baroni, M., Bernardini, S., Ferraresi, A., Zanchetta, E.: The WaCky wide web: A collection of very large linguistically processed web-crawled corpora. Language Resources and Evaluation 43, 209–226 (2009)
  * [4] Bernardini, S.: Collocations in translated language. combining parallel, comparable and reference corpora. In: Proceedings of the Corpus Linguistics Conference. pp. 1–16. Lancaster University (2007)
  * [5] Bestgen, Y., Granger, S.: Collocation et traduction. analyse automatique au moyen d’indices d’association. In: Kauffer, M., Keromnes, Y. (eds.) Theorie und Empirie in der Phraseologie / Approches théoriques et empiriques en phraséologie, pp. 101–113. Stauffenburg Verlag (2019)
  * [6] Bestgen, Y.: Evaluation automatique de textes. Validation interne et externe d’indices phraséologiques pour l’évaluation automatique de textes rédigés en anglais langue étrangère. Traitement automatique des langues 57(3), 91–115 (2016)
  * [7] Bestgen, Y.: Beyond single-word measures: L2 writing assessment, lexical richness and formulaic competence. System 69, 65–78 (2017)
  * [8] Bestgen, Y.: Normalisation en traduction : analyse automatique des collocations dans des corpus. Des mots aux actes 7, 459–469 (2018)
  * [9] Bestgen, Y.: Evaluation de textes en anglais langue ètrangère et séries phraséologiques : comparaison de deux procédures automatiques librement accessibles. Revue française de linguistique appliquée 24, 81–94 (2019)
  * [10] Bestgen, Y., Granger, S.: Quantifying the development of phraseological competence in L2 English writing: An automated approach. Journal of Second Language Writing 26, 28–41 (2014)
  * [11] Cohen, J.: Statistical power analysis for the behavioral sciences. Erlbaum (1988)
  * [12] Constant, M., Eryiǧit, G., Monti, J., van der Plas, L., Ramisch, C., Rosner, M., Todirascu, A.: Survey: Multiword expression processing: A Survey. Computational Linguistics 43(4), 837–892 (dec 2017)
  * [13] Durrant, P., Schmitt, N.: To what extent do native and non-native writers make use of collocations? International Review of Applied Linguistics in Language Teaching 47, 157–177 (2009)
  * [14] Evert, S.: Corpora and collocations. In: Lüdeling, A., Kytö, M. (eds.) Corpus Linguistics. An International Handbook, pp. 1211–1248. Mouton de Gruyter (2009)
  * [15] Koehn, P., Knowles, R.: Six challenges for neural machine translation (2017)
  * [16] Kyle, K., Crossley, S., Berger, C.: The tool for the automatic analysis of lexical sophistication (TAALES): version 2.0. Behavior Research Methods 50, 1030–1046 (2018). https://doi.org/10.3758/s13428-017-0924-4
  * [17] Laufer, B., Waldman, T.: Verb-noun collocations in second language writing: A corpus analysis of learners’ English. Language Learning 61, 647–672 (2011)
  * [18] Lenko-Szymanska, A., Wolk, A.: A corpus-based analysis of the development of phraseological competence in EFL learners using the collgram profile. In: Paper presented at the 7th Conference of the Formulaic Language Research Network (FLaRN) (2016)
  * [19] Li, M., Roller, S., Kulikov, I., Welleck, S., Boureau, Y.L., Cho, K., Weston, J.: Don’t say that! making inconsistent dialogue unlikely with unlikelihood training. In: Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics. pp. 4715–4728. Association for Computational Linguistics (2020)
  * [20] Monti, J., Seretan, V., Pastor, G.C., Mitkov, R.: Multiword units in machine translation and translation technology. In: Mitkov, R., Monti, J., Pastor, G.C., Seretan, V. (eds.) Multiword Units in Machine Translation and Translation Technology, pp. 2–37. John Benjamins (2018)
  * [21] Pawley, A., Syder, F.H.: Two puzzles for linguistic theory: nativelike selection and nativelike fluency. In: Richards, J.C., Schmidt, R.W. (eds.) Language and Communication. Longman (1983)
  * [22] Popel, M., Tomkova, M., Tomek, J., Kaiser, L., Uszkoreit, J., Bojar, O., Zabokrtsky, Z.: Transforming machine translation: A deep learning system reaches news translation quality comparable to human professionals. Nature Communications 11, 1–15 (2020). https://doi.org/10.1038/s41467-020-18073-9
  * [23] Rayson, P.: Matrix: A statistical method and software tool for linguistic analysis through corpus comparison. Ph.D. thesis, Lancaster University (1991)
  * [24] Sinclair, J.: Corpus, Concordance, Collocation. Oxford University Press (1991)
  * [25] Wołk, K., Wołk, A., Marasek, K.: Unsupervised tool for quantification of progress in L2 english phraseological. In: Proceedings of the 2017 Federated Conference on Computer Science and Information Systems. pp. 383–388 (2017)
  * [26] Wu, Y., Schuster, M., Chen, Z., Le, Q.V., Norouzi, M., Macherey, W., Krikun, M., Cao, Y., Gao, Q., Macherey, K., Klingner, J., Shah, A., Johnson, M., Liu, X., Łukasz Kaiser, Gouws, S., Kato, Y., Kudo, T., Kazawa, H., Stevens, K., Kurian, G., Patil, N., Wang, W., Young, C., Smith, J., Riesa, J., Rudnick, A., Vinyals, O., Corrado, G., Hughes, M., Dean, J.: Google’s neural machine translation system: Bridging the gap between human and machine translation (2016)
  * [27] Zaninello, A., Birch, A.: Multiword expression aware neural machine translation. In: Proceedings of the 12th Conference on Language Resources and Evaluation (LREC 2020). pp. 3816–3825. European Language Resources Association (ELRA) (2020)