Pascal Pernot
email: pascal.pernot@universite-paris-saclay.fr

Andreas Savin
email: andreas.savin@lct.jussieu.fr

Using the Gini coefficient to characterize the shape of computational chemistry error distributions

The distribution of errors is a central object in the assesment and benchmarking of computational chemistry methods. The popular and often blind use of the mean unsigned error as a benchmarking statistic leads to ignore distributions features that impact the reliability of the tested methods. We explore how the Gini coefficient offers a global representation of the errors distribution, but, except for extreme values, does not enable an unambiguous diagnostic. We propose to relieve the ambiguity by applying the Gini coefficient to mode-centered error distributions. This version can usefully complement benchmarking statistics and alert on error sets with potentially problematic shapes.

Introduction

The reliability of a computational chemistry method is conditioned by the distribution of its prediction errors. Distributions with heavy tails elicit a risk of large prediction errors. As a benchmarking statistic, the popular mean unsigned error (MUE) bears no information on such a risk. [START_REF] Pernot | Prediction uncertainty of density functional approximations for properties of crystals with cubic symmetry[END_REF][START_REF] Pernot | Probabilistic performance estimators for computational chemistry methods: the empirical cumulative distribution function of absolute errors[END_REF][START_REF] Pernot | Probabilistic performance estimators for computational chemistry methods: Systematic improvement probability and ranking probability matrix[END_REF][START_REF] Pernot | Probabilistic performance estimators for computational chemistry methods: Systematic improvement probability and ranking probability matrix[END_REF] We have recently reported a case where two unbiased error distributions with identical values of the MUE present widely different risks of large errors because of heavy tails in one of them. [START_REF] Pernot | Impact of non-normal error distributions on the benchmarking and ranking of Quantum Machine Learning models[END_REF][START_REF] Pernot | Probabilistic performance estimators for computational chemistry methods: Systematic improvement probability and ranking probability matrix[END_REF] It would therefore be very useful to complement the MUE with a statistic indicating or quantifying the risk of large errors.

We recently proposed alternative statistics such as Q 95 , 2 P η , 2 and systematic improvement probability (SIP). [START_REF] Pernot | Probabilistic performance estimators for computational chemistry methods: Systematic improvement probability and ranking probability matrix[END_REF] In terms of risk, these statistics offer the following interpretations:

• There is a 5 % risk for absolute errors to exceed Q 95 .

• There is a probability P η that absolute errors are larger than a chosen threshold η. P η provides a direct quantification of the risk of large errors, but η has to be defined a priori and might be user-dependent, which complicates its reporting in benchmarking studies.

• For two methods M 1 and M 2 , the SIP provides the system-wise probability that the absolute errors of M 1 are smaller than the absolute errors of M 2 , informing on the risk incurred by switching between two methods. Interestingly, the SIP analysis provides a decomposition of the MUE difference between two methods in terms of gain and loss probabilities. [START_REF] Pernot | Probabilistic performance estimators for computational chemistry methods: Systematic improvement probability and ranking probability matrix[END_REF] The Q 95 and P η partly answer the question, but they are point estimates on the cumulated density function of the absolute errors, and a statistic accounting for the whole distribution might be of interest. Besides, it is well established in econometrics that measures of dispersions such as the variance perform poorly at risk estimation and that higher moments of the distributions have to be considered. [START_REF] Bonato | Robust estimation of skewness and kurtosis in distributions with infinite higher moments[END_REF] This would lead us to such measures as skewness and kurtosis, but none of these alone would be able to cover all the scenarii. The risk of large errors through heavy tails of the errors distribution might be associated with large skewness or large kurtosis or a combination of them.

The Lorenz curve 7 is widely used in econometrics to represent the distribution of wealth in human populations. Its summary statistics, notably the Gini coefficient (noted G), [START_REF] Gini | Variabilità e mutabilità[END_REF][START_REF] Damgaard | Describing inequality in plant size or fecundity[END_REF][START_REF] Eliazar | Measuring statistical heterogeneity: The Pietra index[END_REF] are used to evaluate the degree of inequality within these populations. The Gini coefficient is also used, for instance, in ecology, to estimate the inequality of properties within plant populations, [START_REF] Bendel | Comparison of skewness coefficient, coefficient of variation, and Gini coefficient as inequality measures within populations[END_REF][START_REF] Damgaard | Describing inequality in plant size or fecundity[END_REF] in astronomy, to characterize the morphology of galaxies, [START_REF] Florian | The Gini coefficient as a morphological measurement of strongly lensed galaxies in the image plane[END_REF] or in information theory, to characterize the sparsity of datasets. [START_REF] Hurley | Comparing measures of sparsity[END_REF] The Lorenz curve has many mathematical representations, the most interesting one, for us, being its formulation as an integral of the quantile function, a direct link with our study of probabilistic metrics. [START_REF] Pernot | Probabilistic performance estimators for computational chemistry methods: the empirical cumulative distribution function of absolute errors[END_REF][START_REF] Pernot | Probabilistic performance estimators for computational chemistry methods: Systematic improvement probability and ranking probability matrix[END_REF][START_REF] Pernot | Probabilistic performance estimators for computational chemistry methods: Systematic improvement probability and ranking probability matrix[END_REF] More precisely, we explore here the interest of the Gini coefficient as a complement to the MUE in benchmarking studies.

We introduce the statistical tools and their implementation in Section 2, and apply them to a series of datasets to illustrate the interest and limitations of the Gini coefficient in Section 4. An adaptation of the Gini coefficient is proposed to relieve its main drawbacks when applied to error datasets.

Statistical methods

Definitions

Let us consider a distribution of errors e with probability density function (PDF) f (e). The absolute errors ε = |e| have a folded PDF f F (ε). To avoid ambiguity, statistics based on absolute errors are indexed by F .

CDF and quantile function

The cumulative distribution function (CDF) of the absolute errors is noted

C F (z) = z 0 f F (ε) dε (1)
from which the quantile function is the inverse

q F (p) = C -1 F (p) (2)

Mean unsigned error

The mean unsigned error (MUE), is defined as

µ F = ∞ 0 εf F (ε) dε (3)
Using the change of variable ε = C -1 F (p), p = C F (ε) and dp = f F (ε) dε, the MUE also can be shown to be the integral of the quantile function

µ F = 1 0 C -1 F (p) dp (4) = 1 0 q F (p) dp (5)

The Lorenz curve

The Lorenz curve gives the percentage of cumulated absolute errors due to the 100 × p % smallest values, or equivalently, the portion of the MUE due to the 100 × p % smallest absolute errors:

L F (p) = 1 µ F p 0 q F (t) dt (6)
As shown in Fig. 1(a), it is the ratio between the slanted shaded area and the total slanted area. Its value for p is reported on the corresponding Lorenz curve graph (Fig. 1(b)).

The Lorenz curve provides a scale-invariant representation of the CDF C F (z), [START_REF] Kleiber | The Lorenz curve in economics and econometrics[END_REF] with the following properties: L F (p) is concave, increasing with p, such as 0 ≤ L F (p) ≤ p ≤ 1, L F (0) = 0 and L F (1) = 1. L F (p) lies one the identity line (L F (p) = p) when all the errors are equal, i.e. f F (ε) = δ(ε -c). Note that this case corresponds to a discontinuous CDF, with a jump at ε = c. The deviation of an error distribution from this extreme case can be measured by the Gini coefficient.

The Gini coefficient

It is related to the area between L F (p) and the identity line (Fig. 1(b))

G F = 2 1 0 {p -L(p)} dp (7)
where the factor two is used to scale G F between 0 and 1. The smaller G F , the closer the Lorenz curve to the identity line. The Gini coefficient, usually noted G, is generally used for distributions with positive support. Our notation G F is a reminder that we are working here with distributions of absolute errors f F (ε).

For sets of absolute errors with a normal distribution [START_REF] Bendel | Comparison of skewness coefficient, coefficient of variation, and Gini coefficient as inequality measures within populations[END_REF] where σ F the standard deviation of the absolute errors

N (ε; µ F , σ F ), G F is proportional to the coefficient of variation c v = σ F /µ F ,
G F ∼ c v / √ π ( 8
)
Note that this relationship does hold only when all errors are of the same sign (µ F σ F ), therefore with small c v values.

Two typical values of G F will be useful in the following:

• for any zero-centered normal distribution of errors N (0, σ), the distribution of absolute errors is the half-normal distribution, with value G F N = √ 2 -1 0.41; 15

• for any zero-centered uniform distribution, U (-a, a), or any uniform distribution with a bound at zero, U (-a, 0) or U (0, a), folding produces a uniform distribution with the minimal bound at zero, with value G F U = 1/3. 15

Skewness and kurtosis

