# Multitask methods for predicting molecular properties from heterogeneous
data

K. E. Fisher<EMAIL_ADDRESS>Department of Aeronautics and Astronautics,
Massachusetts Institute of Technology  M. F. Herbst  Mathematics for Materials
Modelling, Institute of Mathematics & Institute of Materials, École
Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland  National
Centre for Computational Design and Discovery of Novel Materials (MARVEL),
École Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland  Y. M.
Marzouk  Department of Aeronautics and Astronautics, Massachusetts Institute
of Technology

###### Abstract

Data generation remains a bottleneck in training surrogate models to predict
molecular properties. We demonstrate that multitask Gaussian process
regression overcomes this limitation by leveraging both expensive and cheap
data sources. In particular, we consider training sets constructed from
coupled-cluster (CC) and density function theory (DFT) data. We report that
multitask surrogates can predict at CC level accuracy with a reduction to data
generation cost by over an order of magnitude. Of note, our approach allows
the training set to include DFT data generated by a heterogeneous mix of
exchange-correlation functionals without imposing any artificial hierarchy on
functional accuracy. More generally, the multitask framework can accommodate a
wider range of training set structures—including full disparity between the
different levels of fidelity—than existing kernel approaches based on
$\Delta$-learning, though we show that the accuracy of the two approaches can
be similar. Consequently, multitask regression can be a tool for reducing data
generation costs even further by opportunistically exploiting existing data
sources.

††preprint: AIP/123-QED

##  I Introduction

Predicting molecular properties based on first-principle approaches drives
innovation across the physical sciences. However, exact prediction is
infeasible due to the exponential scaling of the Schrödinger equation—the
underlying quantum-mechanical description of materials systems. Decades of
research have yielded numerous first-principle models to describe electronic
structures approximately, each representing a compromise between computational
cost and prediction accuracy. Coupled cluster (CC) methods are typical
examples of high fidelity models which offer accuracy in exchange for
considerable computational cost. Specifically, CCSD(T)—the “gold standard” of
quantum-chemical prediction—suffers from a steep $O(N^{7})$ scaling where $N$
is the number of electrons in the target system. Harding  _et al._ (2008) In
contrast, density-functional theory (DFT) methods are popular baseline models
since they provide adequate accuracy for most systems while enjoying
substantially lower cost (scaling only as $O(N^{3})$). Ying, Yu, and Ying
(2019); Perdew and Schmidt (2001); Goerigkab and Grimme (2011); Teale  _et
al._ (2022)

Traditional workflows force users to settle on a single model which is cheap
enough to result in manageable cost but still sufficiently high fidelity to
give meaningful predictions. However, the advent of modern data-driven
modeling has introduced new tools which make it possible to decouple expensive
data generation from prediction. Statistical surrogates are trained with
predictions from first-principle models and, once trained, can be substituted
into workflows, making the prediction step considerably cheaper. In the
context of building potential energy surfaces and interatomic potentials,
these data-driven approaches have become the de facto standard. Smith  _et
al._ (2019); Dral  _et al._ (2020); Goodlett, Turney, and Schaefer (2023);
Zaverkin  _et al._ (2023) Simultaneously, large-scale data generation is
facilitated by sophisticated workflow engines for high-throughput calculations
Curtarolo  _et al._ (2012); Jain  _et al._ (2011); Huber  _et al._ (2020) and
robust DFT methods inspired by mathematical research. Teale  _et al._ (2022);
Cancès  _et al._ (2022); Herbst and Levitt (2020, 2022); Cancès  _et al._
(2023) The increasing amount of openly available data Ruddigkeit  _et al._
(2012); Ramakrishnan  _et al._ (2014); Chanussot*  _et al._ (2021); Smith,
Isayev, and Roitberg (2017) could in principle be readily employed, but
quality is generally not homogeneous as each dataset’s author has made
different modeling choices. Consequently, most surrogate models are single
fidelity—all training data is generated using the same quantum-chemical method
and the same numerical approach. Thus, heterogeneity has been an obstacle to
opportunistically exploiting the already available data. With this work, we
contribute to the promising development of multifidelity surrogates which do
not require a single choice between accurate (but expensive) and fast (but
crude) methods, but instead combine the advantages of both. Specifically, we
consider multitask methods which construct a relationship between multiple
regression problems without requiring precise information about the respective
accuracy of each source. Bonilla, Chai, and Williams (2008); Leen, Peltonen,
and Kaski (2012)

Gaussian process (GP) regression is a popular tool for constructing surrogates
informed by heterogeneous data. Several works Pilania, Gubernatis, and Lookman
(2017); Batra  _et al._ (2019) have explored the use of multifidelity GP
models Kennedy and O’Hagan (2000) trained using two fidelities, namely a high
and a low fidelity DFT approach. Alternatively, $\Delta$-learning may be
employed to model the difference between two levels of theory, then by adding
that difference to low level predictions of the target to produce high level
predictions. Ramakrishnan  _et al._ (2015); Dral, Zubatiuk, and Xue (2023);
Batra  _et al._ (2019) Though GP regression and the related technique of
kernel ridge regression are common choices for fitting such models, in
principle, any method of fitting a surrogate may be employed. Batra  _et al._
(2019); Dral, Zubatiuk, and Xue (2023). $\Delta$-learning can be extended to
leverage more than two training data sets if multiple difference models are
learned and summed. Dral  _et al._ (2020); Dral, Zubatiuk, and Xue (2023)
However, between each pair of training sets, data points must be in some way
aligned to allow for learning a difference model. For example, in a model of
the difference between DFT and CCSD(T), DFT predictions may only contribute to
the training set if a corresponding CCSD(T) prediction is also available.
Further, when either $\Delta$-learning or multifidelity GP regression learn
from more than two heterogeneous training sets, these sets must be assigned
some order—for example, following a hierarchy in accuracy. Recent work Vinod,
Kleinekathöfer, and Zaspel (2023) has proposed an optimized linear combination
of submodels as an alternative, but all coefficients in this linear
combination must be re-optimized whenever any submodel is altered, and
opportunistic data set construction is not fully explored.

Progress towards multifidelity approaches has also been made in the context of
neural network surrogates. For example, transfer learning methods train an
initial neural network on a larger source dataset, then re-optimise part of
the parameters using a target dataset. Towell and Shavlik (1994); Fu (1989);
Smith  _et al._ (2019) For the prediction of atomistic properties, the source
dataset might contain millions of DFT calculations while the target might be a
smaller set of CCSD(T) predictions. In prior studies, it has been observed
that transfer learning can produce higher accuracy than training a network on
only one of these datasets. Dral  _et al._ (2020); Dral, Zubatiuk, and Xue
(2023) While the framework we discuss in this work is not restricted to GP
regression strategies, we will focus exclusively on GP methods. In comparison
to neural network approaches, where a rigorous mathematical explanation for
the low generalization error of well-trained models is still lacking Bartlett,
Montanari, and Rakhlin (2021); Lotfi  _et al._ (2022), GP approaches are more
amenable to mathematical analysis. Rasmussen and Williams (2006) This
comparative transparency simplifies interpretation of the surrogate’s
predictive process in the diverse training data configurations we analyze
here. Additionally, GP methods have the advantage of yielding accurate
predictions even for relatively small training set sizes. While the training
cost of GP regression may become large as the size of the training set
increases, many strategies for efficient scaling have been developed and can
be routinely employed Quinonero-Candela and Rasmussen (2005); Liu  _et al._
(2018); Wilson, Dann, and Nickisch (2015); Cole, Christianson, and Gramacy
(2021).

This work explores a GP regression approach which can accommodate an arbitrary
number of quantum chemical datasets without an assumption of hierarchy. In
particular, we employ the multitask framework Leen, Peltonen, and Kaski
(2012), which has been shown to easily incorporate heterogeneous data sources.
The notion of heterogeneity in fact encompasses many possibilities. Fig. 1
compares the multitask method to other approaches by visualizing mappings
between sets of molecules and the quantum chemistry methods which are used to
generate viable training data sets. The highest level of theory used to
construct these training sets will be called the primary level. In our case,
the primary level will be CCSD(T). All other levels of theory will be
categorized as secondary levels. Our secondary levels will include DFT methods
using different density functional approximations (DFAs). Different levels of
theory may contribute data points corresponding to the same set of molecules,
different sets of molecules, or some union of shared and unique sets of
molecules. Furthermore, the testing data set contains primary level
predictions for target molecules. Thus, the multitask method has the option to
use secondary predictions for the target molecules as part of the training
set. The diagrams in Fig. 1 show several options for constructing a training
set for inference. The impact of these choices on inference is under-explored
in literature, and the multitask method is a useful tool for investigating the
benefits of different training set configurations. Crucially, the multitask
framework enables training with a “dataset of opportunity,” amalgamated from
existing data sets, as an alternative to expending computational time to
generate entirely new data.

Figure 1:  Diagram of training data sets compatible with the methods
considered in this work. Here, we define three sets of molecules—Core,
Additional, and Target systems—and two quantum chemistry methods—primary and
secondary. Our goal is to train a model that will predict our quantity of
interest for the Target molecules with the accuracy of the primary method. By
definition, the Core molecules are all systems for which we can compute
training data with the primary method. Additional molecules are systems for
which we have no primary data but which are not a part of the Target set.
Theoretically, we may train using secondary predictions for molecules of all
sets though not all inference methods support use of all data sets. The
diagram shows the combinations of data sets and quantum chemistry methods
supported by each inference method.

###  I.1  Advantages of the Multitask Method

Our numerical experiments will demonstrate the following advantages of the
multitask modeling approach. Key figures illustrating each claim are
referenced in parentheses.

  * •

Multitask models can achieve target accuracy with less computational cost than
Gaussian process regression (Fig. 3).

  * •

Accuracy of multitask inference improves as more DFT data is added to the
training set. Furthermore, increasing the amount of training data from
different secondary methods for the same molecule improves accuracy (Fig. 6).

  * •

Unlike the $\Delta$ and multifidelity approaches, multitask methods can model
more than two observation data sets without imposing an arbitrary order on the
sets (Figs. 5 and 6)

  * •

Multitask models have flexible training set requirements. The sets are
constructed from multiple levels of theory, and the same molecules do not need
to be included for each level (Fig. 3).

Flexibility is a key advantage of the multitask method compared to a single
task or $\Delta$-learning approach. The method is compatible with many
training set structures. For instance, a multitask model can be trained with
CCSD(T) data for one molecule set and DFT data for a different molecule set,
with CCSD(T) and DFT data for the same molecule set, or with some in
intermediary case. Models also have the option to train with DFT data for the
molecules that we aim to make predictions for at the CCSD(T) level of
accuracy. The flexible multitask framework makes it possible to explore the
impact of different compositions of a training data set on inference as seen
in Fig. 3. Consequently, we can determine the balance of expensive high level
and cost effective low level data which best suits the problem setting and
available data.

