.gifpng.pngconvert gif:#1 png: .gif

# Kernel based quantum machine learning at record rate:
Many-body distribution functionals as compact representations

Danish Khan  Department of Chemistry, University of Toronto, St. George
Campus, Toronto, ON, Canada  Vector Institute for Artificial Intelligence,
Toronto, ON, M5S 1M1, Canada  Stefan Heinen  Vector Institute for Artificial
Intelligence, Toronto, ON, M5S 1M1, Canada  O. Anatole von Lilienfeld
<EMAIL_ADDRESS>Vector Institute for Artificial
Intelligence, Toronto, ON, M5S 1M1, Canada  Departments of Chemistry,
Materials Science and Engineering, and Physics, University of Toronto, St.
George Campus, Toronto, ON, Canada  Machine Learning Group, Technische
Universität Berlin and Institute for the Foundations of Learning and Data,
10587 Berlin, Germany

###### Abstract

The feature vector mapping used to represent chemical systems is a key factor
governing the superior data-efficiency of kernel based quantum machine
learning (QML) models applicable throughout chemical compound space.
Unfortunately, the most accurate representations require a high dimensional
feature mapping, thereby imposing a considerable computational burden on model
training and use. We introduce compact yet accurate, linear scaling QML
representations based on atomic Gaussian many-body distribution functionals
(MBDF), and their derivatives. Weighted density functions (DF) of MBDF values
are used as global representations which are constant in size, i.e. invariant
with respect to the number of atoms. We report predictive performance and
training data efficiency that is competitive with state of the art for two
diverse datasets of organic molecules, QM9 and QMugs. Generalization
capability has been investigated for atomization energies, HOMO-LUMO
eigenvalues and gap, internal energies at 0 K, zero point vibrational
energies, dipole moment norm, static isotropic polarizability, and heat
capacity as encoded in QM9. MBDF based QM9 performance lowers the optimal
Pareto front spanned between sampling and training cost to compute node
minutes, effectively sampling chemical compound space with chemical accuracy
at a sampling rate of $\sim 48$ molecules per core second.

##  I Introduction

Modern data-driven statistical Machine Learning (ML) models have emerged as
powerful tools over the past decade for inferring quantum mechanical
observables throughout chemical compound space, without explicitly solving
electronic Schrödinger equations1; 2; 3. Similar success was obtained for ML
based interatomic potentials and force-fields4; 5; 6; 7; 8; 9 as well as
electronic structure modeling throughout Chemical Compound Space (CCS)10; 11.
For an entire set of extensive in-depth reviews on these and other related ML
applications, we refer the reader to the recent special issue in Chemical
Reviews12; 13. Various aspects in the development of ML model architecture and
training protocols have proven to be essential for data-efficiency. In
particular, the molecular representation is known to have a strong impact on
the performance of similarity based ML models, such as kernel ridge regression
(KRR)14; 15; 16. This is not surprising as the representation controls the
information about the systems, how its weighed and the consistency of these
Quantum Machine Learning (QML) models with ab-initio methods 17. These
representations are non-linear mappings of the atomistic systems to a suitable
Hilbert Space where a statistical regression model can easily be applied. The
Hilbert Space constraint applies due to the requirement of measuring
similarity in terms of an inner product18. These mappings should have some
desirable features of which the most important are i) uniqueness such that
systems with different properties necessarily possess different
representations14and ii) invariance with respect to transformations that leave
the target property invariant, such as global translations, rotations, and
atomic index permutations of the same chemical elements. Other desirable
features include iii) an analytical and continuous form of the representation
function, iv) differentiability with respect to nuclear coordinates, charges,
number of electrons, and number of atoms, v) as general as the Hamiltonian,
vi) computationally efficient evaluation, vii) compact or even constant size,
limiting the computational cost for larger systems19.
Due to their critical role, many representations have been introduced and
investigated within the context of atomistic simulation20; 21; 22; 23; 24; 25;
26; 27; 28; 29; 30; 31; 32. For recent comprehensive reviews, the reader is
referred to Refs. 33; 16. These representations can either describe the
molecule as a whole (global) or each atom (local or atomic) separately. For
the sake of brevity we restricted most of our comparative benchmarks within
this study to the following representations which are commonly used to train
QML models throughout CCS: Faber-Christensen-Huang-Lilienfeld (FCHL19)
representation34; 35, Smooth Overlap of Atomic Positions (SOAP)36, spectrum of
London and Axilrod–Teller–Muto potentials (SLATM)37, atom index sorted Coulomb
matrix (CM)1, and its vectorized form, Bag of Bonds (BOB)38. Other
representations/models tested are mentioned in the data and code section.
While these representations satisfy most of the aforementioned qualities, seen
the immense size of CCS, a more compact and scalable representation would
still be desirable. Formally, the number of degrees of freedom of any material
or molecule would prescribe usage of a $4M$ dimensional feature vector (3
spatial coordinates and one nuclear charge coordinate for each of the atoms).
However, all aforementioned representations when using optimal hyperparameters
require a higher dimensional feature vector mapping in order to be accurate
and training-data-efficient at regression tasks, and some (e.g. CM or BOB)
even scale quadratically with $M$. While verbosity facilitates the inclusion
of invariances, the $4M$ degrees of freedom suggest that the same performance
can be obtained using more compact representations. This is especially an
issue for kernel based ML models, where the size of the representation
directly affects the distance/kernel evaluation time18; 34. Although the
scaling for kernel inversion39 is larger ($\propto\mathcal{O}(n^{3})$ for
Cholesky solvers), for highly data-efficient (i.e. efficient in training data)
QML models it is the kernel generation and evaluation that consumes the most
compute time as demonstrated later. The kernel evaluation pre-factor becomes
even worse when using atomic (or local) representations in conjunction with a
suitable local kernel40. Obvious solutions by reducing finite local cutoffs
within the representation come at the expense of reducing the predictive
power, or, conversely, increasing training data needs. As such, a
computationally efficient yet accurate solution is desirable as was shown in
the discretization of the FCHL18 representation35; 34. Other solutions to this
problem such as sparse kernel models41 and the recently introduced conjugate
gradients based iterative approach for training kernel based models42 would
also be well complemented by a compact molecular representation.

Herein, we propose a methodology for generating representations that minimize
feature size. We use functionals of many-body distributions (MBDFs) and their
derivatives to encode any local chemical environment of any atom in an
integrated compact fashion. The representations thus generated preserve the
system’s various symmetries (translation, rotation, atom index invariance),
and can be tailored to the physical property of interest through scaling
functions. MBDFs are easily extendable to include higher order many-body
interactions with minimal increase in size. In the current formulation, while
including three-body interactions, MBDF scales as $5M$. We further tackle the
issue of storing this information in a manner that remains invariant to the
number of atoms in a molecule. We do this by generating a discretized density
(DF) function of MBDF values, and using it as a global molecular
representation. Using two diverse datasets the performance of MDBF is tested
against aforementioned SOAP and FCHL19, which are commonly used and state-of-
the-art, as well as SLATM, BOB, CM representations and a few other QML models
mentioned later. Lastly, we explore the bottleneck cross-over from kernel
evaluation to inversion.

##  II Theory and Discussion

###  II.1 Many-Body Distribution Functionals

We begin the discussion of our local representation using distribution
functions over the internal coordinates which can be constructed using the
atomic coordinates.

An analytical and continuous distribution over the pair-wise internal
coordinate defined as the inter-atomic distances (pair correlation function),
is easily built using Gaussian probability density functions (PDFs) centered
at each inter-atomic distance with respect to an atom $i$:

$\displaystyle\rho_{i}(r,R_{ij})=\frac{1}{\sqrt{2\pi\sigma_{r}^{2}}}\sum_{j\neq
i}^{M}Z_{j}\exp{\left(-{\frac{(r-R_{ij})^{2}}{\sigma_{r}}}\right)}$ (1)

where $\rho_{i}(r,R_{ij})$ is the normalized distribution with atom $i$ as the
origin, $\sigma_{r}$ is the Gaussian length-scale (or variance) parameter, $M$
denotes the total number of atoms in the system (or within a radial cutoff if
employed), $R_{ij}$ denotes the inter-atomic distance between atoms $i$ and
$j$ and $Z_{j}$ is the nuclear charge. Scaling by nuclear charges defines
elemental identities, but could also be done in other ways e.g. having
different length-scale parameters $\sigma_{r}$ for each unique chemical
element, or multiple dimensions such as period and group specifications 43. In
a similar fashion, a continuous distribution (triplet correlation function)
over the 3-body internal coordinate, inter-atomic angles $\theta$, is defined
as:

$\displaystyle\rho_{i}(\theta,\theta_{jik})=\frac{1}{\sqrt{2\pi\sigma_{\theta}^{2}}}\sum_{j\neq
i}^{M}\sum_{k\neq
j}^{M}Z_{j}Z_{k}~{}\exp{\left(-\frac{(\theta-\theta_{jik})^{2}}{\sigma_{\theta}}\right)},$
$\displaystyle\theta_{jik}=\cos^{-1}\frac{(\mathbf{R}_{i}-\mathbf{R}_{j})^{T}(\mathbf{R}_{i}-\mathbf{R}_{j})}{R_{ij}^{2}}$
(2)

where $\theta_{jik}$ is the inter-atomic angle centered on atom $i$. This can
be generalized to define a continuous distribution, or correlation function,
over any $m$-body internal coordinate $\tau$.
Such analytical distribution functions of the internal coordinates have been
used as descriptors of atomic environments in Atom Centered Symmetry
Functions20; 6 (ACSF), their weighted variants wACSFs23 and other similar
methods. Instead of using these distribution functions as atomic descriptors,
we define the MBDF representation as functionals of these $m$-body
distributions which leads to a more compact descriptor and allows inclusion of
higher order terms with minimal change in size (a single scalar for each
$m$-body term). With the 2- and 3-body distributions defined above, each atom
can then be represented by the two zeroth-order functionals:

$\displaystyle
F_{0}^{(2)}[i]=\int_{0}^{r_{c}}dr~{}g_{0}(r,R_{ij})~{}\rho_{i}(r,R_{ij})$ (3)
$\displaystyle
F_{0}^{(3)}[i]=\int_{0}^{\pi}d\theta~{}g_{0}(\theta,R_{ij},R_{jk},R_{ik})~{}\rho_{i}(\theta,\theta_{jik})$
(4)

where $r_{c}$ denotes the radial cut-off distance, and $g_{0}(r,R_{ij})$ and
$g_{0}(\theta,R_{ij},R_{jk},R_{kj})$ denote 2- and 3-body weighting functions.
Note that when the weighting functions $g_{0}(r,R_{ij})$ and
$g_{0}(\theta,R_{ij},R_{jk},R_{kj})$ correspond to suitable 2 and 3-body
potentials, the functionals $F_{0}^{(2)}$, $F_{0}^{(3)}$ become the average of
the corresponding 2, 3-body inter-atomic interactions weighted by the pair and
triplet correlation functions $\rho_{i}(r,R_{ij})$,
$\rho_{i}(\theta,\theta_{jik})$, respectively. These functionals then form a
coarse approximation to the average 2, 3-body interactions experienced by a
chemical species. Furthermore, we exploit the advantage of using the
infinitely differentiable Gaussian PDFs to define higher order functionals
such as:

$\displaystyle
F_{1}^{(2)}[i]=\int_{0}^{r_{c}}dr~{}g_{1}(r,R_{ij})~{}\frac{\partial\rho_{i}}{\partial
r}(r,R_{ij})$ (5) $\displaystyle
F_{1}^{(3)}[i]=\int_{0}^{\pi}d\theta~{}g_{1}(\theta,R_{ij},R_{jk},R_{ik})~{}\frac{\partial\rho_{i}}{\partial\theta}(\theta,\theta_{jik})$
(6)

with potentially different weighting functions $g_{1}(r)$,
$g_{1}(\theta,R_{ij},R_{jk},R_{kj})$. The derivative functionals are useful
since the functional of any arbitrary distribution is not unique. These
functionals also encode the change in the $m$-body distribution in an atom’s
local neighborhood and have not been used in previous works involving internal
coordinate distribution functions. For any $n$-th derivative of the 2-body
distribution we can define the functional:

$\displaystyle
F_{n}^{(2)}[i]=\int_{0}^{r_{c}}dr~{}g_{n}(r)~{}\partial_{r}^{n}\rho_{i}(r,R_{ij}),$
$\displaystyle\partial_{r}^{n}\rho_{i}(r)=\frac{\partial^{n}\rho_{i}}{\partial
r^{n}}(r,R_{ij})$ (7)

where $g_{n}(r)$ is, again, a suitable radial weighting function. Generalizing
this to all internal coordinates, a functional $F_{n}^{(m)}[i]$ can be defined
over the $n$-th derivative of any $m$-body distribution function centered at
atom $i$:

$F_{n}^{(m)}[i]=\int_{0}^{\tau_{c}}d\tau~{}g_{n}(\tau)~{}\partial_{\tau}^{n}\rho_{i}(\tau,\tau_{i_{1}i_{2}..i_{m}}),\\\
\partial_{\tau}^{n}\rho_{i}(\tau,\tau_{i_{1}i_{2}..i_{m}})=\sum_{i_{1}<i_{2}..<i_{m}}^{M}H_{n}(\tau-\tau_{i_{1}i_{2}..i_{m}})\\\
\times\mathcal{N}(\tau_{i_{1}i_{2}..i_{m}},\sigma_{\tau}^{2})\prod_{j=i_{1}}^{i_{m}}Z_{j}$
(8)

where $\tau$ denotes the $m$-body internal coordinate, $\rho_{i}$ is the
$m$-body distribution function w.r.t atom $i$, $g_{n}$ is the weighting
function for the $n$-th derivative of $\rho_{i}$, $H_{n}$ denotes the Hermite
polynomial of degree $n$ and
$\mathcal{N}(\tau_{i_{1}i_{2}..i_{m}},\sigma_{\tau}^{2})$ denotes the
normalized Gaussian distribution. The Hermite polynomials arise due to the use
of Gaussian PDFs and allow convenient computation of $n$ derivatives of the
distribution function at any point $\tau$.

