{'fraction_non_alphanumeric': 0.02964656497132791, 'fraction_numerical': 0.008620366552977774, 'mean_word_length': 4.857738748627882, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Prediction of Atomization Energy Using Graph Kernel and Active Learning

Data-driven prediction of molecular properties presents unique challenges to the design of machine learning methods concerning data structure/dimensionality, symmetry adaption, and confidence management. In this paper, we present a kernel-based pipeline that can learn and predict the atomization energy of molecules with high accuracy. The framework employs Gaussian process regression to perform predictions based on the similarity between molecules, which is computed using the marginalized graph kernel. We discuss why the graph kernel, paired with a graph representation of the molecules, is particularly useful for predicting extensive properties. We demonstrate that using an active learning procedure, the proposed method can achieve a mean absolute error less than 1.0 kcal/mol on the QM7 data set using as few as 1200 training samples and 1 hour of training time. This is a demonstration, in contrast to common believes, that regression models based on kernel methods can be simultaneously accurate and fast predictors.


Introduction
The chemistry and materials communities have embraced data science and machine learning to bring about revolutionizing solutions to long-standing challenges in molecular modeling, optimal experiment design, and high-throughput structure screening [1][2][3][4][5]. One particularly promising application of machine learning techniques is to train predictive models for molecular properties that are otherwise only available through expensive dynamics simulations and quantum mechanical calculations. This type of regression tasks typically entail quality requirements concerning: 1. accuracy: While the atomization energy of drug-like ligands are on or above the order of 10 3 kcal/mol, it is the binding energy, i.e. the difference between the energy of two ligand-protein complexes, that is sought for applications such as drug screening. It is therefore expected that a good predictive model should carry less than 1 kcal/mol absolute error, or 10 −3 relative error; 2. smoothness: The prediction must be at least twice differentiable with regard to the molecular geometry for a model to be usable as a force field for dynamics simulations; 3. outlier handling: Inevitably, a trained model could encounter outliers which triggers extrapolation far beyond the training domain. Should the predictive model fail to maintain a desired level of accuracy in this situation, a fail-safe mechanism should kick in to let the end user be aware of the potential quality degradation and take remedy actions; 4. symmetry adaptation: The predictive model should be able to recognize configurations and conformations that are merely different by a permutation and/or a rigid body transformation, which only cause covariant changes to the target property. Symmetry adaptation could take advantage of such knowledge to simultaneously improve accuracy, smoothness, and generalization.
Gaussian process regression (GPR) [6] is a robust regression method that features high accuracy, strong smoothness guarantee, and built-in uncertainty estimation. It is a Bayesian inference method that models the similarity among the sample points as covariance between random variables indexed by the input space. The probabilistic model generated by GPR can make not only a point estimate about an unknown sample, but also the associated posterior variance. Since the predictive variance vanishes around training samples and grows in regions that are not supported by the sample data, it can be conveniently used as a measure of predictive uncertainty for outlier detection. As shown in fig. 1, the GPR predictive uncertainty can be readily used to implement an active-learning protocol where test samples with large uncertainties are sought and added to the training set in an iterative process. uncertainty Figure 1: An active learning procedure is carried out to learn the Lennard-Jones potential between two atoms. By iteratively adding the point with the largest predictive uncertainty into the training set, the GPR model quickly converges to the ground truth using a relative small amount of samples.
A critical component in each GPR model is a function that estimates the covariance of target values between pairs of data samples. In the context of molecular property prediction, a covariance function translates structural similarity into correlations between the target properties using the similar property principle [7]. In fact, covariance functions belong to the family of similarity kernels, which are bivariate functions that compute the closeness of two vectors in some normed vector space. A kernel is also a convenient place for encapsulating symmetry adaption mechanisms that are transparent to subsequent GPR computations. The two most wellknown kernels on Euclidean spaces are probably the dot product kernel and the square exponential kernel, which defines similarity based on angle and distance, respectively.
Meanwhile, designing valid kernels on non-Euclidean spaces as those spanned by molecular configurations and conformations is an inherently challenging task, because the kernel must satisfy certain properties such as symmetry and positive-definiteness to ensure convexness of the resulting machine learning model [6]. An approach that is commonly used in existing implementations [8,9] is to apply an Euclidean kernel to fixed length feature vectors that are computed using a molecular fingerprinting algorithm. Examples of this approach include the atom-centered symmetry function fingerprints [10], the Coulomb matrix [11], the densityencoded canonically-aligned fingerprint [12], the encoded bonds descriptor [13] and so on. The smooth overlap of atomic positions kernel [14] used in practice also effectively computes the inner product using a fixed-length vector representation of the spherical power spectra of atomic mass density.
Despite the success of the feature vector approach, it is important to point out that squashing a molecule, which is intrinsically a non-linear structure with a variable number of degrees of freedom, into a fixed-length feature vector could be a detour that compromises the generalizability of machine learning models across the diverse chemical space. This is likely the reason why fingerprints are prevalently used to describe local atomistic neighborhoods that possess bounded complexities. An intrinsically non-linear data structure, such as a graph, could be a powerful and intuitive representation to capture all the geometric and topological information of whole molecules [15,16]. Moreover, there has been extensive work on the design of graph kernels, which compute the similarity between graphs without resorting to a fixed length intermediate presentation [17][18][19][20][21][22][23].
A particularly attractive feature of the graph kernels is their ability to adapt to graphs with arbitrary numbers of vertices and edges.
In this work, we present how a particular type of graph kernel, i.e. the marginalized graph kernel [24], can be integrated into a Gaussian process regression-based active learning pipeline for molecular atomization energy prediction. In section 2, we introduce the mathematical formulations for the techniques used in this study, including the marginalized graph kernel, Gaussian process regression, and active learning. In section 3, we present computational setup and results. We conclude the work in section 4.

Overview
In fig. 2, we present an overview of the machine learning pipeline proposed in this work. Given a set of training molecules, a GPR model can be constructed using the pairwise similarity matrix among the molecules to fit for the assocaited target values. To compute the similarity, the molecules are converted into graphs with vertices labeled by the atoms and edges encoding interatomic distances. The marginalized graph kernel is then applied to average over the similarities of all paths generated from simultaneous random walks on each pair of graphs.
Predictions for the energy of new molecules can be made using the pairwise similarity matrix among the new molecules and the cross-similarity matrix between the new molecule and the training molecules. In the subsections below, we will explain in detail each component of the procedure. Left: given the training molecules and unknown molecules, the pairwise similarity matrices between and within each data set are computed using the marginalized graph kernel. The training set self-similarity matrix and target values are used to construct a GPR model, which makes predictions for the unknown samples using the training-unknown cross-similarity matrix and the unknown-unknown self-similarity matrix. Right: molecules are first converted into labeled graphs with edges weighted by interatomic distances. The marginalized graph kernel defines the similarity between two graphs as the expectation of the similarity between all possible simultaneous random walk paths generated on the product graph.

Preliminaries
We use lower case letters, e.g. a, in bold font to denote vectors, and upper case letters in bold font, e.g. A, for matrices. By default, we assume vectors are column vectors. We use diag(a) to denote a diagonal matrix whose diagonal elements are specified by a, and use I to refer to the identity matrix.

Definition 1 Kronecker product
Given matrices A ∈ R m1×n1 and B ∈ R m2×n2 , the Kronecker product A ⊗ B ∈ R m1m2×n1n2 is defined as: The element-wise product, also known as the Hadamard product, between two matrices of the same size A, B ∈

Definition 3 Undirected graph An undirected graph G is a discrete structure consisting of a set of uniquely-indexed vertices
and a set of undirected edges E ⊂ V × V . The vertices and edges may be labeled using elements from label sets Σ v and Σ e , respectively.

Definition 4 Weighted graph
In a weighted graph, each edge and v j are not connected by an edge. An unweighted graph can be regarded as a specialized weighted graph where w ij = 1 between each pair of (v i , v j ) connected by an edge and 0 elsewhere.

Definition 5 Walk on graph
Two vertices are neighbors if they are connected by an edge. A walk on a graph is a sequence of vertices such that all pairs of adjacent vertices are neighbors.

Definition 6 Adjacency matrix
The adjacency matrix of a graph of n vertices is a matrix A ∈ R n×n with A ij = w ij . The adjacency matrices of undirected graphs are symmetric since w ij ≡ w ji .
The row-normalized adjacency matrix P := diag(d) −1 A, where d i := ∑ j A ij is the vertex degree vector, can be used as the transition matrix of a Markovian random walk process. In other words, P ij could be interpreted as the conditional probability for a random walker at vertex v i to jump to v j during a single step.

Graph Representations for Molecules
The practice of using labeled graphs, with the exemplary ball-and-stick model, to represent molecules has gained popularity well before the era of machine learning. In this work, we represent a molecule of n atoms as an undirected graph Each edge e ij ∈ R between vertices i and j is labeled by the distance between the atoms, while its weight w ij is set by a spatial adjacency rule A(r i , r j ), whose specific form will be discussed in section 3.1. The adjacency matrix of the molecular graph is then Note that the edges are often supersets of the collection of covalent bonds in a molecule.

Marginalized Graph Kernel
The marginalized graph kernel K(G, G ′ ) [24] defines the overall similarity between two graphs as the expectation of the path similarity between all pairs of paths that can be obtained by performing simultaneous random walks on two graphs G and G ′ . The explicit formula for the path similarity expectation takes the form: . (1) Here, ℓ is the length of the path, h and h ′ are paths on the graphs represented by length-ℓ vectors of vertex labels, p s (·) is the starting probability of the random walk on each vertex, p q (·) is the stopping probability of the random walk on each vertex at any given step, p t (·, ·) is the transition probability between a pair of vertices, K v (·, ·) is an elementary kernel that computes the similarity between two vertices, K e (·, ·) is another elementary kernel that computes the similarity between pairs of bonds.
Despite the overwhelming appearance of the infinite summation, eq. (1) can be reformulated under the spirit of dynamic programming as: where R ∞ is the solution to the linear system: with Equations (2) to (4) in fact exhibit a Kronecker product structure, which can be readily recognized in matrix form: with r ∞ being the solution to the linear system where P is the transition probability matrix as defined in definition 6, where the augmented degreẽ ensure that the transition and stopping probabilities on each node sum to 1;

the corresponding vectors and matrices for G ′ ;
Kv ⊗ is the kernelized Kronecker product between v and v ′ where K v (·, ·) replaces the product operation between the vertex labels; Ke ⊗ is the kernelized Kronecker product between E and E ′ where K e (·, ·) replaces the product operation between the edge labels.
For clarity of the discussion below, we denote , To solve eq. (6), first observe that only the product V × r ∞ as a whole is needed to compute K(G, G ′ ). We can thus rearrange eq. (6) to form a symmetric linear system.
The full expression for the marginalized graph kernel in matrix form is then 'Stray' atoms, i.e. atoms that are not neighbors to any other atoms, have d i = 0. This pathological situation can render the linear system singular, and is resolved by setting the corresponding d i = q i = 1. This is equivalent to immediately terminating any random walk path originating from the stray atoms.
with the predictive uncertainty given by the posterior covariance matrix Here, are the training set covariance matrix, training-test cross-covariance matrix, and test set covariance matrix, respectively, that are obtained by iterating the graph kernel over all pairs of molecules from the corresponding datasets.
The marginalized graph kernel is particularly suitable for predicting extensive properties such as molecular energy. To see why, first recall that the kernel possesses a vertex-wise summation structure ∑ h,h ′ · · · as manifested in eq. (2), and that R ∞ (h 1 , h ′ 1 ) can be alternatively expressed as with Equation (2), eq. (18), and eq. (20) together indicate that each of the terms effectively computes the expectation of the similarity between all random walk paths originating from a pair of vertices h 1 and h ′ 1 . We could thus regard it as a vertex-wise atomistic neighborhood similarity kernel, denoted as The overall similarity between two graphs as computed by the marginalized graph kernel is thus merely a summation over the similarity between all pairs of vertices, i.e.
To ensure extensiveness of the energy prediction, we can hypothetically predict the individual contributions of each atom to the total energy of an unknown molecule, and then sum up the contributions. We assume that there exists some energy decomposition algorithm (EDA) [27] that could localize the energy of a molecule among its atoms as {e 1 , e 2 , . . .} = EDA(M, E) subject to ∑ i e i = E. Note that the EDA scheme is purely formal and is never explicitly used, as will be shown shortly. Using κ(·, ·) as a covariance function, we obtain the following GPR formulation: . .
where the superscripts correspond to molecule ids.
Under the multivariate normality condition, the total energy E, being the sum of individual random variables representing atomic energies, is itself a normally distributed random variable, . Note that this exhibits the same pairwise summation structure as in the marginalized graph kernel. Thus, by summing up the cross-covariance submatrices between pairs of molecules in eq. (23) and substituting K(·, ·) for ∑ ij κ(v i , v ′ j ), we could define an alternative GPR model that bypasses energy decomposition: This leads exactly back to eq. (13). Thus, we have shown that when only the total energy is of interest, a GPR prediction for the energy of an entire molecule using the marginalized graph kernel is equivalent to a summation of the GPR predictions for the localized atomic contributions using an atomistic neighborhood kernel.
The former approach, however, avoids the potentially costly and controversial explicit energy decomposition, while still enjoying automatic scaling of the predicted energy with respect to molecule size. This approach can also speed up GPR computation by creating covariance matrices whose size is proportional to the number of molecules, rather than the number of total atoms, in the data set.

Active Learning Protocol
Given training set S, test set T , an acquisition function Q(·, ·) that measures the learning value of a sample given a GPR model, active learning can be performed using the following protocol: 1: repeat 2: G ← TrainGPR(S), 3: α ← argmax i∈T Q(G, i), 4: S ← S ∪ {α},
A natural choice of the acquisition function would be prediction error, which is the absolute difference between the GPR prediction and the truth value of the target function at the sample points. However, since the ground truth may not be unknown or is too difficult to obtain for every point in practice, the GPR predictive uncertainty could be used as an alternative acquisition function. In this case, the 'unsupervised' learning procedure can proceed even without using any target function value by exploiting the fact that the GPR posterior variance does not depend on the target value.

Kernel Specification and Data Set
The vertex elementary kernel that we use for the marginalized graph kernel is an elevated Kronecker delta function on element symbols: The edge elementary kernel is a square exponential function on edge lengths, which evaluates to 1 if two edges are of the same length and smoothly transitions to 0 as the difference in length grows: The adjacency rule that computes the weights for each edge also assumes a square exponential form: where σ ij , as given in the appendix, are element-wise length scale parameters derived from typical bonding lengths. A uniform starting probability p s (·) ≡ s and a uniform stopping probability p q (·) ≡ q are used across all vertices.
In fig. 3, we visualize in matrix form the atomistic neighborhood similarities V × r ∞ between a methoxyethane (CH 3 CH 2 OCH 3 ) molecule and a 2-ethoxyethanol (CH 3 CH 2 OCH 2 CH 2 OH) molecule, computed with ν = 0.25, ζ = 1, λ = 0.02, s = 1, q = 0.01. Note that the marginalized graph kernel can, for example, differentiate between atoms of the same element yet are embedded in different local chemical environments.
The QM7 dataset, which contains 7165 molecules with up to 7 heavy atoms and 23 total atoms along with their atomization energies, is used as the data set for benchmarking purposes [11,28].

Hyperparameter Selection
The complete marginalized graph kernel as specified in section 3.1 is controlled by 5 hyperparameters ν, ζ, λ, s, and q, each of which possesses a unique physical interpretation: • ν is the baseline, or 'prior' , similarity that we assign between atoms of difference elements. ν must be within (0, 1) to ensure that the linear system in eq. (11) is positive-definite.
• ζ controls how quickly the weight of an edge decays with regard to increasing interatomic distance.
It affects the edge density of the molecular graph and is one of the two factors that determines the neighborhood of for each atom.
• λ control the sensitivity of the kernel with respect to differences in edge length. Since the order of a covalent bond is strongly correlated with its length, λ can be effectively used to discriminate different types of bonds.
• s affects the preference among pairs of vertices when sampling random walk paths, and could thus be tor, which is a weighted average (marginalization) of the similarity between all paths generated by simultaneous random walks originating from the two atoms. The matrices effectively visualizes the atomistic neighborhood similarities with darker colors indicate higher similarity. Note for example the difference in similarity between the oxygen atoms in the ether bond to itself and to the oxygen atom in the alcohol group.
thought as a factor of significance. From the perspective of GPR modeling, the magnitude of s is the a priori standard deviation of atomic contributions to molecular energy. In practice, we may abandon the probabilistic interpretation of s and allow it to assume values that sum up to above unity.
• q determines the average length of the random walk paths being sampled, and together with ζ determines the range of the atomistic neighborhoods that are effectively compared between atoms.
Maximum likelihood estimation (MLE) is commonly regarded as the 'default' method for selecting hyper- The extended vertical tail of the error-likelihood joint distribution indicate that GPR models with high log-likelihood could still possess good predictive power. Lower: Despite that the correlation between the predictive uncertainty and error generally decrease as the accuracy of the GPR models improve, it is still always possible to find sets of hyperparameters that lie close to the upper left corner of the plots, which are the regions corresponding to relatively low errors and high uncertainty-error correlation. The latter property is crucial for effective unsupervised active learning. The densities are estimated using a grid of 5400 hyperparameters sets. 300 molecules are randomly selected from the QM7 dataset as the training set, while the remaining 6865 molecules form the test set.
We further examine the conditional distribution of the mean absolute error, log-log-likelihood, and erroruncertainty correlation on specific values of the hyperparameters. As shown by fig. 5, the edge kernel length scale parameter λ appears to be the single most sensitive hyperparameter, followed by ζ, ν, and q, that affects both predictive accuracy and error-uncertainty correlation. A uniform starting probability s seems to only affect the model likelihood, yet it controls the scaling of the magnitude of the predictive uncertainty, which should be at least on the same order of magnitude as the true error. In the subsequent section, we use ζ = 1, ν = 0.3, λ = 0.05, s = 100, and q = 0.05 for all GPR modeling tasks.

Log-log-likelihood
Error-uncertainty correlation  Figure 5: Shown here is the distribution of the mean absolute error, log-log-likelihood, and error-uncertainty correlation of GPR models conditioned on each specific value of the hyperparameters. The accuracy and error-uncertainty correlation of the resulting GPR model are influenced mostly by ζ, ν, λ, and q, while the log-likelihood function is most strongly affected by s and q. The choice of the optimal hyperparameter is apparently a subtle problem due to the non-trial dependency between accuracy and error-uncertainty correlation.

Energy Prediction
As a first benchmark of prediction accuracy, we use a random selection of N molecules from the QM7 dataset to train GPR models and measure the mean absolute error (MAE) and root-mean square error (RMSE) on energy prediction for the rest of the molecules. As shown in the left panel of fig. 6, with a randomly selected training set, the presented method delivers an accuracy as low as 0.92 kcal/mol when trained on 5000 molecules, which is as good as the best of all existing method [13,[29][30][31][32]. It also achieves an MAE as low as 1.   Figure 6: Left: Comparison of the mean absolute error and root-mean square error of the presented method and previous methods on predicting the atomization energy of molecules in the QM7 dataset. Each data point for the presented method are derived using 15 randomly chosen training sets. Reference designation: FHB [32], CGL [13], HMB [29], HBR [30], SAC [31]. Right: Comparison of the performance between GPR models built from random training sets and those built from an active learning process.
We further examined the performance of the method for carring out the active learning protocol as described in section 2.4.2. In each parallel test, the process starts with a seeding set of 50 training samples and iteratively adds the test sample with the highest acquisition score into the training sets until the training set size reaches 1000. As shown in the right panel of fig. 6, supevised active learning, which uses the true predictive error as the acquisition function, leads to GPR models with significantly lower MAE and RMSE than those trained on a randomly selected training set. Meanwhile, unsupervised active learning, which uses the predictive uncertainty as the acquisition function, seems to perform on par with those trained on random training sets in terms of MAE, but does achieve significantly lower RMSE.

Conclusion
In this paper, we presented a new machine learning pipeline, which integrates the marginalized graph kernel, the Gaussian process regression method, and an active learning protocol, for predicting molecular atomization energies. The method achieves excellent accuracy while using significantly smaller training sets as compared to those used in previous studies. We show that the marginalized graph kernel defines a molecular similarity metric that is computed from both topological and geometric information, and that the kernel can naturally adapt to molecules of difference topology and sizes. GPR models created using the marginalized graph kernel can produce energy predictions that scale properly with molecule sizes, while bypassing any explicit energy decomposition procedure. We further demonstrate that the convolutional, or pairwise summation, structure of the kernel is crucial in enabling the prediction of extensive properties that automatically scales with molecule size. It is worth noting that the work is inspired by the localized random walk graph kernel for molecules proposed by Ferré et al [32], which is essentially a localized marginalized graph kernel equipped with trivial vertex and edge kernels.
We expect that the presented method could benefit from a more thorough design and selection of the hyperparameters involving, for example, pairwise priors that are tailored between each type of atoms and bonds. A wider adoption of the graph kernel entail a more efficient computer implementations that can work around the quartic asymptotic computational cost for solving the Kronecker product system induced by large molecules such as polymers, proteins, and polynucleotides.