In Section II, we will review prediction with GP regression as well as the
$\Delta$ and multifidelity methods before describing the multitask framework.
Section III provides detail on our approach to testing different training set
configurations, and Section IV describes our numerical tests. We demonstrate
the performance of the multitask method on two examples: prediction of three-
body energy of water trimer configurations and prediction of ionization
potential of small organic molecules.

##  II Statistical Models

We review the formulation of Gaussian process regression models. Subsections
II.2 through II.4 discuss variants of GP regression which leverage
heterogeneous sets. See references Rasmussen and Williams (2006); Kennedy and
O’Hagan (2000); Leen, Peltonen, and Kaski (2012); Bonilla, Chai, and Williams
(2008); Abrahamsen (1997) for more detail.

###  II.1  Gaussian Process Regression

Our goal is to predict values of some quantity of interest (QoI) $f$, given an
identifying input feature $\bm{X}$. We train our predictive model using
$\\{\bm{X}_{i},\bm{Y}_{i}\\}_{i=1}^{n}$ where each $\bm{Y}_{i}$ is a noisy
observation of $f(\bm{X}_{i})$. In our setting, we interpret imperfect first-
principle predictions of the QoI as noisy observations of the true quantity.
Each $\bm{X}_{i}$ is some representation of the electronic structure of a
molecular system. For our numerical results, we build $\bm{X}_{i}$ using the
well established Smooth Overlap of Atomic Potentials (SOAP). Bartók  _et al._
(2010); Bartók, Kondor, and Csányi (2013); Bartók and Csányi (2015)
$\bm{Y}_{i}$ is the CCSD(T) or DFT prediction of the quantity of interest for
the $i^{th}$ molecular configuration.

We will assume that each observation is independently perturbed by additive
noise which follows identical Gaussian distributions. The corresponding model
is

$\displaystyle\bm{Y}_{i}\ $ $\displaystyle=$ $\displaystyle\ f(\bm{X}_{i})\ +\
\varepsilon_{i},$ (1) $\displaystyle\varepsilon_{i}\ $ $\displaystyle\sim$
$\displaystyle\ N\left(0,\ \sigma^{2}\right).$

Before inference, we choose a mean function, $\mu(\bm{X})$ to describe the
expected value of $f(\bm{X})$ and a kernel function, $k(\cdot,\cdot)$, to
encapsulate our assumptions about the relationship between $f(X_{i})$ and
$f(X_{j})$. Typically, $k(\cdot,\cdot)$ will have hyperparameters that are
selected through optimization, then fixed during inference.

We will employ a Gaussian process prior model for our regression function:

$\displaystyle f(\bm{X})\ \sim\
\text{GP}\Big{(}\mu(\bm{X}),k(\bm{X},\bm{X}^{\prime})\Big{)}.$ (2)

As a consequence of this prior assumption, any finite collection,
$\\{f(\bm{X}_{1}),\dots,f(\bm{X}_{n}\\})$, will follow a multivariate normal
distribution with mean $\bm{\mu}=\\{\mu(\bm{X}_{1}),\dots,\mu(\bm{X}_{n})\\}$
and covariance matrix $K$ with entries at row $i$ and column $j$ given by
$k(\bm{X}_{i},\bm{X}_{j})$. We are interested in inferring the values the
function takes on at target inputs, and we collect these values in the vector,
$\bm{f}_{*}$. Using this notation, we write the joint distribution on the
observations and targets for inference:

$\displaystyle\begin{bmatrix}\bm{Y}\\\ \bm{f}_{*}\end{bmatrix}\ \sim\
N\left(\begin{bmatrix}\bm{\mu}\\\ \bm{\mu}_{*}\end{bmatrix},\
\begin{bmatrix}K_{pp}+\sigma^{2}I&K_{p*}\\\\[3.0pt]
K_{p*}^{T}&K_{**}\end{bmatrix}\right).$ (3)

where $I$ is the identity matrix. The specific subscripts of each covariance
matrix represent the features that the matrix compares: $p$ corresponds to our
primary training data and $*$ indicates the target set. We note that in this
single-fidelity setting all training data is contained in the primary set, but
this will no longer be the case when we consider surrogates informed by
heterogeneous data.

The marginal distribution of $\bm{Y}$ is the model evidence, and conditioning
$\bm{Y}$ on $\bm{f}_{*}$ yields the likelihood. To obtain the posterior, we
employ the well known formula for conditioning part of a multivariate Gaussian
random variable on the rest, Rasmussen and Williams (2006) finding that

$\displaystyle\bm{f}_{*}\ \big{|}\bm{Y}\ \sim\ N\Big{(}\bm{\mu}_{*}\ +\
K_{p*}^{T}(K_{pp}+\sigma^{2}I)^{-1}\big{(}\bm{Y}-\bm{\mu}\big{)},$ (4)
$\displaystyle\quad K_{**}\ -\
K_{p*}^{T}(K_{pp}+\sigma^{2}I)^{-1}K_{p*}\Big{)}.$

Eq. (4) gives the probability distribution of our target values, conditioned
on our observations. The mean of this distribution is a correction to our
prior mean based on the deviations of the observations from that mean. The
predicted covariance indicates how much our observations have reduced our
uncertainty, provided our assumptions are reasonably aligned to reality.

###  II.2 $\Delta$-Learning

Suppose that for each training input, we have data on our QoI from multiple
sources. We thus seek an inference model that leverages
$\left\\{\bm{Y}_{ij}\right\\}_{i=1,j=1}^{n,m}$, where we have data for
molecular configurations $i=1,\dots,n$ and methods $j=1,\dots,m$. A simple
“pooling” approach may put all observations into one vector and proceed with
inference via Eq. (4). Then, $\bm{Y}_{ij_{1}}$ and $\bm{Y}_{ij_{2}}$ are
interpreted as independent observations of the QoI, $f(\bm{X}_{i})$, for
molecular configuration $\bm{X}_{i}$. The distinct behavior of each quantum
chemical method is treated as a perturbation of some mean value by a
realization of the noise distribution. We seek instead a model where the
distribution of each $\bm{Y}_{ij}$ is sensitive to $j$. Ideally, we want to
map differences in electronic structures methods to differences in their
predictions.

Method-specific features could be a solution, but construction of such
features is challenging. One approach involves appending a vector which
indicates the method category to each feature. This strategy puts a large
burden on covariance hyperparameter optimization since one set of
hyperparameters must encode differences between molecular configurations and
between quantum chemical methods. Alternatively, to characterize the
qualitative difference between methods within the features, we would have to
solve the very problem for which we are creating the features in the first
place—for instance, we would need to reliably predict CCSD(T) performance
given CCSD(T) and DFT data from training molecules as well as DFT data for
target molecules. Instead, we use this subsection to discuss
$\Delta$-learning, which directly models the differences between the
predictions of quantum chemistry methods.

$\Delta$-learning performs Gaussian process regression to predict differences,
$\Delta_{i}$, between observations obtained with two different methods for a
given $\bm{X}_{i}$. If we retain our assumption that noise follows an
independent, identically distributed Gaussian distribution, we have the model

$\displaystyle\bm{Y}_{ij_{1}}-\bm{Y}_{ij_{2}}\ $ $\displaystyle\equiv$
$\displaystyle\ \Delta_{i}\ =\ f\big{(}\bm{X}_{i}\big{)}\ +\ \varepsilon_{i},$
(5) $\displaystyle f\ $ $\displaystyle\sim$ $\displaystyle\
\text{GP}\big{(}\bm{\mu}_{\Delta},\
k_{\Delta}(\bm{X},\bm{X}^{\prime})\big{)}.$

where $j_{1}$ and $j_{2}$ indicate two prediction methods. An application of
Eq. (4) yields a distribution on the regression function for $\Delta$. Suppose
our ultimate interest is to predict $\bm{Y}$ at $\bm{X}_{*}$. Within the
$\Delta$ framework, we can supply an observation from one method, say
$\bm{Y}_{*j_{2}}$, in order to make a prediction for the value from the other
method:

$\bm{Y}_{*j_{1}}\ \approx\ \Delta_{i}\ +\ \bm{Y}_{*j_{2}}.$

Thus, $\Delta$-learning is suited to settings where calculations by method
$j_{1}$ are both more accurate and more time consuming than calculations that
use $j_{2}$.

Consider a case where we want to build a model from an existing data set which
includes predictions from CCSD(T) and several DFT methods. If the predictions
all correspond to the same set of molecules, $\Delta$-learning models have the
ability to leverage more of the data set than conventional GP regression
models. Unfortunately, $\Delta$-learning is not-applicable in cases where
predictions by different methods correspond to different molecular
configurations. To use a DFT prediction for a molecule $\bm{X}_{i}$ to train a
$\Delta$ model between CCSD(T) and DFT, we also require a CCSD(T) prediction
for the same molecule. In the following subsections, we consider alternative
approaches which model disparities in data sources and have the ability to use
all available observations.

###  II.3 Multifidelity Method

As a first step toward considering statistical models which train on
heterogeneous data, we describe the multifidelity data fusion model. Kennedy
and O’Hagan (2000) The multifidelity approach defines two regression
functions: $f_{p}$,representing the predictions from some high fidelity model,
and $f_{s}$, representing a low fidelity model. Kennedy and O’Hagan (2000) We
assume the relationships

$\displaystyle\bm{Y}_{ip}\ $ $\displaystyle=$ $\displaystyle\
f_{p}\big{(}\bm{X}_{i}\big{)}\ +\varepsilon_{ip}$ (6) $\displaystyle=$
$\displaystyle\ \rho f_{s}\big{(}\bm{X}_{i}\big{)}\ +\
\delta_{ps}\big{(}\bm{X}_{i}\big{)}\ +\ \varepsilon_{ip},$
$\displaystyle\bm{Y}_{is}\ $ $\displaystyle=$ $\displaystyle\
f_{s}\big{(}\bm{X}_{i}\big{)}\ +\ \varepsilon_{is}.$ (7)

The term $\delta_{ps}$ captures the difference between the high fidelity model
and the low fidelity model, scaled by hyperparameter, $\rho$. This disparity
is endowed with a Gaussian process prior, as is the low fidelity regression
function, $f_{s}$. The correlation, $\rho$, will be optimized before
inference. As before, we assume that the additive noise terms,
$\varepsilon_{ip}$ and $\varepsilon_{is}$, are independent and drawn from
identical, centered Gaussian distributions. Now, suppose we want to predict
our targets, $\bm{f}_{*}$, at the high fidelity level. We can determine that
the joint distribution of our high and low fidelity observations and our
targets has covariance

$\displaystyle\Sigma\ =\
\begin{bmatrix}\rho^{2}K^{s}_{pp}+K^{\delta}_{pp}+\sigma^{2}I&\rho
K^{s}_{ps}&\rho^{2}K^{s}_{p*}+K^{\delta}_{p*}\\\\[5.0pt] \rho
K^{s}_{sp}&K^{s}_{ss}+\sigma^{2}I&\rho K^{s}_{s*}\\\\[5.0pt]
\rho^{2}K^{s}_{*p}+K^{\delta}_{*p}&\rho
K^{s}_{*s}&\rho^{2}K^{s}_{**}+K^{\delta}_{**}\\\\[5.0pt] \end{bmatrix}.$ (8)