We note here that an alternative way to describe a (bounded) distribution in a
compact form is through a moment expansion of the form:

$\displaystyle
G^{(m)}[i]=\int_{0}^{\tau_{c}}d\tau~{}(\tau-\tau_{i_{1}i_{2}..i_{m}})^{m}g_{m}(\tau)\rho_{i}(\tau,\tau_{i_{1}i_{2}..i_{m}})$
(9)

where $G^{(m)}[i]$ denotes the $m$-th moment of the distribution centered at
atom $i$. The set of $m$ moments $G^{(m)}[i]$ would then form the local
representation of the atom $i$. These moments can also be evaluated by placing
a set of Gaussians (or any basis functions) on each atom $i$, and then
evaluating the moments of the atomic density $\mathbf{\rho}_{i}(\mathbf{r})$
w.r.t each atomic position within a radial cutoff $\mathbf{r}_{c}$:

$G^{(m)}[i]=\int_{0}^{\mathbf{r}_{c}}d\mathbf{r}~{}|\mathbf{r}|^{m}g_{m}(|\mathbf{r}|)\rho_{i}(\mathbf{r}),\\\
\rho_{i}(\mathbf{r})=\frac{1}{\sqrt{2\pi\sigma^{2}}}\sum_{j}\exp{\left(-\frac{||\mathbf{r}-\mathbf{R}_{ij}||_{2}^{2}}{2\sigma^{2}}\right)}$
(10)

where $|.|$ is any metric. This form has the advantage of being independent of
many-body orders and the computational cost of evaluating these moments scales
solely with the number of atoms within the cutoff radius $\mathbf{r}_{c}$. The
integral can be simplified in spherical polar coordinates by expanding the
density $\rho_{i}(\mathbf{r})$ in a basis set composed of spherical harmonics
$Y_{l}^{m^{\prime}}$ and orthogonal radial functions $u_{n}$ (which is a
common practice36):

$G^{(m)}[i]=\int_{0}^{r_{c}}\int_{0}^{2\pi}\int_{0}^{\pi}dr~{}d\theta~{}d\phi~{}r^{m+2}~{}g_{m}(r)~{}\sin(\phi)\\\
\times\sum_{nlm^{\prime}}c_{nlm^{\prime}}^{i}~{}u_{n}(r)~{}Y_{l}^{m^{\prime}}(\theta,\phi),\\\
c_{nlm^{\prime}}^{i}=\braket{\rho_{i}(\mathbf{r})}{u_{n}(r)~{}Y_{l}^{m^{\prime}}(\theta,\phi)}$
(11)

Throughout our work we use the derivative formalism from eq. (8) since our
numerical results indicated superior performance than the moments expansion in
the internal coordinates (eq. (9)). However, the moment expansion in eq. (10),
being independent of many-body terms, offers a promising alternative and could
be the subject of a future work.

We have tested multiple weighting functions to identify the best combination.
In particular, since these functionals correspond to correlation function
averages of $m$-body interactions, we have tested the Harmonic, Morse44,
Lennard-Jones45 potentials and simple decaying functions (power laws,
exponential, gaussian decays and their combinations) for 2-body terms. For
3-body terms we have tested Cosine Harmonic, Axilrod-Teller46, Stillinger-
Weber47 potentials and a scaled Fourier series. Through trial and error
combined with cross-validated hyper-parameter optimization, we have identified
the following 5 suitable functionals corresponding to 2 and 3-body
distributions:

$\displaystyle F_{0}^{(2)}[i]=$
$\displaystyle\int_{0}^{r_{c}}dr~{}\left[e^{-\eta
r}-\frac{\sqrt{2}\pi}{10(r+1)^{3}}\right]~{}\rho_{i}(r,R_{ij})$ (12)
$\displaystyle F_{1}^{(2)}[i]=$
$\displaystyle\frac{\sqrt{2}\pi}{10}\int_{0}^{r_{c}}dr~{}\frac{\partial_{r}\rho_{i}(r,R_{ij})}{(r+1)^{6}}$
(13) $\displaystyle F_{2}^{(2)}[i]=$
$\displaystyle\int_{0}^{r_{c}}dr~{}e^{-\alpha
r}~{}\partial_{r}^{2}\rho_{i}(r,R_{ij})$ (14) $\displaystyle F_{0}^{(3)}[i]=$
$\displaystyle\int_{0}^{\pi}d\theta~{}\frac{\sum_{n=0}^{3}a_{n}\cos(n\theta)}{(R_{ij}R_{jk}R_{ik})^{4}}~{}\rho_{i}(\theta,\theta_{jik})$
(15)
$\displaystyle\begin{split}F_{1}^{(3)}[i]=&\int_{0}^{\pi}d\theta~{}\left[\frac{1+\cos(\theta)\cos(\theta_{kji})\cos(\theta_{ikj})}{(R_{ij}R_{jk}R_{ik})^{4}}\right]\\\
&\times\partial_{\theta}\rho_{i}(\theta,\theta_{jik})\end{split}$ (16)

where $\eta$, $\alpha$, $a_{n}$ and the various power laws are all
hyperparameters of the representation. Note the scaling functions of MBDF
contributions $F_{0}^{(2)},F_{1}^{(2)},F_{0}^{(3)},F_{1}^{(3)}$ being
respectively reminiscent of Buckingham type-potential, softened London
dispersion potential, Fourier series scaled by Lennard-Jones repulsion and
Axilrod-Teller-Muto potential scaled by Lennard-Jones repulsion. The specific
reason as to why this particular selection of weighting functions has proven
advantageous will be subject of future research.

Figure 1 shows the effect of each functional on the learning capacity of MBDF
for the task of predicting atomization energies of the QM948 dataset.

Figure 1: MBDF based QML learning curves using concatenated increasingly
higher order and body functionals (Eqs. (12-16)). Mean absolute error (MAE) of
predicting atomization energies of the QM9 48 dataset is shown as a function
of training set size $N$.

It is apparent that both the derivative and many-body terms improve the
learning capacity, albeit by different magnitudes.

Throughout our testing on the QM749, QM948 and QMugs50 datasets, we have found
these 2, 3-body functionals to be sufficient at discriminating between all of
the molecular structures. Cases where the 2-body information does not suffice
include homometric pairs, as already discussed many years ago 14. And even 3-
and 4-body information does not suffice for some cases, as recently discussed
in Ref. 51. We note here that, whenever necessary, arbitrarily higher order
derivative and many-body information could also be included in MBDFs at
minimal increase in size, i.e. one additional term per order (See eq. (8)). In
particular, we believe that the inclusion of the 4-body term as a functional
of the dihedrals could further improve the learning capacity for
conformational isomers. Inclusion of 4-body information has been shown to
result in further improvements of learning curves 52. Also note that the size
of MBDF is invariant to the cutoffs used which can be raised to arbitrarily
higher values while employing a suitable long-range functional (hence
increasing the farsightedness of the representation) without affecting the
kernel evaluation cost. We further note that other weighting functions and
response terms could also be useful for QML models of physical observables
such as dipole moments, vibrational frequencies, heat capacities etc.

###  II.2 Density of functionals

MBDF is a local representation and its size scales linearly with the number of
atoms in the system. In order to eliminate this scaling we can transform to
the frequency space of MBDF functional values. The frequencies can be
evaluated by normalizing the functional values to lie within an arbitrary
range and then using, for e.g., kernel density estimation53. For a finite set
of MBDF functional values $\\{X_{i}\\}_{i=1}^{5M}$ over the range $[a,b]$, a
”smooth histogram” of their frequencies can be constructed by placing a set of
kernel functions $K$ at each point $X_{i}$:

$\displaystyle f(u)=\frac{1}{5M}~{}\sum_{i=1}^{5M}~{}K(u,X_{i})$ (17)

where $f(u)$ gives the density of the samples at any point
$u~{}\epsilon~{}[a,b]$. If the set $\\{X_{i}\\}_{i=1}^{5M}$ are the MBDF
functional values for any molecule, their distribution density can be
evaluated using eq. (17). The density $f(u)$ can then be used as a global
molecular representation whose size is independent of the number of MBDF
functional values, and number of atoms by extension. The (dis-)similarity
between two molecules $A$ and $B$ can be evaluated as, e.g., the l2-distance:

$\displaystyle d(A,B)^{2}=\int_{-c}^{c}du~{}|f_{A}(u)-f_{B}(u)|^{2}$ (18)

where $[-c,c]$ is the normaliztion range chosen as $[-10,10]$ in our work. The
form of the density function used in our work is:

$\displaystyle
f_{A}(u)=\frac{1}{5M}~{}\sum_{i=1}^{5M}~{}\frac{\sqrt{|X_{i}|}}{\sqrt{2\pi}\sigma_{b}}~{}\exp{\left(-\frac{(u-X_{i})^{2}}{2\sigma_{b}^{2}}\right)}$
(19)

where $X_{i}$ are MBDF functionals for molecule $A$, $\sigma_{b}$ is the
variance of the Gaussian function and is a hyperparameter (also called
bandwidth). Comparing with eq. (17), the function $K$ is defined as:

$\displaystyle
K(u,X_{i}):=\frac{\sqrt{|X_{i}|}}{\sqrt{2\pi}\sigma_{b}}~{}\exp{\left(-\frac{(u-X_{i})^{2}}{2\sigma_{b}^{2}}\right)}$
(20)

Note that this is a divergence54 but not a kernel function because it is
asymmetric:

$\displaystyle K(x,y)=$
$\displaystyle\frac{\sqrt{|y|}}{\sqrt{2\pi}\sigma_{b}}~{}\exp{\left(-\frac{(x-y)^{2}}{2\sigma_{b}^{2}}\right)}$
$\displaystyle\begin{split}\neq&\frac{\sqrt{|x|}}{\sqrt{2\pi}\sigma_{b}}~{}\exp{\left(-\frac{(x-y)^{2}}{2\sigma_{b}^{2}}\right)}=K(y,x)\end{split}$
(21)

This function is used because it weights the MBDF functional frequencies by
the functional value itself resulting in the distance measurement (eq. (18))
being weighted by the difference in functional values. Another advantage of
this function is that it eliminates the frequency of null values (or ”ghost
atoms”) within the MBDF representation which might be present due to the
procedure of zero-padding1.
In our work we generate a separate density function $f(u)$ for each of the 5
MBDF functionals in eq. (11-15), and for each unique chemical element present
in the dataset. These are then concatenated to form the global representation
of the molecule. Alternatively, it could be done by using multivariate
Gaussian functions for the density estimation. Let $\mathbf{x}_{i}$ denote the
5-dimensional vector of MBDF functional values (eq. 11-15) for any atom $i$ in
molecule $A$. Then the multivariate density function $f_{A}(\mathbf{u})$ of
this molecule can be evaluated as :

$\displaystyle
f_{A}(\mathbf{u})=\frac{1}{M}~{}\sum_{i=1}^{M}~{}K(\mathbf{u},\mathbf{x}_{i})$
(22) $\displaystyle
f_{A}(\mathbf{u})=\frac{1}{M}\sum_{i=1}^{M}\frac{(\mathbf{x}_{i}^{T}\mathbf{x}_{i})^{1/4}}{\sqrt{2\pi}\sigma_{b}}\exp{\left(-\frac{(\mathbf{u}-\mathbf{x}_{i})^{T}(\mathbf{u}-\mathbf{x}_{i})}{2\sigma_{b}^{2}}\right)}$
(23)

The l2-distance between molecules $A$ and $B$ then takes the form:

$d(A,B)^{2}=\int
d\mathbf{u}~{}(f_{A}(\mathbf{u})-f_{B}(\mathbf{u}))^{T}(f_{A}(\mathbf{u})-f_{B}(\mathbf{u}))$
(24)

where the integral is over the normalization region. Since the former method
generates a more compact representation we have chosen to work with it. The
abbreviation DF (Density of functionals) will be used throughout for this
global representation.

###  II.3 Numerical analysis

The DF method allows generating a representation that does not scale with the
number of atoms in the system. However, in order to use it as a feature vector
the density functions have to be discretized. Through convergence testing we
have set the grid spacing to 0.2 throughout our work. However, we note that
this grid spacing could be changed for a different data-set in order to
achieve the desirable accuracy vs computational cost trade-off. Furthermore,
DF corresponds to a flattened feature vector which can be used with global
kernels (or other ML methods), and which exhibits superior performance when
compared to a straightforward concatenation of all MBDF rows (see Figure 2).
The flattened MBDF representation is generated by sorting the MBDF matrix of
each molecule by the row norm, and then flattening the matrix by concatenation
of the rows to vectorize it1.

Figure 2: QML learning curves for MBDF with local kernel, DF (global), and
flattened MBDF (global and sorted by row norm). Mean absolute error (MAE) of
predicting atomization energies of the QM9 48 data-set as a function of
training set size N.

Figure 3 shows molecular fingerprints generated using the 1 and 5 functional
DF representations for three diverse and relevant organic molecules (glucose,
uric acid, and testosterone) on the same grid. For each molecule, a distinct
fingerprint is obtained, with peak-positions depending on the local chemical
environment of each atom. Consequently, peaks of atoms with chemically similar
environments are located closer to each other. Peak heights encode both number
and type (because of the density estimate being weighted) of chemical
environments [See Eq. 20]. Figure 3 demonstrates that for molecules with
increasing size, corresponding DF based fingerprints will grow in magnitude,
not in size. In the SI, we also show how DF fingerprints distinguish
conformational isomers, as exemplified for the chair and boat conformations of
cyclohexane.

Figure 3: Two-body (left) and three-body (right) DF versions (Eqs. 12, 15)
representations for Glucose, Uric Acid, and Testosterone.

##  III Methods and Data

###  III.1 Kernel ridge regression

The ML method that we focus on, and use throughout this work, is the
supervised learning method called Kernel ridge regression18; 39 (KRR). This
method has been covered extensively earlier1; 35; 34; 3; 2 so we skip the
details here.

