Impact of non-normal error distributions on the benchmarking and ranking of quantum machine learning models

Reliable use of statistical tests of normality require typically at least N = 100 points [8, 9]. For a given sample size, the Shapiro-Wilk W statistics has been shown to have good properties [8], and it is used in this study. The values of W range between 0 and 1, and values of W ≃ 1 are in favor of the normality of the sample. If W lies below a critical value W_c depending on the sample size and the chosen level of type I errors α (typically 0.05), the normality hypothesis cannot be assumed to hold [9].

It might also be useful to assess normality by visual tools: normal quantile-quantile plots (QQ-plots) [4], where the quantiles of the scaled and centered errors sample is plotted against the theoretical quantiles of a standard normal distribution (in the normal case, all points should lie over the unit line); or comparison of the histogram of errors with a gaussian curve having the same mean, estimated by the mean signed error (MSE), and same standard deviation, estimated by the RMSD.

Two other statistics are helpful in characterizing the departure from normality. The skewness (Skew), or third standardized moment of the distribution, quantifies its asymmetry (Skew  =  0 for a symmetric distribution). The kurtosis (Kurt), or fourth standardized moment, quantifies the concentration of data in the tails of the distribution. Kurtosis of a normal distribution is equal to 3; distributions with excess kurtosis (Kurt  3) are called leptokurtic; those with Kurt  3 are named platykurtic.

For error distributions which are non symmetric (Skew ), quantifying the accuracy by a single dispersion-related statistic is not reliable, and one should provide probability intervals or accept to lose information on the sign and use a statistic based on absolute errors, such as Q₉₅ presented below.

For platykurtic and leptokurtic error distribution, the normal probabilistic interpretation of dispersion statistics such as the RMSD is lost ⁴ , and one should rely on dispersion measures providing explicit probability coverage, such as u₉₅ recommended in the thermochemistry literature [10]. u₉₅ is an enlarged uncertainty which enables to define a symmetric 95 % probability interval around the estimated value. If the distribution is also skewed, one cannot rely on a symmetric interval based on u₉₅. A major problem with leptokurtic error distributions is that the excess of data in the tails is related to a risk of large prediction errors, that cannot be appreciated from the usual MAE or RMSD statistics. A solution is again to use Q₉₅ (section 2.2).

For ranking studies, problems might occur when comparing distribution with different shapes, which might lead to conflict between different ranking statistics [2]. A striking example is provided below (section 3).

Pernot and Savin [2] have shown that the non-normality of error distributions prevents the probabilistic interpretation of usual benchmark statistics. However, direct probabilistic information can be extracted from the empirical cumulative distribution function (ECDF) of absolute errors

where η is a threshold value to be chosen, and is the indicator function, with value 1 if x is true and 0 otherwise.

Two statistics based on the ECDF have been proposed as fulfilling the needs of end-users to choose a method suited to their purpose:

• C(η) is the probability that the absolute errors lie below a chosen threshold η, such as the 'chemical accuracy' [2]. A user can define his own needs and accordingly pick methods with a high value of C(η).
• Q₉₅, the 95^th percentile of the absolute error distribution, gives the amplitude of errors that there is a 5% probability to exceed [2]. The choice of this specific percentile was based on several considerations. When the error distribution is normal, Q₉₅ is identical to the enlarged uncertainty u₉₅ recommended for thermochemical calculations and data [10]. Besides, Q₉₅ is much less sensitive to outliers than the maximal error, or even Q₉₉ for small datasets. Thakkar et al [11] proposed a similar statistic (P₉₀) based on the 90^th percentile.

As mentioned previously, the use of such probabilistic scores for non-normal distributions relieves the ambiguity attached to the MAE and RMSD. For instance, in all the cases of error sets that Pernot and Savin [2, 7] have studied previously, the probability for an absolute error to be larger than the MAE lied between 0.2 and 0.45 (it should be about 0.42 for a zero-centered normal distribution). The small values are typically associated with skewed and leptokurtic distributions. Several cases have also been observed where two methods have similar MAE values and very different probabilities of large errors, due to their different error distributions.