where $\sigma^{2}$ is the variance of the noise distribution. The superscript
of each covariance matrix, K, indicates which Gaussian process prior
corresponds to the kernel function used in its construction, and the
subscripts indicate the pairs of features that are compared. Bayesian
inference yields an analytical posterior distribution, analogous to Eq. (4).

###  II.4 Multitask Method

“Multitasking” refers to a category of methods which consider several Gaussian
process regression tasks and assume some relationship between these tasks. As
with $\Delta$-learning or multifidelity fusion, multitasking can allow us to
incorporate a larger dataset into our inference problem without just pooling
the data into one observation vector. For instance, we can define a regression
problem trained on CCSD(T) data as well as one trained on each DFA we
consider. By modeling correlation between regression functions, we can use the
DFA regression tasks to support prediction in the CCSD(T) task.

We can apply multitask regression to problems with $m$ tasks where, for each
task $j\in\\{1,\ldots,m\\}$, we have observational data
$\bm{Y}_{j}\in\mathbb{R}^{n}$. In practice, the method can be generalized, so
that observation sets of different sizes are used for different tasks.
Symmetric variants of multitask modeling treat all tasks equally. For
instance, we may assume that all the data are modeled by a single Gaussian
process whose covariance kernel is the Kronecker product of one kernel which
relates the tasks and another which models relationships between data from the
same task.Bonilla, Chai, and Williams (2008)

In many cases, however, the goal is to perform regression for a primary task,
and all other tasks may be considered secondary tasks which provide data that
may be useful in learning the primary task. For our own setting, CCSD(T) is a
natural choice of data for the primary regression task because we would like
to make predictions which match this method in accuracy. Each set of DFT
predictions made using a different DFA could inform a secondary regression
task which supports the primary task. Thus, we consider an asymmetrical model
Leen, Peltonen, and Kaski (2012) in which secondary tasks are related through
the primary task.

Suppose $\bm{Y}_{p}\in\mathbb{R}^{n_{p}}$ is the observation data for the
primary task and $\bm{Y}_{s_{j}}\in\mathbb{R}^{n_{j}}$ is the data for the
$s_{j}^{th}$ secondary task. As before, we assume that each task has its own
regression function

$\displaystyle\bm{Y}_{ip}\ $ $\displaystyle=$ $\displaystyle\
f_{p}(\bm{X}_{i})\ +\ \varepsilon_{ip},$ (9) $\displaystyle\bm{Y}_{is_{j}}\ $
$\displaystyle=$ $\displaystyle\ f_{s_{j}}(\bm{X}_{i})\ +\
\varepsilon_{is_{j}}\qquad\forall j=1,\dots,m.$ (10)

where each noise term is drawn independently from a centered Gaussian
distribution with variance $\sigma^{2}$. We now model a shared structure in
the regression functions based on the primary function:

$\displaystyle f_{s_{j}}(\bm{X}_{i})\ $ $\displaystyle=\
\rho_{s_{j}}f_{p}(\bm{X}_{i})\ +\ \delta_{s_{j}}(\bm{X}_{i})\qquad\forall
j=1,\dots,m.$ (11)

The specific component, $\delta_{s_{j}}$, of the secondary function serves to
“explain away” behavior that is not captured by the shared component, $f_{p}$.
Leen, Peltonen, and Kaski (2012) Both the correlation parameter,
$\rho_{s_{j}}$ and the specific component aim to mitigate negative
transfer—learning behavior from secondary tasks that is not representative of
the primary task. We make the prior assumption that $f_{p}$ and
$\delta_{s_{j}}$ for $j=1,\dots,m$ follow Gaussian process distributions, each
with a mean and kernel function which we will specify.

As was the case in subsections II.1 through II.3, the first step toward
inference with the multitask method is to determine the joint distribution on
observations and targets. Note that the specific component of each task is
modeled as independent of the others. Then, the covariance between
observations due to the specific component is a block matrix with zeros
corresponding to the covariance of specific functions for different tasks.
Call this block matrix $K^{\delta}$. Let $\bm{R}$ be a diagonal matrix with
dimensions equal to the number of observations, and let each diagonal element
be the correlation hyperparameter of the corresponding task. The covariance of
the joint distribution is

$\displaystyle\Sigma\ =\
\begin{bmatrix}K^{p}_{pp}+\sigma^{2}I&K^{p}_{ps}\bm{R}&K^{p}_{p*}\\\\[5.0pt]
\bm{R}K^{p}_{sp}&\bm{R}K^{p}_{ss}\bm{R}+K^{\delta}+\sigma^{2}I&\bm{R}K^{p}_{s*}\\\\[5.0pt]
K^{p}_{*p}&K^{p}_{*s}\bm{R}&K^{p}_{**}\\\\[5.0pt] \end{bmatrix}.$ (12)

The superscript $p$ indicates that a matrix was created with the primary
kernel, and the matrix subscripts indicate the features that are compared.
Inference proceeds as described in Eq. (4).

###  II.5 Beyond Two Tasks

Asymmetric multitasking assumes a relationship between regression
functions—given by Eq. (11)—which is structurally reminiscent of the
multifidelity fusion model—Eq. (6). Both assume that the regression function
of one dataset is related to the regression function of another by a scaling
parameter, $\rho$, and an additive disparity function, $\delta$.

For only one secondary model, $s_{1}$, the relationships between the
regression functions posited by both methods are mathematically equivalent.
However, asymmetric multitasking diverges more significantly from the
multifidelity method when more than two models are considered. In the
multifidelity approach, model hierarchy is reflected in a nested relationship
between the functions,

$\displaystyle f_{p}(\bm{X})\ $ $\displaystyle=$ $\displaystyle\
\rho_{s_{1}}f_{s_{1}}(\bm{X})\ +\ \delta_{s_{1}}(\bm{X}),$ $\displaystyle
f_{s_{1}}(\bm{X})\ $ $\displaystyle=$ $\displaystyle\
\rho_{s_{2}}f_{s_{2}}(\bm{X})\ +\ \delta_{s_{2}}(\bm{X}),$
$\displaystyle\vdots$ $\displaystyle f_{s_{j-1}}(\bm{X})\ $ $\displaystyle=$
$\displaystyle\ \rho_{s_{j}}f_{s_{j}}(\bm{X})\ +\ \delta_{s_{j}}(\bm{X}),$

whereas in the asymmetric approach, all functions share the same relationship
to the primary function

$\displaystyle f_{s_{1}}(\bm{X})\ $ $\displaystyle=$ $\displaystyle\
\rho_{s_{1}}f_{p}(\bm{X})\ +\ \delta_{s_{1}}(\bm{X}),$ $\displaystyle\vdots$
$\displaystyle f_{s_{j}}(\bm{X})\ $ $\displaystyle=$ $\displaystyle\
\rho_{s_{j}}f_{p}(\bm{X})\ +\ \delta_{s_{j}}(\bm{X}).$

Thus, compared to multifidelity fusion, the asymmetric multitasking approach
more naturally accommodates a reasonably large number of regression tasks and
relates each closely to the primary task of interest. Multitasking also holds
a similar advantage over $\Delta$-learning. Before the notion of a difference
can be imposed it is necessary to choose some order between the source data
sets. Such an order may be arbitrary, for example, if our data sets are
generated by multiple DFAs with similar structure. In Section III, we will
describe different compositions of a training data set which the multitask
framework allows us to test.

##  III  Heterogeneous Training Data

In our numerical examples, we train inference models using data for multiple
molecular configurations computed with multiple quantum chemistry methods. The
multitask method leaves many questions about the particulars of the training
data set up to the researcher. For instance, available data must be
apportioned to different tasks. One of our goals is to demonstrate that by
training with data from multiple quantum chemistry methods, we can construct
an inference model with prediction accuracy comparable to CCSD(T) with less
computational expense than if we trained a single task model on CCSD(T) alone.
Thus, our CCSD(T) data is the primary target of our regression problem and a
natural choice to inform the primary task. We let each secondary task
correspond to a different DFA. Examples with one up to four secondary tasks
will be considered.

Another decision to be made is whether the training data sets for different
tasks should include the same molecular configurations. In other words, we
must determine whether or not to train with predictions for the same quantity
by different methods. If our target for inference is the CCSD(T) prediction of
energy for a system of water molecules, we may also choose to include DFT
predictions of this energy in the training set for a secondary task. The
relative benefits of constructing training sets with different degrees of
intersection between tasks is under-explored in literature. Our work seeks to
understand the impact of training set design on model performance.

We consider a range of possible training set compositions here, and Fig. 2
provides two example diagrams of our system for distinguishing between
different constructions. We assign the molecular configurations to three sub
data sets—Core (C), Additional (A), and Target (T). Dots in the figure
indicate when configurations are available at a given level of theory. Black
indicates availability to the training set, and green indicates membership of
the testing set. The C set contains all configurations for which CCSD(T)
training data is available, the A set contains configurations for which we
have DFT data but no CCSD(T) data, and the T set contains configurations that
we use to test our model’s ability to reproduce CCSD(T) level accuracy.

---
Figure 2:  Example data set structures. C, A, and T are sets of molecular
configurations. C contains all configurations that the primary task, CCSD(T),
is trained on. T contains the Target configurations; our goal is to estimate
their properties with the accuracy of CCSD(T). The testing set is colored
green. All secondary tasks are highlighted in gray. Configurations in set A
can only be used for training secondary tasks. The vertical alignment of dots
indicates that multiple tasks have access to the same molecular
configurations.

In our numerical results, we always train models with the CCSD(T) data from
the C set and consider various compositions of the secondary sets. All
secondary tasks are highlighted in gray. We identify different training set
compositions by the letters of the sets of molecule which are included in our
secondary training data. For example, the top diagram in Fig. 2 shows an A
training set because the only DFT predictions available to train the model
come from the Additional set. By contrast, the lower diagram shows a CAT
training set because dots are present in the slots for DFT predictions for
Core, Additional, and Target molecular configurations.

Furthermore, the vertical alignment of dots indicate when some of the same
molecular configurations are present in the training sets for multiple tasks.
In the CAT example in Fig. 2, the dots indicate that the secondary tasks share
$20\%$ of the C and A sets while the rest of the molecules in these sets are
assigned to only one task each. In Section IV, we will show results for the C,
CT, A, AT, CA, and CAT training set configurations. Fig. 6 shows the impact of
varying the fraction of molecules that the secondary tasks have in common.