Skewness measures the asymmetry of a distribution, while kurtosis is used as a measure either of its "peakedness" or "tailedness" [START_REF] Ruppert | What is kurtosis?: An influence function approach[END_REF] The moments-based formulae for skewness and kurtosis are not robust to outliers, and more robust quantile-based formulae have been proposed by several authors. [START_REF] Groeneveld | Measuring skewness and kurtosis[END_REF][START_REF] Ruppert | What is kurtosis?: An influence function approach[END_REF][START_REF] Bonato | Robust estimation of skewness and kurtosis in distributions with infinite higher moments[END_REF][START_REF] Suaray | On the asymptotic distribution of an alternative measure of kurtosis[END_REF] For the skewness, we use a measure using the difference between the mean and median

β GM = µ -q(0.5) < |e -q(0.5)| > (9)
where the brackets indicate the mean value, q(0.5) is the median of signed error, e, and the GM subscript refers to the authors of this definition, Goeneveld and Meeden. 17 β GM takes its values between -1 and 1, and is 0 for symmetric distributions.

For kurtosis, an estimate based on quantiles is used 6 (originating from a similar form proposed by Crown and Siddiqui, [START_REF] Crow | Robust estimation of location[END_REF] hence the CS subscript) κ CS = q(0.975) -q(0.025) q(0.75) -q(0.25) -2.91 (10) where q(.) is the quantile function for signed errors. The correction factor for the normal distribution (2.91) makes that κ CS measures an excess kurtosis. Datasets with κ CS > 0 have heavier tails than a normal distribution, and the opposite for negative values. Specific notations will be introduced below when these statistics are applied to sets of absolute errors.

Estimation

We consider in this section the application of the previous statistics to finite error samples, and the corresponding formulae. Let us consider a set of errors

(E = {e i } N i=1 ), derived from a set of N calculated values (C = {c i } N i=1 ) and reference data (R = {r i } N i=1 ), by e i = r i -c i . The absolute errors are noted Ξ = {ε i = |e i |} N i=1 .

MSE, MUE and mode. The mean signed error (MSE) is estimated as

µ = 1 N N i=1 e i ,

and the mean unsigned error (MUE) as µ

F = 1 N N i=1 ε i .
As one is not dealing with necessarily symmetric distributions, the mode is an interesting location statistic, notably in correspondence to the tails of a distribution. The mode locates the part of the population with the highest density, which is expected to bring a large contribution to µ F , and therefore influence the Lorenz curve and Gini coefficient. As one cannot assume the unimodality of the underlying distributions, one will consider the main mode. A non-parametric robust method, Bickel's half-range mode (HRM) estimator, [START_REF] Bickel | Robust estimators of the mode and skewness of continuous data[END_REF][START_REF] Hedges | Comparison of mode estimation methods and application in molecular clock analysis[END_REF] is used to estimate the location of the error samples main mode. This methods proceeds by iterative bipartition of modal intervals (intervals with highest density).

L F (p) and G F . Let us introduce the cumulated sum of the n ≤ N smallest absolute errors

S n = n i=1 ε [i] (11)
where ε [i] is the i th order statistic (i.e., the value with rank i ) of the sample of absolute errors. For consistency, one sets S 0 = 0. The Lorenz curve is estimated as

L F (p) = S p×N S N ( 12
)
where 0 ≤ p ≤ 1. Note that, due to the use of finite samples, p takes its values in {i/N } N i=0 . Using a fast sorting of the sample of absolute errors, an efficient estimation of G F uses the formula [START_REF] Damgaard | Describing inequality in plant size or fecundity[END_REF][START_REF] Glasser | Variance formulas for the mean difference and coefficient of concentration[END_REF][START_REF] Dixon | Bootstrapping the Gini coefficient of inequality[END_REF][START_REF] Zeileis | ineq: Measuring Inequality, Concentration, and Poverty, 2014. R package version 0[END_REF]

G F = N i=1 (2i -N -1)ε [i] N N i=1 ε [i] (13)
A slower, but equivalent expression in terms of mean values is 10

G F = 1 M U E < max{0, ε 1 -ε 2 } > ( 14
)
where ε 1 and ε 2 are two elements of Ξ and the mean is taken on all pairs.

β GM and κ CS . For skewness and kurtosis, Eq. 9 and Eq. 10 are applied directly, with the robust method to estimate quantiles due to Harrell and Davis. [START_REF] Harrell | A new distribution-free quantile estimator[END_REF][START_REF] Wilcox | Comparing two dependent groups via quantiles[END_REF][START_REF] Pernot | Probabilistic performance estimators for computational chemistry methods: Systematic improvement probability and ranking probability matrix[END_REF] Uncertainty. Uncertainty on any statistic X, noted u(X), is estimated by bootstrapping 26 with 1000 samples. Note that there is a known risk of underestimation of G F for small datasets (N < 100). 15

Implementation