This non-ambiguity is essential if one wants to assess a risk of large prediction errors. Note that the transition from descriptive statistics to predictive statistics for risk assesment requires further constraints on the reference dataset and error distributions. Ideally, the reference dataset has to be representative enough to cover future prediction cases and there should be no notable trend in the errors with respect to the calculated values by the method of interest (in the statistical prediction framework, the predictor variable is the calculated value of a property, from which one wants to infer a credible interval for its true value).

Several statistical trend correction methods have been proposed in the computational chemistry literature, from the simple scaling of the calculated values [12–14], or linear corrections [1, 4, 15–17], to more complex, ML-based corrections, such as Δ-ML [6, 18, 19] or Gaussian Processes [20].

There is no unique method to identify outliers for a non-normal distribution. One might, for instance, use visual tools, such as QQ-plots [4], or automatic selection tools, such as selecting points for which the absolute error is larger than the 95th percentile (Q₉₅), or another percentile corresponding to a predefined error threshold.

In cases where the errors distribution E(x) seems heterogeneous, one can attempt to analyze it as a mixture of normal distributions. One might for instance use a bi-normal fit

where is a normal distribution of x with mean value µ and standard deviation σ. Outliers are then defined as the points which lie outside of a interval, where µ₁ and σ₁ are the parameters of the most concentrated component. The enlargement factor n should be chosen large enough to exclude as much points from this component as possible, but small enough to capture enough points of the wider component to enable tests of chemical or physical hypotheses.

All calculations have been made in the R language [21], using several packages, notably for bootstrap(boot [22]), skewness and kurtosis (moments [23]) and mixture analysis (mixtools [24, 25]). Bootstrap estimates are based on 1000 draws. The R implementation of the Shapiro-Wilk test (function shapiro.test) has a limit of dataset size at 5000 [8]. For the large datasets used in the application, 5000 points are randomly chosen to evaluate W.