Our testing framework is designed so that models trained with DFT Target data
can readily be compared with equivalent models that do not train on this data.
Since the inference model can make Target predictions one at a time, for each
Target molecule we train a new model using the C set, A set, and the DFT
prediction for only that Target molecule. While in practice we could take the
much more efficient approach of training on the entire DFT Target set at once,
our approach enables us to contrast each CAT model to a CA model with a
training set which only differs by one molecule. Consequently, we can isolate
the effect of the DFT level Target data on prediction performance.

In Section IV, we will examine the performance of the $\Delta$ and Multitask
methods, informed by CCSD(T) and DFT predictions, and compare them to basic
Gaussian process regression, informed only by CCSD(T).

##  IV Numerical Examples

###  IV.1 Numerical Setup

We consider two case studies. In the first, the quantity of interest is the
three-body (3-b) interaction energy of randomly selected water trimer
configurations. The 3-b interaction term is isolated by subtraction of all 1-b
and 2-b energies from the full energy of the trimer. We selected water because
it is a simple molecule with unique properties which make it crucial for life
and a popular subject for research. Chen  _et al._ (2017); Dasgupta  _et al._
(2021); Gillan, Alfè, and Michaelides (2016) The second case study considers
the prediction of the ionization potential of small organic molecules.
Ionization potential is known to be challenging for DFT to capture accurately,
so it is insightful to test the ability of inference models to predict this
quantity.

####  IV.1.1 Data Generation

For our first test case—inference of 3-b interaction energy—5986 water trimer
configurations and corresponding CCSD(T) calculations were randomly selected
from 45,332 configurations available from  https://github.com/jmbowma/q-AQUA.
Yu  _et al._ (2022) Density functional theory data was generated for 2992 of
these configurations using Psi4. Smith  _et al._ (2020) We employed the
PBEPerdew, Burke, and Ernzerhof (1996) and SCANSun, Ruzsinszky, and Perdew
(2015) functional approximations. The former was selected for its popularity
in materials science research and the latter because its description of
intermediate range dispersion has been found to contribute to accurate
prediction of water energy differences Chen  _et al._ (2017); Dasgupta  _et
al._ (2021) Following recent work on dispersion corrections, Price, Bryenton,
and Johnson (2021) we paired the PBE approximation with a D3 correction using
finite damping and the SCAN approximation with a D3 correction using zero
damping. Both DFT and CCSD(T) calculations employed a standard counterpoise
correction to alleviate basis set superposition error. Huang, Braams, and
Bowman (2006)

The second test case is informed by the ionization potential of 3,165
molecular configurations, combinations of 479 different molecules and 7
configuration options. Duan, Fisher, and Kulik (2023) Configurations are
subselected Duan  _et al._ (2020) from the ANI-1 dataset Smith, Isayev, and
Roitberg (2017) and contain only elements from the set $\\{H,C,N,O\\}$. We
compute Duan  _et al._ (2021) ionization potential for these configurations
using CCSD(T) as well as secondary tasks corresponding to data from four
density functionals: PBE0_DH, Brémond and Adamo (2011) PBE0, Adamo and Barone
(1999) PBE, Perdew, Burke, and Ernzerhof (1996) and BLYP. Lee, Yang, and Parr
(1988)

Quantification of cost in our numerical results refers to the cost of data set
generation. The cost of inference is negligible next to quantum chemistry
calculations. We base our cost model for each method on the average runtime of
Psi4 computations for ten representative systems. For more details about the
cost model, see Appendix A.

####  IV.1.2 Features

We require features which can distinguish between molecular systems to serve
as inputs, $\bm{X}_{i}$, to our Gaussian process models. Kernel functions use
features to model the relationships between quantities of interest. As noted
in Section II.1, in our experiments we consider Smooth Overlap of Atomic
Positions (SOAP) features, which have been successfully used to fit
interatomic potentials in a Gaussian process based approach. Bartók  _et al._
(2010); Bartók, Kondor, and Csányi (2013); Bartók and Csányi (2015) SOAP
features are constructed based on the neighborhood of each atom in a system,
so challenges arise when two systems contain different numbers of constituent
atoms of elements. Multiple strategies have been proposed for constructing
“global” features to capture entire systems. De  _et al._ (2016) The simplest
approach is to take an average over the SOAP features for each atom in the
system. Since this approach is computationally efficient and works well in
practice, we apply it in our test cases. We chose SOAP parameters following
experimentation and the advice of literature. Deringer  _et al._ (2021); Musil
_et al._ (2021) More details are available in Appendices B and C.

####  IV.1.3 Hyperparameter Optimization

Squared exponential kernels are used within all Gaussian process priors in our
experiments. Abrahamsen (1997) These kernels introduce mean, variance ($v$),
and length scale ($\ell$) hyperparameters in addition to the multitask
correlation ($\rho$). We use the mean of our training data for each task as
that task’s prior mean, calculated at training time. To estimate the remaining
parameters, we construct an optimization data set containing $150$ randomly
selected molecular configurations which we subsequently exclude from our
training and testing data sets. We assume here that we have access to
predictions of our quantities of interest from each quantum chemistry method
for each of these $150$ molecular configurations, but when this assumption
does not hold, a conventional GP regression model may be used to fill in
missing data.

We estimate correlation parameters, $\rho$, using Pearson’s correlation
coefficient between each secondary task optimization data set and the
corresponding primary task data. Popular belief in literature holds that
correlation parameters require a detailed optimization procedure, Dral,
Zubatiuk, and Xue (2023) but we find that these easy to obtain estimates work
reasonably well in practice. A more detailed optimization is left for future
work. We estimate $\ell\in\mathbb{R}$ by maximizing its log-marginal
likelihood. This approach is easily extended to multifidelity and multitask
cases by treating $Y_{ij}-\rho_{j}Y_{i1}$ as the $i^{th}$ observation.
Forrester and A. Sóbester (2007) For the organic molecules example, we also
use a maximum likelihood estimate for $v$ while for the water trimer example,
we find better results from using an estimate of $v$ which minimizes mean
absolute prediction error.

###  IV.2 Three-Body Interaction Energy of Water Trimer

####  IV.2.1 Multitask Performance

---
Figure 3:  Water Trimer Case. Comparison of the MAE of the multitask method
for a range of training set constructions (scatter points) to the MAE of a GP
regression model trained only with CCSD(T) data (black line). Predictions by
the multitask method are informed by CCSD(T) and SCAN. The legends of the
plots indicate the structure of the secondary training set according to the
system presented by Fig. 2. All MAE reported are the average results obtained
from six randomly chosen training sets.

We first compare the performance of multitask models with only one secondary
task to the performance of a traditional GP—or single task—model. Fig. 3 shows
the error incurred by several implementations of each model when predicting
the three-body (3-b) energy of water trimer configurations. We used the SCAN
functional to generate the secondary task data, and Pearson’s correlation
coefficient between the primary and secondary tasks is $0.997$, indicating
strong positive correlation. The three subplots of Fig. 3 cover six secondary
training set configurations (C, CT, A, AT, CA, and CAT), labeled within the
legends.

All models were tested on a Target set of 320 molecular configurations. We let
the size of the Core set, $N_{c}$, take on values in
$\\{5,10,20,40,80,160,320\\}$ for both inference methods and additionally
consider $N_{c}\in\\{640,1280,2560\\}$ for the single task method. The
subplots of Fig. 3 also show seven clear groupings of data points with steep
slope. This effect is a consequence of the high cost of CCSD(T) relative to
DFT. Each of the seven groupings corresponds to a different $N_{c}$. Given the
size of the Core set, we let the number of configurations in the A set range
from $0.5N_{c}$ to $6N_{c}$. Each point on the plot represents a model with a
different training set size. The x axis indicates the cost of the training set
while the y axis gives the average mean absolute error for six random
assignments of molecular configurations to the C, A, and T sets. We provide
details of our cost model in Appendix A.

|
---|---
(a) | (b)
Figure 4:  Water Trimer Case. Effect of including Core and Additional data in
secondary training sets. Each indigo point represents a comparison between two
multitask models where one model uses a training set which is a subset of the
other. In (a), this comparison is between A and CA training configurations,
and in (b), we compare AT and CAT configurations. For each pair of models
compared, both the primary training data and the secondary data from the A set
are identical. The $y=x$ line is plotted in green for reference.

For every training set configuration except C, the multitask method
outperforms the single task reference, shown by a black line. Including T data
from DFT in the training set leads to better results than corresponding cases
without T in every example. In particular, training with secondary T data can
enable the multitask method to exhibit an accuracy comparable to a GP model
that is an order of magnitude more expensive.

While the C results demonstrate no significant improvement on the single task
model, both A and CA results show steady improvement as the size of the A set
increases. This behavior suggests that for the multitask method to be
beneficial, secondary training data must cover molecular space that is not
included in the CCSD(T) training data. We mention that it is possible that the
behavior of C models would improve if we consider task dependent features or
different relationships between task regression functions. We do not
investigate these alternatives further here; instead, we focus on the usual
setting where secondary DFT data for additional systems is cheap to obtain.

CT and CAT models are the best-performing multitask models, demonstrating very
similar MAE per cost to each other. Furthermore, unlike the AT model, the CAT
model does not show improvement as the size of the A set increases. It seems
that for our setting, the combination of secondary C and T training data is so
valuable to the model that any benefit from A data is dwarfed. We can
hypothesize that when the model has access to data for the Core molecules from
the primary and secondary tasks, it implicitly learns information about the
difference between the tasks. Given secondary data on the Target, it can
leverage the captured relationship between tasks to make an accurate
prediction for the Target on the primary level.

Fig. 4 provides more insight to the relative roles of the training data sets
and to the inherent flexibility of the multitask framework. While Fig. 3
demonstrates that multitask models can be successful when trained on many
different training data configurations from the CAT framework, it does not
directly compare the performance of the A and CA models. Such a comparison
will illuminate the resilience of the multitask model to the loss of points in
the training set which are common to primary and secondary tasks. We emphasize
that this resilience is an advantage of the multitask method in comparison
with $\Delta$-learning—where no training is feasible without access to
secondary C data. Each point on subplot (a) of Fig. 4 compares the average MAE
of two multitask models: one with a CA secondary training set, the other with
an identical A set but no secondary C data. As before, we are inferring the
3-b energy of water trimer configurations with only one secondary task. The
closeness of the points to the reference $y=x$ line shows that there is
minimal effect on accuracy when the secondary C data is excluded from the
training set. In fact, the median improvement in accuracy from using the CA
training model instead of an A training model model is $0.0007$ eV. By
comparison, the median improvement in accuracy from combining DFT data from
the A set with CCSD(T) training data—that is, switching from a single task
approach to a multitask method trained on an A set—is $0.004$ eV.

Subplot (b) of Fig. 4 also compares models trained with and without C
secondary data in the case where all of these models also train on secondary
Target data. While the AT models perform well, there is a penalty to losing
access to secondary C training data when secondary Target data is also
available. The median improvement in accuracy when switching from an AT model
to a CAT model is $0.002$ eV, a third of the $0.006$ eV median improvement
incurred from switching from a single task approach to an AT multitask model
with the same CCSD(T) training cost. As we remarked previously, when the model
trains on secondary C data, it likely gains information about the difference
between the primary and secondary tasks. This information is most useful to
the model when it also has access to secondary predictions for the Target
molecules.