The kernel functions55; 56 we use in our work along with global
representations are the Gaussian kernel,

$\displaystyle
k(\mathbf{x}_{I},\mathbf{x}_{J})=\exp{\left(-\frac{||\mathbf{x}_{I}-\mathbf{x}_{J}||^{2}_{2}}{2\sigma^{2}}\right)}$
(25)

and Laplacian kernel,

$\displaystyle
k(\mathbf{x}_{I},\mathbf{x}_{J})=\exp{\left(-\frac{||\mathbf{x}_{I}-\mathbf{x}_{J}||_{1}}{\sigma}\right)}$
(26)

where $\mathbf{x}_{I}$ denotes the representation vector of molecule $I$.
The kernel function used for the local representations FCHL19, SOAP and MBDF
is a summation of atomic kernels:

$\displaystyle k(\mathbf{M}_{I},\mathbf{M}_{J})=\sum_{a\epsilon
i}\sum_{b\epsilon j}k^{l}(\mathbf{x}_{Ia},\mathbf{x}_{Jb})$ (27)

with the local Gaussian kernel:

$\displaystyle
k^{l}(\mathbf{x}_{Ia},\mathbf{x}_{Jb})=\mathbf{\delta}_{Z_{a},Z_{b}}\exp{\left(-\frac{||\mathbf{x}_{Ia}-\mathbf{x}_{Jb}||^{2}_{2}}{2\sigma^{2}}\right)}~{}$
(28)

the local Laplacian kernel:

$\displaystyle
k^{l}(\mathbf{x}_{Ia},\mathbf{x}_{Jb})=\mathbf{\delta}_{Z_{a},Z_{b}}\exp{\left(-\frac{||\mathbf{x}_{Ia}-\mathbf{x}_{Jb}||_{1}}{\sigma}\right)}~{}$
(29)

or the local exponential kernel with the Euclidean norm:

$\displaystyle
k^{l}(\mathbf{x}_{Ia},\mathbf{x}_{Jb})=\mathbf{\delta}_{Z_{a},Z_{b}}\exp{\left(-\frac{||\mathbf{x}_{Ia}-\mathbf{x}_{Jb}||_{2}}{\sigma}\right)}~{}$
(30)

where $\mathbf{M}_{I}$ denotes the representation matrix of molecule $I$,
$\mathbf{x}_{Ia}$ denotes the representation vector of atom $a$ within
molecule $I$ and $\mathbf{\delta}_{Z_{a},Z_{b}}$ denotes a Kronecker Delta
over the nuclear charges $Z_{a},Z_{b}$ which restricts the similarity
measurement between atoms of the same chemical element34. Other kernel
functions will also be tested in the future57; 58; 40.
Throughout this study, we evaluate performance of ML methods through learning
curves for the task of predicting physical properties of molecular systems.
Learning curves quantify the model prediction error $\varepsilon$ (often
measured as mean absolute error (MAE)) against the number of training samples
$N$ and are key to understand the efficiency of ML models. It is generally
known18; 59; 60 that they are linear on a log-log scale,

$\displaystyle\log{\left(\varepsilon\right)}\approx I+S\log{(N)}~{}$ (31)

where $I$ is the initial error and $S$ is the slope indicating model
improvement given more training data. We also note that according to the
central limit theorem the distribution of the errors $\epsilon$ approaches the
normal distribution with standard deviation $\frac{\sigma}{\sqrt{N}}$, and
mean 0 as $N\xrightarrow{}\infty$. Hence eq. (31) becomes :

$\displaystyle\log{\left(\varepsilon\right)}\approx\log(\sigma)-\frac{1}{2}\log{(N)}$
(32)

###  III.2 Hyperparameter Optimization

The current form of the representations has been optimized for Kernel based
learning models. It depends on the weighting functions used and a number of
hyperparameters which include variances ($\sigma_{r},\sigma_{\theta}$) of the
Gaussian PDFs, weighting function hyperparameters mentioned in eq. (12)-(16),
and bandwidth $\sigma_{b}$ for the DF representation. The hyperparameter
optimization was done on a random subset of two thousand molecules from the
QM7 dataset49, and then kept fixed for all other data-sets. We note that
further improvements might be possible if they had been optimized
simultaneously on all data-sets. The weighting functions $g_{n}(\tau)$ in eq.
(8) were chosen by straightforward screening of the functions mentioned
earlier. The optimization minimized the atomization energy prediction errors
on the QM7 subset using Gaussian Process (GP) based Bayesian Optimization
(BOpt)39. Starting with a Gaussian prior, the method fits a posterior
distribution over the objective function using successive function
evaluations. The posterior distribution is then used to construct an efficient
acquisition function which can be optimized using, for instance, a quasi-
Newton method to determine the next query point. Table 1 shows the optimized
hyperparameter values used throughout this work.

Parameter | Value
---|---
$\sigma_{r}$ | 1
$\sigma_{\theta}$ | 2
$\eta$ | 10.8
$\alpha$ | 1.5
$a_{0},a_{1},a_{2},a_{3}$ | 3, 100, -200, -164
$\sigma_{b}$ | 0.07
$r_{c}$ (Å) | 6
Table 1: MBDF and DF hyperparameters after optimization on atomization
energies of QM7 49 subset.

We used the scikit-optimize61 implementation of GP based BOpt using the
default Matérn kernel with unit variance and the limited memory BFGS
optimizer62 for the acquisition function. In order to enable comparison to
other representations (such as FCHL or SOAP) that rely on distance cutoffs, we
have chosen to set inter-atomic distance cutoff $r_{c}$ for MBDF to 6 Å
throughout this work. We note that larger cutoffs for MBDF would not change
its size.

All MBDF functionals were evaluated using the trapezoidal numerical
integration method. The grid spacing for discretizing DF densities has been
set to 0.2 throughout our work as noted earlier. The bandwidth $\sigma_{b}$ =
0.07 was found to work well on the QM7 subset however it is recommended to be
screened once (in the range $\sim$ [0.01,1]) along with the grid-spacing when
using with other datasets. The Numpy63 and Numba64 libraries are used in the
representation generation code.

###  III.3 Data and Code

The QM9 dataset48 consists of $\sim$134k small organic molecules with up to 9
heavy atoms (C, N, O, F). The calculations were performed at the
B3LYP/6-31G(2df,p)65; 66; 67 level of theory.
QMugs is a dataset containing $\sim$665k biologically and pharmacologically
relevant drug-like molecules. It consists of larger molecules than QM9 with up
to 100 heavy atoms (C, N, O, F, P, S, Cl, Br, or I) per molecule. The training
and predictions were performed on the DFT ($\omega$B97X-D/def2-SVP)68; 69
values reported in the dataset. The QMugs subsets we used for Figure 5 were
drawn at random and consist of 20k molecules. Throughout, we used zero-padding
for all representations studied in order to accommodate training and test
molecules smaller than the maximum present in the data.

In order to keep the FCHL19 and SOAP kernel evaluations computationally
tractable, we have (a) restricted ourselves to a maximum 100 atoms per QMugs-
molecule, and (b) reduced the default hyperparameters of the FCHL19 and SOAP
representations to nRs2 = 12, nRs3 = 10, $r_{cut}$ = 6$\mathrm{\AA}$ and
$n_{max}$ = $l_{max}$ = 3, $\sigma$ = 0.1, $r_{cut}$ = 6$\mathrm{\AA}$,
respectively. For consistency, we used the same parameters for all other
results reported in this article. These two versions of the representations
with the reduced basis sets are denoted as FCHL19* and SOAP* in all reported
figures. Note that the latter choice of hyperparameters negligibly
deteriorates the predictive accuracy for QM9 (as assessed below and when
comparing to the published prediction errors on QM9 for FCHL19 and SOAP2013).
For FCHL19 and SOAP based prediction errors reported here within for QMugs, it
could still be that the accuracy could improve further if these parameters
were optimized.

Figure 4 also includes results for QML models based on the 2- (k = 2) and
3-body (k = 3) Many-Body Tensor Representations (MBTR)22 and a variant70 of
the Atom Centered Symmetry Functions20 (ACSF) as implemented in the QMLcode
library71. The MBTR representations were generated with the same
hyperparameters as those used in Ref.33 for the 10k training point on the QM9
dataset.

Throughout our work, the FCHL19, SLATM and ACSF representations were generated
using the QMLcode library71, SOAP was generated using the Dscribe library72
with the default gaussian type radial basis functions, MBTR was generated
using the qmmlpack library73 and SchNet74, PaiNN75 were built and trained
using the SchNetPack76 libray. Allegro77 was built and trained using the
nequip78 code79 based on the E(3)-NN framework80; 81.
For FCHL19, SOAP and ACSF we employed the local Gaussian (eq. (28)), for
SLATM, MBTR and DF we use the global Gaussian (eq. (25)) and for CM, BOB we
use the global Laplacian (eq. (26)) kernels respectively. For MBDF we used the
local Gaussian kernel on the MD17 dataset and the local exponential (eq. (30))
kernel for all other models. These choices were made based on the best
performances for each representation. All kernel evaluations were performed
using the QMLcode library.

Figure 4: MBDF/DF performance and comparison on atomization energies from QM9
data set (drawn at random from $\sim$134k organic molecules) 48. a) (Left)
Mean absolute error (MAE) of prediction as a function of training set size for
representations CM 1, BOB 38, MBTR 22, SLATM 37, PI 82, ACSF 20; 70, FCHL19
34, SOAP 36. Numbers in legend denote representation size (feature vector
dimensions), G and L denote Global and Local kernels, respectively. b) (Right)
Timing for training and testing as a function of training set size required
for making chemically accurate (MAE = 1 kcal/mol) predictions on 100k
molecules. Blue, red, and green points indicate local kernels, global kernels,
and neural network, respectively. Dashed gray line corresponds to the optimal
Pareto front. For SchNet 74, PaiNN75, Allegro77 an Nvidia RTX 3070 GPU has
been used. All other timings were evaluated on a compute node equipped with a
24-core AMD EPYC 7402P @ 3.3 GHz CPU and 512 GB RAM. Timings for FCHL18 35,
BOB, CM are estimated using smaller kernels (not taking into account kernel
inversion). Asterisk denotes representations with reduced hyperparameters used
in this work. $N$ values for ACSF, MBTR, BOB, CM estimated via extrapolation.
Numbers in brackets denote year of publication.  Figure 5: MBDF/DF performance
and comparison to CM 1, BOB 38, SLATM 37, FCHL19 34, and SOAP 36
representation based atomization energy estimates using the QMugs ($\sim$667k
organic molecules with up to one hundred heavy atoms) data set 50. Training
and testing data drawn at random. Prediction mean absolute errors (MAE) on
holdout test set of 4k molecules shown as a function of training set size.
Numbers in legend denote representation size (feature vector dimensions), G
and L denote Global and Local kernels respectively. Shaded region denotes the
standard deviation across 4 different learning curves (except for FCHL19 and
SOAP for which only one learning curve was tractable). Asterisk denotes
representations with reduced hyperparameters used in this work.

##  IV Results

###  IV.1 Atomization Energies

#### QM9

Figure 4 a) shows learning curves for QM9 and the size (dimension of feature
vector) of the representation arrays in the legend. For the task of predicting
atomization energies, local representations have previously been shown to be
more efficient34; 35. Results for the local representations FCHL19 and SOAP
are closely reproduced, and reach chemical accuracy after training on less
than 10k molecules 34. Among the global representations, SLATM has previously
also been shown34; 37 to perform quite well reaching chemical accuracy after
training on $\sim$9k molecules although it shows a smaller slope at larger
training sizes. This is closely followed by MBDF which reaches chemical
accuracy after training on $\sim$10k molecules (less than 10% of the dataset).
The global DF representation also performs decently well reaching chemical
accuracy at $\sim$32k training data. The local ACSF representation shows a
larger offset but a better slope and it reaches chemical accuracy at $\sim$50k
training set size. We note here that for consistency with the other
representations used in this work, we did not optimize the hyperparameters of
the MBTR representation on every training point, but rather kept them fixed
throughout. Only the KRR hyperparameters were optimized at each training set
size as with all other representations used here. The 3-body MBTR reaches
chemical accuracy at $\sim$60k training set size while the 2-body MBTR
performs better than the other 2-body representations, BOB and CM. We have
also included the recently introduced constant size Persistence Images82(PI)
representation for comparison.

Note that MBDF has the smallest size, requiring only 5 numbers per atom (145
dimensions for QM9 molecules). By contrast, other local representations such
as FCHL19, SOAP require $\sim$400 numbers per atom, while ACSF uses a 150
dimensional feature vector per atom. Encouragingly, and despite its compact
size, MBDF outperforms most of the other larger representations with the
exception of SOAP and FCHL. We note here that while the size of the global DF
representation is larger than MBDF, it utilises a global kernel implying
training and prediction cost invariance with respect to system size.

This compactness of the representation translates into faster ML model
evaluation timings. This is shown in Figure 4 b) which plots the trade-off
between training and prediction timings vs. training data needs for reaching
mean absolute prediction errors of atomization energies of 1 kcal/mol. We note
that there are only two representations located on the Pareto-front, FCHL18 35
and MBDF [this work]. We also point out that currently the best performing
model on the QM9 dataset is the recently proposed Wigner Kernels method83
which is not included in this study.
As noted earlier, local kernels based on representations such as FCHL18/19,
SOAP, or ACSF exhibit very good training data efficiency, but this comes at
the expense of a larger computational overhead. The exception is the local
MBDF based kernel which achieves the fastest training timing of $\sim$0.07
compute node minutes (14k training molecules) due to its compact size.
Predictions on 100k QM9 molecules using the local MBDF kernel are made in
$\sim$1.46 compute node minutes which translates to an unprecedented
chemically accurate navigation rate of $\sim$1140 molecules/second. SLATM
being clost to the Pareto front, and the DF representation, both represent
fast global kernel based KRR models. While requiring more training data than
SLATM in order to reach chemical accuracy, DF has the advantage that it is
largely invariant with respect to system size (see below). For the sake of
completeness, Fig. 4 b) also includes results for the deep learning based
models SchNet74, PaiNN75 and the current state-of-the-art, Allegro77 which
were all trained on a GPU. The reported timing for SchNet refers to 3000
epochs of training on 20k molecules and predictions of 100k molecules taking
about 3 h: 27 min and 7 sec respectively. For PaiNN the reported timing
corresponds to 1000 epochs of training on 7k molecules which took 0 h: 56 mins
and a prediction time of 9 sec. Allegro reached chemical accuracy after
training on 4500 molecules for 728 epochs, with early stoppping, which took 17
h: 8 min while the prediction on 100k molecules took $\sim$7 mins (using a
maximum possible batch size of 5 on the used GPU).
MBDF achieves the fastest training time out of all models and the fastest
prediction rate among the kernel based models. Numerical results for porting
KRR results based on FCHL19 to GPUs 84 would suggest that it seems likely for
the prediction rate to increase significantly once MBDF is reimplemented in
CUDA.