The code and datasets to reproduce the tables and figures of this article are available at Zenodo (https://doi.org/10.528 1/zenodo.3733 272). The datasets can also be analyzed with the ErrView code (https://doi.org/10.528 1/zenodo.3628 489), or its web interface (http://upsa.shinyapps.io/ErrView).

Summary statistics and their uncertainty have been estimated by bootstrap [27]. They are reported in table 1. In this part, one concentrates on the unmodified methods, HF, MP2, CM-L2, and SLATM-L2. Considering MAE alone, one might conclude that SLATM-L2 produces slightly smaller errors than MP2 and a significant improvement over all other methods. MSE tells us that there is no important bias when compared to the MAE or RMSD (all distributions are nearly centered on zero). It is striking that MP2 has a smaller RMSD than SLATM-L2. As the RMSD is often taken as an estimator of prediction uncertainty [1], this contradicts the MAE ranking.

This is also a strong clue that not all the error distributions are normal. In the case of normal distributions, and considering that MSE≃ 0, the RMSD should confirm the MAE ranking ( for a zero-centered normal distribution). Indeed, the Shapiro-Wilk normality test rejects the normality of all the error sets, with a critical value W_c  = 0.999 3 for N = 5000 and a type I error rate of 0.05. However the W statistic shows that MP2 is much closer to normal (W = 0.996 1) than SLATM-L2 (W = 0.71). This is confirmed by the kurtosis, with a very large value for SLATM-L2 (Kurt = 29.0) whereas for MP2 the value (3.33) is much closer to the value for a normal distribution (3.0). HF and CM-L2 are also leptokurtic (Kurt > 3.0), but much less than SLATM-L2.

The large kurtosis value indicates that the SLATM-L2 error distribution presents tails much heavier than for a normal distribution. This excess in large errors can be appreciated with Q₉₅, which informs us that 5 % of the absolute errors are above 3.3 kcal mol⁻¹ for MP2 and 4.7 kcal mol⁻¹ for SLATM-L2. Despite a smaller MAE, one has thus a non-negligible risk to get larger errors by SLATM-L2 than by MP2. Note that in parallel with large errors, there is an excess in small errors: we estimated that the probability to have absolute errors smaller than the RMSD, C(2.44), is 0.9 for SLATM-L2, when it should be 0.67 for a normal distribution.

In order to substantiate the interpretation of the error statistics, it is important (1) to check that the error distributions are not dominated by systematic errors (trends), and (2) to assess their normality or lack thereof. This is essential if any quantification of the predictive ability of the methods is of interest.

3.1.1. Trends

Plots of the error sets (figure 1(a-d)) show that the ML errors present weak (linear) trends, whereas HF and MP2 are more affected. For comparison of their performances, the data sets have been linearly corrected to remove the trends, according to Pernot et al[1]. Summary statistics for corrected methods (with a 'lc-' prefix) are reported in table 1, showing a small impact on all statistics. The normality Shapiro-Wilk statistic is improved for HF and MP2, and the error distribution for lc-MP2 is practically normal. After linear correction, and considering the small uncertainty on the linear correction parameters due to the large size of the dataset, the prediction uncertainty can be estimated by the RMSD [1]. This enables to assess the better performance of MP2 in terms of prediction uncertainty.

Zoom In Zoom Out Reset image size

3.1.2. Normality

3.1.3. ECDFs

Considering figure 2, one see that the absolute errors ECDFs for SLATM-L2 and MP2 intersect around 3 kcal mol⁻¹. Below this value, SLATM-L2 presents smaller absolute errors, with the opposite above. It is clearly visible here why SLATM-L2 has a larger Q₉₅ than MP2. Indeed, one sees also on figures 1 (b), (d) that SLATM-L2 has a non-negligible set of errors much larger than the largest MP2 errors. A similar effect occurs even between HF and SLATM-L2, the latter presenting larger absolute errors above approximately 8 kcal mol⁻¹.

Zoom In Zoom Out Reset image size

Computational time aside, the choice between MP2 and SLATM-L2 as a predictor for CCSD(T) energies depends on the acceptance by the user of a small percentage of large errors. SLATM-L2 would largely benefit from a reduction of such large errors. In the following, we focus on the SLATM-L2 error set to identify underlying features linked to large prediction errors.

In order to get a chemical insight on the error trends in the present dataset, the absolute errors have been analyzed as functions of several compositional data, such as the number of individual species in the composition (H, C, N, O, S, Cl) and the Double Bond Equivalent (DBE), which estimates the degree of unsaturation (number of double bonds and rings) [28]

where is the number of atoms in the molecule. Note that in this formula divalent atoms do not contribute to the DBE.

Sensitivity of the errors to the various variables has been estimated by using rank correlation coefficients. The largest sensitivity has been found with respect to the number of hydrogen atoms in the molecule () and the DBE. Note that both indicators are not independent (DBE is anti-correlated with ), and that is strongly correlated with the number of atoms in the molecule.

The distribution of the absolute errors sorted by DBE class is shown in figure 3(d). Distributions of the absolute errors according to and DBE are presented in (figure 3(a),(b)), where each datum is represented by a circle of diameter proportional to its absolute error. The scale is identical for all the series of plots. The imprint of the points cloud reflects the correlation between and DBE.

Zoom In Zoom Out Reset image size

For SLATM-L2, some systems with DBE  =  3-5 present very large errors, and the maximal absolute error increases almost linearly from DBE  =  0 to 4. All compositions with DBE  =  6 are correctly predicted, although, except for , none of those compositions are found in the learning set. For molecules with DBE  =  0, all heavy atoms are sp³-hybridized, being the most local and thus transferable and easy to learn by QML models. For molecules with DBE  =  6, bond environments are non-local, but could be easily covered by the training set. Therefore, low prediction errors are expected in these cases.

For MP2, the points cloud in figure 3(b) presents no strong feature, but the errors seem smaller for DBE  =  0, and for each DBE class, there seems to be a maximum around the mid-range of admissible values.

The non-normality of the errors distribution for SLATM-L2 is clearly visible on the QQ-plot in figure 4(a). There is a near linear section at the center if the distribution, seen by comparison with the interquartile line, and also linear sections in the tails of the distribution. The core of the distribution is therefore dominated by a normal component, as well as the tails, albeit with different widths (and centers). Indeed, a bi-normal fit (equation 3) provides a good representation of the distribution. Note that to get a perfect fit, as much as four components are needed, but the bi-normal representation is sufficient for the present purpose of outliers analysis.

Zoom In Zoom Out Reset image size

The graphical results and optimal parameters are reported in figure 4(b). This decomposition provides a strong component (weight 0.83) with a sub-kcal/mol standard deviation (σ₁ = 0.87 kcal/mol), and a weaker component (weight: 0.17) with a much larger dispersion (σ₂ = 5.5 kcal/mol). Both components are nearly zero-centered.

In order to identify the systems involved in the widest contribution, one selects outliers with respect to the most concentrated component, i.e. points lying at more than 5*σ₁ from the center of the concentrated component (µ₁).

This filter selects 350 molecules (about 6 % of the test set population). The chemical formulae of these systems are analyzed in terms of and DBE (figure 3(c)). Systems with DBE  =  0 or 6 are practically absent from the list of outliers (at the exception of a few DBE  =  0 ones). The outliers with the largest absolute errors are molecules with some unsaturation (DBE ≥ 3).

Figure 5 reports the DBE distributions in the leaning set figure 5(a) and test set figure 5(b). In both cases, the distributions are peaked around DBE = 2. The distributions are similar because the learning set is a random selection in the full dataset. The ratio of the learning set DBE distribution over the full dataset DBE distribution is practically constant, with slight increasing trend that might be a random effect (figure 5(c)). In contrast, figure 5(d) shows that the DBE classes around DBE = 4 contain a larger proportion of outliers than the lower DBE classes, up to 12 % for DBE  =  4. Outliers are not uniformly distributed amongst the DBE classes.

Zoom In Zoom Out Reset image size

Compositions with a maximum absolute error larger than 10 kcal mol⁻¹ are listed in table 2. The composition with the largest absolute error is (27.9 kcal mol⁻¹ for the molecule SMILES:C#CC( = N#N)C = O [29]). This molecule, with DBE = 5, has 3 isomers/configurations in the learning set and 12 in the test set, 4 of which are outliers. The second one is (27.4 kcal mol⁻¹, SMILES:c1cnn2c1CC2), with 17 configurations in the learning set, 60 in the test set, 10 of which are outliers.

Table 2. Composition of the 5σ outliers with maximum absolute error above 10 kcal mol⁻¹. For each composition, one defines medAE: median Absolute error in outlier isomers; maxAE: max Absolute error in outlier isomers; nOutl: number of isomers in the outlier set; nTest: number of isomers in the test set; nLearn: number of isomers in the learning set.

Composition	DBE	medAE	maxAE	nOutl/nTest	nLearn
	5	10.28	27.92	4/12	3
	4	6.28	27.40	10/60	17
	4	8.13	25.79	14/64	13
	3	9.14	24.46	9/152	25
	4	10.93	21.66	4/25	3
	3	6.22	21.36	37/223	39
	3	6.07	21.29	11/115	22
	3	6.02	19.79	10/203	33
	2	12.71	17.90	6/207	29
	2	9.86	17.62	5/6	1
	4	13.13	16.76	2/6	2
	4	16.60	16.60	1/3	0
	3	8.44	15.36	20/109	18
	2	5.32	15.09	4/242	51
	4	9.48	15.08	6/51	6
	3	6.71	14.91	14/40	8
	1	13.63	14.58	5/5	0
	2	12.10	14.18	12/12	0
	4	9.98	13.71	2/43	7
	4	8.67	13.71	5/66	6
	1	6.88	12.07	9/15	1
	3	8.57	12.04	2/2	1
	4	11.99	11.99	1/5	5
	2	5.72	11.87	6/229	26
	2	7.39	11.31	5/5	2
	0	6.49	11.30	6/8	0
	4	6.08	11.20	7/83	10
	1	7.77	10.96	9/14	4
	2	5.21	10.83	4/194	28
	2	5.66	10.20	20/245	50
	1	5.56	10.14	9/151	25

One can see in this list that most outliers compositions present some isomers in the learning set. Unsurprisingly, it contains also a few compositions absent from the learning set, e.g., (SMILES:c1c2c(cs1)OC2), but as seen above for the DBE = 6 compounds, this is not systematically associated with poor predictions.

The impact on the error statistics of the removal of the 350 outliers from the test set is reported in table 1, as method SLATM-L2(-o). The improvement is noticeable on all aspects: better normality of the distribution, reduced bias (MSE ≃ 0), and a notable reduction, as expected, of MAE, RMSD and Q₉₅.

Having identified 350 systems for which SLATM-L2 provides exceedingly large errors, we want to check if their transfer to the learning set has a positive impact on prediction performances. To validate the impact of this learning set augmentation, a list of 350 systems has been transferred from the test set to the learning set according to two selections: (1) the outliers identified above (SLATM-L2(+o)), and (2) a random choice (SLATM-L2(+r)), both resulting in a learning set of size 1350. The second scenario is a control to assess the size effect of the learning set. Note that the present random selection contains 22 of the 350 outliers (≃ 6 %), in proportion with the random selection of 350 points among 6211.

The SLATM-L2 model has been retrained with both learning sets, and the error distribution analysis has been performed for these two new datasets. The results are reported in table 1 and figure 1(e), (f), 2(b), figure 6.

Zoom In Zoom Out Reset image size

A first aspect to consider is how the inclusion of the outliers in the learning set (SLATM-L2(+o)) affects the prediction errors with respect to the pruned error set (SLATM-L2(-o)) presented in the previous Section. Globally, the improvement is minor, but consistent. The MSE and RMSD are identical, the MAE is barely decreased, and the Q₉₅ is slightly improved from 2.70 to 2.49 kcal mol⁻¹. The learning on the outliers has been transferred to some of the predictions.

Then, we compare the gain in prediction quality due to the augmentation of the learning set by outliers versus random systems. The MAE values of SLATM-L2(+o): 0.83 and SLATM-L2(+r): 0.95 kcal mol⁻¹ indicate a statistically significant, but rather weak effect of the outliers selection, whereas the sheer augmentation of the learning set size improves notably the MAE of the original method by 0.3 kcal mol⁻¹. The MSE shows that the slight bias of the original SLATM-L2 has been essentially corrected in both cases. Impact on the RMSD is larger than on the MAE (1.20 for '+o' vs. 1.89 kcal mol⁻¹ for '+r'), with a notable reduction of 1.24 kcal mol⁻¹ from the initial RMSD. Similarly, the decrease in Q₉₅ is important for SLATM-L2(+o), from 4.7 to 2.49 kcal mol⁻¹, and much less for SLATM-L2(+r), to 3.15 kcal mol⁻¹, which reaches the level of MP2 (3.34 kcal mol⁻¹). One sees also that the normality of the distribution has been improved for SLATM-L2(+o) (Kurt from 29 to 8.1; W from 0.71 to 0.935), which is not the case for the random selection (Kurt = 46; W  =  0.66).

Without looking at the plots of the distributions, this set of statistics points to a notable improvement of the error distribution on all criteria for the outliers-augmented learning set. By contrast, the improvement by random selection of the 350 systems, although significant, does not improve the quality of the error distribution, which stays strongly leptokurtic. This can be appreciated on the QQ-plots for both error distributions in figure 6(a),(b).

The ECDFs of the error distributions illustrate clearly the analysis of the summary statistics (figure 2(b)). The augmented SLATM-L2 versions have more concentrated absolute errors than SLATM-L2. In the case of SLATM-L2(+o) the large errors have essentially vanished and the overall performance is better than SLATM-L2 and even than MP2. Albeit SLATM-L2(+r) is also more concentrated, it still presents a set of systems with large errors, exceeding those of MP2. The comparison of the /DBE plots in figure 6(c),(d) and depicts precisely this point.

The QM7b dataset has a nice feature of diversity, albeit it is small in total size. A natural yet undesired consequence is that the dataset is highly inhomogeneous in chemical space, the conclusions drawn above from the results on QM7b may therefore not be extensible to other datasets, such as QM9 [30, 31] (less diverse in composition, yet larger in total set size). Considering these properties of QM9, one might expect the prediction error distribution for SLATM-L2 to be less leptokurtic than for QM7b, and a less problematic use of MAE for performance assesment. This has however to be checked, but at the moment, one has no access to hierarchical properties for the full QM9, as for QM7b.

In this paper we assessed the prediction of energy, which is extensive in nature. There exists a set of intensive properties as well and which are typically difficult to learn. The generalizability of our observations to these properties should also be confirmed in future studies.

A major feature of the SLATM-L2 error distribution is the presence of a strong peak of small errors and wide tails of large errors. This leptokurtic distribution has been successfully modeled as a mixture of two normal distributions, both practically centered on zero, with standard deviations 0.87 and 5.5 kcal mol⁻¹ and a mixture ratio 83:17. Based on this analysis, we explore the chemical identity of the systems with large errors.

It is not possible to unscramble both populations in the area where the most concentrated population is dominant, but if one goes far enough in the tails (beyond 5 times 0.87 kcal mol⁻¹) one finds a large majority of systems belonging to the population with large errors. Such 350 systems were identified and tagged as outliers. A chemical analysis of these systems based on their DBE revealed that compositions with high levels of unsaturation (DBE= 3 − 5) are overrepresented with respect to their abundance in the test set. It seems therefore that the large prediction errors are not randomly distributed over the QM7b molecules, which might provide insights for improved molecular descriptors or a better design of the learning set.

In fact, if the learning set is augmented with these outliers (SLATM-L2(+o)), one gets better prediction performances than if one augments the learning set with 350 randomly chosen systems (SLATM-L2(+r)). This suggests an iterative method for the design of learning sets that might deserve further consideration.

For ML methods to be used trustfully as replacements of quantum chemistry methods, one needs to have a reliable measure of their accuracy. Ideally, this should be quantified by a single number characterizing the predictive ability of a method. This requires a finite variance, homogeneous, error distribution (with a dispersion that does not vary strongly within the calculated property range) and a zero-centered symmetric error distribution. Moreover, we have seen that leptokurtic error distributions with a strong core of small errors and tails of large errors, as described above by a bi-normal distribution might not be properly summarized by a single dispersion statistic. The MAE is not a good choice in such cases, and Q₉₅ presents a more reliable alternative, which quantifies the risk of large prediction errors, even in the case of non-zero-centered symmetric error distributions.

We have shown that inclusion of outliers in the learning set of the SLATM-L2 method could result in a notable improvement of the shape of its prediction errors distribution. One might also presume that the use of much larger learning sets should improve the structure of the prediction error distribution.

Whatever the pertinence of this question, it cannot be answered by considering the MAE, contrary to a still too common practice in the benchmarking literature to use it as an 'accuracy' measure. The example of the SLATM-L2(+o) method shows that a sub-chemical-accuracy MAE (0.86 kcal mol⁻¹) can hide a non-negligible risk (about 30 %) to exceed this limit.

A direct answer to the question is provided by the ECDF of the errors and C(η) with η = 1 kcal mol⁻¹. The values reported in table 1 tell us that none of the studied methods achieves the 1 kcal mol⁻¹ threshold with a high probability. The best performers would be SLATM-L2 methods with augmented learning sets, with C(1) = 0.743 for the (+r) version. Even in this case, one has no firm guarantee to reach chemical accuracy for any new prediction, the risk being higher than 25 % to overpass the threshold. Note that this is a notable improvement over MP2, for which the risk is about 55 %. Besides, this qualifies chemical accuracy with respect to CCSD(T), not with respect to experimental values.

In section 3.1, we outlined that the basic summary statistics (MAE, RMSD) provided contradictory information about the SLATM-L2 prediction performances, notably when compared to MP2. SLATM-L2 has a smaller slightly smaller MAE than MP2 (1.26 vs. 1.31 kcal mol⁻¹) but a much larger RMSD (2.44 vs. 1.67). This is a consequence of the peculiar shape of the error distribution for SLATM-L2. In line with the RMSD values, MP2 achieves a notably smaller Q₉₅ value (3.34 vs. 4.7 kcal mol⁻¹), meaning that SLATM-L2 is susceptible of larger errors than MP2.

Given this set of information and computation time aside, one would most certainly pick MP2 over SLATM-L2 as a reasonable predictor of CCSD(T) energies, with a prediction uncertainty of approximately 1.7 kcal mol⁻¹ for E^* (RMSD in table 1), and with added confidence from the near-normality of the distribution, reducing the risk of rogue predictions that might be expected from SLATM-L2. Finally, the modified SLATM-L2(+o) method, for which the learning set was augmented with SLATM-L2 outliers, clearly outperforms both SLATM-L2 and MP2 for all statistics.

This illustrates once more that the MAE is not sensitive to major differences in the shape of the errors distributions, as was shown by Pernot and Savin [2]: even for normal distributions, many combinations of MSE and RMSD can produce the same MAE. The case is even worse for non-normal distributions, where the MAE loses its ability to quantify prediction uncertainty. It is clear that the MAE should not be used as a unique scoring statistic, and that it should at least be complemented (or replaced) by probabilistic indicators such as Q₉₅ .