####  IV.2.2 Comparison with $\Delta$-learning

| CCSD(T) [eV] | SCAN [eV] | PBE [eV]
---|---|---|---
absolute mean | $0.0205$ | $0.0223$ | $0.0225$
absolute median | $0.0126$ | $0.0144$ | $0.0155$
variance | $0.0009$ | $0.0010$ | $0.0011$

| CCSD(T)$-$SCAN [eV] | CCSD(T)$-$PBE [eV]
---|---|---
absolute mean | $0.0021$ | $0.0037$
absolute median | $0.0013$ | $0.0018$
variance | $1.1\times 10^{-5}$ | $3.2\times 10^{-5}$

Table 1:  Water Trimer case. Summary statistics for absolute quantum chemistry
predictions and for differences.
---
(a) Estimation of CCSD(T)-PBE
(b) Estimation of CCSD(T)-SCAN
Figure 5:  Water Trimer Case. The MAE of different inference methods: GP
regression (line), multitask regression trained with a CAT secondary set
(orange), two level $\Delta$-learning (indigo), three level $\Delta$-learning
(light blue), and multitask inference on $\Delta$s (gold). Reported
performance is the average from six draws of each training set.

$\Delta$-learning is trained on differences between predictions of our
quantity of interest while the single and multitask methods are trained on
absolute predictions of the quantity of interest. Table 1 reports summary
statistics for these data sets. In this case, the mean and variance of the
differences between quantum chemistry methods are smaller by at least an order
of magnitude than those of any set of quantum chemistry predictions. Since
$\Delta$-learning is trained on differences to predict differences, it is
expected that this method will produce smaller mean absolute error than a
method trained on absolutes to predict absolutes.

By modifying our choice of tasks, we can design a multitask model which trains
on and predicts differences. We call this the multitask $\Delta$ option. Fig.
5 compares this option to the multitask method on absolutes as well
$\Delta$-learning. For each inference method, the average mean absolute error
obtained for six random assignments of the training data is plotted for each
of the same data set sizes which were considered in our previous experiments.
The CAT secondary training set structure is used for the multitask absolute
models as well as the multitask difference models. As expected, the multitask
absolute model is not as successful as those which train on and learn
differences.

In subplot (a) of Fig. 5, the object of estimation is the difference between
CCSD(T) and PBE computations for each molecular system. For the multitask
$\Delta$ option, we let the primary task be exactly this difference. The
secondary task is the difference between SCAN and PBE computations which has
correlation coefficient $0.782$ with the primary task data. Results of the
multitask $\Delta$ option are plotted in gold. The basic $\Delta$-learning
model is plotted in indigo and uses only two levels of theory, CCSD(T) and
PBE, to predict their difference. Note that these models are the single task
variant of the multitask $\Delta$ approach—they are trained only on the
primary task. Therefore, Subplot (a) provides another example where the
multitask method outperforms its single task counterpart, in this case with
correlation parameter $\rho=0.782$. The turquoise points represent a
$\Delta$-learning model trained with three levels of theory: CCSD(T), SCAN,
and PBE. This model is the sum of a pair of two level models where the top
level of the first difference—between SCAN and PBE—is the bottom level of the
second—between CCSD(T) and SCAN. Given sufficiently large training sets for
the lower levels of theory, this method is comparable with the multitask
$\Delta$ option. The multitask $\Delta$ option outperforms most other models,
particularly when either C or A is relatively small.

Subplot (b) of Fig. 5 shows the performance of the same set of methods when
the primary and secondary tasks demonstrate a slight negative correlation. In
this case, the extra data from the secondary task does not contribute to
improved inference of the primary task. Here, we estimate the difference
between CCSD(T) and SCAN predictions. For all descriptions accompanying
Subplot (a), we reverse the roles of SCAN and PBE for Subplot (b). The major
difference between the cases is the correlation coefficient between data for
CCSD(T)-SCAN and PBE-SCAN: $\rho=-0.271$. The multitask models shown in
Subplot (b) demonstrate the limitations of choosing Pearson’s coefficient as
the correlation parameter. The multitask $\Delta$ option and three level
$\Delta$-learning do not perform as well as two level $\Delta$-learning in
this setting. This behavior suggests that the multitask method requires a
stronger correlation between tasks to outperform the single task method.
Further, the subplot shows that the multitask $\Delta$ option outperforms
three level $\Delta$-learning for models trained on larger amounts of DFT
data. Thus, three level $\Delta$-learning seems to suffer more from
problematic data than the Multitask $\Delta$ option. Importantly,
$\Delta$-learning has no correlation parameter and cannot adjust the
contribution of each difference model. Though we have left a detailed
optimization of the Multitask $\Delta$ option’s correlation parameter for
future work, we observe that $\rho=0$ can always be chosen. Therefore, the
fully optimized multitask $\Delta$ option will always perform at least as well
as two level $\Delta$-learning. Consequently, the multitask framework makes it
possible to benefit from useful secondary datasets which two level
$\Delta$-learning cannot accommodate while offering more protection against
uncorrelated data than three level $\Delta$ models.

In practice, we may have more than two tasks available to support inference,
so we can further contrast the methods discussed in this section by
considering how they scale to more tasks. Suppose that alongside PBE and
CCSD(T) we want to include $n$ _additional_ DFT functionals
$y_{1},\dots,y_{n}$. Then, $\Delta$-learning requires us to put all quantum
chemistry methods in some order—with CCSD(T) first and PBE last and the
$y_{j}$ in between. The difference between neighboring pairs of methods is
then learned. However, hierarchies in DFA accuracy can be difficult to predict
and differ between problem settings. This necessary ordering can thus become
rather arbitrary. By contrast, the multitask $\Delta$ option would incorporate
the $n$ new secondary tasks based on the differences to PBE, i.e.
$y_{j}-\text{PBE}$. Each of these secondary tasks can be directly related to
the primary task through a correlation parameter—which thus eliminates the
need to impose an arbitrary order.

###  IV.3 Ionization Potential of Small Organic Molecules

| PBE0_DH | PBE0 | PBE | BLYP
---|---|---|---|---
$\ \rho\ $ | $0.970$ | $0.967$ | $0.954$ | $0.955$
Table 2:  Organic Molecules Case. Pearson’s correlation coefficients between each secondary task and the primary task data set. The column headers give the DFA used to generate the secondary set. Fraction of shared molecules in CA: $0.0$ |
---|---
|
Fraction of shared molecules in CA: $0.5$ |
|
Fraction of shared molecules in CA: $1.0$ |
|
(a) Trained without DFT data from Target Set | (b) Trained with DFT data from Target Set
Figure 6:  Organic Molecules Case. For different numbers of tasks (indicated
by color), plots show MAE versus cost. Note that $10^{-0.3}\textrm{eV}$ is
approximately half an electronvolt. For each row, a given fraction of all
molecules in the C and A sets are used to generate training data for all
secondary levels, and the remaining molecules are split among the secondary
levels. The fractions are given at the top of each row. Models represented in
column (a) were trained without DFT data from the T set while models in column
(b) had access to this data.

We now turn our attention to the prediction of the ionization potential of
small organic molecules made up of elements in $\\{$C,N,O,H$\\}$. Our
secondary data for this example includes up to four tasks, generated with the
PBE, PBE0, PBE_DH, and BLYP functionals, respectively. Table 2 shows the
correlation parameter used to relate each of these secondary tasks to the
primary task data, generated with CCSD(T). As in the water trimer example, we
consider Core sets of sizes $\\{5,10,20,40,80,160,320\\}$. The Target set
contains 320 molecular configurations. We let the number of molecules in the A
set range from $1$ to $6$ times the number of molecules in the C set. As
Section III visualizes, once the C and A sets have been fixed, there is still
a design choice in allotting the molecular configurations to secondary tasks.
We consider three cases: (1) the molecules in both C and A are randomly
distributed among tasks, so no two tasks have access to training data for the
same molecule, (2) half of the molecules in these sets are common to all tasks
while the other half are distributed among the tasks with no redundancy, and
(3) all tasks have access to training data for all molecules in both C and A.
These cases correspond—in order, from the top down—to the three rows of Fig.
6.

In Fig. 6, as for previous results, the reported mean absolute error is the
average over six tests produced by randomly assigning molecules to the C, A,
and T sets. The left column of subplots compares results for multitask models
trained on CA data sets to single task models trained only on CCSD(T) while
the right column considers multitask models trained on CAT data. Translucent
stripes encompass the results for models with different total numbers of
tasks. For the two task case, the secondary task data is generated with PBE,
and for the three task case the secondary data comes from PBE and PBE0. In all
cases, secondary Target data is only generated with PBE. Recall that in our
testing framework, when we train with the T data set, we train a custom model
for each Target molecule. Specifically, the only T molecule in the training
set is the molecule which we will use the model to predict. We remark that in
practice it would not be necessary to train a new model for each molecule in
the Target set.

As we observed in the water trimer case, the multitask models consistently
perform better than single task models of comparable cost. Especially for
relatively small training sets, the multitask method shows an improvement in
performance when secondary Target data is included in the training set. Fig. 6
also informs us about the utility of including more than one secondary task in
a model. Note that for each number of tasks we consider, we implement models
with the same combinations of C, A, and T set sizes. In the top row, molecules
in the C and A sets are each assigned to exactly one secondary task,
regardless of the total number of secondary tasks. Therefore, in this row, for
each model implemented with two total tasks, we consider models of the same
cost implemented with three tasks and with five tasks. The results show that
for sufficiently large training sets, there is a benefit to training with
multiple tasks for a given budget. For smaller datasets, two tasks may lead to
performance as good or better than a larger number of tasks. In particular,
the top subplot in column (b) suggests that two task models perform
particularly well when given a relatively inexpensive training set of the CAT
structure. Note that, in each model, only one secondary task is given access
to T data. Thus, in this setting, the decrease in performance when new tasks
are added likely comes from siphoning off part of the computational budget to
train tasks which do not have access to T data.

We can also consider whether it is useful to train on more than two
predictions for the same molecular configuration. The middle row of Fig. 6
shows cases where each task has access to training data for the same 50% of
molecules in the C set as well as 50% of the molecules in the A set while the
remaining molecules are distributed one per task, and the final row shows
cases where all tasks have access to all data in the C and A sets. From these
subplots, we conclude that there is a benefit to training our model on
multiple secondary predictions for the same molecule.

##  V Conclusion

The present work has demonstrated that inference models see a general benefit
from leveraging all available data. Furthermore, use of the multitask
framework can reduce the cost of generating entirely new training data by an
order of magnitude and can also facilitate the opportunistic incorporation of
existing datasets of heterogeneous quality into the training set. Our work is
the first to explore a multitask GP regression model in a materials science
setting. Significantly, the multitask framework is not restricted to the
choice of GP inference, and future work can investigate training models with
neural networks and other tools for regression.