#### QMugs

We have tested the generalizability of our method to larger systems and more
diverse chemistries using the QMugs 50 data set. Figure 5 shows the
atomization energy learning curves. Due to the large variety in the dataset,
the predictive error is larger for all representations than their QM9
counterparts even when predicting on a much smaller test set. MBDF reaches
$~{}\sim$2 kcal/mol prediction error after training on 16k molecules. This is
better than the QML based neural network predictions published in Ref. 85, and
similar to the $\Delta$-QML numbers they also report. In terms of speed,
generating the local MBDF kernel for training and testing on 20k molecules on
this dataset takes $\sim$1.8 compute node mins (see below) which corresponds
to a navigation rate of $\sim$185 molecules/second. By comparison, this is
substantially faster than the GPU based prediction rates of approximately 50
and 5 molecules per second for the direct and $\Delta$-learning (using
GFN2-xTB 86) based ML models, respectively using the convolutional neural
network reported in Ref. 85. Only SLATM and FCHL19 exhibit lower off-set than
MBDF, while the performance for SOAP and DF is similar, albeit slightly worse
than MBDF. As mentioned before, however, in order to make FCHL19 and SOAP
tractable, we have dramatically reduced the hyper parameters. In particular,
we believe that the learning efficiency of SOAP for QMugs is being reduced due
to the use of small basis sets ($n_{max}$ = $l_{max}$ = 3). Note that no
representation reaches chemical accuracy within 16k training molecules,
indicating that QMugs possesses substantially more chemical diversity than
QM9.

In terms of representation sizes, MBDF again remains the smallest
representation since it still requires only 5 dimensions per atom regardless
of the chemical composition. However, being a local representation, on average
its size increases $\sim$ 3.4 times compared to QM9. FCHL19 and SOAP, on the
other hand, now require more than 1000 dimensions to represent each atom for
this larger dataset. CM, BOB, FCHL19, and SOAP show larger than 10 fold
increase in the representation size compared to QM9, followed by SLATM which
shows an increase of $\sim$ 6.6 times. This results in a considerable increase
in the train/test time (vide infra) which precludes the straightforward
application of these representations to the entire QMugs dataset. The DF
representation changes the least in size since it does not formally scale with
number of atoms but only with number of different chemical elements.
Consequently, its size doubles compared to QM9 since a separate density
function is generated for each unique chemical element in the dataset using
eq. (19) as mentioned earlier.

#### QM7b and MD17

Figure 2 in the SI shows learning curves for the QM7b87; 88 atomization
energies and the size (dimension of feature vector) of the representation
arrays in the legend. Similar trends in performance and representation size
noted so far are observed on this dataset as well.
Figure 3 in the SI shows learning curves of energies for a few molecules from
the revised89; 90 MD1791; 92 molecular dynamics dataset. Although the
comparative performance trend on this dataset is similar to the others, we
note that the chosen functionals for MBDF and the current representation
hyperparameters are not optimized for potential energy surface learning.
Furthermore, the implementation of gradients for MBDF to enable force based
learning should significantly enhance the performance for such tasks as was
shown89 for the FCHL19 representation. However, the relatively larger
difference in performance between MBDF and FCHL19 on this dataset might
indicate that compact representations are sufficient for regression through
equilibrium structures across CCS while more verbose representations are
required for potential energy surface learning.

###  IV.2 Timings

Figure 6 shows scaling plots of kernel evaluation timings across both QM9 and
QMugs datasets for various representations and as a function of training set
size. As one would expect, the off-set increases systematically with
increasing representation and training molecule size. More specifically, for
the larger molecules of QMugs, the FCHL19 and SOAP kernel evaluations become
rapidly intractable very quickly with the 16k kernel (Fig. 5) already taking
an entire compute node day. Encouragingly and in stark contrast, DF, being a
size invariant representation, shows hardly any change in computational
overhead when moving from the small QM9 molecules to the much larger QMugs
molecules.

For context, the time required in the kernel inversion step for each training
kernel is shown as well. The bottleneck crossover from kernel generation
($\mathcal{O}(n^{2})$) to the inversion step ($\mathcal{O}(n^{3})$) occurs
rather late. When using Cholesky decomposition it occurs for MBDF at $N\sim$
64k training molecules. For less compact representations (for SLATM and
larger) the same cross-over occurs at training set sizes that exceed $\sim$1
M. As demonstrated above, chemical accuracy is already achieved at
substantially smaller training set sizes. Consequently and contrary to popular
belief, for any of the more modern and accurate representations, kernel
inversion does not constitute a bottleneck.

Figure 6: Compute node time required for kernel inversion and evaluation as a
function training set sizes $N$ drawn from QM9 (squares) and QMugs (diamonds)
datasets. QML results shown for SOAP 36, FCHL19 34, SLATM 37, and MBDF and DF.
Dotted lines indicate extrapolation (using quadratic (kernel evaluation) and
cubic (kernel inversion) polynomial fits). G and L denote Global and Local
kernels, respectively. Compute node: 24-core AMD EPYC 7402P @ 3.3 GHz CPU with
512 GB RAM.

Table 2 reports the kernel evaluation and representation generation timings
for both the QM9 dataset (130k molecules) and the QMugs subset (20k molecules)
used in our work. It can be seen that MBDF reduces the local kernel evaluation
timings from days to a few minutes for both small and large molecules. For the
representation generation step we note that our code is currently written in
Python and uses the Numba64 library and could be further optimized with a low-
level implementation. However, the current timings as well do not affect the
overall QML model cost much given the kernel evaluation bottleneck.

| QM9 130k Molecules | QMugs 20k Molecules
---|---|---
Representation |  $t_{\mathrm{rep}}$ [min] |  $t_{\mathrm{kernel}}$ [min] | Dimension |  $t_{\mathrm{rep}}$ [min] |  $t_{\mathrm{kernel}}$ [min] | Dimension
CMa (G) | $0.186$ | $2.862$ | 435 | $0.012$ | $1.146$ | 5050
BoBa (G) | $0.216$ | $7.296$ | 1128 | $1.362$ | $3.396$ | 18528
SLATMa (G) | $18.60$ | $86.32$ | 11960 | $15.76$ | $14.58$ | 79045
FCHL19a (L) | $0.846$ | $1071$ | 10440 | $1.764$ | $1566$ | 122000
SOAPb (L) | $0.216$ | $1873$ | 13920 | $0.246$ | $2925$ | 186000
MBDF (L) | $1.626$ | $11.81$ | 145 | $4.182$ | $1.848$ | 500
DF (G) | $2.262$ | $12.16$ | 2500 | $2.442$ | $0.996$ | 5000
Table 2: Compute node times for generating representations
($t_{\mathrm{rep}}$) and kernel matrices ($t_{\mathrm{kernel}}$) for 130k
molecules from the QM9 dataset and 20k molecules from the QMugs dataset.
Global and Local kernels are again denoted by (G) and (L) respectively.
Representations with superscript $a$ were generated with QMLcode71 library and
$b$ with the Dscribe72 library. Compute node: 24-core AMD EPYC 7402P @ 3.3 GHz
CPU and 512 GB RAM.

###  IV.3 Performance for molecular quantum properties

Figure 7: Learning curves for various representations in QML models of highest
occupied molecular orbital (HOMO), lowest unoccupied molecular orbital (LUMO)
eigenvalues, HOMO-LUMO gap ($\Delta\epsilon$), internal energy at 0 K
($U_{0}$), dipole moment norm ($|\boldsymbol{\mu}|$), static isotropic
polarizability ($\alpha$), zero point vibrational energy (ZPVE), and heat
capacity ($C_{v}$) from the QM9 dataset. G and L denote Global and Local
kernels, respectively. Asterisk denotes representations with reduced
hyperparameters used in this work.

We assessed the generalization capacities of MBDF/DF on physical properties
other than atomization energies. Figure 7 shows the learning curves for the
task of predicting 8 important molecular quantum properties from the QM9
dataset. These include the highest occupied molecular orbital (HOMO), lowest
unoccupied molecular orbital (LUMO) eigenvalues and the HOMO-LUMO gap,
internal energy at 0 K ($U_{0}$), dipole moment norm ($|\boldsymbol{\mu}|$),
static isotropic polarizability ($\alpha$), zero point vibrational energy
(ZPVE) and heat capacity ($C_{v}$). Due to substantial computational costs,
the KRR hyper-parameters were not optimized at each training size for the
FCHL19 and SOAP representations, and we picked the same parameters as those
used for learning atomization energies in figure 4. We reproduce earlier
trends among intensive and extensive properties when using local/global
kernels 93; 3. MBDF and DF match this trend: They perform better on extensive
and on intensive properties, respectively.

Again, note that the performance on these properties of MBDF/DF could be
further improved by including different functionals suited to the property, or
by augmenting them with response terms as was done for the FCHL19
representation93. It would also be interesting to see how the learning
capacities across these different physical properties is affected by the
inclusion of higher order functional terms.

##  V Conclusion

We have introduced ultra-compact atomic local many body distribution
functional (MBDF) and global density of functionals (DF) representations for
use within kernel ridge regression based Quantum Machine Learning (QML) models
for rapid sampling of chemical compound space. MBDF/DF can accurately describe
any atom/molecule using a minimal number of discrete elements and thereby
reduce QML training and prediction times by multiple orders of magnitude for
small and large molecules. MBDF and DF correspond to functionals of analytical
weighted many body distributions in interatomic distances and angles, as well
as their derivatives. Chemical identity is encoded as a prefactor to the
atomic functionals. DF is a weighted density estimations of MBDF, i.e. a
global molecular fingerprint that is invariant with respect to number of atoms
(not number of chemical elements).

We have demonstrated predictive power competitive with the state-of-the-art
for a variety of quantum physical properties, as encoded in the QM9 dataset
48. On the QM9 dataset MBDF reaches a MAE of atomization energies of only 0.69
kcal/mol after training on 32k molecules while using only 5 dimensions to
describe an atom. Regarding different molecular properties, it is beneficial
to use DF representation along with a global kernel for intensive properties,
and MBDF along with a local kernel for the extensive properties. MBDF and DF
generalize also well to other compound spaces, as evinced for the chemically
more diverse QMugs dataset 50 consisting of substantially larger molecules:
After training on 16k molecules, MBDF reaches a peak MAE for atomization
energies of 1.97 kcal/mol while still using only 5 dimensions per atom.
Corresponding training and prediction times for both data-sets is $\sim$2
compute node time minutes for MBDF based models. Furthermore, our results
indicate that using MBDF/DF brings the train/test timings of kernel based ML
models to their ”lower-bound” as imposed by the kernel inversion bottleneck
for both small and large molecules.

We have analyzed the comparative performance for the sampling cost vs.
training set size needs to reach chemical accuracy for predicting atomization
energies in the QM9 data-set. MBDF has extended the corresponding optimal
Pareto front towards minimal time requirements.

While the numerical evidence collected is indicative of a surprising
effectiveness of MBDF, it is clear that the truncations used in this work may
lead to lack of uniqueness. However, neither within QM9 nor within the QMugs
subset we sample, have we encountered a case where using this small set of
functionals maps two different molecular structures to the same feature
vector. Furthermore, the likelihood of uniqueness can easily be increased by
inclusion of higher order derivatives and many-body terms. More rigorous
numerical analysis would be desirable to quantify this trend.

Thanks to its numerical efficiency, we believe that this approach holds great
promise for further accelerating the virtual discovery of novel molecules and
materials. Furthermore, this framework provides a possible solution to the
general problem of unfavourable scaling due to i) inclusion of higher order
many body and long range physics and ii) applying these ML models to larger
molecules with greater chemical diversity. Future work will deal with
extension of the representations as described to deal with various types of
conformers, and an assessment of its sensitivity to changes in molecular
structures.

##  VI Acknowledgments

D.K. is thankful for discussions with B. Huang, S. Krug, D. Lemm, M. Sahre,
and J. Weinreich. O.A.v.L. has received funding from the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme (grant agreement No. 772834). O.A.v.L. has received support as the
Ed Clark Chair of Advanced Materials and as a Canada CIFAR AI Chair.

##  VII Supplementary material

See supplementary material for QM7b87, MD1791; 92; 90; 89 learning curves and
a figure showing the transition of cyclohexane from chair to boat conformation
alongside the response by its DF fingerprint. It also contains heat maps
showing ideal KRR hyperparameters to be used with the MBDF representation and
kernel PCAs of some molecules from the QM7b dataset.

##  VIII Code and Data Availability

Python script for generating the MBDF and DF representations is available at
https://github.com/dkhan42/MBDF. Train/test split of the QM7 dataset used for
optimizing representation hyperparameters, other data and code for reproducing
the results in this study and some tools for performing KRR are openly
available at https://github.com/dkhan42/QML.