All calculations have been made in the R language, 27 using several packages, notably for the Gini coefficient (package ineq [START_REF] Zeileis | ineq: Measuring Inequality, Concentration, and Poverty, 2014. R package version 0[END_REF] ), the HRM estimator (package genefilter [START_REF] Gentleman | genefilter: methods for filtering genes from high-throughput experiments[END_REF] ) and bootstrapping (package boot [START_REF] Canty | boot: Bootstrap R (S-Plus) Functions[END_REF] ). The Gini coefficient, Lorenz curves, G M CF , β GM , κ CS and mode estimator have been implemented in the freely available R [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF] package ErrViewLib (v1.3, https://doi.org/10.5281/ zenodo.3628475). The datasets can be analyzed with the ErrView graphical interface (source: https://doi.org/10.5281/zenodo.3628489; web interface: http://upsa.shinyapps.io/ErrView).

Datasets

Model datasets

Before applying the Gini coefficient to literature datasets, one explores its properties on error sets generated from distributions with controlled properties: uniform, normal, Student's-t, lognormal [START_REF] Evans | Statistical Distributions[END_REF] and g-and-h. [START_REF] Hoaglin | Exploring data tables, trends, and shapes, chapter Summarizing shape numerically: The g-and-h distributions[END_REF][START_REF] Wilcox | Comparing two dependent groups via quantiles[END_REF][START_REF] Pernot | Probabilistic performance estimators for computational chemistry methods: Systematic improvement probability and ranking probability matrix[END_REF] It is important to note that we explore only distributions with a single dominant, more or less structured peak, such as the ones encountered in most computational chemistry error datasets. In the list of analytical distributions above, the uniform is an exception because of its undefined mode. We use it as an extreme case of single peaked, continuous distribution, with negative excess kurtosis. Besides, it is easy to design multi-peaked distributions for which our conclusions on the Gini coefficient would not be valid. In fact, none of the usual summary statistics (MUE, MSE, skewness, kurtosis...) would describe properly such distributions.

Case

Literature datasets

The statistical tools described above are applied to datasets gathered in the computational chemistry literature. These are summarized in Table 1, and strongly overlap with those studied in more details in a previous article, 4 from which we removed small datasets (N < 100). The statistics, empirical cumulative density functions and Lorenz curves corresponding to these datasets are provided as Supplementary Information.

Applications

In its usual application fields, the Gini coefficient is applied to distributions with positive support. Our application to computational chemistry error sets involves the intermediate folding operation, which is not reversible. In a first part, we show how this limits the information on the errors distribution that can be inferred from the Gini coefficient. To relieve this difficulty, we propose a mode-centering operation before folding, which better preserves some of the tail properties of the original distributions. The Gini coefficient is then compared to other tails statistics, notably at the level of statistical uncertainty.

Gini coefficient vs. bias

The link between the Gini coefficient and the coefficient of variation (Eq. 8) tells us that, for a normal distribution of given standard deviation, an decreasing bias should result in increasing values of G F . At some point, this relation is broken by the folding operation: as noted earlier, for a centered normal distribution, one has G F = 0.41, for which Eq. 8 does not hold. The dependence between a distribution shift and G F has been plotted in Fig. 2(a) for uniform U (-1, 1), normal N (0, 1), Student's-t(ν = 2), lognormal LN (1, 0.5) and g-and-h GH(g = 1, h = 0) distributions.

The zero-centered unimodal symmetric distributions (normal and Student's-t) have their maximal G F value when the bias is null. The G F curve for the uniform distribution reaches G F U = 1/3 when the distribution is centered or shifted by ±1. For intermediate absolute values of the bias, the folded distribution is not uniform and presents higher values of G F . The value decreases when the bias is large enough to exclude zero from the range of non-null densities. Having heavier tails, the Student's-t distribution has a larger Gini coefficient than the normal. Any added bias leads to a decrease of G F .

The decay curves are symmetric with respect to the sign of the bias. This is no longer the case for asymmetric distributions (lognormal and g-and-h), for which the peak is reached for non-null values of the bias and the decay curves are non-symmetric.

In Fig. 2(b) we plotted similar curves for distributions centered on their mode before adding a bias (the symmetric distributions have been left unchanged). This shows that the maximal value of G F is reached when the mode of the distribution is at, or near, the origin. This assertion is validated in the next section (Sect. 4.2).

Another important point illustrated by these curves is that the level of information that can be recovered from G F is not uniform over the range of G F values. For instance, G F = 0.41 might as well occur for a centered normal distribution as for positively or negatively biased Student'st, lognormal of g-and-h distributions, whereas values above 0.41 exclude normal and uniform distributions. Within the restrictions on the distributions shapes we considered above, one might infer that a value smaller than G F U = 1/3 is likely to reveal a biased error distribution, while a value above G F N = 0.41 is likely to signal distributions with one or two extended tail(s) (compared to normal tails), and with a possible bias. Between these bounds lies a blind zone where compensation between bias and shape factors prevent any inference on either of them.

Does mode-centering maximize G F ?

In order to avoid the blind zone effect observed above and be able to characterize the shape of a distribution from its G F value, mode-centering seems to offer an interesting way to relieve the bias/shape compensation. Mode-centering a distribution results in a folded distribution where both tails overlap and mix, but it ensures that the contribution of the most-extended tail will prevail. In absence of mode-centering, when a biased distribution has a large tail encompassing zero, the folding around zero might considerably reduce this tail.

The assertion that mode-centering maximizes G F is tested here by comparing the results for mode-centering with those obtained by explicitly maximizing the Gini coefficient with respect to a bias value. We define b max as the value of the bias which maximizes

G F b max = max b G F (|E -b|) ( 15
)
and note

G F max = G F (|E -b max |).
This equation is solved numerically by the Nelder and Mead optimizer. [START_REF] Nelder | A Simplex Method for Function Minimization[END_REF] The values of b max and G F max were computed for the literature datasets and compared to the mode m(E) and G M CF respectively, through z-scores

z b = m(E) -b max u(m(E)) 2 + u(b max ) 2 (16)
and

z G = G M CF -G F max u(G M CF ) 2 + u(G F max ) 2 (17)
where the uncertainties are estimated by bootstrap.

In the hypothesis of a normal distribution of z-scores, a test threshold of 2 is generally chosen for the absolute value of the z-score. [START_REF] Kacker | Assessing differences between results determined according to the guide to the expression of uncertainty in measurement[END_REF] For absolute values above 2, there is less than 5 percent of probability that the difference is due to sampling effects. For values below, one does not reject the hypothesis that the compared values are equal. [START_REF] Pernot | Probabilistic performance estimators for computational chemistry methods: Systematic improvement probability and ranking probability matrix[END_REF] Histograms for the z-scores are shown in Fig. 3. At the exception of one point, the absolute value of all z-score values are smaller than 2, and we have therefore no reason to reject the hypothesis that these values are equal considering their uncertainty. The outlying point, with z b = -3.8 and z G = -2.5 corresponds to the MP2 method in dataset ZAS2019, which has a practically normal distribution. [START_REF] Pernot | Impact of non-normal error distributions on the benchmarking and ranking of Quantum Machine Learning models[END_REF] One has G M CF = 0.418 (6) and G F max = 0.436(3) for a distance of 1.07 between the mode and b max , to be compared to the standard deviation of the distribution, 1.7. As seen in Fig. 2(a), there is a flat area near the top of the G F curve as a function of bias for a normal distribution: very small deviations from a perfect normal distribution (as hinted to by the value of G F max being larger than 0.41) can deviate the optimal point over a wide range.

For all practical purposes in the present study, one can therefore estimate that mode-centering maximizes the value of the Gini coefficient, at a fraction of the computing cost for the search for b max . We further note that for distributions having distinct mean and mode, centering on the mean would not maximize G F and therefore preserve some of the ambiguity due to bias/shape compensation.

G M CF versus G F

We note the Gini coefficient of mode-centered folded distributions G M CF . Fig. 4(a) displays G M CF versus G F for the literature datasets. It is clear that mode-centering increases all G values, i.e. G M CF ≥ G F for all datasets, within the estimation uncertainties. The G-scale is now reduced to The uncertainties are reported in Fig. 4(b), showing that for some datasets the uncertainty on G M CF is larger than the uncertainty on G F , up to a factor two. This extra uncertainty is due to the uncertainty on the mode value. We note also a size effect, the smallest datasets (THA2015, N = 145) having the largest uncertainty, and the largest dataset (ZAS2019, N = 6211), the smallest one.

G M CF as a shape statistic

In parallel with the Gini coefficient, the skewness of the distribution has also been considered as an estimator of inequality. [START_REF] Bendel | Comparison of skewness coefficient, coefficient of variation, and Gini coefficient as inequality measures within populations[END_REF] One is interested here in comparing G M CF with β M CF , which is the skewness β GM (Eq. 9) of the mode-centered folded distribution.

The values for our selection of literature datasets is shown in Fig. 4(c). There is an excellent correlation between those statistics, considering the uncertainties reported in Fig. 4(d). Using model datasets of large size (N = 10 6 ) for Student's-t and g-and-h distributions with a range of shape parameters, one observes a nearly perfect non-linear correlation (dashed line, resulting from a quadratic fit of the sampled values). The dispersion of the points for the literature datasets about this curve is mostly due to statistical uncertainty (the points for the largest dataset (ZAS2019, N = 6211) are very close to the curve). It is important to note that the uncertainty on G M CF is a factor two to five smaller than the uncertainty on β M CF , and therefore performs better for smaller datasets.

We can therefore conclude that G M CF is apt at estimating the heaviness of the errors distribution tail after mode-centering and folding. In order to appreciate the information about the signed error distribution that can be extracted from G M CF , we plotted it against the skewness β GM and excess kurtosis κ CS (Fig. 5). Considering skewness (Fig. 5(a)), all the points seem to lie within an angular sector, indicating that distributions with large skewness have necessarily large G M CF values. For instance, if the absolute value of the skewness is above 0.3, G M CF is larger than 0.5. Reciprocally, G M CF provides only an upper limit to |β GM | (for G M CF = 0.55, the absolute value of the skewness cannot be above 0.4). For kurtosis (Fig. 5(b)), there is a lax positive trend between both statistics, and, globally, large values of G M CF are associated with large excess kurtosis, which might be due to heavy tails or outliers. In both graphs, values of G M CF below 0.45 are associated with low skewness and excess kurtosis. Although information is lost because of folding, G M CF can still provide some information about the shape of the distribution of signed errors, and notably about the kurtosis.

Let us consider a few examples to illustrate this point. We see in Fig. 5 that most points fall between 0.4 and 0.55, but a few methods reach higher values. The largest G M CF value in this study is 0.66 (orange dot) for method CAM-B3LYP in the PER2018 set 2 (cf. Table 2). This corresponds to large values of both β GM and κ CS . The authors discussed how this DFA is in the head group of two methods with similar MUE values, but does not minimize the risk of large errors (Sect. II.A 2 ). From the same set, B3LYP has the second largest G M CF value (0.61) and presents the same tail features than CAM-B3LYP. The third largest G M CF value (0.61, violet dot) belongs to the ZAS2019 dataset, and it presents a null skewness and a large kurtosis. An in-depth study has been published for this case, [START_REF] Pernot | Impact of non-normal error distributions on the benchmarking and ranking of Quantum Machine Learning models[END_REF] where the errors distribution for the SLATM-L2 method was shown to have large tails, despite having the smallest MUE among the compared methods. In the same set, the MP2 method has a practically normal errors distribution and can be found in the lower part of the Gini scale (0.42). This analysis can be repeated for the ten methods with largest G M CF values (Table 2), showing that large G M CF values point indeed towards error sets with high kurtosis and/or skewness.

Application of G M CF to ranking

To evaluate the interest of G M CF in a practical scenario, one might consider it as an alert mechanism to complement a MUE ranking. But a question remains: "What is the threshold for G M CF one should use to detect problematic error distributions ?" We have seen above that the ten largest G M CF values, above 0.55, point to distributions with notable tails. We propose for the present study to adopt an alert value based on the median of the G M CF values for our full dataset (0.51) and round it to 0.5. This might be reevaluated when more data are analyzed. Using this threshold, one might flag distributions suspected of having tails unsuitable for reliable predictions. Fig. 6(a) shows, for each dataset, the flagging of the methods with best ranking. If one considers the first rank, five methods are flagged, but it is striking that three datasets have all their 10 lowest MUE-ranked methods flagged (BOR2019, PER2018 and THA2015). For BOR2019, all methods present some excess kurtosis and variable levels of skewness. We can relate this to an increasing trend of the errors with the bandgap value. [START_REF] Pernot | Probabilistic performance estimators for computational chemistry methods: Systematic improvement probability and ranking probability matrix[END_REF] In the case of PER2018, the error distributions present also large skewness and kurtosis, that can be associated with the chemical heterogeneity of the dataset. [START_REF] Pernot | Probabilistic performance estimators for computational chemistry methods: the empirical cumulative distribution function of absolute errors[END_REF] For THA2015, it was noted previously [START_REF] Thakkar | How well do static electronic dipole polarizabilities from gasphase experiments compare with density functional and MP2 computations?[END_REF][START_REF] Pernot | Probabilistic performance estimators for computational chemistry methods: Systematic improvement probability and ranking probability matrix[END_REF] that some experimental reference data with large measurement uncertainty could not be reproduced by any method in the studied set. These outliers contribute to the tails of all the error distributions (so-called global outliers) and affect G M CF values. Note that, more generally, reference data are not necessarily the origin of global outliers, as a missing physical contribution in the tested methods could produce similar effects.

To explore the role of global outliers, Fig. of global outliers, defined as systems lying out of the [q(0.025), q(0.975)] interval for all methods of a dataset. The removal of 6 systems affects strongly the case THA2015, confirming the previous analysis. The results for BOR2019 are mostly unchanged, except for the best MUE-ranked method (mBJ) which benefits from the removal of a single global outlier. No effect is observed for the PER2018 dataset, confirming the intrinsic heavy-tailed shape of these heterogeneous atomization energy error sets. [START_REF] Pernot | Probabilistic performance estimators for computational chemistry methods: the empirical cumulative distribution function of absolute errors[END_REF] The other datasets with leading G M CF -flagged methods are ZAS2019 and ZHA2018. The former has already been discussed (Sect. 4.4), and the removal of several global outliers has no impact. In the case ZHA2018, the first MUE-ranked method is SCAN, which has a G M CF value just above the threshold (0.503), presents no skewness and a slight level of kurtosis. Removal of 4 global outliers does not improve the shape of the errors distribution. However, the errors on the formation enthalpies present a linear trend as a function of the calculated values. Correcting this trend 1 improves slightly the performance and the shape of the SCAN distribution, but most notably of the PBE distribution, which performances get indistinguishable from those of SCAN (Table 3).

Limits of the G M CF coefficient

We have shown above that G M CF might be a useful complement to the usual ranking statistics, in order to detect error distributions with shapes that might reveal problem in prediction reliability. However there remains cases where the G M CF index is insufficient to reveal underlying problems. Fig. 7 proposes a scenario of two normal distributions (G M CF = 0.41) with the same value of the MUE (1.0), and yet very different risks of large prediction errors. This is clearly a case showing that a quantile-based statistic, such as Q 95 , is an essential complement to the MUE.

Conclusion

The Gini coefficient presents an interesting addition to the computational chemistry benchmarking statistical toolbox. We focused here on its properties in relation with features of the error distributions, such as bias and shape (skewness and kurtosis). The interest of the Gini coefficient is that it correlates with these features, and offers a one-number summary. This is also one of its weaknesses, as there is no unique mapping from the Gini coefficient to these features.

To unscramble this situation, we propose to use the Gini coefficient of the mode-centered distributions, G M CF , which offers a simpler to interpret, shape-based, measure of tailedness. Large G M CF values, e.g. above 0.5, alert us about large tails that might be due to large skewness, kurtosis and/or the presence of outliers. For high ranking methods, this is an incentive to inspect for some tail statistics. The probability to have absolute errors larger than 1.0 is P 1 = 0.50 for the blue curve and 0.42 for the red curve, hiding the fact that the red distribution contains much worse results that the blue one. In this case, the problem is solved by the values of Q 95 , giving 1.16 for the blue curve, and 2.46 for the red one. Shape statistics, such as the kurtosis, would not enable to discriminate between both normal distributions. closely the error distributions and check if the selected methods might have problems of reliability in their predictions. It might then be worth to investigate if the distorted shape of the distribution is due to systematic trends in the errors, as they can often be corrected by simple linear transformations. [START_REF] Lejaeghere | Ab initio based thermal property predictions at a low cost: An error analysis[END_REF][START_REF] Pernot | Prediction uncertainty of density functional approximations for properties of crystals with cubic symmetry[END_REF][START_REF] Lejaeghere | Is the error on first-principles volume predictions absolute or relative?[END_REF][START_REF] Proppe | Uncertainty quantification for quantum chemical models of complex reaction networks[END_REF][START_REF] Proppe | Reliable estimation of prediction uncertainty for physicochemical property models[END_REF] The impact of such corrections on the shape of error distributions is the prospect of further studies.

Figure 1 :

 1 Figure 1: Schematic cumulative density function (a) and Lorenz curve (b) for a folded standard normal density function of absolute errors. The area above the CDF (slanted) is the mean unsigned error (MUE). For a given probability p , the ratio of the shaded area to the total slanted area gives the value of the Lorenz curve L p = L(p ). Q p is the quantile for probability p . The area between the Lorenz curve and the identity axis (vertical bars) is half the Gini coefficient.

Figure 2 :

 2 Figure 2: Variation of G F with a bias value added to several model distributions: (a) no centering applied; (b) mode-centering applied before adding bias (Uniform, Student and Normal are unchanged).

Figure 3 :

 3 Figure 3: Histograms of z-scores for (a) the position of the mode versus the G F maximizer; and (b) the values of the corresponding Gini coefficients.

Figure 4 :

 4 Figure 4: Comparison of G M CF with other statistics: (a) correlation of G F and G M CF for the literature datasets, and (b) comparison of their uncertainties (the dashed line has a slope of 2); (c) correlation of G M CF and β M CF for the literature datasets (points) and a series of large samples (N = 10 6 ) of Student's-t and g-and-h distributions with a range of shape parameters (dashed line), and (d) comparison of their uncertainties (the points are bracketed by lines of slope 2 and 5).

Figure 5 :

 5 Figure 5: Comparison of G M CF with shape statistics of error distributions: (a) absolute value of skewness |β GM |; (b) excess kurtosis κ CS .

Figure 6 :

 6 Figure 6: Flagging of large G M CF values after mode-centering of at most 10 lowest MUE-ranked methods in each dataset.

Figure 7 :

 7 Figure 7: Example of different error distributions having the same MUE (1.0) and offering contradictory resultsfor some tail statistics. The probability to have absolute errors larger than 1.0 is P 1 = 0.50 for the blue curve and 0.42 for the red curve, hiding the fact that the red distribution contains much worse results that the blue one. In this case, the problem is solved by the values of Q 95 , giving 1.16 for the blue curve, and 2.46 for the red one. Shape statistics, such as the kurtosis, would not enable to discriminate between both normal distributions.

Table 1 :

 1 Literature datasets: N is the number of systems in the dataset and K is the number of methods.

	Property	N	K Source
	BOR2019 Band gaps (eV)	471 15 32
	NAR2019 Enthalpies of formation (kcal/mol)	469 4 33
	PER2018 Intensive atomization energies (kcal/mol) 222 9 2
	SCH2018 Adsorption energies (eV)	195 7 34
	THA2015 Polarizability (relative errors, in %)	135 7 35
	WU2015 Polarizability (relative errors, in %)	145 36 36
	ZAS2019 Effective atomization energies (kcal/mol) 6211 3 37
	ZHA2018 Solid formation enthalpies (kcal/mol)	196 2 38

Table 2 :

 2 The ten methods with the largest G M CF values, and the corresponding skewness and kurtosis.

Table 3 :

 3 Statistics for the methods of the ZHA2018 dataset, before and after linear correction ('lc-' prefix).

	Methods	MUE	MSE	Q 95	β GM	κ CS	G M CF
		(kcal/mol) (kcal/mol) (kcal/mol)
	PBE	0.2106(98) -0.205(10) 0.467(34)	-0.043(64) -0.17(28) 0.408(18)
	SCAN	0.1024(69) -0.0165(97) 0.291(21)	-0.007(67) 1.49(53) 0.503(19)
	lc-PBE	0.0923(61) 0.0	0.287(38)	-0.138(71) 1.02(43) 0.479(24)
	lc-SCAN 0.0917(63) 0.0	0.276(26)	-0.082(67) 1.27(41) 0.475(20)

6(b)

reports the same analysis after search and removal

Data availability statement

The data and codes that support the findings of this study are openly available at the following URL: https://doi.org/10.5281/zenodo.4333217

Supplementary Information

Statistics, ECDFs and Lorenz curves for the literature datasets are openly available at the following URL: https://doi.org/10.5281/zenodo.4333217

Dedication

This work is dedicated to Ramon Carbó-Dorca for his 80th birthday. It reflects his interest for cross-disciplinary aspects of science, especially the role of mathematics in chemistry.


domain: stat info

text:
Malik Tiomoko
email: malik.tiomoko@u-psud.fr

Tiomoko Hafiz
email: hafiz.tiomoko.ali@huawei.com

Ali

Romain Couillet
email: romain.couillet@gipsa-lab.grenoble-inp.fr

DECIPHERING AND OPTIMIZING MULTI-TASK LEARNING: A RANDOM MATRIX APPROACH

This article provides theoretical insights into the inner workings of multi-task and transfer learning methods, by studying the tractable least-square support vector machine multi-task learning (LS-SVM MTL) method, in the limit of large (p) and numerous (n) data. By a random matrix analysis applied to a Gaussian mixture data model, the performance of MTL LS-SVM is shown to converge, as n, p → ∞, to a deterministic limit involving simple (small-dimensional) statistics of the data.

We prove (i) that the standard MTL LS-SVM algorithm is in general strongly biased and may dramatically fail (to the point that individual single-task LS-SVMs may outperform the MTL approach, even for quite resembling tasks): our analysis provides a simple method to correct these biases, and that we reveal (ii) the sufficient statistics at play in the method, which can be efficiently estimated, even for quite small datasets. The latter result is exploited to automatically optimize the hyperparameters without resorting to any cross-validation procedure.

Experiments on popular datasets demonstrate that our improved MTL LS-SVM method is computationally-efficient and outperforms sometimes much more elaborate state-of-the-art multi-task and transfer learning techniques.

INTRODUCTION

The advent of elaborate learning machines capable to surpass human performances on dedicated tasks has reopened past challenges in machine learning. Transfer learning, and multitask learning (MTL) in general, by which known tasks are used to help a machine learn other related tasks, is one of them. The particularly interesting aspects of multi-task learning lie in the possibility (i) to exploit the resemblance between the datasets associated to each task so the tasks "help each other" and (ii) to train a machine on a specific target dataset comprised of few labelled data by exploiting much larger labelled datasets, however composed of different data. Practical applications are numerous, ranging from the prediction of student test results for a collection of schools [START_REF] Aitkin | Statistical modelling issues in school effectiveness studies[END_REF], to survival of patients in different clinics, to the value of many possibly related financial indicators [START_REF] Greg | Marketing models of consumer heterogeneity[END_REF], to the preference modelling of individuals in a marketing context, etc.

Since MTL seeks to improve the performance of a task with the help of related tasks, a central issue to (i) understand the functioning of MTL, (ii) adequately adapt its hyperparameters and eventually (iii) improve its performances consists in characterizing how MTL relates tasks to one another and in identifying which features are "transferred". The article aims to decipher these fundamental aspects for sufficiently general data models.

Several data models may be accounted for to enforce relatedness between tasks. A common assumption is that the data lie close to each other in a geometrical sense [START_REF] Evgeniou | Regularized multi-task learning[END_REF], live in a low dimensional manifold [START_REF] Agarwal | Learning multiple tasks using manifold regularization[END_REF], or share a common prior [START_REF] Daumé | Bayesian multitask learning with latent hierarchies[END_REF]. We follow here the latter assumption in assuming that, for each task, the data arise from a 2-class Gaussian mixture. 1Methodologically, in its simplest approach, MTL algorithms can be obtained from a mere extension of support vector machines (SVM), accounting for more than one task. That is, instead of finding the hyperplane (through its normal vector ω) best separating the two classes of a unique dataset, [START_REF] Evgeniou | Regularized multi-task learning[END_REF] proposes to produce best separating hyperplanes (or normal vectors) ω 1 , . . . , ω k for each pair of data classes of k tasks, with the additional constraint that the normal vectors take the form ω i = ω 0 + v i for some common vector w 0 and dedicated vectors v i . The amplitude of the vectors v i is controlled (through an additional hyperparameter) to enforce or relax task relatedness. We study this approach here. Yet, to obtain explicit and thus more insightful results, we specifically resort to a least-square SVM (as proposed e.g., in [START_REF] Xu | Multi-output least-squares support vector regression machines[END_REF]) rather than a margin-based SVM. This only marginally alters the overall behavior of the MTL algorithm and has no impact on the main insights drawn in the article. Moreover, by a now well-established universality argument in large dimensional statistics, [START_REF] Mai | A large scale analysis of logistic regression: Asymptotic performance and new insights[END_REF] show that quadratic (leastsquare) cost functions are asymptotically optimal and uniformly outperform alternative costs (such as margin-based methods or logistic approaches), even in a classification setting. This argument further motivates the choice of considering first and foremost the LS-SVM version of MTL-SVM.

Technically, the article exploits the powerful random matrix theory to study the performance of the MTL least-square SVM algorithm (MTL LS-SVM) for data arising from a Gaussian mixture model, assuming the total number n and dimension p of the data are both large, i.e., as n, p → ∞ with p/n → c ∈ (0, ∞). As such, our work follows after the recent wave of interest into the asymptotics of machine learning algorithms, such as studied lately in e.g., (Liao & Couillet, 2019;[START_REF] Deng | A model of double descent for highdimensional binary linear classification[END_REF][START_REF] Mai | A random matrix analysis and improvement of semi-supervised learning for large dimensional data[END_REF][START_REF] El | On information plus noise kernel random matrices[END_REF]. Our analysis reveals the following major conclusions:

• we exhibit the sufficient statistics, which concretely enable task comparison in the MTL LS-SVM algorithm; we show that, even when data are of large dimensions (p 1), these statistics remain small dimensional (they only scale with the number k of tasks); • while it is conventional to manually set labels associated to each dataset within {-1, 1}, we prove that this choice is largely suboptimal and may even cause MTL to severely fail (causing "negative transfer"); we instead provide the optimal values for the labels of each dataset, which depend on the sought-for objective: these optimal values are furthermore easily estimated from very few training data (i.e., no cross-validation is needed); • for unknown new data x, the MTL LS-SVM algorithm allocates a class based on the comparison of a score g(x) to a threshold ζ, usually set to zero. We prove that, depending on the statistics and number of elements of the training dataset, a bias is naturally induced that makes ζ = 0 a largely suboptimal choice in general. We provide a correction for this bias, which again can be estimated from the training data alone; • we demonstrate on popular real datasets that our proposed optimized MTL LS-SVM is both resilient to real data and also manages, despite its not being a best-in-class MTL algorithm, to rival and sometimes largely outperform competing state-of-the-art algorithms.

These conclusions thus allow for an optimal use of MTL LS-SVM with performance-maximizing hyperparameters and strong theoretical guarantees. As such, the present article offers through MTL LS-SVM a viable fully-controlled (even better performing) alternative to state-of-the-art MTL.

Reproducibility. Matlab and Julia codes for reproducing the results of the article are available in the supplementary materials.

Notation. e

[n]

m ∈ R n is the canonical vector of R n with [e [n] m ] i = δ mi . Moreover, e [2k] ij = e [2k] 2(i-1)+j . Similarly, E [n] ij ∈ R n×n is the matrix with [E [n] ij ] ab = δ ia δ jb . The notations A ⊗ B and A B
for matrices or vectors A, B are respectively the Kronecker and Hadamard products. D x is the diagonal matrix containing on its diagonal the elements of the vector x and A i• is the i-th row of A. The notation Å is used when a centering operation is performed on the matrix or vector A. Uppercase calligraphic letters (A, K, Γ, M, V,...) are used for deterministic small dimensional matrices. Finally, 1 m and I m are respectively the vector of all one's of dimension m and the identity matrix of dimension m × m. The index pair i, j generally refers to Class j in Task i.

RELATED WORKS

Let us first point out the difference between MTL and transfer learning: while MTL makes no distinction between tasks and aims to improve the performance of all tasks, transfer learning aims to maximize the performance of a target task with the help of all source tasks. Yet, both methods mostly sharing the same learning process, in this section, we mainly focus on the MTL literature, which is divided into parameter-based versus feature-based MTL.

In the parameter-based MTL approach, the tasks are assumed to share some parameters (e.g., the hyperplanes best separating each class) or their hyperparameters have a common prior distribution. Existing learning methods (SVM, logistic regression, etc.) can then be appropriately modified to incorporate these relatedness assumptions. In this context, [START_REF] Evgeniou | Regularized multi-task learning[END_REF][START_REF] Xu | Multi-output least-squares support vector regression machines[END_REF][START_REF] Parameswaran | Large margin multi-task metric learning[END_REF] respectively adapt the SVM, LS-SVM, and Large Margin Nearest Neighbor (LMNN) algorithms to the MTL paradigm. The present article borrows ideas from [START_REF] Evgeniou | Regularized multi-task learning[END_REF]; [START_REF] Xu | Multi-output least-squares support vector regression machines[END_REF].

In the feature-based MTL approach, the tasks data are instead assumed to share a low-dimensional common representation. In this context, most of the works aim to determine a mapping of the ambient data space into a low-dimensional subspace (through sparse coding, deep neural networks, principal component analysis, etc.) in which the tasks have high similarity [START_REF] Argyriou | Multi-task feature learning[END_REF][START_REF] Maurer | Sparse coding for multitask and transfer learning[END_REF][START_REF] Zhang | Deep model based transfer and multi-task learning for biological image analysis[END_REF][START_REF] Sinno Jialin Pan | Domain adaptation via transfer component analysis[END_REF]; other works simply use a feature selection method by merely extracting a subset of the original feature space [START_REF] Obozinski | Multi-task feature selection[END_REF][START_REF] Wang | Safe screening for multi-task feature learning with multiple data matrices[END_REF][START_REF] Gong | Multi-stage multi-task feature learning[END_REF]. We must insist that, in the present work, our ultimate objective is to study and improve "data-generic" MTL mechanisms under no structural assumption on the data; this approach is quite unlike recent works exploiting convolutive techniques in deep neural nets or low dimensional feature-based methods to perform transfer or multi-task learning mostly for computer vision.

From a theoretical standpoint though, few works have provided a proper understanding of the various MTL algorithms. To our knowledge, the only such results arise from elementary learning theory (Rademacher complexity, VC dimension, covering number, stability) and only provide loose performance bounds [START_REF] Baxter | A model of inductive bias learning[END_REF][START_REF] Ben | Exploiting task relatedness for multiple task learning[END_REF][START_REF] Baxter | A bayesian/information theoretic model of learning to learn via multiple task sampling[END_REF]. As such, the present work fills a long-standing gap in the MTL research.

THE MULTI-TASK LEARNING SETTING

Let X ∈ R p×n be a collection of n independent data vectors of dimension p. The data are divided into k subsets attached to individual "tasks". Specifically, letting

X = [X 1 , . . . , X k ], Task i is a binary classification problem from the training samples X i = [X (1) i , X (2) i ] ∈ R p×ni with X (j) i = [x (j) i1 , . . . , x (j) inij ] ∈ R p×nij the n ij vectors of class j ∈ {1, 2} for Task i. In particular, n = k i=1 n i and n i = n i1 + n i2 for each i ∈ {1, . . . , k}.
To each x il ∈ R p of the training set is attached a corresponding "label" (or score) y il ∈ R. We denote y i = [y i1 , . . . , y ini ] T ∈ R ni the vector of all labels for Task i, and y = [y T 1 , . . . , y T k ] T ∈ R n the vector of all labels. These labels are generally chosen to be ±1 but, for reasons that will become clear in the course of the article, we voluntarily do not enforce binary labels here.

Before detailing the multitask classification scheme, a preliminary task-wise centering operation is performed on the data, i.e., we consider in the following the datasets

Xi = X i I ni - 1 n i 1 ni 1 T ni , ∀i ∈ {1, . . . , k}.
As such, we systematically work on the labeled datasets ( X1 , y 1 ), . . . , ( Xk , y k ). Remark 1 in the supplementary material motivates this choice, which avoids extra biases produced by the algorithm.

THE OPTIMIZATION FRAMEWORK

The multitask learning least square support vector machine (MTL LS-SVM) aims to predict, for input vectors x ∈ R p not belonging to the training set, their associated score y upon which a decision on the class allocation of x is taken, for a given target task. To this end, based on the labeled sets ( X1 , y 1 ), . . . , ( Xk , y k ), MTL LS-SVM determines the normal vectors

W = [ω 1 , ω 2 , . . . , ω k ] ∈ R p×k and intercepts b = [b 1 , b 2 , . . . , b k ] T ∈ R k
defining k separating hyperplanes for the corresponding k binary classification tasks. In order to account for task relatedness, each ω i assumes the form

ω i = ω 0 + v i for some common ω 0 ∈ R p and task-dedicated v i ∈ R p . Formally, writing V = [v 1 , . . . , v k ] ∈ R p×k (so that W = ω 0 1 k + V )
and following the work of [START_REF] Evgeniou | Regularized multi-task learning[END_REF][START_REF] Xu | Multi-output least-squares support vector regression machines[END_REF], the optimization function is given by min

(ω0,V,b)∈ R p ×R p×k ×R k 1 2λ ω 0 2 + 1 2 k i=1 v i 2 γ i + 1 2 k i=1 ξ i 2 , ξ i = y i -( XT i ω i + b i 1 ni ), 1 ≤ i ≤ k.
In this expression, the parameter λ enforces more task relatedness while the parameters γ 1 , . . . , γ k enforce better classification of the data in their respective classes.

Being a quadratic optimization problem under linear equality constraints, ω 0 , V, b are obtained explicitly (see details in Section 1 of the supplementary material). The solution is best described through the expression of the hyperplanes ω 1 , . . . , ω k ∈ R p which take the form:

ω i = e [k]T i ⊗ I p AZα,
and b = (P T QP ) -1 P T Qy, where α = Q(y -P b) = Q 1 2 (I n -Q 1 2 P (P T QP ) -1 P T Q 1 2 )Q 1 2 y ∈ R n is the Lagrangian dual and Q = 1 kp Z T AZ + I n -1 ∈ R n×n , Z = k i=1 E [k] ii ⊗ Xi ∈ R pk×n A = D γ + λ11 T ⊗ I p ∈ R kp×kp , P = k i=1 E [k] ii ⊗ 1 ni ∈ R n×k with γ = [γ 1 , . . . , γ k ] T and D γ = diag(γ).
MTL LS-SVM differs from a single-task joint LS-SVM for all data in that the data X1 , . . . , Xk are not treated simultaneously but through k distinct filters: this explains why Z ∈ R kp×n is not the mere concatenation [ X1 , . . . , Xk ] but a block-diagonal structure isolating each Xi . As such, the Xi 's-relating matrix A plays an important role in the MTL learning process.

With this formulation for the solution (W, b), the prediction of the class of any new data point x ∈ R p for the target Task i is then obtained from the classification score

g i (x) = 1 kp e [k] i ⊗ x T AZα + b i (1)
where x = x -1 ni X i 1 ni is a centered version of x with respect to the training dataset for Task i.

LARGE DIMENSIONAL STATISTICAL MODELLING

The first objective of the article is to quantify the MTL performance, and thus of the (a priori intricate) statistics of g i (x), under a sufficiently simple but telling Gaussian mixture model for training and test data.

Assumption 1 (Distribution of X and x). The columns of [X, x] are independent Gaussian random variables. Specifically, the n ij samples x

(j) i1 , . . . , x (j)
inij of class j for Task i are independent N (µ ij , I p ) vectors, and we let ∆µ i ≡ µ i1 -µ i2 . As for x, it follows an independent N (µ x , I p ) vector.

In the supplementary material, Assumption 1 is relaxed to [X, x] arising from a generative model of the type x

(j) il = h ij (z (j) il ) for z (j)
il ∼ N (0, I p ) and h ij : R p → R p a 1-Lipschitz function. This model encompasses extremely realistic data models, including data arising from generative networks (e.g., GANs [START_REF] Goodfellow | Generative adversarial nets[END_REF]) and is shown in the supplementary material to be universal in the sense that, as n, p → ∞, the asymptotic performances of MTL LS-SVM only depend on the statistical means and covariances of the x (j) il : the performances under complex mixtures are thus proved to coincide with those under an elementary Gaussian mixture. This generalized study however comes at the expense of more complex definitions and formulas, which impedes readability; hence the simpler isotropic Gaussian mixture model here.

Our central technical approach for the performance evaluation of the MTL LS-SVM algorithm consists in placing ourselves under the large p, n regime of random matrix theory. Assumption 2 (Growth Rate). As n → ∞, n/p → c 0 > 0 and, for 1

≤ i ≤ k, 1 ≤ j ≤ 2, nij n → c ij > 0. We let c i = c i1 + c i2 , c = [c 1 , . . . , c k ] T ∈ R k .
With these notations and assumptions, we are in position to present our main theoretical results.

THE MULTI-TASK LEARNING ANALYSIS

TECHNICAL STRATEGY AND NOTATIONS

To evaluate the statistics of g i (x) (equation 1), we resort to finding so-called deterministic equivalents for the matrices Q, AZQ, etc., which appear at the core of the formulation of g i (x). Those are provided in Lemma 1 of the supplementary material. Our strategy then consists in "decoupling" the effect of the data statistics from those of the MTL hyperparameters λ, γ 1 , . . . , γ k . Specifically, we extract two fundamental quantities for our analysis: a data-related matrix M ∈ R 2k×2k and a hyperparameter matrix A ∈ R k×k :

M = k i,j=1 ∆µ T i ∆µ j E [k] ij ⊗ c i c T j , c i =   ci2 ci ci1 ci -ci1 ci ci2 ci   , A = I k + D -1 2 ∆ D γ + λ1 k 1 T k -1 D -1 2 ∆ -1
where ∆ = [ ∆1 , . . . , ∆k ] T are the unique positive solutions to the implicit system ∆i = ci c0 -A ii (this is implicit because A is a function of the ∆j 's). In passing, it will appear convenient to use the shortcut notation ∆ = [ ∆11 , . . . , ∆k2 ] T ∈ R 2k where ∆ij = cij ci c 0 ∆i . We will see that M plays the role, in the limit of large p, n, of a sufficient statistic for the performance of the MTL LS-SVM algorithm only involving (i) the data statistics µ 11 , . . . , µ k2 and (ii) the (limiting) relative number c 11 /c 1 , . . . , c k2 /c k of elements per class in each task. As for A, it captures the information about the impact of the hyperparameters λ, γ 1 , . . . , γ k and of the dimension ratios c 1 , . . . , c k and c 0 . These two matrices will be combined in the core matrix Γ ∈ R 2k×2k of the upcoming MTL LS-SVM performance analysis, defined as

Γ = I 2k + A ⊗ 1 2 1 T 2 M -1
where we recall that ' ' is the Hadamard (element-wise) matrix product.

We raised in the introduction of Section 3 that we purposely relax the binary "labels" y ij associated to each datum x ij in each task to become "scores" y ij ∈ R. This will have fundamental consequences to the MTL performance. Yet, since x i1 , . . . , x ini are i.i.d. data vectors, we impose equal scores y i1 = . . . = y ini within each class. As such, we may reduce the complete score vector y ∈ R n under the form y = ỹ11 1 T n11 , . . . , ỹk2 1 T n k2

T for ỹ = [ỹ 11 , . . . , ỹk2 ] T ∈ R 2k . From Remark 1 of the supplementary material, the performances of MTL are insensitive to a constant shift in the scores y i1 and y i2 of every given Task i: as such, the recentered version ẙ = [ ẙ11 , . . . , ẙk2 ] T of ỹ, where ẙij = ỹij -( ni1 ni ỹi1 + ni2 ni ỹi2 ) will be central in the upcoming results.

MAIN RESULTS

Theorem 1 (Asymptotics of g i (x)). Under Assumptions 1-2, for x ∼ N (µ x , I p ) with µ x = µ ij ,

g i (x) -G ij → 0, G ij ∼ N (m ij , σ 2 i ) in distribution where, letting m = [m 11 , . . . , m k2 ] T , m = ỹ -D -1 2 ∆ ΓD 1 2 ∆ẙ , σ 2 i = 1 ∆i ẙT D 1 2 ∆Γ D K i. c ⊗12 + V i ΓD 1 2 ∆ẙ with V i = 1 c0 (AD c 0 K i• c +e [k] i A ⊗ 1 2 1 T 2 ) M and K = [A A](I k -D c 0 kc [A A]) -1 .
As anticipated, the statistics of the classification scores g i (x) mainly depend on the data statistics µ i j and on the hyperparameters λ and γ 1 , . . . , γ k through the matrix Γ (and more marginally through V i and K for the variances).

Since g i (x) has a Gaussian limit centered about m ij , the (asymptotic) standard decision for x to be allocated to Class 1 (x → C 1 ) or Class 2 (x → C 2 ) for Task i is obtained by the "averaged-mean" test

g i (x) C1 ≷ C2 1 2 (m i1 + m i2 ) (2)
the classification error rate i1 ≡ P (x → C 2 |x ∈ C 1 ) of which is then

i1 ≡ P g i (x) ≥ m i1 + m i2 2 x ∈ C 1 = Q m i1 -m i2 2σ i + o(1) (3)
with m ij , σ i as in Theorem 1 and

Q(t) = ∞ t e -u 2 2 du.
Further comment on Γ is due before moving to practical consequences of Theorem 1. From the expression (A ⊗ 1 2 1 T 2 ) -1 M, we observe that: (i) if λ 1, then A is diagonal dominant and thus "filters out" in the Hadamard product all off-diagonal entries of M, i.e., all cross-terms ∆µ T i ∆µ j for i = j, therefore refusing to exploit the correlation between tasks; (ii) if instead λ is not small, A may be developed (using the Sherman-Morrison matrix inverse formulas) as the sum of a diagonal matrix, which again filters out the ∆µ T i ∆µ j for i = j, and of a rank-one matrix which instead performs a weighted sum (through the γ i 's and the ∆i 's) of the entries of M. Specifically, letting γ -1 = (γ -1 1 , . . . , γ -1 k ) T in the expression of A, we have

D γ + λ1 k 1 T k -1 = D -1 γ -λγ -1 γ -1T 1+λ 1 k k i=1 γ -1 i .
As such, disregarding the regularization effect of the ∆i 's, the off-diagonal ∆µ T i ∆µ j entry of M is weighted with coefficient (γ i γ j ) -1 : the impact of the γ i 's is thus strongly associated to the relevance of the correlation between tasks.

A fundamental aspect of Theorem 1 is that it concludes that the performances of the large dimensional (n, p 1) classification problem at hand merely boils down to 2k-dimensional statistics, as all objects defined in the theorem statement are at most of size 2k. More importantly from a practical perspective, these "sufficient statistics" are easily amenable to fast and efficient estimation: it only requires a few training samples to estimate all quantities involved in the theorem. This, as a corollary, lets one envision the possibility of efficient transfer learning methods based on very scarce data samples as discussed in Remark 1.

Estimating m ij and σ i not only allows one to anticipate theoretical performances but also enables the actual estimation of the decision threshold 1 2 (m i1 + m i2 ) in equation 2 and opens the possibility to largely optimize MTL LS-SVM through an (asymptotically) optimal choice of the training scores ỹ. Indeed, the asymptotics in Theorem 1 depend in an elegant manner on the training data labels (scores) ỹ. Since the variance σ2 i is independent of the classes, we easily determine the vector ỹ = ỹ minimizing the misclassification probability for Task i as

ỹ = arg max ỹ∈R 2k (m i1 -m i2 ) 2 σ 2 i = arg max ỹ∈R 2k ỹT (I 2k -D 1 2 ∆ΓD -1 2 ∆ )(e [2k] i1 -e [2k] i2 ) 2 ỹT D 1 2 ∆Γ(D K i. c ⊗12 + V i )ΓD 1 2 ∆ ỹ the solution of which is explicit ỹ = D -1 2 ∆ Γ -1 H[(A⊗1 2 1 T 2 ) M]D -1 2 ∆ (e [2k] i1 -e [2k] i2 ), H ≡ (D K i. c ⊗12 + V i ) -1 (4)
with corresponding (asymptotically) optimal classification error i1 (equation 3) given by

i1 = Q 1 2 (e [2k]
i1 -e

[2k] i2 ) T G(e [2k] i1 -e [2k] i2 ) , G = D 1 2 ∆[(A⊗1 2 1 T 2 ) M ]V -1 i [(A⊗1 2 1 T 2 ) M ]D
The only non-diagonal matrices in equation 4 are Γ and V i in which M plays the role of a "variance profile" matrix. In particular, assuming ∆µ T i ∆µ = 0 for all = i (i.e., the statistical means of all tasks are orthogonal to those of Task i), then the two rows and columns of M associated to Task i are all zero but on the 2 × 2 diagonal block. Therefore, ỹ must be filled with zero entries but on its Task i two elements. All other values at the zero entry locations of ỹ (such as the usual ỹ = [1, -1, . . . , 1, -1] T ) would be suboptimal and possibly severely detrimental to the classification performance of Task i (not by altering the means m i1 , m i2 but by increasing the variance σ 2 i ). This extreme example strongly suggests that, in order to maximize the MTL performance on Task i, one must impose low scores ỹjl to all Tasks j strongly different from Task i.

The choice ỹ = [1, -1, . . . , 1, -1] T can also be very detrimental when ∆µ T i ∆µ j < 0 for some i, j: that is, when the mapping of the two classes within each task is reversed (e.g., if Class 1 in Task 1 is closer to Class 2 than Class 1 in Task 2). In this setting, it is easily seen that ỹ = [1, -1, . . . , 1, -1] T works against the classification and performs much worse than a single-task LS-SVM.

Another interesting conclusion arises from the simplified setting of equal number of samples per task and per class, i.e., n 11 = . . . = n k2 . In this case,

ỹ = Γ -1 H (A ⊗ 1 2 1 T 2 ) M (e [2k]
i1 -e

[2k]
i2 ) in which all matrices are organized in 2 × 2 blocks of equal entries. This immediately implies that ỹ j1 = -ỹ j2 for all j. So in particular the detection threshold 1 2 (m i1 + m i2 ) of the averaged-mean test (equation 2) is zero, as conventionally assumed. In all other settings for the n jl 's, it is very unlikely that ỹ i1 = -ỹ i2 and the optimal decision threshold must also be estimated.

These various conclusions give rise to an improved MTL LS-SVM algorithm. A pseudo-code (Algorithm 1) along with discussions on (i) the estimation of the statistics met in Theorem 1 (Remark 1) and on (ii) its multi-class extension (Remark 2) are covered next. Matlab and Julia implementations of Algorithm 1 and its extensions are proposed in the supplementary material.

EXPERIMENTS

Our theoretical results (the data-driven optimal tuning of MTL LS-SVM, as well as the anticipation of classification performance) find various practical applications and consequences. We exploit them here in the context of transfer learning, first on a binary decision on synthetic data, and then on a multiclass classification on real data.

To this end, and before proceeding to our experimental setup, a few key remarks to practically exploit Theorem 1 are in order. Remark 1 (On the estimation of m ij and σ i ). All quantities defined in Theorem 1 are a priori known, apart from the k 2 inner products ∆µ T i ∆µ j for 1 ≤ i, j ≤ k. For these, define S il ⊂ {1, . . . , n il } (l = 1, 2) and the corresponding indicator vector j il ∈ R ni with [j ij ] a = δ a∈Sij . For i = j, further let S il ⊂ {1, . . . , n il } with S il ∩ S il = ∅ and the corresponding indicator vector j il ∈ R ni . Then, the following estimates hold:

∆µ T i ∆µ j - j i1 |S i1 | - j i2 |S i2 | T X T i X j j j1 |S j1 | - j j2 |S j2 | = O (p min l∈{1,2} {|S il |, |S jl |}) -1 2 ∆µ T i ∆µ i - j i1 |S i1 | - j i2 |S i2 | T X T i X i j i1 |S i1 | - j i2 |S i2 | = O (p min l∈{1,2} {|S il |, |S il |}) -1 2 .
Observe in particular that a single sample (two when i = j) per task and per class (|S il | = 1) is sufficient to obtain a consistent estimate for all quantities, so long that p is large.

Remark 2 (From binary to multiclass MTL). Accessing the vector m in Theorem 1 allows for the extension of the MTL framework to a multiclass-per-task MTL by discarding the well known inherent biases of multiclass SVM. In the context of L i classes for Task i, using a one-versusall approach (for each ∈ {1, . . . , L i }, one MTL LS-SVM algorithms with Class "1" being the target class and Class "2" all other classes at once), one needs to access L i pairs of values (m i1 ( ), m i2 ( )) and, for a new x, decide on the genuine class of x based on the largest value among g i (x; 1) -m i1 (1), . . . , g i (x; L i ) -m i1 (L i ) with g i (x; ) the output score for "Class versus all".

For simplicity, from Remark 1 in the supplementary material, one may choose smart shift vectors ȳ( ) ∈ R k for the scores ỹ( ) ∈ R 2k (i.e., replace ỹ( ) by ỹ( ) + ȳ( ) ⊗ 1 2 ) so that m i1 ( ) = 0 for each . Under this shift, the selected class is the one for which g i (x; ) is maximum.

EFFECT OF INPUT SCORE (LABEL) AND THRESHOLD DECISION CHOICES

In order to support the theoretical insights drawn in the article, our first experiment illustrates the effects of the bias in the decision threshold for g i (x) (in general not centered on zero) and of the input score (label) optimization ỹ on synthetic data. Specifically, MTL-LSSVM is first applied to the following two-task (k = 2) setting: for Task 1, x (j) 1l ∼ N ((-1) j µ 1 , I p ) and for Task 2, x (j) 2l ∼ N ((-1) j µ 2 , I p ), where µ 2 = βµ 1 + 1 -β 2 µ ⊥ 1 and µ ⊥ 1 is any vector orthogonal to µ 1 and β ∈ [0, 1]. This setting allows us to tune, through β, the similarity between tasks. For four different values of β, Figure 1 depicts the distribution of the binary output scores g i (x) both for the classical MTL-LSSVM (top displays) and for our proposed random matrix improved scheme, with optimized input labels (bottom displays).

As a first remark, note that both theoretical prediction and empirical outputs closely fit for all values of β, thereby corroborating our theoretical findings. In practical terms, the figure supports (i) the importance to estimate the threshold decision which is non-trivial (not always close to zero) and (ii) the relevance of an appropriate choice of the input labels to improve the discrimination performance between both classes, especially when the two tasks are not quite related as shown by the classification error presented in red in the figure.

MULTICLASS TRANSFER LEARNING

We next turn to the classical Office+Caltech256 [START_REF] Saenko | Adapting visual category models to new domains[END_REF][START_REF] Griffin | Caltech-256 object category dataset[END_REF] real data (images) benchmark for transfer learning, consisting of the 10 categories shared by both datasets. For fair comparison with previous works, we compare images using p = 4096 VGG features. Half of the samples of the target is randomly selected for the test data and the accuracy is evaluated over 20 trials. We use here Algorithm 1, the results of which (Proposed) are reported in Table 1 against the non-optimized LS-SVM [START_REF] Xu | Multi-output least-squares support vector regression machines[END_REF] and alternative state-of-the-art algorithms: MMDT, CDLS and ILS. 3Algorithm 1 Proposed Multi Task Learning algorithm.

Input: Training samples X = [X 1 , . . . , X k ] with X i = [X (1) i , . . . , X (Li) i
], X i ∈ R p×n i and test data x. Output: Estimated class ˆ ∈ {1, . . . , L t } of x for target Task t. for j = 1 to L t do Center and normalize data per task: for all i ∈ {1, . . . , k},

• Xi ← X i I n i -1 n i 1 n i 1 T n i
• Xi ← Xi / 1 n i p tr( Xi XT i ) Estimate: Matrix M from Remark 1 and ∆, using X

ỹ (j) = D -1 2 ∆ Γ -1 H[(A⊗1 2 1 T 2 ) M]D -1 2 ∆ (e [2k]
t1 -e

[2k] t2 ). Estimate and center m = m(j) from Theorem 1 as per Remark 1. (Optional) Estimate the theoretical classification error t (λ, γ) minimize over (λ, γ) 2 . Compute classification scores g t (x; j) according to equation 1. end for Output: ˆ = arg max ∈{1,...,Lt} g t (x; ). 1 demonstrates that our proposed improved MTL LS-SVM, despite its simplicity and unlike the competing methods used for comparison, has stable performances and is quite competitive. We further recall that, in additions to these high levels of performance, the method comes along with theoretical guarantees, which none of the competing works are able to provide.

CONCLUSION: BEYOND MTL

Through the example of transfer learning (and more generally multitask learning), we have demonstrated the capability of random matrix theory to (i) predict and improve the performance of machine learning algorithms and, most importantly, to (ii) turn simplistic (and in theory largely suboptimal) methods, such as here LS-SVM, into competitive state-of-the-art algorithms. As Gaussian mixtures are quite "universal" and thus already appropriate to handle real data (as shown in supplementary material), one may surmise the optimality of the least square approach, thereby opening the possibility to prove that MTL LS-SVM is likely close-to-optimal even in real data settings.

This is yet merely a first step into a generalized use of random matrix theory and large dimensional statistics to devise much-needed low computational cost and explainable, yet highly competitive, machine learning methods from elementary optimization schemes.

Figure 1 :

 1 Figure 1: Scores g 2 (x) [empirical histogram vs. theory in solid lines] for x of Class C 1 (red) or Class C 2 (blue) for Task 2 in a 2-task (k = 2) setting of isotropic Gaussian mixtures for: (top) classical MTL-LSSVM with y ∈ {±1} and threshold ζ = 0; (bottom) proposed optimized MTL-LSSVM with ỹ and estimated threshold ζ; decision thresholds ζ represented in dashed vertical lines; red numbers are misclassification rates; task relatedness with β = 0 for orthogonal tasks, β > 0 for positively correlated tasks, β < 0 for negatively correlated tasks; p = 100, [c 11 , c 12 , c 21 , c 22 ] = [0.3, 0.4, 0.1, 0.2], γ = 1 2 , λ = 10. Histograms drawn from 1 000 test samples of each class.

  k } as data of class 2. Create scores

Table 1 :

 1 Classification accuracy over Office+Caltech256 database. c(Caltech), w(Webcam), a(Amazon), d(dslr), for different "Source to target" task pairs (S → T ) based on VGG features. Best score in boldface, second-to-best in italic. Mean score LSSVM 96.69 89.90 92.90 90.00 93.80 78.70 93.50 95.00 85.00 90.20 94.70 100 91.70 MMDT 93.90 87.05 90.83 84.40 94.17 86.25 94.58 97.50 86.25 87.23 92.05 97.35 90.96 ILS 77.89 73.55 86.85 76.22 86.22 71.34 74.53 82.80 68.15 63.49 78.98 92.88 77.74 CDLS 97.60 88.30 93.54 88.30 93.54 92.50 93.54 93.75 93.75 88.30 97.35 96.70 93.10 Ours 98.68 89.90 94.40 90.60 94.40 93.80 94.20 100 92.50 89.90 98.70 99.30 94.70

	S/T	c →	w →	c →	a →	w →	a →	d →	w →	c →	d →	a →	d →
		w	c	a	c	a	d	a	d	d	c	w	w
		Table

In the supplementary material, we extend this setting to a much broader and more realistic scope, and justify in passing the relevance of a Gaussian mixture modelling to address multi-task learning with real data.

∆.

MMDT: Max Margin Domain Transform, proposed in (Hoffman et al., 2013), applies a weighted SVM on a learnt transformation of the source and target. CDLS: Cross-Domain Landmark Selection, proposed in (Hubert[START_REF] Tsai | Learning cross-domain landmarks for heterogeneous domain adaptation[END_REF], derives a domain invariant feature space. ILS: Invariant Latent Space, proposed in[START_REF] Herath | Learning an invariant hilbert space for domain adaptation[END_REF], learns an invariant latent space to reduce the discrepancy between source and target.

ACKNOWLEDGMENTS

Couillet's work is partially supported by MIAI at University Grenoble-Alpes (ANR-19-P3IA-0003) and the HUAWEI LarDist project.