Unlike other machine learning methods for constructing surrogates from ab
initio data, the multitask model does not require all data sources be ordered
in a hierarchy of accuracy. In this work, the data which informed our primary
task—generated with CCSD(T)—was the most accurate in our training set, but we
were also able to leverage an arbitrary number of secondary data sets without
any information on their accuracy relative to each other. This feature of the
multitask method is particularly useful when employing datasets constructed
from different density functional approximations. We found that multitask
model performance improves as we increase the number of secondary sets in the
training data which include predictions for a molecule of interest. The
influence of each secondary dataset is controlled by a correlation parameter,
and we found that in most settings an estimate of Pearson’s correlation
coefficient from training data works well. Use of this estimate simplifies
implementation of the method. The multitask method also supports training on
differences, and with appropriate choice of correlation parameter, multitask
difference models can perform at least as well as and often better than
$\Delta$ models. To explore the framework’s flexibility, we demonstrated the
success of the multitask method on diverse training data set configurations.

###### Acknowledgements.

The authors would like to thank Yeongsu Cho, Chenru Duan, Dallas Foster, and
Heather Kulik for insightful discussions. MIT SuperCloud and Lincoln
Laboratory Supercomputing Center are acknowledged for providing the HPC
resources that have contributed to the research results reported within this
paper. Reuther  _et al._ (2018) This material is based upon work supported by
the Department of Energy, National Nuclear Security Administration PSAAP-III
program, under Award Number DE-NA0003965 as well as work supported by the
National Science Foundation Graduate Research Fellowship under Grant No.
1745302. This research was supported by the NCCR MARVEL, a National Centre of
Competence in Research, funded by the Swiss National Science Foundation (grant
number 205602).

## Author Contributions Statement

Katharine Fisher: Conceptualization (equal); Data Curation (equal);
Methodology (equal); Software (lead); Writing – original draft (lead); Writing
– review and editing (equal); Visualization (equal). Michael Herbst:
Conceptualization (equal); Data Curation (equal); Methodology (equal);
Supervision (equal); Writing – review and editing (equal); Visualization
(equal). Youssef Marzouk: Conceptualization (equal); Methodology (equal);
Supervision (equal); Writing – review and editing (equal); Visualization
(equal).

## Conflict of Interest Statement

The authors have no conflicts to disclose.

## Data Availability Statement

All data supporting this work are openly available. All code to reproduce
figures as well as all DFT predictions for the three-body energy of water
trimers can be found at https://doi.org/10.5281/zenodo.10387647 or
https://github.com/kefisher98/Multitask_GP/tree/Multitask. Water trimer
configurations and CCSD(T) predictions for three-body energy are available at
https://github.com/jmbowma/q-AQUA. The data that support the findings of this
study for the small organic molecules example are at
https://doi.org/10.5281/zenodo.10215421.