## References

  * Rupp  _et al._ (2012a) M. Rupp, A. Tkatchenko, K.-R. Müller, and O. A. von Lilienfeld, Phys. Rev. Lett. 108, 058301 (2012a).
  * von Lilienfeld (2018) O. A. von Lilienfeld, Angewandte Chemie International Edition 57, 4164 (2018), http://dx.doi.org/10.1002/anie.201709686.
  * Faber  _et al._ (2017) F. A. Faber, L. Hutchison, B. Huang, J. Gilmer, S. S. Schoenholz, G. E. Dahl, O. Vinyals, S. Kearnes, P. F. Riley, and O. A. von Lilienfeld, Journal of Chemical Theory and Computation 13, 5255 (2017), pMID: 28926232, https://doi.org/10.1021/acs.jctc.7b00577 .
  * Ho and Rabitz (1996) T. Ho and H. Rabitz, J. Chem. Phys. 104, 2584 (1996).
  * Manzhos and Carrington, Jr. (2006) S. Manzhos and T. Carrington, Jr., J. Chem. Phys. 125, 084109 (2006).
  * Behler and Parrinello (2007) J. Behler and M. Parrinello, Phys. Rev. Lett. 98, 146401 (2007).
  * Bartók  _et al._ (2010) A. P. Bartók, M. C. Payne, R. Kondor, and G. Csányi, Phys. Rev. Lett. 104, 136403 (2010).
  * Handley and Popelier (2010) C. M. Handley and P. L. A. Popelier, The Journal of Physical Chemistry A 114, 3371 (2010), pMID: 20131763, https://doi.org/10.1021/jp9105585 .
  * Smith  _et al._ (2017) J. S. Smith, O. Isayev, and A. E. Roitberg, Chem. Sci. 8, 3192 (2017).
  * Snyder  _et al._ (2012) J. C. Snyder, M. Rupp, K. Hansen, K.-R. Müller, and K. Burke, Phys. Rev. Lett. 108, 253002 (2012).
  * Pozun  _et al._ (2012) Z. D. Pozun, K. Hansen, D. Sheppard, M. Rupp, K.-R. Müller, and G. Henkelman, J. Chem. Phys. 136, 174101 (2012).
  * Ceriotti  _et al._ (2021a) M. Ceriotti, C. Clementi, and O. Anatole von Lilienfeld, Chem. Rev. 121, 9719 (2021a).
  * Ceriotti  _et al._ (2021b) M. Ceriotti, C. Clementi, and O. Anatole von Lilienfeld, The Journal of Chemical Physics 154, 160401 (2021b), https://doi.org/10.1063/5.0051418 .
  * von Lilienfeld  _et al._ (2015) O. A. von Lilienfeld, R. Ramakrishnan, M. Rupp, and A. Knoll, International Journal of Quantum Chemistry 115, 1084 (2015), https://onlinelibrary.wiley.com/doi/pdf/10.1002/qua.24912 .
  * Ghiringhelli  _et al._ (2015) L. M. Ghiringhelli, J. Vybiral, S. V. Levchenko, C. Draxl, and M. Scheffler, Phys. Rev. Lett. 114, 105503 (2015).
  * Musil  _et al._ (2021) F. Musil, A. Grisafi, A. P. Bartók, C. Ortner, G. Csányi, and M. Ceriotti, Chemical Reviews 121, 9759 (2021), pMID: 34310133, https://doi.org/10.1021/acs.chemrev.1c00021 .
  * Ramakrishnan and von Lilienfeld (2015) R. Ramakrishnan and O. A. von Lilienfeld, CHIMIA 69, 182 (2015).
  * Vapnik (1998) V. N. Vapnik,  _Statistical Learning Theory_ (Wiley-Interscience, 1998).
  * Collins  _et al._ (2018) C. R. Collins, G. J. Gordon, O. A. von Lilienfeld, and D. J. Yaron, The Journal of Chemical Physics 148, 241718 (2018), https://doi.org/10.1063/1.5020441 .
  * Behler (2011) J. Behler, The Journal of Chemical Physics 134, 074106 (2011), https://doi.org/10.1063/1.3553717 .
  * Zhu  _et al._ (2016) L. Zhu, M. Amsler, T. Fuhrer, B. Schaefer, S. Faraji, S. Rostami, S. A. Ghasemi, A. Sadeghi, M. Grauzinyte, C. Wolverton, and S. Goedecker, The Journal of Chemical Physics 144, 034203 (2016), https://doi.org/10.1063/1.4940026 .
  * Huo and Rupp (2022) H. Huo and M. Rupp, Machine Learning: Science and Technology 3, 045017 (2022).
  * Gastegger  _et al._ (2018) M. Gastegger, L. Schwiedrzik, M. Bittermann, F. Berzsenyi, and P. Marquetand, The Journal of Chemical Physics 148, 241709 (2018), https://doi.org/10.1063/1.5019667 .
  * Drautz (2019) R. Drautz, Phys. Rev. B 99, 014104 (2019).
  * Braams and Bowman (2009) B. J. Braams and J. M. Bowman, International Reviews in Physical Chemistry 28, 577 (2009), https://doi.org/10.1080/01442350903234923 .
  * Hirn  _et al._ (2015) M. J. Hirn, N. Poilvert, and S. Mallat, CoRR abs/1502.02077 (2015), 1502.02077 .
  * Shapeev (2016) A. V. Shapeev, Multiscale Modeling & Simulation 14, 1153 (2016), https://doi.org/10.1137/15M1054183 .
  * Nigam  _et al._ (2020) J. Nigam, S. Pozdnyakov, and M. Ceriotti, The Journal of Chemical Physics 153, 121101 (2020), https://doi.org/10.1063/5.0021116 .
  * Zaverkin and Kästner (2020) V. Zaverkin and J. Kästner, Journal of Chemical Theory and Computation 16, 5410 (2020), pMID: 32672968, https://doi.org/10.1021/acs.jctc.0c00347 .
  * Jinnouchi  _et al._ (2020) R. Jinnouchi, F. Karsai, C. Verdi, R. Asahi, and G. Kresse, The Journal of Chemical Physics 152, 234102 (2020), https://doi.org/10.1063/5.0009491 .
  * Thompson  _et al._ (2015) A. Thompson, L. Swiler, C. Trott, S. Foiles, and G. Tucker, Journal of Computational Physics 285, 316 (2015).
  * Willatt  _et al._ (2018) M. J. Willatt, F. Musil, and M. Ceriotti, Phys. Chem. Chem. Phys. 20, 29661 (2018).
  * Langer  _et al._ (2022) M. F. Langer, A. Goeßmann, and M. Rupp, npj Computational Materials 8, 41 (2022).
  * Christensen  _et al._ (2020) A. S. Christensen, L. A. Bratholm, F. A. Faber, and O. Anatole von Lilienfeld, The Journal of Chemical Physics 152, 044107 (2020), https://doi.org/10.1063/1.5126701 .
  * Faber  _et al._ (2018) F. A. Faber, A. S. Christensen, B. Huang, and O. A. von Lilienfeld, The Journal of Chemical Physics 148, 241717 (2018), https://doi.org/10.1063/1.5020710 .
  * Bartók  _et al._ (2013) A. P. Bartók, R. Kondor, and G. Csányi, Phys. Rev. B 87, 184115 (2013).
  * Huang and von Lilienfeld (2017) B. Huang and O. A. von Lilienfeld, “Quantum machine learning using atom-in-molecule-based fragments selected on-the-fly,”  (2017).
  * Hansen  _et al._ (2015) K. Hansen, F. Biegler, R. Ramakrishnan, W. Pronobis, O. A. von Lilienfeld, K.-R. Müller, and A. Tkatchenko, The Journal of Physical Chemistry Letters 6, 2326 (2015), pMID: 26113956, https://doi.org/10.1021/acs.jpclett.5b00831 .
  * Rasmussen and Williams (2006) C. E. Rasmussen and C. K. I. Williams,  _Gaussian processes for machine learning._ , Adaptive computation and machine learning (MIT Press, 2006) pp. I–XVIII, 1–248.
  * De  _et al._ (2016) S. De, A. P. Bartók, G. Csányi, and M. Ceriotti, Phys. Chem. Chem. Phys. 18, 13754 (2016).
  * Bartók and Csányi (2015) A. P. Bartók and G. Csányi, International Journal of Quantum Chemistry 115, 1051 (2015).
  * Chmiela  _et al._ (2023) S. Chmiela, V. Vassilev-Galindo, O. T. Unke, A. Kabylda, H. E. Sauceda, A. Tkatchenko, and K.-R. Müller, Science Advances 9, eadf0873 (2023), https://www.science.org/doi/pdf/10.1126/sciadv.adf0873 .
  * Faber  _et al._ (2016) F. A. Faber, A. Lindmaa, O. A. von Lilienfeld, and R. Armiento, Phys. Rev. Lett. 117, 135502 (2016).
  * Morse (1929) P. M. Morse, Phys. Rev. 34, 57 (1929).
  * Jones (1924) J. E. Jones, Proceedings of the Royal Society of London. Series A (1924).
  * Axilrod and Teller (1943) B. M. Axilrod and E. Teller, The Journal of Chemical Physics 11, 299 (1943), https://doi.org/10.1063/1.1723844 .
  * Stillinger and Weber (1985) F. H. Stillinger and T. A. Weber, Phys. Rev. B 31, 5262 (1985).
  * Ramakrishnan  _et al._ (2014) R. Ramakrishnan, P. O. Dral, M. Rupp, and O. A. von Lilienfeld, Scientific Data 1 (2014).
  * Rupp  _et al._ (2012b) M. Rupp, A. Tkatchenko, K.-R. Müller, and O. A. von Lilienfeld, Phys. Rev. Lett. 108, 058301 (2012b).
  * Isert  _et al._ (2022) C. Isert, K. Atz, J. Jiménez-Luna, and G. Schneider, Scientific Data 9, 273 (2022).
  * Pozdnyakov  _et al._ (2020) S. N. Pozdnyakov, M. J. Willatt, A. P. Bartók, C. Ortner, G. Csányi, and M. Ceriotti, Phys. Rev. Lett. 125, 166001 (2020).
  * Huang and von Lilienfeld (2016) B. Huang and O. A. von Lilienfeld, The Journal of Chemical Physics 145, 161102 (2016), https://doi.org/10.1063/1.4964627 .
  * Parzen (1962) E. Parzen, The Annals of Mathematical Statistics 33, 1065 (1962).
  * Amari and Nagaoka (2000) S.-i. Amari and H. Nagaoka,  _Methods of information geometry_ , Vol. 191 (American Mathematical Soc., 2000).
  * Deisenroth  _et al._ (2020) M. P. Deisenroth, A. A. Faisal, and C. S. Ong,  _Mathematics for Machine Learning_ (Cambridge University Press, 2020).
  * Schütt  _et al._ (2020) K. T. Schütt, S. Chmiela, O. A. von Lilienfeld, A. Tkatchenko, K. Tsuda, and K.-R. Müller,  _Machine Learning Meets Quantum Physics_ (Springer, 2020).
  * Çaylak  _et al._ (2020) O. Çaylak, O. A. von Lilienfeld, and B. Baumeier, Machine Learning: Science and Technology 1, 03LT01 (2020).
  * Fabregat  _et al._ (2022) R. Fabregat, P. van Gerwen, M. Haeberle, F. Eisenbrand, and C. Corminboeuf, Machine Learning: Science and Technology 3, 035015 (2022).
  * Cortes  _et al._ (1994) C. Cortes, L. D. Jackel, S. A. Solla, V. Vapnik, and J. S. Denker, in  _Advances in Neural Information Processing Systems_ (1994) pp. 327–334.
  * Müller  _et al._ (1996) K. R. Müller, M. Finke, N. Murata, K. Schulten, and S. Amari, Neural Comp. 8, 1085 (1996).
  * Head  _et al._ (2021) T. Head, M. Kumar, H. Nahrstaedt, G. Louppe, and I. Shcherbatyi, “scikit-optimize/scikit-optimize,”  (2021).
  * Liu and Nocedal (1989) D. C. Liu and J. Nocedal, Mathematical Programming 45, 503 (1989).
  * Harris  _et al._ (2020) C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Nature 585, 357 (2020).
  * Lam  _et al._ (2015) S. K. Lam, A. Pitrou, and S. Seibert, in  _Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC_, LLVM ’15 (Association for Computing Machinery, New York, NY, USA, 2015).
  * Becke (1993) A. D. Becke, The Journal of Chemical Physics 98, 5648 (1993), https://doi.org/10.1063/1.464913 .
  * Lee  _et al._ (1988) C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B 37, 785 (1988).
  * Petersson  _et al._ (1988) G. A. Petersson, A. Bennett, T. G. Tensfeldt, M. A. Al‐Laham, W. A. Shirley, and J. Mantzaris, The Journal of Chemical Physics 89, 2193 (1988), https://doi.org/10.1063/1.455064 .
  * Chai and Head-Gordon (2008) J.-D. Chai and M. Head-Gordon, Phys. Chem. Chem. Phys. 10, 6615 (2008).
  * Weigend and Ahlrichs (2005) F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys. 7, 3297 (2005).
  * Yao  _et al._ (2018) K. Yao, J. E. Herr, D. W. Toth, R. Mckintyre, and J. Parkhill, Chemical Science 9, 2261 (2018).
  * Christensen  _et al._ (2017) A. S. Christensen, L. A. Bratholm, F. A. Faber, B. Huang, A. Tkatchenko, K. R. Müller, and O. A. von Lilienfeld, “Qml: A python toolkit for quantum machine learning,”  (2017).
  * 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, Computer Physics Communications 247, 106949 (2020).
  * Rupp (2015) M. Rupp, International Journal of Quantum Chemistry 115, 1058 (2015).
  * Schütt  _et al._ (2018) K. T. Schütt, P. Kessel, M. Gastegger, K. A. Nicoli, A. Tkatchenko, and K.-R. Müller, Journal of Chemical Theory and Computation 15, 448 (2018).
  * Schütt  _et al._ (2021) K. T. Schütt, O. T. Unke, and M. Gastegger, CoRR abs/2102.03150 (2021), 2102.03150 .
  * Schütt  _et al._ (2017) K. Schütt, P.-J. Kindermans, H. E. Sauceda Felix, S. Chmiela, A. Tkatchenko, and K.-R. Müller, in  _Advances in Neural Information Processing Systems_, Vol. 30, edited by I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (Curran Associates, Inc., 2017).
  * Musaelian  _et al._ (2023) A. Musaelian, S. Batzner, A. Johansson, L. Sun, C. J. Owen, M. Kornbluth, and B. Kozinsky, Nature Communications 14 (2023), 10.1038/s41467-023-36329-y.
  * Batzner  _et al._ (2022) S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth, N. Molinari, T. E. Smidt, and B. Kozinsky, Nature Communications 13 (2022), 10.1038/s41467-022-29939-5.
  * Geiger  _et al._ (2022a) M. Geiger, T. Smidt, M. Alby, B. K. Miller, W. Boomsma, B. Dice, K. Lapchevskyi, M. Weiler, Michał Tyszkiewicz, M. Uhrin, S. Batzner, D. Madisetti, J. Frellsen, N. Jung, S. Sanborn, Jkh, Mingjian Wen, J. Rackers, M. Rød, and M. Bailey, “e3nn/e3nn: 2022-12-12,”  (2022a).
  * Geiger  _et al._ (2022b) M. Geiger, T. Smidt, A. M., B. K. Miller, W. Boomsma, B. Dice, K. Lapchevskyi, M. Weiler, M. Tyszkiewicz, S. Batzner, D. Madisetti, M. Uhrin, J. Frellsen, N. Jung, S. Sanborn, M. Wen, J. Rackers, M. Rød, and M. Bailey, “Euclidean neural networks: e3nn,”  (2022b).
  * Geiger and Smidt (2022) M. Geiger and T. Smidt, “e3nn: Euclidean neural networks,”  (2022).
  * Townsend  _et al._ (2020) J. Townsend, C. P. Micucci, J. H. Hymel, V. Maroulas, and K. D. Vogiatzis, Nature Communications 11, 3230 (2020).
  * Bigi  _et al._ (2023) F. Bigi, S. N. Pozdnyakov, and M. Ceriotti, “Wigner kernels: body-ordered equivariant machine learning without a basis,” (2023), arXiv:2303.04124 [physics.chem-ph] .
  * Browning  _et al._ (2022) N. J. Browning, F. A. Faber, and O. Anatole von Lilienfeld, The Journal of Chemical Physics 157, 214801 (2022).
  * Atz  _et al._ (2022) K. Atz, C. Isert, M. N. Böcker, J. Jiménez-Luna, and G. Schneider, Physical Chemistry Chemical Physics 24, 10775 (2022).
  * Bannwarth  _et al._ (2019) C. Bannwarth, S. Ehlert, and S. Grimme, Journal of Chemical Theory and Computation 15, 1652 (2019).
  * Montavon  _et al._ (2013) G. Montavon, M. Rupp, V. Gobre, A. Vazquez-Mayagoitia, K. Hansen, A. Tkatchenko, K.-R. Müller, and O. A. von Lilienfeld, New Journal of Physics 15, 095003 (2013).
  * Blum and Reymond (2009) L. C. Blum and J.-L. Reymond, J. Am. Chem. Soc. 131, 8732 (2009).
  * Christensen and von Lilienfeld (2020) A. S. Christensen and O. A. von Lilienfeld, Machine Learning: Science and Technology 1, 045018 (2020).
  * Christensen and lilienfeld (2020) A. S. Christensen and A. V. lilienfeld,  (2020), 10.6084/m9.figshare.12672038.v3.
  * Chmiela  _et al._ (2017) S. Chmiela, A. Tkatchenko, H. E. Sauceda, I. Poltavsky, K. T. Schütt, and K.-R. Müller, Science Advances 3 (2017), 10.1126/sciadv.1603015.
  * Chmiela  _et al._ (2018) S. Chmiela, H. E. Sauceda, K.-R. Müller, and A. Tkatchenko, Nature Communications 9 (2018), 10.1038/s41467-018-06169-2.
  * Christensen  _et al._ (2019) A. S. Christensen, F. A. Faber, and O. A. von Lilienfeld, The Journal of Chemical Physics 150, 064105 (2019), https://doi.org/10.1063/1.5053562 .