## References

  * Harding  _et al._ (2008) M. E. Harding, T. Metzroth, J. Gauss, and A. A. Auer, “Parallel calculation of ccsd and ccsd(t) analytic first and second derivatives,” Journal of Chemical Theory and Computation 4, 64–74 (2008).
  * Ying, Yu, and Ying (2019) L. Ying, J. Yu, and L. Ying, “Numerical methods for kohn–sham density functional theory,” Acta Numerica 28, 405–539 (2019).
  * Perdew and Schmidt (2001) J. Perdew and K. Schmidt, “Jacob’s ladder of density functional approximations for the exchange-correlation energy,” in  _AIP Conference Proceedings_ (2001) presented at AIP Conference Proceedings 577.
  * Goerigkab and Grimme (2011) L. Goerigkab and S. Grimme, “A thorough benchmark of density functional methods for general main group thermochemistry, kinetics, and noncovalent interactions,” Physical Chemistry Chemical Physics 13, 6670–6688 (2011).
  * Teale  _et al._ (2022) A. M. Teale, T. Helgaker, A. Savin, C. Adano, B. Aradi, A. V. Arbuznikov, P. Ayers, E. J. Baerends, V. Barone, P. Calaminici, E. Cances, E. A. Carter, P. K. Chattaraj, H. Chermette, I. Ciofini, T. D. Crawford, F. De Proft, J. Dobson, C. Draxl, T. Frauenheim, E. Fromager, P. Fuentealba, L. Gagliardi, G. Galli, J. Gao, P. Geerlings, N. Gidopoulos, P. M. W. Gill, P. Gori-Giorgi, A. Görling, T. Gould, S. Grimme, O. Gritsenko, H. J. A. Jensen, E. R. Johnson, R. O. Jones, M. Kaupp, A. Koster, L. Kronik, A. I. Krylov, S. Kvaal, A. Laestadius, M. P. Levy, M. Lewin, S. Liu, P.-F. Loos, N. T. Maitra, F. Neese, J. Perdew, K. Pernal, P. Pernot, P. Piecuch, E. Rebolini, L. Reining, P. Romaniello, A. Ruzsinszky, D. Salahub, M. Scheffler, P. Schwerdtfeger, V. N. Staroverov, J. Sun, E. Tellgren, D. J. Tozer, S. Trickey, C. A. Ullrich, A. Vela, G. Vignale, T. A. Wesolowski, X. Xu, and W. Yang, “Dft exchange: Sharing perspectives on the workhorse of quantum chemistry and materials science,” Physical Chemistry Chemical Physics (2022), 10.1039/d2cp02827a.
  * Smith  _et al._ (2019) J. S. Smith, B. T. Nebgen, R. Zubatyuk, N. Lubbers, C. Devereux, K. Barros, S. Tretiak, O. Isayev, and A. E. Roitberg, “Approaching coupled cluster accuracy with a general-purpose neural network potential through transfer learning,” Nature Communications 10 (2019), 10.1038/s41467-019-10827-4.
  * Dral  _et al._ (2020) P. O. Dral, A. Owens, A. Dral, and G. Csányi, “Hierarchical machine learning of potential energy surfaces,” The Journal of Chemical Physics 152 (2020), 10.1063/5.0006498.
  * Goodlett, Turney, and Schaefer (2023) S. M. Goodlett, J. M. Turney, and H. F. Schaefer, “Comparison of multifidelity machine learning models for potential energy surfaces,” The Journal of Chemical Physics 159 (2023), 10.1063/5.0158919.
  * Zaverkin  _et al._ (2023) V. Zaverkin, D. Holzmüller, L. Bonfirraro, and J. Kästner, “Transfer learning for chemically accurate interatomic neural network potentials,” Physical Chemistry Chemical Physics 25, 5383–5396 (2023).
  * Curtarolo  _et al._ (2012) S. Curtarolo, W. Setyawan, G. L. Hart, M. Jahnatek, R. V. Chepulskii, R. H. Taylor, S. Wang, J. Xue, K. Yang, O. Levy, M. J. Mehl, H. T. Stokes, D. O. Demchenko, and D. Morgan, “Aflow: An automatic framework for high-throughput materials discovery,” Computational Materials Science 58, 218–226 (2012).
  * Jain  _et al._ (2011) A. Jain, G. Hautier, C. J. Moore, S. Ping Ong, C. C. Fischer, T. Mueller, K. A. Persson, and G. Ceder, “A high-throughput infrastructure for density functional theory calculations,” Computational Materials Science 50, 2295–2310 (2011).
  * Huber  _et al._ (2020) S. P. Huber, S. Zoupanos, M. Uhrin, L. Talirz, L. Kahle, R. Häuselmann, D. Gresch, T. Müller, A. V. Yakutovich, C. W. Andersen, F. F. Ramirez, C. S. Adorf, F. Gargiulo, S. Kumbhar, E. Passaro, C. Johnston, A. Merkys, A. Cepellotti, N. Mounet, N. Marzari, B. Kozinsky, and G. Pizzi, “Aiida 1.0, a scalable computational infrastructure for automated reproducible workflows and data provenance,” Scientific Data 7 (2020), 10.1038/s41597-020-00638-4.
  * Cancès  _et al._ (2022) E. Cancès, A. Levitt, Y. Maday, and C. Yang, “Numerical methods for kohn–sham models: Discretization, algorithms, and error analysis,” in  _Density Functional Theory: Modeling, Mathematical Analysis, Computational Methods, and Applications_ (Springer, 2022) pp. 333–400.
  * Herbst and Levitt (2020) M. F. Herbst and A. Levitt, “Black-box inhomogeneous preconditioning for self-consistent field iterations in density functional theory,” Journal of Physics: Condensed Matter (2020), 10.1088/1361-648x/abcbdb.
  * Herbst and Levitt (2022) M. F. Herbst and A. Levitt, “A robust and efficient line search for self-consistent field iterations,” Journal of Computational Physics 459, 111127 (2022).
  * Cancès  _et al._ (2023) E. Cancès, M. F. Herbst, G. Kemlin, A. Levitt, and B. Stamm, “Numerical stability and efficiency of response property calculations in density functional theory,” Letters in Mathematical Physics 113 (2023), 10.1007/s11005-023-01645-3.
  * Ruddigkeit  _et al._ (2012) L. Ruddigkeit, R. van Deursen, L. C. Blum, and J.-L. Reymond, “Enumeration of 166 billion organic small molecules in the chemical universe database gdb-17,” Journal of Chemical Information and Modeling 52, 2864–2875 (2012), pMID: 23088335, https://doi.org/10.1021/ci300415d .
  * Ramakrishnan  _et al._ (2014) R. Ramakrishnan, P. O. Dral, M. Rupp, and O. A. von Lilienfeld, “Quantum chemistry structures and properties of 134 kilo molecules,” Scientific Data 1 (2014).
  * Chanussot*  _et al._ (2021) L. Chanussot*, A. Das*, S. Goyal*, T. Lavril*, M. Shuaibi*, M. Riviere, K. Tran, J. Heras-Domingo, C. Ho, W. Hu, A. Palizhati, A. Sriram, B. Wood, J. Yoon, D. Parikh, C. L. Zitnick, and Z. Ulissi, “Open catalyst 2020 (oc20) dataset and community challenges,” ACS Catalysis (2021), 10.1021/acscatal.0c04525.
  * Smith, Isayev, and Roitberg (2017) J. Smith, O. Isayev, and A. Roitberg, “A data set of 20 million calculated off-equilibrium conformations for organic molecules,” Scientific Data 4 (2017).
  * Bonilla, Chai, and Williams (2008) E. Bonilla, K. Chai, and C. Williams, “Multi-task gaussian process prediction,” in  _Advances in neural information processing systems_ , edited by J. Platt, D. Koller, Y. Singer, and S. Roweis (MIT Press, Cambridge, Massachusetts, 2008) pp. 153–160.
  * Leen, Peltonen, and Kaski (2012) G. Leen, J. Peltonen, and S. Kaski, “Focused multi-task learning in a gaussian process framework,” Machine Learning 1-2, 157–182 (2012).
  * Pilania, Gubernatis, and Lookman (2017) G. Pilania, J. Gubernatis, and T. Lookman, “Multi-fidelity machine learning models for accurate bandgap predictions of solids,” Computational Materials Science 129, 156–162 (2017).
  * Batra  _et al._ (2019) R. Batra, G. Pilania, B. Uberuaga, and R. Ramprasad, “Multifidelity information fusion with machine learning: A case study of dopant formation energies in hafnia,” ACS Applied Materials & Inference 11, 24906–24918 (2019).
  * Kennedy and O’Hagan (2000) M. Kennedy and A. O’Hagan, “Predicting the output from a complex computer code when fast approximations are available,” Biometrika 87, 1–13 (2000).
  * Ramakrishnan  _et al._ (2015) R. Ramakrishnan, P. O. Dral, M. Rupp, and O. A. von Lilienfeld, “Big data meets quantum chemistry approximations: The $\Delta$-machine learning approach,” Journal of Chemical Theory and Computation 11, 2087–2096 (2015).
  * Dral, Zubatiuk, and Xue (2023) P. O. Dral, T. Zubatiuk, and B.-X. Xue, “Learning from multiple quantum chemical methods: $\delta$-learning, transfer learning, co-kriging, and beyond,” in  _Quantum Chemistry in the Age of Machine Learning_ (Elsevier, 2023) pp. 491–507.
  * Vinod, Kleinekathöfer, and Zaspel (2023) V. Vinod, U. Kleinekathöfer, and P. Zaspel, “Optimized multifidelity machine learning for quantum chemistry,” (2023), arXiv:2312.05661 [physics.chem-ph] .
  * Towell and Shavlik (1994) G. G. Towell and J. W. Shavlik, “Knowledge-based artificial neural networks,” Artificial Intelligence 70, 119–165 (1994).
  * Fu (1989) L.-M. Fu, “Integration of neural heuristics into knowledge-based inference,” in  _International 1989 Joint Conference on Neural Networks_ (1989) pp. 606 vol.2–.
  * Bartlett, Montanari, and Rakhlin (2021) P. L. Bartlett, A. Montanari, and A. Rakhlin, “Deep learning: a statistical viewpoint,” Acta Numerica 30, 87–201 (2021).
  * Lotfi  _et al._ (2022) S. Lotfi, M. Finzi, S. Kapoor, A. Potapczynski, M. Goldblum, and A. G. Wilson, “Pac-bayes compression bounds so tight that they can explain generalization,” (2022), arXiv:2211.13609 [cs.LG] .
  * Rasmussen and Williams (2006) C. E. Rasmussen and C. K. I. Williams,  _Gaussian Processes for Machine Learning_ (the MIT Press, 2006).
  * Quinonero-Candela and Rasmussen (2005) J. Quinonero-Candela and C. Rasmussen, “A unifying view of sparse approximate gaussian process regression,” Journal of Machine Learning Research 6, 1939–1959 (2005).
  * Liu  _et al._ (2018) H. Liu, Y.-S. Ong, X. Shen, and J. Cai, “When gaussian process meets big data: A review of scalable gps,”  (2018).
  * Wilson, Dann, and Nickisch (2015) A. G. Wilson, C. Dann, and H. Nickisch, “Thoughts on massively scalable gaussian processes,” CoRR abs/1511.01870 (2015), 1511.01870 .
  * Cole, Christianson, and Gramacy (2021) D. A. Cole, R. B. Christianson, and R. B. Gramacy, “Locally induced gaussian processes for large-scale simulation experiments,” Statistics and Computing 31, 1573–1375 (2021).
  * Abrahamsen (1997) P. Abrahamsen, “A review of gaussian random fields and correlation functions,”  (1997), 10.13140/RG.2.2.23937.20325.
  * Bartók  _et al._ (2010) 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), 10.1103/physrevlett.104.136403.
  * Bartók, Kondor, and Csányi (2013) A. P. Bartók, R. Kondor, and G. Csányi, “On representing chemical environments,” Physical Review B 87, 184115 (2013).
  * Bartók and Csányi (2015) A. P. Bartók and G. Csányi, “Gaussian approximation potentials: A brief tutorial introduction,” International Journal of Quantum Chemistry 115, 1051–1057 (2015).
  * Chen  _et al._ (2017) M. Chen, H.-Y. Ko, R. Remsing, M. Andrade, B. Santra, Z. Sun, A. Selloni, R. Car, M. Klein, J. Perdew, and W. Xifan, “Ab initio theory and modeling of water,” Proceedings of the National Academy of Sciences 114, 201712499 (2017).
  * Dasgupta  _et al._ (2021) S. Dasgupta, E. Lambros, J. Perdew, and F. Paesani, “Elevating density functional theory to chemical accuracy for water simulations through a density-corrected many-body formalism,”  (2021), 10.33774/chemrxiv-2021-hstgf.
  * Gillan, Alfè, and Michaelides (2016) M. J. Gillan, D. Alfè, and A. Michaelides, “Perspective: How good is DFT for water?” The Journal of Chemical Physics 144, 130901 (2016), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.4944633/13696839/130901_1_online.pdf .
  * Yu  _et al._ (2022) Q. Yu, C. Qu, P. L. Houston, R. Conte, A. Nandi, and J. M. Bowman, “q-aqua: A many-body ccsd(t) water potential, including four-body interactions, demonstrates the quantum nature of water from clusters to the liquid phase,” The Journal of Physical Chemistry Letters 13, 5068–5074 (2022), pMID: 35652912.
  * Smith  _et al._ (2020) D. G. A. Smith, L. A. Burns, A. C. Simmonett, R. M. Parrish, M. C. Schieber, R. Galvelis, P. Kraus, H. Kruse, R. Di Remigio, A. Alenaizan, A. M. James, S. Lehtola, J. P. Misiewicz, M. Scheurer, R. A. Shaw, J. B. Schriber, Y. Xie, Z. L. Glick, D. A. Sirianni, J. S. O’Brien, J. M. Waldrop, A. Kumar, E. G. Hohenstein, B. P. Pritchard, B. R. Brooks, I. Schaefer, Henry F., A. Y. Sokolov, K. Patkowski, I. DePrince, A. Eugene, U. Bozkaya, R. A. King, F. A. Evangelista, J. M. Turney, T. D. Crawford, and C. D. Sherrill, “PSI4 1.4: Open-source software for high-throughput quantum chemistry,” The Journal of Chemical Physics 152, 184108 (2020), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0006002/16684807/184108_1_online.pdf .
  * Perdew, Burke, and Ernzerhof (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized gradient approximation made simple,” Phys. Rev. Lett. 77, 3865–3868 (1996).
  * Sun, Ruzsinszky, and Perdew (2015) J. Sun, A. Ruzsinszky, and J. P. Perdew, “Strongly constrained and appropriately normed semilocal density functional,” Phys. Rev. Lett. 115, 036402 (2015).
  * Price, Bryenton, and Johnson (2021) A. J. A. Price, K. R. Bryenton, and E. R. Johnson, “Requirements for an accurate dispersion-corrected density functional,” The Journal of Chemical Physics 154, 230902 (2021).
  * Huang, Braams, and Bowman (2006) X. Huang, B. J. Braams, and J. M. Bowman, “Ab initio potential energy and dipole moment surfaces of (h2o)2,” The Journal of Physical Chemistry A 110, 445–451 (2006), pMID: 16405316, https://doi.org/10.1021/jp053583d .
  * Duan, Fisher, and Kulik (2023) C. Duan, K. Fisher, and H. Kulik, “kefisher98/IP_EA_deltaSCF: Ionization Potential, Electron Affinity, and Delta SCF for Small Organic Molecules,”  (2023).
  * Duan  _et al._ (2020) C. Duan, F. Liu, A. Nandy, and H. J. Kulik, “Data-driven approaches can overcome the cost–accuracy trade-off in multireference diagnostics,” Journal of Chemical Theory and Computation 16, 4373–4387 (2020).
  * Duan  _et al._ (2021) C. Duan, S. Chen, M. G. Taylor, F. Liu, and H. J. Kulik, “Machine learning to tame divergent density functional approximations: a new path to consensus materials design principles,” Chem. Sci. 12, 13021–13036 (2021).
  * Brémond and Adamo (2011) E. Brémond and C. Adamo, “Seeking for parameter-free double-hybrid functionals: The PBE0-DH model,” The Journal of Chemical Physics 135, 024106 (2011).
  * Adamo and Barone (1999) C. Adamo and V. Barone, “Toward reliable density functional methods without adjustable parameters: The PBE0 model,” J. Chem. Phys. 110, 6158–6170 (1999).
  * Lee, Yang, and Parr (1988) C. Lee, W. Yang, and R. G. Parr, “Development of the colle-salvetti correlation-energy formula into a functional of the electron density,” Phys. Rev. B 37, 785–789 (1988).
  * De  _et al._ (2016) 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, 13754–13769 (2016).
  * Deringer  _et al._ (2021) V. Deringer, A. Bartók, N. Bernstein, D. Wilkins, M. Ceriotti, and G. Csányi, “Gaussian process regression for materials and modelling,” Chemical Reviews 121, 10073–10041 (2021).
  * Musil  _et al._ (2021) F. Musil, A. Grisafi, A. P. Bartók, C. Ortner, G. Csányi, and M. Ceriotti, “Physics-inspired structural representations for molecules and materials,” Chemical Reviews 121, 9759–9815 (2021).
  * Forrester and A. Sóbester (2007) A. I. J. Forrester and A. J. K. A. Sóbester, “Multi-fidelity optimization via surrogate modelling,” Proceedings of Royal Society A 463, 3251–3269 (2007).
  * Reuther  _et al._ (2018) A. Reuther, J. Kepner, C. Byun, S. Samsi, W. Arcand, D. Bestor, B. Bergeron, V. Gadepally, M. Houle, M. Hubbell, M. Jones, A. Klein, L. Milechin, J. Mullen, A. Prout, A. Rosa, C. Yee, and P. Michaleas, “Interactive supercomputing on 40,000 cores for machine learning and data analysis,” in  _2018 IEEE High Performance extreme Computing Conference (HPEC)_ (IEEE, 2018) pp. 1–6.
  * Civalleri  _et al._ (2012) B. Civalleri, D. Presti, R. Dovesi, and A. Savin, “On choosing the best density functional approximation,” in  _Uncertainty Quantification in Multiscale Materials Modeling_, edited by M. Springborg (RSC Publishing, 2012) Chap. 6, pp. 168–185.
  * Himanen  _et al._ (2020) 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, 106949 (2020).

##  Appendix A Training Data Cost Model

In our numerical results, we estimate the training cost as

$\displaystyle\text{model training cost}\ =\ \sum_{j}^{m}n_{j}\hat{c}_{j}$

where $n_{j}$ is the number of training systems from the $j^{th}$ level of
theory and $\hat{c}_{j}$ is an estimate of prediction cost for that level of
theory. We produce one set of estimates for $\hat{c}_{1},\dots,\hat{c}_{m}$
for the water trimer example and a different set of estimates for the small
organic molecules example.

For the water trimers example, we randomly selected ten molecular
configurations from our data set as the basis for our cost model. Using Psi4,
we computed the 3-b interaction energy at the level of SCAN and CCSD(T) for
each configuration. We applied a counterpoise correction to ameliorate basis
set superposition error. All computations for the water cost model where
performed on physically the same Intel Xeon Platinum 8260 node of MIT
SuperCloud (https://supercloud.mit.edu/systems-and-software). All computations
for a given molecular configuration where performed on the same basis set
subasis set sucore. For this case study, we estimate that a CCSD(T)
calculation incurs $37$ times the cost of a DFT calculation based on the
average prediction time of these methods on our example molecular
configurations.

Several configurations in the small organic molecules data set are larger than
the water trimer configurations. The median number of electrons in this data
set is $40$, the minimum is $14$, and the maximum is $50$. Due to the expense
of ionization potential calculations, we base our cost model on energy
calculations for ten randomly sampled configurations with $36$ electrons. Note
that CCSD(T) scales as $N^{7}$ where $N$ is the number of electrons in a
system while DFT scales as $N^{3}$. Therefore, an estimate of the ratio
between the cost of CCSD(T) and the cost of DFT based on $36$ electrons will
be smaller than the actual ratio for most molecular configurations in our data
set. Furthermore, the prediction of ionization potential was performed using
${\Delta}CC$ (delta coupled cluster). Thus, while our reported results (Fig.
6) estimate that CCSD(T) calculations demand $244$ times the expense of DFT,
the true ratio between the cost of CCSD(T) and DFT is even larger.

##  Appendix B SOAP parameters

To fix reasonable values for parameters required by SOAP ($r_{cut}$,
$\sigma_{atom}$, $n_{max}$, $l_{max}$), we refer to established conventions as
well as experimentation. The cutoff radius, $r_{cut}$ controls the size of
each local neighborhood, and a small value may lead to lost geometric insight.
Unfortunately, surveys on feature construction Deringer  _et al._ (2021);
Musil  _et al._ (2021) hold that a large cutoff radius does not necessarily
provide proportionate insight: increasing radii larger than 6-8 Å is rarely if
ever useful. Our choice of $\sigma_{atom}$ also influences our model of each
atom’s neighborhood: it determines the lengthscale of the Gaussian shapes
located on each of the surrounding atoms. The larger our choice of
$\sigma_{atom}$, the more chance that Gaussian tails will slip passed the
$r_{cut}$ border. Current literature Deringer  _et al._ (2021) indicates that
the best practice when working with the first three rows of the periodic table
is to use $\sigma_{atom}=0.3\ $\mathrm{\SIUnitSymbolAngstrom}$$ for systems
with hydrogen and $\sigma_{atom}=0.5\ $\mathrm{\SIUnitSymbolAngstrom}$$ for
systems without. Finally, we choose the parameters $n_{max}$ and $l_{max}$ to
control the size of our expansions of the local neighborhoods. Generally
Deringer  _et al._ (2021), choosing $n_{max}=12$ with $l_{max}=6$ is
sufficient for high accuracy. In practice, the best values for $n_{max}$ and
$l_{max}$ will be sensitive to the choice of radial basis set.

Figure 7: Example of inference results for a range of SOAP parameters. The
colorbar gives the mean absolute error of predictions made for the ionization
potential of small organic molecules. Both $r_{cut}$ and $\sigma_{atom}$ have
units of $\mathrm{\SIUnitSymbolAngstrom}$.

We selected reasonable candidate values for each SOAP parameter and tested the
performance of the resultant features in Gaussian process regression. All
features for this work were calculated using the DScribe Python package.
Himanen  _et al._ (2020) Fig. 7 shows the mean absolute error (MAE) obtained
with different SOAP parameters when GP regression trained with PBE level data
is used to predict ionization potential (IP). Tests on CCSD(T) and PBE0 data
produced similar results. Our tests include all combinations of
$r_{cut}\in\\{3,4,5\\}$ , $\sigma_{atom}\in\\{0.3,0.4,0.5\\}$,
$l_{max}\in\\{2,4,6,8\\}$, and $n_{max}\in\\{6,8,10,12\\}$. Units for
$r_{cut}$ and $\sigma_{atom}$ are $\mathrm{\SIUnitSymbolAngstrom}$. We find
that our accuracy is comparable to probabilistic predictions of IP found in
literature. De  _et al._ (2016) Of the SOAP parameters, the choice of
$\sigma_{atom}$ has the largest impact on predictive performance.

For our numerical experiments in this work, we fix $\sigma_{atom}=0.4\
$\mathrm{\SIUnitSymbolAngstrom}$$ because this choice demonstrated the best
overall performance across our tests. The other parameters are chosen to
balance reasonable accuracy and cost. We fix both $l_{max}$ and $n_{max}$ at
8, and for the small organic molecules example, we use $r_{cut}=4\
$\mathrm{\SIUnitSymbolAngstrom}$$. For water trimers, we use $r_{cut}=10\
$\mathrm{\SIUnitSymbolAngstrom}$$ to capture all three component water
molecules.

##  Appendix C Global Features

Figure 8: Correlation between different methods of constructing global
features when calculating the distance between molecule pairs. De  _et al._
(2016)

When using SOAP to compare two molecular systems, we must make modeling
choices to standardize our representation of entire systems. SOAP features are
constructed based on the neighborhood of each atom in a system, so challenges
arise when two systems contain different numbers of constituent atoms of
elements. There are multiple strategies for constructing “global” features to
capture entire systems. De  _et al._ (2016) The simplest approach is to
average the SOAP features for each atom in the system. It is possible the loss
of information from averaging may make it challenging to distinguish between
the features of similar molecular systems. A variation to this approach
averages the local features corresponding to each element contained in a
molecular system. De  _et al._ (2016) If we wish to compare systems A and B,
and system B contains more elements than A, we insert SOAP features
corresponding to isolated atoms of those excess elements in the representation
of A. This model suggests that an atom of such an element does not interact
with the rest of the system. A Regularized entropy match (REMatch) strategy
has also been proposed for constructing global features. De  _et al._ (2016)
This approach solves a regularized optimization problem to find the match
between the sets of local features for each molecular system which maximizes
information entropy.

With the same set of small organic molecules which we use to test SOAP
parameters, we investigate how well the global features constructed by
different approaches agree. For pairs of molecules, we construct global
features using each of the three approaches outlined in the previous
paragraph. We represent difference between each pair by squaring an inner
product of their features. This procedure is the same as using a polynomial
kernel of degree 2. By comparing the output of the kernel for different
featurization strategies, we determine whether the strategies generally agree
about which molecules are similar.

Fig. 8 shows the correlation in kernel outputs corresponding to each pair of
molecules for different global featurization strategies. Consistently, in our
tests, we find strong correlation (Pearson’s coefficient $>0.99)$ between the
REMatch features and features computed by averaging all local representations.
These two strategies are also both positively correlated with the “Average by
Species” approach which produces an averaged feature for each element in the
system and inserts isolated atoms to represent elements not included in the
system, but these correlations are not as strong as that between the REMatch
and totally averaged features. We may see this result because the REMatch and
totally averaged features share the same dimension, $n$, whereas features that
are constructed by species specific averaging have dimension of $2n$ to $4n$
in the cases we considered. Because averaging features is much more
computationally efficient than constructing REMatch features, we use this
approach in our numerical experiments.


* [14] Higson, N., and Kasparov, G. E-theory and KK-theory for groups which act properly and isometrically on hilbert space.  Inventiones mathematicae 144, 1 (2001), 23–74.
  * [15] Holm, T., Jorgensen, P., and Rouquier, R. Triangulated Categories.  Cambridge University Press, 2010.
  * [16] Lafforgue, V. La conjecture de Baum-Connes à coefficients pour les groupes hyperboliques.  J. Noncommutative Geometry (2012), 1–197.
  * [17] Lemeux, F., and Tarrago, P. Free wreath product quantum groups: the monoidal category, approximation properties and free probability.  Journal of Functional Analysis 270, 10 (2016), 3828–3883.
  * [18] Martos, R. The Baum-Connes conjecture for Quantum Groups. Stability properties and ${K}$-theory computations.  PhD thesis, Université Paris VII, 2018.
  * [19] Martos, R. The Baum-Connes property for a quantum (semi-)direct product.  J. Noncommut. Geom. 13, 4 (2020), 1295–1357.
  * [20] Martos, R. On the torsion-freeness property for divisible discrete quantum subgroups.  Submitted (2021).
  * [21] Meyer, R. Homological algebra in bivariant K-theory and other triangulated categories II.  Tbilisi Mathematical Journal 1 (2008), 165–210.
  * [22] Meyer, R., and Nest, R. The Baum–Connes conjecture via localisation of categories.  Topology 45, 2 (2006), 209–259.
  * [23] Meyer, R., and Nest, R. An analogue of the Baum-Connes isomorphism for coactions of compact groups.  Mathematica Scandinavica 100, 2 (2007), 301–316.
  * [24] Nest, R., and Voigt, C. Equivariant Poincaré duality for quantum group actions.  Journal of Functional Analysis 258, 5 (2010), 1466–1503.
  * [25] Nikshych, D., and Vainerman, L. Finite quantum groupoids and their applications.  In New directions in Hopf algebras, vol. 43 of Math. Sci. Res. Inst. Publ. Cambridge Univ. Press, Cambridge, 2002, pp. 211–262.
  * [26] Rosenberg, J., and Schochet, C. The Künneth theorem and the universal coefficient theorem for equivariant ${K}$-theory and ${K}{K}$-theory.  American Mathematical Soc. 348 (1986).
  * [27] Tahara, K.-i. On the Second Cohomology Groups of Semidirect Products.  Math. Z. 129 (1972), 365–379.
  * [28] Thomsen, K. The universal property of equivariant ${K}{K}$-theory.  J. reine angew. Math. 504 (1998), 55–71.
  * [29] Vergnioux, R., and Voigt, C. The ${K}$-theory of free quantum groups.  Mathematische Annalen 357, 1 (2013), 355–400.
  * [30] Voigt, C. The Baum–Connes conjecture for free orthogonal quantum groups.  Advances in Mathematics 227, 5 (2011), 1873–1913.
  * [31] Voigt, C. On the structure of quantum automorphism groups.  Journal für die reine und angewandte Mathematik (Crelles Journal) (2015).
  * [32] Wang, S. Tensor products and crossed products of compact quantum groups.  Proc. London Math. Soc. 71 (1995), 695–720.

R. Martos, Department of Mathematical Sciences, University of Copenhagen,
Denmark.

E-mail address<EMAIL_ADDRESS>