Z^{({\rm adj})} \underset{w_{N-1} \rightarrow \infty}{\longrightarrow } \left( \sum_Y \prod_{(i,j)\in Y} \frac{\left(h_{(i,j)}-\frac{M}{\tilde{\epsilon}}\right)\left(h_{(i,j)}+\frac{M}{\tilde{\epsilon}})\right)}{h_{(i,j)}^2}\right)^N =  \left( q^{-\frac{1}{24}}\eta(q)\right)^{-N\left(1-\frac{M^2}{\tilde{\epsilon}^2} \right)}.
\end{eqnarray}

For $\epsilon_1 \neq -\epsilon_2$ the product in (<ref>) was not rigorously computed yet, but in [2], [6] a suggestion has been made for $SU( 2)$ and $SU(3)$ theories which appears to be correct. Embracing this conjecture we get

\begin{eqnarray}
Z^{({\rm adj})}\underset{w_{N-1} \rightarrow \infty}{\longrightarrow }\left( \sum_Y \prod_{(i,j)\in Y}\frac{f_{(i,j)}(M)}{f_{(i,j)}(0)} \right)^N= \left( q^{-\frac{1}{24}}\eta(q)\right)^{-N\left(1+\frac{M}{\epsilon_1} \right)\left(1+\frac{M}{\epsilon_2} \right)}.
\end{eqnarray}

* Fundamental and anti-fundamental matter. In the case of fundamental and anti-fundamental hypermultiplets we see from (<ref>) that the leading term of $Z_k^{({\rm fund})}$ is

\begin{equation}
Z^{({\rm fund})}_k \underset{w_{N-1} \rightarrow \infty}{\longrightarrow } \frac{a_N^{k(N_f+N_a)}}{a_N^{2k(N-1)}}c_k,
\end{equation}

where $c_k$ is some constant. Therefore if $N_f+N_a<2(N-1)$, then again the only non-vanishing contribution is $Z_0^{({\rm fund})}$ and

\begin{equation}
Z^{({\rm fund})} \underset{w_{N-1} \rightarrow \infty}{\longrightarrow } 1, \quad N_f+N_a<2(N-1).
\end{equation}

In the case of the critical number of matters $N_f+N_a=2(N-1)$ the limit of $Z^{({\rm fund})}$ is a constant not depending on ${\bf a}$. Let us find this constant.
Only the leading term of $Z_k$ matters in this case, so

\begin{eqnarray} \label{Zkdunflim}
Z_k^{({\rm fund})} \underset{w_{N-1} \rightarrow \infty}{\longrightarrow }&&\sum_{\underset{|\vec{Y}|=k}{\vec{Y}}}\frac{\varepsilon^{(2N-2)\sum_{u=1}^{N}u |Y_u|}}{\prod_{u=1}^N(\prod_{v \neq u}(\varepsilon^v-\varepsilon^u)^{|Y_u|})\cdot\prod_{v=1}^N(\prod_{u \neq v}(\varepsilon^v-\varepsilon^u)^{|Y_v|})} \frac{1}{\prod_{u=1}^N  \prod_{(i,j)\in Y_u} \epsilon_1\epsilon_2 f_{(i,j)}} \nonumber \\
&=&\sum_{\underset{|\vec{Y}|=k}{\vec{Y}}}\prod_{u=1}^N\frac{1}{\prod_{v \neq N}(-(1-\varepsilon^v))^{2|Y_u|}} \frac{1}{ \prod_{(i,j)\in Y_u} \epsilon_1\epsilon_2 f_{(i,j)}}.
\end{eqnarray}

Note that

\begin{equation}
\prod_{v \neq N}(1-\varepsilon^v)=\underset{x \rightarrow 1}{\rm lim}\frac{x^N-1}{x-1}=\frac{{\rm d} x^N}{{\rm d}x}\bigg|_{x=1}=N.
\end{equation}

\begin{equation} \label{Zkfactorised}
Z_k^{({\rm fund})}
\underset{w_{N-1} \rightarrow \infty}{\longrightarrow }
\sum_{\underset{|\vec{Y}|=k}{\vec{Y}}}\prod_{u=1}^N \left(\frac{(-1)^{N-1}}{N^2\epsilon_1 \epsilon_2}\right)^{|Y_u|}\prod_{(i,j)\in Y_u}\frac{1}{f_{(i,j)}}.
\end{equation}

We see again that there is no interaction between the Young diagrams in (<ref>), and hence we can write

\begin{equation}
Z^{({\rm fund})}
\underset{w_{N-1} \rightarrow \infty}{\longrightarrow } \sum_{k} q^k \sum_{\underset{|\vec{Y}|=k}{\vec{Y}}}\prod_{u=1}^N\left(\frac{(-1)^{N-1}}{N^2\epsilon_1 \epsilon_2}\right)^{|Y_u|}\prod_{(i,j)\in Y_u}\frac{1}{f_{(i,j)}} =\left( \sum_{Y} q^{|Y|}\left(\frac{(-1)^{N-1}}{N^2\epsilon_1 \epsilon_2}\right)^{|Y|}\prod_{(i,j)\in Y_u}\frac{1}{f_{(i,j)}}\right)^N.
\end{equation}

Using a relation provided[To see (<ref>) from (4.5) of [9] one has to replace $t_1 \rightarrow e^{\epsilon_1\delta}$, $t_2 \rightarrow e^{\epsilon_2\delta}$, $q \rightarrow \delta^2 q$ and send $\delta$ to zero.] in [9]

\begin{equation} \label{knownrelation}
\sum_{Y}x^{|Y|}\prod_{I \in Y}\frac{1}{f_{(i,j)}}=e^x
\end{equation}

we immediately get

\begin{equation} \label{Zfundasympt}
Z^{({\rm fund})}\underset{w_{N-1} \rightarrow \infty}{\longrightarrow }{\rm exp}\left( (-1)^{N-1}\frac{q}{N \epsilon_1 \epsilon_2}\right), \quad N_f+N_a=2(N-1).
\end{equation}

Finally if the number of fundamental and anti-fundamental hypermultiplets is above critical $N_f+N_a>2(N-1)$ the asymptotic behaviour of $Z^{({\rm fund})}({\bf a})$ can be a nontrivial function of ${\bf a}$ and finding it goes beyond this paper. We refer an interested reader to [2], [6], where this behaviour was studied with the AGT approach in $SU(2)$ and $SU(3)$ theories. Although our residue formula (<ref>) is valid also in this case, the recurrence relation which we will find below does not describe it.

The recurrence relations will be written for a partition function $\bar{Z}({\bf a})$ normalised to the constants discussed above, so the asymptotic behaviour of the normalised partition function is

\begin{equation}
\bar{Z}^{({\rm R})} \underset{w_{N-1} \rightarrow \infty}{\longrightarrow } 1.
\end{equation}

The recurrent relation via symmetric variables. Putting all together and taking into account the behaviour at infinity we get the recurrent Zamolodchikov-like relations in terms of the symmetric variables.

\begin{equation}
\bar{Z}^{\rm (R)}({\bm \omega},w_{N-1})=1+\sum_{k=1}^{N-1}\sum_{m,n=1}^{\infty} \frac{q^{mn} J^{(m,n)}({\bm \omega},\bar{w}_{N-1}^{(k|m,n)}) }{(w_{N-1}-\bar{w}^{(k|m,n)}_{N-1})}\frac{\mathcal{P}^{(m,n)}_{N, {\rm R}}({\bm \omega},\bar{w}_{N-1}^{(k|m,n)})}{\mathcal{P}^{(m,n)}_{N}({\bm \omega},\bar{w}_{N-1}^{(k|m,n)})}\bar{Z}^{\rm (R)}({\bf \hat{w}}^{(k|m,n)}).
\end{equation}

Explicit examples.
Let us now explicitly write the recurrence relations in several simplest cases.

* SU(2) case.
For completeness let us formulate the $SU(2)$ case in terms of symmetric variables, although the Weyl symmetry is evident even in terms of the variables ${\bf a}$ in this case.
The only symmetric variable in this case is

\begin{equation*}
w_1=a_1 a_2 = - \frac{1}{4} a_{12}^2.
\end{equation*}

The polynomial $Q^{({\bf w})}(x)$ is quadratic

\begin{equation*}
Q^{({\bf w})}(x)=x^2+w_1
\end{equation*}

and the resultant gives us the denominator

\begin{equation*}
\Delta^{(m,n)}({\bf w})=-4 w_1 - \epsilon_{m,n}^2,
\end{equation*}

so the pole is located at

\begin{equation*}
\bar{w}_1=-\frac{1}{4}\epsilon_{m,n}^2
\end{equation*}

as it should.
The dual point is

\begin{equation*}
\hat{w}_1^{(m,n)}=\bar{w}_1-\Delta w_1^{(m,n)} =-\frac{1}{4}\epsilon_{m,n}^2+m n \, \epsilon_1 \epsilon_2=-\frac{1}{4} \epsilon_{m,-n}^2,
\end{equation*}

which obviously corresponds to the partial Weyl permutation in $a_1$, $a_2$.
The polynomial $\Delta Q(x)$ in this case is just a constant  $\Delta Q(x)=1$, so the most part of   $\mathcal{P}^{(m,n)}$ and $J^{(m,n)}(\bar{w}_1)$ disappears, and the recurrent relation for the pure theory is just

\begin{equation}
\bar{Z}^{({\rm R})}(w_1)=1+\sum_{m,n=1}^{\infty}\frac{ q^{mn} 2 \epsilon_{m,n} }{  \left(-4w_1-\epsilon_{m,n}^2 \right)}\frac{\mathcal{P}^{(m,n)}_{2,{\rm R}}}{\mathcal{P}^{(m,n)}_{2}}\bar{Z}^{({\rm R})}\left(-\frac{1}{4}\epsilon_{m,-n}^2\right),
\end{equation}

\begin{equation}
\mathcal{P}^{(m,n)}_{2}=  \prod_{i=-m}^{m-1}\sideset{}{'} \prod_{j=-n}^{n-1} \epsilon_{i,j},
\end{equation}

\begin{equation}
\mathcal{P}^{(m,n)}_{2,{\rm adj}}= \prod_{i=-m}^{m-1} \prod_{j=-n}^{n-1} (\epsilon_{i,j}-M),
\end{equation}

\begin{equation}
\mathcal{P}^{(m,n)}_{2,{\rm fund}}=\prod_{i=1}^{m}\prod_{j=1}^{n} \Bigg[ \prod_{t=1}^{N_f}   \left(-\frac{1}{2}\epsilon_{m,n}+\epsilon_{i,j}-m_t  \right)\cdot  \prod_{t=1}^{N_a}  \left(-\frac{1}{2}\epsilon_{m,n}+\epsilon_{i,j}+m_t \right) \Bigg].
\end{equation}

* SU(3) case.
Let us now compare our result in the $SU(3)$ case with the one obtained in [6].
The variables $(u,v)$ used in  [6] differ from ours by a numerical factor

\begin{eqnarray}
w_1&=&\sum_{i<j} a_i a_j=-\frac{1}{3} (a_{12}^2+a_{12}a_{23}+a_{23}^2)=
-\frac{1}{3} u, \nonumber \\
w_2&=& a_1 a_2 a_3=-\frac{1}{27} (a_{12}-a_{23})(2a_{12}+a_{23})(a_{12}+2a_{23})=-\frac{1}{27} v. \nonumber
\end{eqnarray}

The polynomial $Q(x)$ in this case is

\begin{equation*}
Q^{({\bf w})}(x)=x^3+x w_1 - w_2,
\end{equation*}

which gives us a square equation for the poles

\begin{equation}
27 {w}_2^2+4 w_1^3+9w_1^2\epsilon_{m,n}^2+6w_1\epsilon_{m,n}^4+\epsilon_{m,n}^6=0
\end{equation}

and hence the positions of the poles are

\begin{equation} \label{wmn}
\bar{w}^{(k|m,n)}_2 = \pm \frac{1}{3}(w_1+\epsilon_{m,n}^2)\sqrt{-\frac{1}{3}(4 w_1+\epsilon_{m,n}^2)}\triangleq\pm w_{(m,n)}(w_1).
\end{equation}

System of recurrent equations (<ref>) boils down to only one equation

\begin{equation} \label{recn3}
w_1+\frac{3}{4}\frac{(\Delta w_2^{(k|m,n)} )^2}{(\Delta w_1^{(m,n)})^2}+\frac{1}{4}\epsilon_{m,n}^2=0.
\end{equation}

It has two roots and we have to pick one for each of $\bar{w}_2^{(k|m,n)}$ using the Euclidean algorithm to divide $\bar{Q}(x+\epsilon_{m,n})$ by $\bar{Q}(x)$. After two steps of the division we get a linear remainder

\begin{equation}
r_2(x)=x\left(\frac{2 w_1}{3}+\frac{2 \epsilon_{m,n}^2}{3} \right)-\bar{w}_2^{(k|m,n)}+\frac{w_1 \epsilon_{m,n}}{3}+\frac{\epsilon_{m,n}^3}{3}
\end{equation}

and demanding that $r_2(a_2(\Delta w_2^{(k|m,n)} ))=0$ we find the shift $\Delta w_2^{(k|m,n)}$

\begin{equation} \label{dw2}
\Delta w_2^{(k|m,n)} =-\frac{3 \bar{w}_2^{(k|m,n)} \Delta w_1^{(m,n)} }{w_1+\epsilon_{m,n}}.
\end{equation}

Substituting two roots $w_{2}^{(k|m,n)}$ given by (<ref>) we see that (<ref>) is indeed the two roots of (<ref>).
Therefore the dual points $(\hat{w}_1^{(m,n)},\hat{w}_2^{(k|m,n)})$ are

\begin{eqnarray}
\hat{w}_1^{(m,n)} &=&w_1+ m n \,\epsilon_1 \epsilon_2, \nonumber \\
\hat{w}_2^{(k|m,n)}&=& \pm w_{(m,n)}(w_1) \left(1+\frac{3\Delta w_1^{(m,n)} }{w_1+\epsilon_{m,n}} \right)=\pm  w_{(m,-n)}(\hat{w}_1^{(m,n)} )
\end{eqnarray}

in agreement with (<ref>).
The polynomial $\Delta Q(x)$ is linear

\begin{equation}
\Delta Q(x)=x-\frac{\Delta w_2^{(k|m,n)} }{\Delta w_1^{(m,n)} }.
\end{equation}

The recurrence relation in $SU(3)$ theory is the following

\begin{eqnarray} \label{recsymn3}
\bar{Z}^{({\rm R})}(w_1,w_2)=1+\sum_{k=1}^{2}\sum_{m,n=1}^{\infty}\bigg( &&\frac{q^{mn} (3w_1+\epsilon_{m,n}^2)(-w_1-\epsilon_{m,n}^2)\epsilon_{mn}}{9\bar{w}_2^{(k|m,n)} (w_2-\hat{w}_2^{(k|m,n)})} \nonumber \\ && \cdot\frac{\mathcal{P}^{(m,n)}_{3,{\rm R}}(w_1,\bar{w}_2^{(k|m,n)})}{\mathcal{P}^{(m,n)}_{3}(w_1,\bar{w}_2^{(k|m,n)})}\bar{Z}^{({\rm R})}(\hat{w}_1^{(m,n)} ,\hat{w}_2^{(k|m,n)}) \bigg),
\end{eqnarray}

\begin{equation}
\mathcal{P}^{(m,n)}_{3}(w_1,\bar{w}_2^{(k|m,n)})=\prod_{i=-m}^{m-1}\sideset{}{'} \prod_{j=-n}^{n-1} \epsilon_{i,j} \cdot \prod_{i=1}^{m} \prod_{j=1}^{n}(3 w_1+\epsilon_{m,n}^2-\epsilon_{i,j}\epsilon_{m-i,n-j}),
\end{equation}

\begin{equation}
\mathcal{P}^{(m,n)}_{3,{\rm adj}}(w_1,\bar{w}_2^{(k|m,n)})= \prod_{i=-m}^{m-1} \prod_{j=-n}^{n-1} (\epsilon_{i,j}-M) \cdot \prod_{i=1}^{m} \prod_{j=1}^{n}(3 w_1+\epsilon_{m,n}^2-(\epsilon_{i,j}+M)(\epsilon_{m-i,n-j}-M)),
\end{equation}

and $\mathcal{P}_{3, {\rm fund}}^{(m,n)}({\bm \omega}, \bar{w}_{N-1}^{(k|m,n)})$ is always the same and given by (<ref>).
One can recurrently see that for a pure theory and a theory with the adjoint matter the partition function actually depends only on $(\Delta w_2^{(k|m,n)})^2$, i.e. only on  $w_{(m,n)}^2$, so we can write (<ref>) as

\begin{equation}
\bar{Z}^{({\rm R})}(w_1,w_2^2)=1+\sum_{m,n=1}^{\infty}\frac{q^{mn}(3w_1+\epsilon_{m,n}^2)(-w_1-\epsilon_{m,n}^2)2 \epsilon_{m,n}}{9(w_2^2-w_{(m,n)}^2)}\frac{\mathcal{P}^{(m,n)}_{3,{\rm R}}(w_1)}{\mathcal{P}^{(m,n)}_{3}(w_1)}\bar{Z}^{({\rm R})}(\hat{w}_1^{(m,n)} ,w_{(m,-n)}^2(\hat{w}_1^{(m,n)} )).
\end{equation}

The results for pure theory and theory with adjoint hypermultiplet coincide with the ones found in [6] up to the sign of the mass of adjoint multiplet $M$. To compare also the case with fundamental and anti-fundamental matter we should consider a particular case of $N_f=N_a=N$ and redefine the masses as

\begin{eqnarray*}
m_t \rightarrow \epsilon-m_{t} \quad {\rm fundamental} \\
m_{t} \rightarrow -m_{t}  \quad {\rm anti-fundamental}\\
\end{eqnarray*}

In this case we do not know the behaviour of the partition function at infinity since $2N$ is above the critical number $2(N-1)$, but if we embrace the asymptotic behaviour provided in [6], we will recover exactly the recurrence relation found in there.

* SU(4) case.
Finally we are to write the recurrence relations for a non considered before case of $SU(4)$ theory.
The polynomial $Q^{({\bf w})}(x)$ in this case is

\begin{equation}
Q^{({\bf w})}(x)=x^4+w_1 x^2-w_2 x+w_3,
\end{equation}

and the equation for the poles

\begin{equation} \label{polesSU4}
{\rm \bf res}(Q^{({\bm \omega},{w}_{3})}(x), Q^{({\bm \omega},{w}_{3})}(x+\epsilon_{m,n}))=0
\end{equation}

is a cubic equation on ${w}_{3}$ of general form with three roots $\bar{w}^{(k|m,n)}_{3}({\bm \omega})$.
The system of recurrence equations (<ref>) consists of two equations

\begin{equation}
\begin{cases}
w_2+\frac{\Delta w_3^{(k|m,n)} \Delta w_2^{(k|m,n)} }{(\Delta w_1^{(m,n)} )^2}+\frac{1}{4}\frac{\Delta w_2^{(k|m,n)} }{\Delta w_1^{(m,n)} } \left( \epsilon_{m,n}^2-\frac{(\Delta w_2^{(k|m,n)} )^2}{(\Delta w_1^{(m,n)} )^2}\right)=0 \\
w_1-\frac{\Delta w_3^{(k|m,n)} }{\Delta w_1^{(m,n)} }+\frac{3}{4}\frac{(\Delta w_2^{(k|m,n)} )^2}{(\Delta w_1^{(m,n)} )^2}+\frac{1}{4}\epsilon_{m,n}^2=0
\end{cases}
\end{equation}

For the shift $\Delta w_3\textbf{}$ we get

\begin{equation}
\Delta w_{3}\textbf{}=\frac{4 \bar{w}_{3}^{(k|m,n)}}{\frac{(\Delta w_2^{(k|m,n)} )^2}{(\Delta w_1^{(m,n)} )^2}-\epsilon_{m,n}^2}.
\end{equation}

After three steps of the Euclidean algorithm of the division of $\bar{Q}(x+\epsilon_{m,n})$ by $\bar{Q}(x)$ we find the remainder

\begin{eqnarray}
r_3(x)=\frac{2\epsilon_{m,n}}{(2w_1+5\epsilon_{m,n}^2)^2}&\cdot&\bigg( \big( 2 x+\epsilon_{m,n})(-8w_1 \bar{w}_3^{(k|m,n)}-20 \bar{w}_3^{(k|m,n)}\epsilon_{m,n}^2 \nonumber \\ &+& 2 w_1^3+9w_2^2+9w_1^2\epsilon_{m,n}^2+12w_1\epsilon_{m,n}^4+5\epsilon_{m,n}^6 \big) \\ \nonumber
&-&2 w_2 (12\bar{w}_3^{(k|m,n)}+w_1^2+8w_1 \epsilon_{m,n}^2+7\epsilon_{m,n}^4 ) \bigg)
\end{eqnarray}

and since $r_3(a_2(\Delta w_2^{(k|m,n)} ))=0$ we get

\begin{equation}
\frac{\Delta w_2^{(k|m,n)} }{\Delta w_1^{(m,n)} }=\frac{-2 w_2 (12\bar{w}_3^{(k|m,n)}+w_1^2+8w_1 \epsilon_{m,n}^2+7\epsilon_{m,n}^4 )}{(-8w_1 \bar{w}_3^{(k|m,n)}-20 \bar{w}_3^{(k|m,n)}\epsilon_{m,n}^2+2 w_1^3+9w_2^2+9w_1^2\epsilon_{m,n}^2+12w_1\epsilon_{m,n}^4+5\epsilon_{m,n}^6) }.
\end{equation}

The polynomial $\Delta Q(x)$ in this case is

\begin{equation}
\Delta Q(x)=x^2-\frac{\Delta w_2^{(k|m,n)}}{\Delta w_1^{(m,n)} }x+\frac{\Delta w_3^{(k|m,n)}}{\Delta w_1^{(m,n)} }.
\end{equation}

The recurrence relation is

\begin{equation}
\bar{Z}^{({\rm R})}({\bm \omega},w_3)=1+\sum_{k=1}^{3}\sum_{m,n=1}^{\infty} \frac{q^{mn}J^{(m,n)}({\bm \omega},\bar{w}_{3}^{(k|m,n)})}{(w_3-\bar{w}^{(k|m,n)}_3)}\frac{ \mathcal{P}^{(m,n)}_{4,{\rm R}} ({\bm \omega},\bar{w}_3^{(k|m,n)})}{\mathcal{P}^{(m,n)}_{4}({\bm \omega},\bar{w}_3^{(k|m,n)})}\bar{Z}^{({\rm R})}({\hat{\bm \omega}^{(k|m,n)} ,\hat{w}_3^{(k|m,n)})}),
\end{equation}

where $\bar{w}^{(k|m,n)}_3$ are the roots of equation (<ref>),

\begin{eqnarray}
J^{(m,n)}({\bm \omega},\bar{w}_{3}^{(k|m,n)})&=&(-1)\frac{1}{3}\epsilon_{m,n}(w_1+2\frac{\Delta w_3^{(k|m,n)} }{\Delta w_1^{(m,n)} }+\epsilon^2) \nonumber \\ &\cdot& (w_1^2-4w_1\frac{\Delta w_3^{(k|m,n)} }{\Delta w_1^{(m,n)} }+4\frac{(\Delta w_3^{(k|m,n)})^2}{(\Delta w_1^{(m,n)})^2}+8w_1 \epsilon_{m,n}^2-4\frac{\Delta w_3^{(k|m,n)}}{\Delta w_1^{(m,n)}}\epsilon_{m,n}^2+7\epsilon_{m,n}^4) \nonumber \\&\cdot&(w_1^2-4w_1\frac{\Delta w_3^{(k|m,n)} }{\Delta w_1^{(m,n)} }+4\frac{(\Delta w_3^{(k|m,n)})^2}{(\Delta w_1^{(m,n)})^2}+\frac{4}{3}w_1 \epsilon_{m,n}^2-\frac{4}{3}\frac{\Delta w_3^{(k|m,n)} }{\Delta w_1^{(m,n)} }\epsilon_{m,n}^2+\frac{1}{3}\epsilon_{m,n}^4) \nonumber\\
&\cdot&\bigg( 2 w_1^4+18 w_1 w_2^2-32 w_1^2 \bar{w}^{(k|m,n)}_3+96 (\bar{w}^{(k|m,n)}_3)^2+  4w_1^3\epsilon_{m,n}^2+27w_2^2\epsilon_{m,n}^2 \nonumber \\ &-&48w_1\bar{w}^{(k|m,n)}_3\epsilon_{m,n}^2+3 w_1^2 \epsilon_{m,n}^4-28\bar{w}^{(k|m,n)}_3 \epsilon_{m,n}^4+2w_1\epsilon_{m,n}^6+\epsilon_{m,n}^8 \bigg)^{-1}
\end{eqnarray}

\begin{equation}
\mathcal{P}^{(m,n)}_{4}({\bm \omega},\bar{w}_3^{(k|m,n)})=\prod_{i=-m}^{m-1} \sideset{}{'} \prod_{j=-n}^{n-1} \epsilon_{i,j} \cdot \prod_{i=1}^{m} \prod_{j=1}^{n}\left( (w_1-2\frac{\Delta w_3^{(k|m,n)} }{\Delta w_1^{(m,n)} }+\epsilon_{i,j}\epsilon_{m-i,n-j})^2 -\frac{(\Delta w_2^{(k|m,n)} )^2}{(\Delta w_1^{(m,n)} )^2}\epsilon_{m-2i,n-2j}^2 \right),
\end{equation}

\begin{eqnarray}
\mathcal{P}^{(m,n)}_{4,{\rm adj}}({\bm \omega},\bar{w}_3^{(k|m,n)})&=& \prod_{i=-m}^{m-1}\prod_{j=-n}^{n-1} (\epsilon_{i,j}-M)  \\
&\cdot& \prod_{i=1}^{m} \prod_{j=1}^{n}\bigg( (w_1-2\frac{\Delta w_3^{(k|m,n)} }{\Delta w_1^{(m,n)} }+(\epsilon_{i,j}+M)(\epsilon_{m-i,n-j}-M))^2 \nonumber \\ &-& \frac{(\Delta w_2^{(k|m,n)} )^2}{(\Delta w_1^{(m,n)} )^2}(\epsilon_{m-2i,n-2j}-2M)^2 \bigg), \nonumber
\end{eqnarray}

and $\mathcal{P}_{4, {\rm fund}}^{(m,n)}({\bm \omega}, \bar{w}_{N-1}^{(k|m,n)})$ is given by (<ref>).

§ SUMMARY OF RESULTS

In this Section we collect the main relations obtained throughout the paper.
We showed that the instanton partition function has poles only at the points $a_{uv}=\epsilon_{m,n}$ with $m\cdot n>0$, the poles are simple and for the positive $m$, $n$ the residue of the instanton partition function with respect to the variable $a_{uv}$ can be expressed via its value at the point distinguished by the partial Weyl permutation between $a_u$ and $a_v$.

\begin{equation}
{\rm Res}_{a_{uv}=\epsilon_{m,n}} Z^{({\rm R})}( {\bf a})=q^{m n}\frac{\mathcal{P}^{(uv)}_{N,{\rm R}}(m,n|{\bf a})}{\mathcal{P}^{(uv)}_{N}(m,n|{\bf a}) }Z^{({\rm R})}(\hat{{\bf a}}^{(uv)}),
\end{equation}

\begin{eqnarray}
\mathcal{P}^{(uv)}_{N}(m,n|{\bf a})=
\prod_{i=-m}^{m-1}\sideset{}{'} \prod_{j=-n}^{n-1} \epsilon_{i,j} \cdot\underset{w \neq u, \, v}{\prod_{w=1}^{N}} \prod_{i=1}^m\prod_{j=1}^n\left[  (a_{vw}+\epsilon_{i,j})(-a_{uw}+\epsilon_{i,j}) \right],
\end{eqnarray}

for the pure theory we found a trivial numerator

\begin{equation}
\mathcal{P}^{(uv)}_{N,{\rm 0}}=1,
\end{equation}

for a theory with the adjoint matter there is a polynomial

\begin{eqnarray}
\mathcal{P}^{(uv)}_{N,{\rm adj}}(m,n|{\bf a})=    \prod_{i=-m}^{m-1}\prod_{j=-n}^{n-1}  \left(\epsilon_{i,j}-M \right)  \cdot\underset{w \neq u, \, v}{\prod_{w=1}^{N}} \prod_{i=1}^m\prod_{j=1}^n\left[  (a_{vw}+\epsilon_{i,j}+M)(-a_{uw}+\epsilon_{i,j}+M) \right],
\end{eqnarray}

and finally for the theory with fundamental and anti-fundamental multiplets the polynomial is

\begin{eqnarray}
\mathcal{P}^{(uv)}_{N,{\rm fund}}(m,n|{\bf a})=\prod_{i=1}^{m}\prod_{j=1}^{n} \Bigg[ \prod_{t=1}^{N_f}   \left(-\frac{1}{2}\epsilon_{m,n}+\epsilon_{i,j}-m_t-\frac{1}{2}\sum_{\underset{w \neq u,v}{w=1}}^N a_w \right) \cdot \\
\prod_{t=1}^{N_a}  \left(-\frac{1}{2}\epsilon_{m,n}+\epsilon_{i,j}+m_t+\frac{1}{2}\sum_{\underset{w \neq u,v}{w=1}}^N a_w \right) \Bigg]. \nonumber
\end{eqnarray}

One can easily see from the proof of these relations that these hypermultiplets can be considered together, and the polynomial in the numerator will be simply the product of the polynomials above[It is actually not difficult to see that the residue formula holds for the matter in any representation, but in the general case the polynomial in the numerator is too long to write in a paper.].
For the residue at the point $a_{uv}=\epsilon_{-m,-n}$ there is an additional minus sign.
There is an equivalent form in terms of the full partition function

\begin{equation}
\lim_{a_{uv} \rightarrow \epsilon_{m,n}}\frac{\mathcal{Z}^{({\rm R})}({\bf a})}{\mathcal{Z}^{({\rm R})}({\bf \hat{a}}^{(uv)})}=-{\rm Sign}(\epsilon_1), \quad m,n \in \mathbb{Z}\setminus \{0\}
\end{equation}

for the permutation of $\epsilon_2$-coefficients (and $-{\rm Sign}(\epsilon_2)$ for the permutation of $\epsilon_1$-coefficients).
For the pure theory we found a recurrence relation for the instanton partition function in terms of the variables $a_{uv}$.

\begin{eqnarray}
{Z}^{({\rm 0})}({\bf a})=  1+
\sum_{w=2}^{N}\sum_{m,\, n=1}^{\infty} \frac{q^{mn} {Z}^{({\rm 0})}(\hat{\bf a}^{(1w)})}{(a_{1N}-a_{wN}+\epsilon_{m,n})(a_{1N}-a_{wN}-\epsilon_{m,n})} \frac{2 \epsilon_{m,n}\mathcal{P}^{(1w)}_{N,{\rm 0}}(m,n|{\bf a})}{\mathcal{P}^{(1w)}_{N}(m,n|{\bf a})}.
\end{eqnarray}

To write the recurrence relation for theories with matter hypermultiplets we switched to the symmetrical variables defined as

\begin{equation}
\omega_{l}=\sum_{i_1<\ldots<i_k} a_{i_{1}}\cdot \ldots \cdot a_{i_{l+1}} \, \quad l=1, \ldots, N-1.
\end{equation}

The recurrence relation for the normalised instanton partition function $\bar{Z}^{\rm (R)}$ is

\begin{equation}
\bar{Z}^{\rm (R)}({\bm \omega},w_{N-1})=1+\sum_{k=1}^{N-1}\sum_{m,n=1}^{\infty} \frac{q^{mn} J^{(m,n)}({\bm \omega},\bar{w}_{N-1}^{(k|m,n)}) }{(w_{N-1}-\bar{w}^{(k|m,n)}_{N-1})}\frac{\mathcal{P}^{(m,n)}_{N, {\rm R}}({\bm \omega},\bar{w}_{N-1}^{(k|m,n)})}{\mathcal{P}^{(m,n)}_{N}({\bm \omega},\bar{w}_{N-1}^{(k|m,n)})}\bar{Z}^{\rm (R)}({\bf \hat{w}}^{(k|m,n)}),
\end{equation}

where $\bar{w}^{(k|m,n)}_{N-1}$ are the roots of equation (<ref>), the polynomials $\mathcal{P}^{(m,n)}_{N}$, $\mathcal{P}^{(m,n)}_{N,{\rm R}}$ and Jacobian $J^{(m,n)}$ are defined in (<ref>), (<ref>), (<ref>) and (<ref>).
The normalisation constants are defined by the behaviour at infinity, which is the following

\begin{equation}
Z^{({\rm 0})} \underset{w_{N-1} \rightarrow \infty}{\longrightarrow } 1.
\end{equation}

\begin{equation}
Z^{({\rm fund})} \underset{w_{N-1} \rightarrow \infty}{\longrightarrow } 1, \quad N_f+N_a<2(N-1).
\end{equation}

\begin{equation}
Z^{({\rm fund})}\underset{w_{N-1} \rightarrow \infty}{\longrightarrow }{\rm exp}\left( (-1)^{N-1}\frac{q}{N \epsilon_1 \epsilon_2}\right), \quad N_f+N_a=2(N-1).
\end{equation}

\begin{eqnarray}
Z^{({\rm adj})}\underset{w_{N-1} \rightarrow \infty}{\longrightarrow }\left( q^{-\frac{1}{24}}\eta(q)\right)^{-N\left(1+\frac{M}{\epsilon_1} \right)\left(1+\frac{M}{\epsilon_2} \right)}.
\end{eqnarray}

[1] Nikita Nekrasov. (2003). Seiberg-Witten Prepotential from Instanton Counting. Advances in Theoretical and Mathematical Physics. 7. 10.4310/ATMP.2003.v7.n5.a4.

[2] Rubik Poghossian. (2009). Recursion relations in CFT and N=2 SYM theory. Journal of High Energy Physics. 2009. 10.1088/1126-6708/2009/12/038.

[3] A.B. Zamolodchikov. (1984). Conformal symmetry in two dimensions: An explicit recurrence formula for the conformal partial wave amplitude. Commun.Math. Phys. 96, 419–422. 10.1007/BF01214585

[4] M. Bershtein., G. Bonelli, M. Ronzani and A. Tanzini. (2016). Gauge theories on compact toric surfaces, conformal field theories and equivariant Donaldson invariants. J.Geom.Phys. 118 (2017) 40-50. 10.1016/j.geomphys.2017.01.012

[5] G. Bonelli, F. Fucito, J. F. Morales, M. Ronzani, E. Sysoeva, A. Tanzini. (2021). Gauge theories on compact toric manifolds. Letters in Mathematical Physics. 111. 10.1007/s11005-021-01419-9.

[6] Rubik Poghossian. (2017). Recurrence relations for the W3 conformal blocks and N=2 SYM partition functions. J. High Energ. Phys. 2017, 53. 10.1007/JHEP11(2017)053

[7] Nikita Nekrasov, Andrei Okounkov. (2006). Seiberg-Witten Theory and Random Partitions. Etingof, P., Retakh, V., Singer, I.M. (eds) The Unity of Mathematics. Progress in Mathematics, vol 244. Birkhäuser Boston. 10.1007/0-8176-4467-9_15

[8] U. Bruzzo, F. Fucito, J. F. Morales, A. Tanzini. (2002). Multi-Instanton Calculus and Equivariant Cohomology. Journal of High Energy Physics. 2003. 10.1088/1126-6708/2003/05/054.

[9] Hiraku Nakajima, Kota Yoshioka. (2003). Instanton counting on blowup. I. 4-dimensional pure gauge theory. Inventiones mathematicae. 162. 10.1007/s00222-005-0444-1.

[10] Hiraku Nakajima, Kota Yoshioka. (2003). Lectures on instanton counting. CRM Workshop on Algebraic Structures and Moduli Spaces. e-Print: math/0311058 [math.AG]

[11] Hiraku Nakajima. (1999). Lectures on Hilbert Schemes of Points on Surfaces. AMS University Lecture Series, ISBN 0-8218-1956-9.

[12] L. Göttsche, H. Nakajima, K. Yoshioka. (2008). Instanton counting and Donaldson invariants. J.Diff.Geom. 80 (2008) 3, 343-390. 10.4310/jdg/1226090481.

[13] M.F. Atiyah, R.Bott. (1984). The moment map and equivariant cohomology. Topology. Volume 23, Issue 1, Pages 1-28. https://doi.org/10.1016/0040-9383(84)90021-1

[14]
Alexander Klyachko. (1991). Vector bundles and torsion free sheaves on the projective plane. Preprint MPI/91-59.

[15]
Allen Knutson, Eric Sharp. (1998). Sheaves on toric varieties for physics.
Advances in Theoretical and Mathematical Physics 2 873. 10.4310/ATMP.1998.v2.n4.a6.

[16] Guo-Niu Han. (2008). The Nekrasov-Okounkov hook length formula: Refinement, elementary proof, extension and applications. Annales de l’institut Fourier. 60. 10.5802/aif.2515.

[18] I. M. Gelfand, M. M. Kapranov, A. V. Zelevinsky. (1994). Discriminants, resultants, and multidimensional determinants. Boston: Birkhäuser, ISBN 978-0-8176-3660-9.