Classification of Heavy-tailed Features in High Dimensions: a Superstatistical Approach
Authors: Urte Adomaityte
Gabriele Sicuro
Pierpaolo Vivo
Venue: NeurIPS 2023 poster (2023)
Keywords: Classification
Gaussian Mixture Model
Superstatistics
Empirical Risk Minimization
Replica theory
Power-law distribution

Abstract: We characterise the learning of a mixture of two clouds of data points with generic centroids via empirical risk minimisation in the high dimensional regime, under the assumptions of generic convex loss and convex regularisation. Each cloud of data points is obtained via a double-stochastic process, where the sample is obtained from a Gaussian distribution whose variance is itself a random parameter sampled from a scalar distribution $\varrho$. As a result, our analysis covers a large family of data distributions, including the case of power-law-tailed distributions with no covariance, and allows us to test recent ''Gaussian universality'' claims. We study the generalisation performance of the obtained estimator, we analyse the role of regularisation, and we analytically characterise the separability transition.

Classification of Heavy-tailed Features in High Dimensions: a Superstatistical Approach
Urte Adomaityte Department of Mathematics King's 's College London
Gabriele Sicuro Department of Mathematics
King's College London
urte. . wdw.aidttttteekcllccoo ac. uk
gabriele. sicuro@kcl. ac. uk
Pierpaolo Vivo
Department of Mathematics King's College London
pierpaolo. vivo@kcl. ac. uk
Abstract
We characterise the learning of mixture of clouds of data points with generic centroids via empirical risk minimisation in the high dimensional regime, under the assumptions of generic convex loss and convex regularisation. Each cloud of data points is obtained via double-stochanee drawn from a Gaussian distribution whose variance itself random parameter sampled from scalar distribution P. As result, our analysis covers large family data distributions, including the case of power-law-tailed distributions with no covariance, and allows to test recent "Gaussian universality" claims. We study the generalisation performance of the obtained estimator, we analyse the role regularisation, and analytically characterise the separability transition.
process, where the sample is
Introduction
Generalised linear models (GLMs) are still ubiquitous in the theoretical research on machine learning, despite their simplicity. Their nontrivial phenomenology often amenable to complete analytical treatment and has been key to understanding the unexpected behavior of large and complex archi- tectures. Random features models [61, 47, 24], for example, allowed clarify many aspects of the well-known double-descent phenomenon in neural networks [54. 8, 9]. A line of research spanning more than three decades [68, 74] has considered a variety of GLMs to investigate a number of aspects of learning in high dimensions. Yet, crucial assumption adopted in many such theoretical models that the covariates are obtained from a Gaussian distribution, or from a mixture of Gaussian distributions [42, 49, 3, 37, 43]. Although such Gaussian design has served
convenient working hypothesis (in some cases experimentally and theoretically justified [67, 44]), is however not obvious how much it limits the scope of the results. It is reasonable to expect structure, heavy tails, and large fluctuations in real (often non-Gaussian) data [1, 62) to play an important role in the learning process, and it would be therefore desirable to include structured and beavy-tailedddistribuut covariates in our theoretical toolbox.
This paper presents, to the best of our knowledge for the first time, the exact asymptotics for classification tasks on covariates obtained from a mixture of heavy-tailed distributions. We focus supervised binary classification assuming that the sample size and the dimensionality d of the space where the covariates live are both sent to infinity, keeping their ratio "/d a fixed. The paper fits therefore in the line of works exact high-dimensional asymptoties for classification learning via a GLM [49, 43], but, crucially, we relax the usual Gaussian data hypothesis by including in our
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

analysis power-law-tare law-t-la distributions with possibly no covariance. The mixture is obtained from two distributions, each centered around a centroid j E Rd and resulting from a double stochastic process. Namely, each sample point is obtained from a distribution N(H, whose covariance I Ald E Rdxd is itself random variable so that A has density e. supported on R Using, for example, an appropriately parametrised rnverse-Gamma distribution for e. such a superstatistical construction (as known in the physics literature [5, 7J) can provide, for example, a heavy-tailed, Cauchy-like data distribution with infinite covariance. The replica method [48], an analytical tool widely adopted in statistical physics, provides asymptotic formulas for a generic density e, allowing us to study the effects of non- non-Gaussianity on the performance curves, and test Gaussian universality hypotheses [59] in such a classification task.
Learning a rule to classify data points clustered in clouds in high dimension is a classical problem in statistics [32]. Its ubiquity is exemplified by the recently observed neural collapse phenomenon in deep neural networks, which the last layer classifier was found to operate on features clustered in clouds around the vertices of a simplex [38, 56]; more recently, Seddik et al. [67] showed that Gaussian mixtures suitably describe the deep learning representation of GAN data. In theoretical models, data points are typically assumed to be organised in K Gaussian clouds, each corresponding to a class. Each class k, k E {l..... .. K]. is centered around a mean vector Mg and has covariance Σg. The binary classification case, K 2. is the simplest, and possibly most studied, setting. For this case, Mai and Liao [46] considered an ERM task with generie convex loss and ridge regularisation in the high- -igh--imensional regime n, d too with n/a E (0, +00), In this setting, they gave a precise prediction of the classification error, showing the optimality of the square loss in the unregularised case. Their results have been extended by Mignacco et al. [49], who showed that the presence of a regularisation can actually drastically affect the performance, improving it. In the same setting, the Bayes optimal estimator has been studied in both the supervised and the cemi-supervised setting [49, 41]; mamnost-Byyes--. cottimal solutions have been put in relation to wide flat landscape regions [3]. In the context of binary classification, the maximum number of samples that can be perfectly fitted by a linear model in dimension d [49, 16, 37] has been the topic of investigation since the seminal work of Cover [12] and ís related to the classical storage capacity problems on networks [20, 39]. The corresponding separability transition value a n/a remarkably associated with the existence transition of the maximum likelihood estimator [70. 76]. The precise asymptotic characterisation of the test error in learning to classify K 2 2 Gaussian clouds with generic means and covariances has been recently obtained by Loureiro al. [43]. Within this line of research, rigorous results have been obtained by variety of methods, such as Gordon's 's inequality technique [31, 71, 49] or mapping to approximate message passing schemes [4, 36, 10, 43, 27].
Motivation and related works
As previously mentioned, the working hypothesis of Gaussian design is widely adopted in high- dimensional classification problems. On top of being convenient technical hypothesis, this assump- tion has been justified in terms of a "Gaussian universality" principle. In other words, in a number of cireumstances, non-Gaussian features distributions are found to be effectively described by Gaussian ones with matching first and second moments as far as the asymptotic properties of the estimators are concerned [44]. Such a "universality property" has been rigorously proven for example in compressed sensing [50] or in the case of LASSO with non-Gaussian dictionaries [55]. A Gaussian equivalence principle introduced by Goldt et al. [30] has been proven to hold for wide class of generalised linear estimation problems [47]. Extending work by Hu and Lu [34], Montanari and Saeed [51] recently proved such a principle in a GLM under the assumption of pointwise normality of the distribution of the features. This crucial assumption can be intended, roughly speaking, the assumption of sub-Gaussian decay of a marginal of the feature distribution in any direction (see also [25, 13) for further developments).
The Gaussian universality principle, however, can break down by relaxing some of the aforementioned assumptions on the feature distribution. Montanari and Saeed [5 for example showed that pointwise normality is a necessary hypothesis for universality. Studying regression tasks on an elliptically distributed dataset, El Karoui [18] showed that claims of "universality" obtained in the Gaussian setting require serious scrutiny, as the statistics of the estimators might strongly depend on the non- Gaussianity of the covariates [18, 69, 72]. More recently. Pesce et al. [59] showed that a structured Gaussian mixture cannot be effectively described by a single Gaussian cloud. Building upon series of contributions related to the asymptotic performance of models in the proportional high-dimensional limit [33, 47, 28, 24, 30, 44], we aim precisely to explore the validity of Gaussian universality within classification problems.
2

---

In this paper. we work in a "superstatistical" data setting [5. 7], meaning that we superpose the statistics of Gaussian data distribution with an assumed distribution its variance. Such construction is adopted in number of disciplines and contexts to go beyond Gaussianity, albeit is known under different names. In statistical physics, it is known as "superstatistics" and is employed the analysis of non-equilibruum and non-linear systems [6]. In Bayesian modeling, it is common to refer to hierarchical priors and models [22, 23], while in probability, statistics, and actuarial sciences such distributions are known as compound probability distributions [63], doubly-stochastic models [60, 66]. Crucially, this construction allows us to consider a very large family of non-Gaussian distributions, which include, but are not limited to, any power-law decay and Cauchy-like with possible ciinite-ovvarance parametrisations.
Our contributions In this manuscript provide the following results.
We study a classification task on a non-Gaussian mixture model (see Eq. (1) below) via a generalised
linear model (GLM) and analytically derive, using the replica method [14, 19, 48), an asymptotic characterisation of the statistics of the empirical risk minimisation (ERM) estimator. Our results go therefore beyond the usual Gaussian assumption for the dataset, and include for example the case of covariates obtained from a mixture of distributions with infinite variance. The analysis is performed in the high-dimensional, proportional limit and for any convex loss and convex regularisation. By using this result, we provide asymptotic formulas for the generalisation, training errors and training loss.
We analyse the performance of this ERM task on a specific family of dataset distributions by using different convex loss functions (quadratic and logistic) and ridge regularisation. We show in particular that, in the case of two balanced clusters with a specific non- non-Gaussian distribution, the optimal ridge regularisation strength 1 is finite, at odds with the Gaussian case, for which x - co [49]. this respect, by considering distributions with matching first and second moments, we analytically show that the performances the analysed GLM do depend in general on higher moments, and therefore the "Gaussian universality principle" breaks down when heavy-tailed distributions are considered.
We derive the separability threshold on large family of non-Gaussian dataset distributions, possibly with unbounded covariance, generalising the known asymptotics for the separability of Gaussian clouds [49]. The result of Cover [12] is recovered in the case of infinite distribution width.
Under some moment conditions, we derive the Bayes-optimal performance in the case of binary classifications with symmetric centroids, generalising the argument in Ref. [49].
We finally extend recent results on Gaussian universality in the Gaussian mixture model with random labels [25, 599 to the case of non-Gaussian distributions. We show that the universal formula
for the square training loss at zero regularisation found in Ref. [25] holds in more general setting as well.
2 Main result
Dataset construction We consider the task of classifying two data clusters in the d-dimensional space R.. The dataset D = {(x", v"))rr[[n] is obtained by extracting n independent datapoints x", each associated with a label y" E 1) (a more general setting, involving a multiclass classification task, is discussed in Appendix A). The data points are correlated with the labels via a law Pxx,y) which we assume to have the form
P(x, y) 6(x1))(x())))++++++))))))))))))) (1-1
P PE(O.1) 1),
H_ EER.
(la)
Here P(x| (4) is distribution with mean μl, and the mean vectors μ. E Rd are distributed according to some density, such that B [|] = O(1). and correspond to the center of the two clusters. The scalar quantity o weighs the relative contribution of the two clusters: in the following, we will denote P+ = === p-..
p.. Each cluster distribution P(x|μt around a vector ş is assumed to have the form
(1b) where N(x|u, Σ) is a Gaussian distribution with mean f E Id and covariance E E Rdxd and A is randomly distributed with some density e with support on Rt == (0, +oo). The family of "ellipticllike" distributions in Eq. (1b) has been extensively studied, for instance, by the physics community, in the context of superstatistics [5, 7J. Mixtures of normals in the form of Eq. 1b are a central tool in Bayesian statistics [65] due to their ability to approximate any distribution given sufficient number
P(x|μ) = E. [N ATa)]
3

---

of components [53, 2, 29]. Although the family in Eq. (1b) is not the most general of such mixtures, it is sufficient include large family of power-law-tailes densities and allow us go beyond the usual Gaussian approximation. El Karoui [18] considered the statistical properties, in the same asymptotic proportional regime as here, of ridge-regularised regression estimators on datasets with distribution as in Eq. 1b, under the assumption (here relaxed) that B[A²] +oo. This elliptic family includes a large class of distributions with properties markedly different from the ones of Gaussians. For example, the inverse-Gamin with density ((A) (2m e leads to d-dimensional
Cauchy-like distribution p(xlμ) (1 - || x=1²)= _ having as marginals Cauchy distributions and [[||x ==l²] = +00. We will perform our classification task by searching for a set of parameters (w*,b b*). called weights and bias, respectively, that will allow to construct an estimator via certain classifier y: IR 1)
x Tw Vd
x
-
+ b*
By means of the above law, we will predict the label for new, unseen data point x sampled from the same law P(x,y) y). Our analysis will be performed in the high-dimensional limit where both the sample size n and dimensionality d are sent to infinity, with "/a a kept constant.
Learning task In the most general setting, the parameters are estimated by minimising an empirical risk function in the form
w TTxx
(w*,b*)= arg min R(w, b) where R(w,b) = weR"
b + +r(ww.
(3)
Vd
beR
v=1
Here is a strictly convex loss function with respect to its second argument, and r is a strictly convex regularisation function with the parameter à 2 0 tuning its strength.
State evolution equations. Let us now present our main result, namely the exact asymptotic characterisation of the distribution of the estimator w* Rd and of 1/vaw*T X R", where X E Rd is the concatenation the n dataset column vectors x E Rd, v E [n]. Such an asymptotic characterisation is performed via of order parameters satisfying system of self-consistent "state- evolution" equations, which we will solve numerically. The asymptotic expressions for generalisation and training errors, in particular, are written in terms of their solutions, the order parameters. In the following, we use the shorthand l- (u) = ((=1, u), and given a function Φ: [- - R, we write P. ((=1) and B. [+2]= = p.P+ + P_ Q_..
Result 2.1 Les 5 N(O, 1) and A o. both independent from other quantities. Let also be & ~ N(0, I)). In the setting described above, given 61 Rd - I and Φ2: R² - R, the estimator w* and the vector z* = X IR" verify:
(1(w*) n.d-+0o
Be (g)).
22(2*) n,d-seoo Bg [d2(h))
(4)
where we have introduced the proximal g Rd, defined as
(D-uml² Ex= mew Mk V Ar(w))
g arg minw
(5)
and where we also have introduced the proximal for the loss h RR" obtained by concutenating p+m quantities h4 with p_n quantities h-, with b2 given by
h = arg ming
(1=-0r)22 -2v 24v
+ +c4(w)
where w := m + b vqAC.
(6)
The collection of parameters (q. m.m v. q, m= V, appearing in the equations above is given by the fixed-point solution of the following self-consistent equations:
m+ - [g #_].
q = a...c.s |
fz := h-
q y [llg
9 aq B..... VAf+
vA
(7)
== (
m = Bp.. [f2].
b E [12 m²).

---

The derivation of the result above is given, in a multiclass setting, in Appendix A using the replica method, which can be put on rigorous ground by mapping of the risk minimisation problem into an approximate casssage-passing iteration scheme [17, 43, 26]. As in the purely Gaussian case, the state evolution equations naturally split the dependence on the loss function E, which is relevant in the computation of h2g and on the regularisation function entering in the computation of In the most general setting, we are required to solve two convex minimisation problems, consisting of the computation of the proximals h and g [57]; the convexity of the problem guarantees the uniqueness of the fixed point solution. Result 2.1 implies the following.
Result 2.2 Under the assumptions of Result 1.1, the training loss Et, the training error 6, and the test (or generalisation) error Eg are given, in the proportional asymptotic limit, by
1 n | w"ix"
+b*
- Bf
val
- 1(
w-'x" #y"
- (+))
(8)
vid
B(y
+b*
#y)
vil
Result 2.3 Under the hypothesis that B[A] < too and B[A- < too, the Bayes optimal test error for binary classification is
Ao
:
RO RD -pp.B(0cc)))pp_BB00.
K1
A/
BLAJELAI11 -3
1
2
aBIA
A
In
P-
(9)
where the expectation is over the two identically distributed random variables Do, A o, and (xx) is the complementary error function.
Quadratic loss with ridge regularisation. Explicit choices for or r can highly simplify the fixed- point equations above. As an example, let us consider the case of quadratic loss ((y,x) = (y x)) with ridge regularisation r(w) +|lw] 2) In this setting, an explicit formula for the proximal h2 in Eq. (6) and for the proximal in Eq. (5) can be easily found. As result, we can obtain explicit expressions for the corresponding set of state evolution equations which become particularly suitable for fast numerical solution and explicit the dependence on higher moments of A:
m+ Σx=- M.
q = aB (+1 m2 - -b)
+ taqBs
d+D
(I+rA)²
(I+vA)²
q=
E. Hhau P4g (2+0)²
= aBA 1AV
(10)
m ap= (+1 m b)BA J.
2+0
3 Application to synthetic datasets
In this section, we compare our theoretical predictions with the results of numerical experiments for a large family of data distributions. The results have been obtained using ridge (2-regularisation, and both quadratic and logistic losses, with various data cluster balances p. We will also assume, without loss of generality, that μ, = + u. where N(, I.). The synthetic data sets will be produced using, an inverse-Gammen Gaamaa-disrrbuute variance A, with density parametrised
((A) Qa.e (A) Flajaat e f
(11a)
depending on the shape parameter a > 0 and on the scale parameter c > 0. As result, each class will be distributed as
P(x|μ) (2<=x=((
(2c)*T(a+*))
(11b)
which produces a cloud centered in μ.. The reason for our choice is that the distribution Ca.e allows us easily explore different regimes by changing two parameters only. Indeed, the distribution has, for a > 1, covariance matrix Σ = B - (4) @ (x 4] = where 2² = - By fixing 2 and taking the limit a - +oo, p(x|μ) - et al. [49]. For (0, 1], instead, 2² = +00. Note that the choice of ivveseeGaama distribution is not that unusual. It has been adopted, for example, to describe Gaus-aa data in quantitative finance [15, 40] or econometrics models [52]. Finally, we will construct our label estimator in Eq. (2) with (xx) sign(x).
i.e., the Gaussian case discussed by Mignacco
5

---

3.1 Finite-covariance case
Let us start by considering a data distribution as in Eq. (11) with shape parameter a = c+ c++1> 1, which is therefore given by
04
0.3
0.2
P(x(p) = I(a) 41 (2a-2+
2 (-(0-1)T(a+421)
(12)
0.10
and decays as ||x|||-a-1 in the radial direction for ||x - Mlll >> O. As a consequence, the distribution has lima //B[llx- re|l*] +oo for k 2 2a. With this choice, all elements of the family Eg. (12) have the same covariance I B[ (x μ1) @ (x - 4 La as a is varied, including the Gaussian limit P(x - N( La) obtained for
0.05
0.00
0.2
+o0.
0.1-
Fig. we present the results of our numerical ex- periments using the square loss and small regularisa- tion. An excellent agreement between the theoretical predictions and the results of numerical experiments is found for a range of values of a and sample com- plexity both for balanced, i.e., equally sized, and unbalanced clusters of data (the plot for this case can be found in Appendix A.3). The test error Eg presents the classical interpolation peak at a = I, smoothened by the presence of non-zero regularisation strength d, and the typical double-descent behavior [54]. As a is increased, the results of Mignacco al. [49] for Gaussian clouds with P(x|μ) = N(x La) are approached. The a 400 curves correspond, there- fore, to the Gaussian mixture model we would have constructed by simply fitting the first and second mo- ments of each class in the finite-a case. The plot shows that, at given population covariance the classification of power-law distributed clouds is associated with different, and in particular smaller, er- rors than the corresponding task on Gaussian clouds with the same covariance: this sense, no Gaussian universality principle holds. This (perhaps) counter-intuitive effect stems from the fact that, although for all classes Σ = Id, in our example the mean absolute deviation smaller for finite
Test error Eg (solid line, top), train- ing error (center) and training loss Eç (bot- tom) as predicted by Eq. (8) in the balanced 1/2 case. The dataset distribution is parametrised as Eq. (12). The classifica- tion task is solved using quadratic loss with ridge regularisation with A 10-5. In the top figure, the dashed line corresponds to the Bayes optimal bound. Dots correspond to the average outcome of 50 numerical experiments in dimension d = 10². In our parametrisation, the population covariance is Σ La for all values of a and moreover, for a - +oo, the case of Gaussian clouds with the same cen- troids and covariance is recovered. For fur- ther details on the numerical solutions, see Appendix B.
Figure 1:
ld.
lim
B[ -ilx-=llll
F(a- Va- a a-+000
(13)
d=0400
[(a)
The learning process benefits from the fact that as a decreases while keeping the variance fixed, points move on average closer their means μ to compensate heavier tails. Finally, we observed that test errors are systematically smaller for unbalanced clusters than for balanced clusters (see Appendix A.3).
The same numerical experiment was repeated using the logistic loss for training. Once again, in Fig. 2, we focus on p = 1/2 and show that the theoretical predictions of the generalisation and training emrors agree with the results of numerical experiments for a range of sample complexity values a and various values of a. Just as for the square loss, the test error Eg of the Gaussian model (a - +00) is larger than the one observed for plwer-law distributed clouds at finite values of a. In the A - 0 limit, the typical interpolation cusp the generalisation error observed: the cusp occurs at different values of a as a is changed; interpolation occurs at smaller values of for larger (i.e., for "more Gaussian" distributions).
comparison with the test error 6, confirms that this cusp occurs the value of a where the training eттor becomes non-zero and the two training data clouds become non-separan non-separable [70, 16). This sharp transition in a also corresponds to the existence transition of maximum- maximum--likelihoo estimator in high-dimensional logistic regression, analysed for Gaussian data in Refs. [70, 49] and, in our setting, in Section 4 below. For larger regularisation strength, the cusp smoothens, and the training error becomes non-zero at smaller values of

---

3.2 Infinite-covariance case
0.4
The family of distributions specified by Eq. (11) allows to also consider the case of infinite covariance, i.e., o² +00, To test our formulas in this setting, we considered the case in which each cloud obtained by "contaminating"" standard Gaussian with an infinite-covariance distribution as in Eq. (11b) with c and < [35]. In other words, we use the den- sity @(A) = ( + r)6(A D), with E [0, I] interpolation parameter, for A in Eq. (1). Each class has therefore infinite co- variance for 0 < S I and the Gaussian case recovered for r = Fig. collects our re- sults for two balanced clouds, analysed using square loss and logistic loss. Good agreement between the theoretical predictions and the sults of numerical experiments is found for range of values of sample complexity and of the ratio r. In this case, the varianvariance case surprisingly, to the low- est test error with both square loss and logistic loss. which grows with
0.2
0.1
0.0
0.0
2
8 10
Figure 2: Test error Eg (top), training error 6, (cen- ter) and training loss Eg (bottom) via logistic loss training on balanced clusters parametrised as in Eq. (12) (I = Ia). Ridge regularisation with A = 10 4 is adopted. Dots correspond to the aver- age over 20 numerical experiments with d 10 The Gaussian limit recovered for a - +00, Fur- ther details the numerical solutions can be found in Appendix B.
corresponds,
3.3 The role of regularisation in the classification of non-Gaussian clouds
One of the main results in the work of Mignacco et al. [49] is related the effect of regularisation on the performances in a classification task on a Gaussian mixture. They observed that, for all values of the sample complexity a, the optimal ridge classification performance obtained for an infinite regularisation strength A, using both hinge loss and square loss. We tested the robustness of this result beyond the purely Gaussian setting. In Fig. 4 (left), we plot the test error obtained using square loss and different regularisation strengths à on the dataset distribution as in Eq. (12) with a 2. We observe that, both in the balanced case and in the unbalanced case, the optimal regularisation strength A is finite and, moreover, a-dependent. On the other hand, if we fix and let grow towards +00, the optimal 2 grows as well to recover the result of Mignacco et al. [49] of diverging optimal regularisation for balanced Gaussian clouds, as shown in Fig. 4 (center)
pair of balanced clouds is
0.50
0.45
0.40
0.40
0.35
0.35
0.30
7 0.8
0.30
0.25
0.25
0.4
0.2
0.15
0,10
Figure 3: Test error Eg in the classification of two balanced clouds, via quadratic loss (left) and logistic loss (right). In both cases, ridge regularisation is adopted (A 10-4 for the square loss case, 1 for the logistic loss case). Each cloud is superposition of power-law distribution with infinite variance and Gaussian with covariance Σ = Ia- The parameter r allows us to contaminate the purely Gaussian case 0) with an infinite-varian variance contribution (0 r 1) as Eq. (11) with = and 1/2 (leff) or 3/4 (right). Dots correspond to the average test error of 20 numerical experiments in dimension associated with the lowest test error for both losses.
10 Note that, at a given sample complexity, Gaussian clouds are

---

for a = 2. This result is represented also in Fig. 4 (right), where it is shown that, in the unbalanced case, the optimal regularisation strength 2* instead saturates finite value for - +00, i.c. for Gaussian clouds [49].
0.5
Se4tS
0,215
0.5 p
0.1
0.50
0.25
3
10y
0.25-
0.2
0,200
1o4t
0.150
0.1
0.30
0.125
0.25
0.100
100
10
0.15
0.75
e-30
2-
20
Figure 4: (Leff) Test error for ridge regularised quadratic loss for various regularisation strengths. The data points of each cloud in the training set are distributed as in Eq. (12), with shape parameter 2, for balanced clusters (top) and unbalanced clusters (p = 1/4, bottom). Points are the results of 50 numerical experiments, and the dashed lines are Bayes-optimal bounds. (Center) Test error for different regularisation strengths 2 for two balanced clusters with quadratic loss at sample complexity 2 using the data distribution (12). The optimal regularisation strength value A obtained from averaging 5 runs for each a is marked with a cross. (Right) Optimal regularisation strength 2 at 2 for different values of E [1.5, 10²] for both balanced and unbalanced clusters, obtained from averaging runs. Note that, for p 1/2, 2 - +00 as t00.
The separability threshold
By studying the training error e, with logistic loss at zero regularisation can obtain information on the boundary between the regime in which the training data are perfectly separable and the regime in which they are not. In other words, we can extract the value of the so-called separability transition complexity a* [70, 49, 43]. Once again, a complete characterisation of this transition point available in the Gaussian mixture case, where a* can be explicitly given function of 2 [49].
not trivial to extend this result to the general case of non-Gaussian mixtures. It is, however, possible to derive a* within the large family of "superstatistical" models we considering here under the assumption that Result 2. holds. In this case, in Appendix A.3 we prove the following.
Result 4.1 In the considered double-stochastic model, data are linearly separable for a*, where
a = max 1-02 8(0,y)=
/ PBB P.N 0+y 10.1 (ppN)))N f+ -- 101101) /o.1]ar. (14) VA
0e(0.11.78(0,y)
Fig. 5 shows the values of a for different choices of the shape parameters and c in the case of two balanced clouds of datapoints distributed in Eq. (lIb) as predicted by Eq. (14). In Fig. 5 (left), we fixed 2² q- +co, and plotted the separability threshold a* for a range of values of the shape parameter a 1. As grows, the known threshold value for the Gaussian mixture case is recovered [49]. The double limit a - co and 2 + co, on the other hand, provides the Cover [12] transition value a²
Fig. 5 (right), instead, we analyse the case of infinite-variance clusters, by fixing a e (0. 1) and by varying the scale parameter in Eq. (1)b), which controls the spread (width) of the distribution. The Cover transition therefore correctly recovered as e diverges for all values of a.
as expected.
8

---

10
10
10a
10a
3/4
109
109
10
10
10-2
10-b
10-4
10-3
10
10"
10
10-2
10
10
10-1
10-2
10
10"
10'
10a
Figure 5: Separability threshold a* obtained by solving the equations in Eq. (7) with logistic loss, ridge regularisation strength = 10 and 1/2. The data points of each cloud are distributed around their mean μ with two-parameter distribution in Eq. (1lb). (Left). Finite covariance ==02, case, x² = -- as function of a. Dashed lines are the threshold values of the Gaussian case derived by Mignacco et al. [49]. At large and large variance o², the Cover's transition a = 2 for balanced clusters recovered. (Right) udfnite-varianee data clouds obtained by fixing < a < in Eq. (1lb).
Random labels in the teacher-student scenario and Gaussian universality
Our analysis can be extended to the case in which point x is labeled not according to its class but by a "teacher", represented by a probabilistic law Po(yTT) and a vector 0, e Rd so that each dataset element independently generated with the joint distribution
P(x,y)= Po
10 Vd
Σ PeENN Ad)) k=*
(15)
Let us assume now that On such that lim 2//70J = 00, i.e., the teacher is "uncorrelated" with the structure of the mixture. In Appendix A.5 we show that, if CCy, v) = ((-y, -7) and Po(ylr) single cloud centered in the origin,
Po(-yl r), the task is equivalent to classification problem i... to the case μ = 0, so that (xx,y) Po (o(v) (/50, BB[N Al))]. This mean- inn-iinaria,ce, recently proven in Ref. [25] in the pure Gaussian setting, holds therefore our non-Gaussian setting as well.
special, relevant case of "uncorrelated teacher" is the one of random labels, where Po(ylz) 1/100,++1 + 6y.-1). -)) This apparently vacuous model relevant in the study of the storage capacity problem [21, 39, 11] and of the separability threshold [12], and it has also been effective tool to better understand deep learning [45, 75]. In this case, the following universal result holds with respect our setup.
Result 5.1 In the random label setting studied via ridge regression case with zero regularisation, the asymptotic training loss is independent of the distribution of A, and given by the universal formula
n a-1
0(a-1)=EE.
(16)
2n
Vd
a/Ea 2a
The proof of the previous result is given in Appendix A.5. Note that this result does not require 2² : B[A] to be finite, and it includes, as a special case, the Gaussian mixture model with random labels considered in Ref. [25, 59]. Therein, Eq. (16) has been obtained under the Gaussian design hypothesis. turns out therefore that Eq. (16) holds for a much larger ensemble, which includes heavy-tailed distributions with possibly infinite covariance. On the other hand, universality does not hold, e.g., for the training error.
In Fig. 6 we exemplify these results by using andge-regularssed square loss. The numerical experiments correspond a classification task on a dataset of two clouds centered in p., , -H_ N(0. a/l obtained using the distribution
P(x,y) y) 6y.11+6y.-1 Σ P(x)pg).
(17)

---

where P(xlu.) is the distribution in Eq. (11b) depending on the shape parameters and c. These results are compared with the theoretical prediction obtained for single cloud, with the same parameters a and c, but centered in the origin, so that P(x. = _u.utt P(x jr = 0): the agreement confirms the mmea--nvvariance" property stated above. As anticipated, as we are considering random labels, the training loss follows Eq. (16) in all cases, whereas the test and training errors are not universal. Note that, in place of the test error (which trivially Eg = 1/2) we plot the mean square error on the ve-activations &g = x(y,x) !/vax w J. Also, observe that this quantity divergent in the case of infinite covariance, hence is absent from the plot for E (0, 1J.
2
Conclusions and perspectives
The framework introduced in this work has allowed us to perform an exact asymptotic analysis of the learning of mixture of two data point clouds with possibly power-law tails. our dataset, each sample point x is labeled by y {-1,1} and obtained as yu Vāz, where R N(0, La) and ~ e. We have derived exact expressions for the test and training errors and training loss achieved by a GLM in the limit of large size n the dataset and large dimensionality d under the assumption that a = "/a is kept fixed. We have shown that the perfor- mance crucially depends higher-order moments of the data distribution, so that a Gaussian universal- ity assumption is feasible. Moreover, we have shown the effects of non- Con-Gaussiantty on the optimal regularisation strength and derived exact formula for the critical a* at which the two clouds become linearly inseparable. Finally, we have considered the so-called random labels setting (in which sample la- bel assignment randomised) and shown that, in this case, the training square loss takes a universal asymptotic form independent of the distribution of this result further supports the crucial role of labels' randomness in some recent universality results [25].
10
0.15
a=
0.10
a 5
0.05
0.00
0.2
0.1-
0.0
0.25
0.5
0.75
1.0
1.25
1.5
1.75
2.0
Figure 6: Mean square error &g (top), train- ing error E, (center) and training loss Er (bot- tom) obtained using square loss dataset distributed as Eq (17) with random labels (dots) compared with the theoretical predic- tions for the single cluster case (continuous lines). The distributions ing Eq. (12) for a and Eq. (11b) for 0
In conclusion, we have provided an exactly solvable toy model exhibiting power-law tails and possibly unbounded covariance, highlighting the fact that a dataset deviation from Gaussianity high dimen- sions cannot, in general, be neglected. The applica- tion of the introduced theoretical framework real datasets is an interesting direction to explore in the future: the main difficulty, this case, the choice of the best distribution given the observed dataset, a problem that has a long tradition the context Bayesian estimation [2, 2
parametrised us-
ridge regularisation with A = 10 is adopted. Note that & is absent in the top plot for 1/2 and 1 being divergent quantity in this case.
Acknowledgements
The authors would like to thank Nikolas Nüsken for stimulating discussions. Also, we would like to thank Bruno Loureiro for his feedback and for his careful reading of the manuscript.
I In more simplistic approach, it can be observed that in the case in which the square loss is adopted, equa- tions (10) depend (I v1)1) the dataset, for example, by empirically evaluating the identity B[(1 vA)- vvä So" (v u) 11/2 ax? ]d when it is assumed that x ([N( A)].
-11 and 3[(1 + =2] only. Such quantities can be numerically estimated from
10

---

References
[1] R.J. Adler, R.E. Feldman, and M.S. Taqqu, editors. A Practical Guide to Heavy Tails: Statistical Techniques and Applications. Birkhauser Boston Inc., USA, 1998. ISBN 0817639519.
[2] D. Alspach and H. Sorenson. Nonlinear Bayesian estimation using Gaussian sum approxima- tions. IEEE Trans. Autom. Control, 17(4):439-48,, 1972.
[3] C. Baldassi, E.M. Malatesta, M. Negri, and R. Zecchina. Wide flat minima and optimal generalization classifying high-dimensional gaussian mixtures. J. Stat. Mech.: Theory Exp., 2020(12) 2020.
[4] M. Bayati and A. Montanari. The dynamics of message passing dense graphs, with applica- tions to compressed sensing. IEEE Trans. Inf. Theory, 57(2):764-785, 2011.
[5] C. Beck. Superstatistics: Theory and applications. Contin. Mech. Thermodyn, 16, 2003.
[6] C. Beck. Recent developments in superstatistics. Braz. J. Phys., 39:3573363,,
[7] C. Beck and E.G.D. Cohen. Superstatistics. Physica A, 322:267-277 2003.
[8] M. Belkin, D. Hsu, S. Ma, and S. Mandal. Reconeiling modern machine-leaming, practice and the classical bias-variance trade-off. Proc. Natl. Acad. Sci. U.S.A., 16632)):1589--15544 2019.
[9] M. Belkin, A. Rakhlin, and A.B. Tsybakov. Does data interpolation contradict statistical optimality? K. Chaudhuri and M. Sugiyama, editors, Proceedings of the 22nd International Conference on Artificial Intelligence and Statistics, volume 89 of Proc. Mach. Learn. Res., pages 1611-1619. PMLR, 2019.
[10] R. Berthier, A. Montanari, and P.-M. Nguyen. State evolution for approximate message passing with non-separable functions. Inf. Inference, 9(1):33-77) 2020.
[11] N. Brunel, J.-P. Nadal, and G. Toulouse. Information capacity of a perceptron. JJ. Phys. A: Math. Gen., 25(19)50117, 1992.
[12] T.M. Cover. Geometrical and statistical properties of systems of linear inequalities with
applications in pattern recognition. IEEE Trans. Comput., EC-14(3)326-334, 1965.
[13] Y. Dandi, L. Stephan, F. Krzakala, B. Loureiro, and L. Zdeborová'. Universality laws for gaussian mixtures in generalized linear models. mrXii:23302.08933 2023.
[14] del Giudice, S. Franz, and M. Virasoro. Perceptron beyond the limit of capacity. Physique, 50:121-134, 1989.
[15] D. Delpini and G. Bormetti. Minimal model of financial stylized facts. Phys. Rev. E, 83:041111, 2011.
[16] A. Deng, Z.and Kammoun and C. Thrampoulidis. A model of double descent for high- dimensional binary linear classification. Inf. Inference, 11, 2021.
[17] D.L., Donoho, A. Maleki, and A. Montanari. Message-passing algorithms for compressed sensing. Proc. Natl. Acad. Sci. U.S.A., 106(45) Nov 2009.
[18] N. El Karoui. On the impact of predictor geometry the performance on high-dimensional ridge-regularized generalized robust regression estimators. Probab. Theory Relat. Fields, 170, 2018.
[19] S. Franz, D.J. Amit, and M.A. Virasoro. Prosopagnosia in high capacity neural networks storing uncorrelated classes. J. Physique, 51(5):387 1990.
[20] E Gardner. The space of interactions in neural network models. J. Phys. A: Math. Gen., 21(1): 257, 1988.
[21] E. Gardner and B. Derrida. Three unfinished works on the optimal storage capacity of networks.
J. Phys. A: Math. Gen., 22(12):998,, 1989.

---

[22] A. Gelman and J.1/ Hill. Data Analysis Using Regression and Multileve//Hierachiicd Models. Analytical Methods for Social Research. Cambridge University Press, 2006.
[23] A. Gelman, J.B. Carlin, H.S. Stern, D.B. Dunson, A. Vehtari, and D.B. Rubin. Bayesian Data Analysis, Third Edition. Chapman & Hall/CRC Texts in Statistical Science. Taylor & Francis, 2013.
[24] F. Gerace, B. Loureiro, F. Krzakala, M. Mézard, and L. Zdeborová. Generalisation error in learning with random features and the hidden manifold model. In ICML, pages 3452-3462, 2020.
[25] F. Gerace, F. Krzakala, B. Loureiro, L. Stephan, and L. Zdeborová. Gaussian universality of perceptrons with random labels. arXiv:2205.13300 2023.
[26] C. Gerbelot and R. Berthier. Graph-based approximate message passing iterations. arXiv preprint arXiv:2109.1 1995, 2021.
[27] C. Gerbelot, A. Abbara, and F. Krzakala. Asymptotic errors for teacher-student convex general- ized linear models (or: How to prove kabashima's replica formula). IEEE Trans. Inf. Theory, 69 (8):1824-1852, 2023.
[28] B. Ghorbani. S. Mei, T. Misiakiewicz. and A. Montanari. Limitations of lazy training of two-layers neural network. In H. Wallach, H. Larochelle, A. Beygelzimer, d'Alché-Bust E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32, 2019.
[29] J.K. Ghosh and R.V. Ramamoorthi. Bayesian Nonparametrics. Springer Series in Statistics. Springer New York, 2006.
[30] S. Goldt, M. Mézard, F. Krzakala, and L. Zdeborová. Modeling the influence of data structure on learning in neural networks: The hidden manifold model. Phys. Rev. X, 10:041044, 2020.
[31] Y. Gordon. Some inequalities for gaussian processes and applications. Isr. J. Math., 50:265-289, 1985.
[32] T. Hastie, R. Tibshirani, and J.H. Friedman. The elements of statistical learning: data mining, inference, and prediction, volume 2. Springer, 2009.
[33] T. Hastie, A. Montanari, S. Rosset, and R.J. Tibshirani. Surprises in high-dimensional ridgeless least squares interpolation. Ann. Stat., 50(2):949 986, 2022.
[34] H. Hu and Y.M. Lu. Universality laws for high-dimensional learning with random features. 20XY:20009006669 2022.
[35] P.J. Huber. Robust estimation of location parameter. Ann. Math. Stat., 35(1):73-101, 1964.
[36] A. Javanmard and A. Montanari. State evolution for general approximate message passing algorithms, with applications to spatial coupling. Inf. Inference, 2(2):115-144, 2013.
[37] G.R. Kini and C. Thrampoulidis. Phase transitions for one-vs-one and one-vs-all linear separa- bility in multiclass gaussian mixtures. ICASSP 2021, pages 4020-4024, 2021.
[38] V. Kothapalli, E. Rasromani, and V. Awatramani. Neural collapse: A review on modelling principles and generalization. arXiv:2206. 04041, 2022.
[39] W. Krauth and Mézard. Storage capacity of memory networks with binary couplings. Physique, 50000)3057-33666 1989.
[40] N. Langrené, G. Lee, and 7.. Zhu. Switching non-affine stochastic volatility: A closed-form expansion for the inverse gamma model. Int. J. Theor. Appl. Finance, 19, 2015.
[41] M. Lelarge and L.. Miolane. Asymptotic bayes risk for gaussian mixture in semi-supervised setting. In 2019 IEFE 8th International Workshop on Computational Advances in Senso- Adaptive Processing, pages 639-643, 2019.
12

---

[42] M. Lelarge and L. Miolane. Asymptotic bayes risk for gaussian mixture in semi- sume-supee setting. In 2019 IEEE 8th International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), pages 639 633-6433,3 2019.
[43] B. Loureiro, G. Sicuro, C. Gerbelot, A. Pacco, F. Krzakala, and L. Zdeborová. Learning Gaussian Mixtures with Generalized Linear Models: Precise Asymptotics in High-dimensions Advances in Neural Information Processing Systems, 14:10144-00157, 2021.
[44] B. Loureiro, C. Gerbelot, H. Cui, S. Goldt, F. Krzakala, M. Mézard, and Zdeborová. Learning curves of generic features maps for realistic datasets with a teacher-student model. J. Stat. Mech.:
Theory Exp., 2022(11): 14001, 2022.
[45] H. Maennel, L.M. Alabdulmohsin, I.O. Tolstikhin, R. Baldock, O. Bousquet, S. Gelly, and D. Keysers. What do neural networks learn when trained with random labels? In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, editors, Advances Neural Information Processing Systems, volume 33, pages 19693-19704. Curran Associates, Inc., 2020.
[46] X. Mai and Z. Liao. High dimensional classification via regularized and unregularized empirical risk minimization: Precise error and optimal loss. 1905. 2020.
[47] S. Mei and A. Montanari. The Generalization Error of Random Features Regression: Precise Asymptotics and the Double Descent Curve. Commun. Pure Appl. Math., 75, 2019.
[48] M. Mézard, G. Parisi, and M.A. Virasoro. Spin glass theory and beyond: An Introduction the Replica Method and Its Applications, volume 9, World Scientific Publishing Company, 1987.
[49] F. Mignacco. F. Krzakala, Y. Lu, P. Urbani, and L. Zdeborová. The role of regularization in classification high-dimensional noisy Gaussian mixture. In Hal Daumé III and Aarti Singh, editors, Proceedings of the 37th International Conference Machine Learning, volume 1 119 of Proc. Mach. Learn. Res., pages 6874-6883. PMLR, 2020.
[50] A. Montanari and P.-M. Nguyen. Universality of the elastic net error. In 2017 IEEE ISIT, pages 2338-2342 IEEE, 2017.
[51] A. Montanari and B.N. Saced. Universality of empirical risk minimization. In Po-Ling Loh and Maxim Raginsky, editors, Proceedings of 35th Conference on Learning Theory, volume 178 of Proceedings of Machine Learning Research, pages 4310 4312. PMLR, 2022.
[52] D.B. Nelson. Arch models as diffusion approximations. Econom., 45(1-2):7-38, 1990.
[53] V. Nestoridis, S. Schmutzhard, and V. Stefanopoulos. Universal series induced by approximate identities and some relevant applications. J. Approx. Theory, (63(1)):1783-119) 2011.
[54] M. Opper, W. Kinzel, J. Kleinz, and R. Nehl. On the ability of the optimal perceptron to generalise. J. Phys. A: Math. Gen., 23(1 238111, 1990.
[55] A. Panahi and B. Hassibi. A universal analysis of large-scale regularized least squares solutions. In Guyon, U. Von Luxburg, S. Bengio. H. Wallach, R. Fergus, Vishwanathan, and Garnett, editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
[56] V. Papyan, X.Y. Han, and D.L. Donoho. Prevalence of neural collapse during the terminal phase of deep learning training. Proc. Natl. Acad. Sci. U.S.A., 117(40))24652-24466 2020.
[57] N. Parikh and S. Boyd. Proximal algorithms. Found. Trends Optim., 1(3):127-23,, 2014.
[58] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python. J. Mach. Learn. Res.. 12:2825-288 2830, 2011.
[59] I. Pesce, F. Krzakala, B. Loureiro, and L. Stephan. Are Gaussian data all you need? Extents and limits of universality in high-dimensional generalized linear estimation. arXiv:2302.08892. 2023.
13

---

[60] M.A. Pinsky and S. Karlin. An Introduction to Stochastic Modeling: Fourth Edition, pages 1-563. Elsevier Inc, December 2010. ISBN 9780123814166
[61] A. Rahimi and B. Recht. Random features for large-scale kernel machines. In J. Platt, D. Koller, Y. Singer, and S. Roweis, editors, Advances in Neural Information Processing Systems, vol- ume 20. Curran Associates, Inc., 2007.
[62] W.J. Reed and B.D. Hughes. From gene families and genera to incomes and internet file sizes: Why power laws are so common in nature. Phys. Rev. E, 66:067103, 2002.
[63] H.E. Robbins. Asymptotically subminimax solutions of compound statistical decision problems. 1985.
[64] S. Rosset, J. Zhu, and T. Hastie. Margin maximizing loss functions. In S. Thrun, L. Saul, and B. Schölkopf, editors, Advances Neural Information Processing Systems, volume 16. MIT Press, 2003.
[65] P.E. Rossi. Bayesian Non- and Semi- mamiipprammtricc Methods and Applications. Princeton University Press, 2014.
[66] D. Schnoerr, R. Grima, and G. Sanguinetti. Cox process representation and inference for stochastic reaction-diffusion processes. Nat. Comm., 7:11729, 05 2016.
[67] M.E.A. Seddik, C. Louart, M. Tamaazousti, and Couillet. Random matrix theory proves that deep learning representations of GAN-data behave as Gaussian mixtures. In H.D. III and A. Singh, editors, Proceedings of the 37th International Conference on Machine Learning. volume 119 of Proc. Mach. Learn. Res., pages
PMLR, 2020.
[68] H. S. Seung, H. Sompolinsky, and N. Tishby. Statistical mechanics of learning from examples. Phys. Rev. A, 45:6056-6091, 1992.
[69] L. Steinberger and H. Leeb. Conditional predictive inference for stable algorithms. Ann. Stat., 51(1):290- 311, 2018.
[70] Sur and E.J. Candès. A modern maximum-likelihood theory for high-dimensional logistic regression. Proc. Natl. Acad. Sci. U.S.A., 116629)145166-152) 2019.
[71] C. Thrampoulidis, S. Oymak, and B. Hassibi. Regularized linear regression: A precise analysis of the estimation error. In Peter Grünwald, Elad Hazan, and Satyen Kale, editors, Proceedings of The 28th Conference on Learning Theory, volume 40 of Proc. Mach. Learn. Res., pages 1683-1709, Paris, France, 2015. PMLR
[72] C. Thrampoulidis, Abbasi, and B. Hassibi. Precise error analysis of regularized m -estimators in high dimensions. IEEE Trans. Inf. Theory, 64(8):5592-5628 2018.
[73] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. Nelson, Eric Jones, Robert Kern, Eric Larson, C Carey, ilhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul Mulbregt, and SeiPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17:261-272, 2020. doi: 10.
[74] T.L.H. Watkin, A. Rau, and M. Biehl. The statistical mechanics of learning rule. Rev. Mod. Phys., 65:499-556, 1993.
[75] C. Zhang. S. Bengio, M. Hardt, B. Recht, and O. Vinyals. Understanding deep learning (still) requires rethinking generalization. Commun. ACM, 64(3):107-115, 2021.
[76] Q. Zhao, P. Sur, and E.J. Candès. The asymptotic distribution of the MLE in high-dimensional logistic models: Arbitrary covariance, Bernoulli, 28(3):1835 1861, 2022.
14


Reviews:
[{"review_id": "fkTFOOpAiA", "reviewer": "Reviewer_K5kC", "summary": "The paper derive a theory for training and generalization error when classifying a large number of points from a non-Gaussian high-dimensional data distribution. The data model is a double-stochastic process where a parameter is sampled from a scalar distribution and then a sample is taken from a Gaussian distribution with this parameter as variance. A self-consistent mean-field theory is provided for the case where the number of points is large and proportional to the dimensionality and the equations can be numerically solved for logistic and square losses. The theory is applied to data-sets with finite and infinite covariance, to study the role of regularization, and to estimate the separability threshold for such data. This highlighting both cases of “Gaussian universality” where the results coincide with previous “Gaussian” literature and deviation from such universality.", "questions": " * Why bother with the classifier estimator phi? Is there any reasonable choice beyond sign?\n * Why do you refer to z* as “the matrix”?\n * Are all the 8 (or 10) order parameters scalars?\n * Can you clarify the interpretation of the proximal h and g? Do their distribution is a mean-field version of some real-worlds quantity?\n", "limitations": "The authors adequately addressed the limitations of their work. Those include the diagonal structure of the covariance matrix (conditional on the value of delta), the use of K=2 which leads to lack of discussion about the mean of the distribution (because they do not affect anything for K=2 beyond their norm), and the resulting theory solvable only numerically.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": " *\tOriginality: the tasks and methods are not new, and previous contributions are very well introduced. The originality of the work lies in the successful calculation of the theory for the non-Gaussian case, which is a valuable contribution. Furthermore, the work provides a basis for analyzing when an extrapolation beyond the Gaussian case is justified.\n *\tSignificance: the paper is important in highlighting where non-Gaussian data may diverge from the Gaussian case discussed in the literature, and as such it opens a venue for future work to use non-Gaussian analysis of real-worlds data, which is an important direction. The conclusions about test error in Gaussian vs non-Gaussian cases is non-trivial (the inversion between figure 1+2 and 3) and as is the finding or optimal finite regularization value for the non-Gaussian case. \n *\tClarity: the paper is in general well-written and can be served as exemplar for providing complicated theoretical results without sacrificing the clarity of the ideas.\n *\tQuality: the paper seems technically very sound, with impressive combination of theory and simulations.\n", "weaknesses": " *\tClarity: some of the notations and ideas presented are only hastily introduced, with two prominent examples being “Superstatistical Features” (from the title) and “uncountable superposition” (from the abstract, which seem overly complicated. To me, a presentation through “double stochastic” process (as in my summary) is straight-forward and require no extra jargon.\nAnother avenue for improving the understanding of the reader may lie in “Quadratic loss with ridge regularisation” where results are more amendable to interpretation. The authors should have provided more intuition for those results and furthermore point out where does non-Gaussian enters in the self-consistent equations (i.e., what part is shared with the Gaussian case).\n *\tOriginality: the main part of the work focus on reproducing known results from Gaussian literature and exploring the deviation from them for the non-Gaussian case. In that sense, there is no originality in this work beyond the (impressive) achievement of providing a theory which describe this non-Gaussian case.\n *\tSignificance: the work would have been more influential if it provided new tools which can be applied to datasets, where a small number of shape parameters is fit to non-Gaussian data and the ability to classify this data can be predicted from theory and then compared to actual classification of the data.\nCases where the Gaussian case predicts the behavior for the non-Gaussian case might deserve a fuller theoretical analysis, perhaps through the analysis suggested for “Quadratic loss with ridge regularisation” above.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for her/his remarks and positive evaluation of our paper, and for capturing the spirit of our contribution very well. We are grateful to the referee for her/his suggestions about improving the clarity of the manuscript, which we implemented in the new version. We also thank her/him for suggestions regarding possible applications of the theory, in particular, in its simplest form given in the case of ridge regression: we have added a comment, in this sense, in a newly introduced *Conclusions and perspectives* section.\n\nWe list our answers to her/his questions below, hoping that they will be satisfactory.\n* The referee is right about the fact that, in the given setting, the sign function is indeed the only reasonable choice, and the one we indeed adopted through the text. We decided to leave the classifier more general as, in the new version, we have extended the replica calculation to the multiclass case in the Appendix, where the classifier’s choice is less obvious. \n* Thanks for pointing out the typo regarding $\\boldsymbol z^\\star$, which is, indeed, a vector!\n* All order parameters are scalars: in the paper, we tried to consistently use bold fonts for vectorial/matricial objects (e.g., $\\boldsymbol g$) and normal fonts for scalars.\n* The proximals $h_\\pm$ can be seen as expressing the statistics of the preactivation $\\frac{\\boldsymbol w^\\star{}^\\intercal \\boldsymbol x^\\nu}{\\sqrt d}+b^\\star$ when $\\boldsymbol x^\\nu$ is a training set datapoint with label $y^\\nu=\\pm1$. The proximal $\\boldsymbol g$, instead, captures the statistics of $\\boldsymbol w^\\star$ itself (as expressed by Eq. 8)."}, {"review_id": "ku1WBCPcL0", "reviewer": "Reviewer_26qY", "summary": "The paper is focused on the non-Gaussian mixture model and asymptotical investigation of the asymptotic characterization of the statistics of the empirical risk minimization estimator. The paper takes under consideration the models with two clusters and applies their analysis to the convex loss functions and regularizers. The empirical evaluations investigate the theoretical aspects in practice. ", "questions": "Please refer to the remarks from the Weaknesses section. ", "limitations": "The paper is not organized properly, it is problematic to identify contributions and a strong plot in the work. Empirical evaluation on real cases is missing. ", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 1, "contribution": 3, "strengths": "- The theoretical analysis of generalized linear models is an exciting direction, especially in the contents of high-dimensional data. \n- The empirical evaluation of synthetic use cases seems to confirm the theoretical investigations. \n  ", "weaknesses": "- The paper is chaotic and very difficult to follow. It makes the work very difficult to understand. The listed contusions are not defined to the point. They are focused on studying and analyzing, not on practical outcomes. The work is not even summarized (taking into account that there is still some space in the manuscript, it seems to be strange) in the conclusions section, and limitations are not discussed well. Some explanations of crucial symbols are missing in the paper, and the motivations behind some steps are not explained well.\n\n- The proposed theoretical investigation is limited to the models with two classes. How can the results scale to multiclass scenarios?\n\n- The empirical evaluation is limited only to artificial cases. The problem investigated by the authors is very practical, and it is crucial to provide some empirical evaluation using real datasets. \n\n- There are many ways to go beyond the Gaussian distribution. Normalizing flows may be used to model the distributions for each of the considered clusters as an alternative to this approach. It would be beneficial to discuss this issue in the paper and even provide some empirical comparison to the approach. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the referee for her/his time in reading and evaluating the manuscript. As a general comment, we would like to stress that the goal of our work was to provide, for the first time, a theoretical model to analytically handle the asymptotic properties of classification estimators on non-Gaussian mixtures, in such a way that the effect of non-Gaussianity can be kept *fully* under control and compared with equivalent performances on Gaussian models. \n\nAlthough the framework of the analysis is therefore, at least at this stage, purely theoretical, it provides a number of important insights: it sheds light on the validity of recent Gaussian universality claims going beyond a large literature which was, up to now, limited to the purely Gaussian case as stated in our *Introduction*.  The insights include, but are not limited to, analytical expressions for the generalisation, training errors, and training loss; analytical formulas for the Bayes-optimal performance and data separability threshold (maximum number of samples possible to perfectly interpolate for a given dimension); role of the regularisation strength; validity of Gaussian universality principle on structured and random-labeled datasets. All these results were determined on a very large class of non-Gaussian distributions, with full control over the fatness of tails and even (non-)existence of data covariance.\n* We would like to thank the referee for the feedback: we revisited the manuscript and improved the readability, taking into account all comments. In particular, we have added a *Conclusions and perspectives* section in which the results have been summarised. We have also discussed more clearly the limitations with respect to the treatment of real datasets in that section. There, we highlight that the main difficulty, in this case, is the choice of the best distribution $\\varrho$ given the observed dataset, a problem that however has a long tradition in the context of Bayesian estimation [Alspach and Sorenson, *Nonlinear Bayesian estimation using Gaussian sum approximations*, 1972; Gelman et al., *Bayesian Data Analysis*, 2013]. In a more simplistic approach, it can be observed that, in the case in which the square loss is adopted, the self-consistent equations depend on $\\mathbb E[(1+v\\Delta)^{-1}]$ and $\\mathbb E[(1+v\\Delta)^{-2}]$ only, and these quantities can be numerically estimated from the dataset. The exact evaluation of the quality and limitations of such an approximation on a real dataset are left for future investigation. We would like to stress, however, that taking parameters of distributions themselves as random variables, resulting in superpositions of distributions are readily used in direct applications to describe non-Gaussian data in quantitative finance [Delpini and Bormetti, *Minimal model of financial stylized facts*, 2011; Langrene et al., *Switching to non-affine stochastic volatility: A closed-form expansion for the Inverse Gamma model*, 2015] or econometrics models [Nelson, *ARCH models as diffusion approximations*, 1990], so there already exist schemes for convenient choices in for some types of datasets. \n* The theoretical model can be easily generalised to the case of $K$ classes: in the main text, we limited ourselves to the case of 2 classes for the sake of simplicity and clarity. However, to support this answer, we have modified, in Appendix A, the derivation of our results to include the case of $K$ clusters with scalar labels. Other variations (eg, one-hot-encoding labeling) can be obtained following the formalism in the cited reference [Loureiro et al., 2021].\n* The restriction to synthetic datasets is due to the following reason: the provided theoretical prediction relies on the knowledge of the distribution $\\varrho$ which, in the case of empirical datasets, can be estimated but is in general not known (see also above). One experiment we would like to perform in the future is indeed the comparison of the results of numerical experiments on real datasets with theoretical predictions obtained after an empirical estimation of $\\varrho$, in the spirit of the analysis in [Loureiro et al., 2021] for the case of Gaussian mixtures. On the other hand, as mentioned above, our focus in this contribution was the fact that it is indeed possible to have a simple model within a non-Gaussian setting fully under control and observe therefore a breakdown of Gaussian universality results in it.\n* Although normalising flows are definitely a way to produce non-Gaussian distributions and possibly model a non-Gaussian dataset, as mentioned above the goal of our paper is actually not to provide a tool to model datasets but rather an exactly solvable model, such that the presence of non-Gaussian features can be taken into account and exactly treated. Note that the superposition of the variance distribution with a Gaussian data distribution specifically, while generating a very large family of distributions including ones with any power-law tail, is necessary to employ our analytic method, where Gaussian integration can be conveniently performed by virtue of this data construction starting from Gaussian. We are not aware of references where high-dimensional asymptotics on non-Gaussian distributions is analytically obtained via normalising flows, and we would be extremely grateful to the referee if she/he could point us to some pertinent references, which we will gladly add to the manuscript. \n"}, {"review_id": "mRwKGKcBz7", "reviewer": "Reviewer_X846", "summary": "This paper investigates the asymptotic behavior of Generalized Linear Models (GLM) when the number of training samples $n$, and the dimension of the feature-space $d$ both go to infinity, but the ratio $n/d$ is fixed to some known bounded value $\\alpha$. Moreover, authors assume the training data points are drawn from a mixture of two heavy-tailed distributions, which is different from the usual Gaussian assumption in most of the existing works (the mentioned heavy-tailed distributions are constructed by combining uncountably infinitely many Gaussian distributions).\n\nPaper claims to achieve a non-trivial asymptotic characterization of this problem setting, and also try to validate it via a number of experiments on synthetic data. Paper has a number of shortcomings, therefore my current vote is borderline reject. Presentation of the main results needs to be significantly improved, and also I would like to see the comments from other reviewers with more expertise in this particular field to assess the level of technical contribution in this work.", "questions": "L.130 to L.158: Results are not properly presented. I suggest using a formal theorem and a set of lemmas.\n\nL.139, Eq (4): What are $\\boldsymbol{g}$ and $\\boldsymbol{h}$?\n\nL.140, Eq (5): What are $h_{\\pm}$, $\\omega_{\\pm}$, $q$ and ... Actually this list can go on.\n\nThe main theoretical contributions are presented in Eq (8) and Eq (9). However, the vague explanation preceding them, would impose a huge negative impact on the potential reader.", "limitations": "-", "rating": 5, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- Paper is well-written (at least in most parts), and the literature review part in the introduction section is very informative.\n- I have not completely checked the proofs, however, I have not noticed any mathematical mistakes. The technical validity of the theoretical part looks fine. I have not checked the experimental parts.", "weaknesses": "- All the theoretical derivations are based on asymptotics, while any result in the non-asymptotic case would be far more interesting.\n- I suggest presenting the results more formally, i.e., in the form of Theorems, Lemmas, and etc. Otherwise, the actual level of technical contribution in this work becomes hard to assess. Right now, there are no theorems inside the manuscript. Also, the explanations from L.130 to L.158 are vague (please see the questions section).\n- The main motivation behind this work is to assume non-Gaussian distributions (with a possibly infinite covariance) as the components of the mixture model which generates the data. I am concerned with how much this setting would look imporntant and/or interesting to the community. Due to the Gaussian universality principle, the analysis based on the Gaussian assumptions applies (more or less) to all \"Gaussian-like\" distributions (which covers almost all distributions with bounded moments) as well. Heavy-tailed distributions with power-law tails which do not have a bounded covariance are of course excluded from this list, but how important are they? IMO, authors have not given enough motivation regarding this issue.\n- The process which is used to generate the above-mentioned heavy-tailed distributions is very specific: superposition of an uncountably many Gaussians, or equivalently assuming that the covariance matrix of the Gaussian itself is a R.V. with an inverse-Gamma distribution. Authors have not discussed the limitation of this process. How general is it? does it include almost all heavy-tailed distributions?\n- I have not completely checked the proofs in supp. However, the mathematical tools used for deriving the results are not sophisticated. Not using sophisticated math or not relying on existing elegant theorems is fine, as long as an important problem has been solved or an interesting discovery has been made. This again takes us to a previous comment, on the importance level of this problem setting. I am not familiar with this particular line of research, so I have to wait for other reviewers to comment on that.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the referee for her/his time and comments. Here are some observations concerning the points raised in the report:\n* We agree on the fact that the non-asymptotic behavior is also an interesting problem to consider. In this paper, we have worked in line with a large body of literature that focuses on the asymptotic behavior, however, finite-size *corrections* to the asymptotic behavior are also accessible with the adopted approach (although quite cumbersome) by taking into account corrections to the saddle-point approximation [see, for a review, Lucibello, arXiv:1502.02471]. To our knowledge, this is actually an interesting problem that has not yet been investigated even within the (a priori much simpler) pure Gaussian case. We will keep this suggestion in mind for a future investigation and we thank the referee for raising this point.\n* We have re-organized the material highlighting our main results throughout the text. The set of fixed point equations (7) has been obtained using the heuristic replica method in Appendix A, which we are confident could be made fully rigorous via standard techniques such as Approximate Message Passing (AMP) [Loureiro et al., 2021] or Gordon Minimax techniques [Mignacco et al., 2020]. We will clarify this point in the text with references to the original works. Note that rigorous proof may well require finer technical assumptions. All the following results in our work are instead obtained through fully rigorous derivations from the result in Section 2. We hope that the new form improves the readability of the manuscript.\n* We thank the referee for giving us the opportunity to clarify this important point. The goal of the paper is to “challenge” the validity of a Gaussian universality principle by proposing a model that is under full analytical control and includes the possibility of non-Gaussian datasets. On the one hand, our contribution covers the case of distributions with infinite covariance: such distributions are clearly out of the applicability range of a possible Gaussian universality principle, as correctly stated by the referee. We present this result in this case as a powerful byproduct of our analysis, as, typically, the available theoretical works impose some constraints on the finiteness of the moments, whereas we have no such condition. On the other hand, with respect to the Gaussian universality literature, our most surprising result concerns the case of distributions with *finite* second moment. In the examples of Section 3.1 and Section 3.5, in particular, the existence of moments is controlled via the parameter $a$ of the variance distribution, so that the data distribution has an unbounded k-th moment if  $k\\geq 2a$. Crucially, in such cases Gaussian universality breaks down due to the presence of fat tails, i.e., a Gaussian approximation of the dataset (obtained by matching first and second moment) *does not* reproduce the correct asymptotics. We thus *analytically* show that the performances in terms of generalisation do depend on higher moments. We have further clarified this point in our *Introduction*.\n* We would like to thank the referee for raising this important point. It is true in general that any distribution can be approximated (in the L1 sense) by a possible uncountable superposition of Gaussians, whose means and covariances are given according to some law (see, e.g. [Alspach and Sorenson, IEEE Trans. Autom. Control, 17(4), 439-448, 1972]). We leverage this powerful result by considering the case of a law on the second moment of a Gaussian distribution (as a result, the generated functions are even in $\\boldsymbol x$); the resulting family is large enough to allow the analysis of *any* tail behavior which was a central goal of the paper. We have added a comment with respect to this in the *Introduction*, with new due references. Moreover, the idea of constructing distributions by taking their parameters as random variables themselves appeared in various disciplines, albeit is known under different names: if in statistical physics it is known as superstatistics [Beck, *Recent developments in superstatistics* 2008], there is also a considerable line of work in Bayesian modeling regarding hierarchical priors and models [Gelman and Hill, *Data Analysis Using Regression and Multilevel/Hierarchical Models*, 2006; Gelman et al., *Bayesian Data Analysis*, 2013], while in probability and statistics such distributions are known as compound probability distributions [Robbins, *Asymptotically Subminimax Solutions of Compound Statistical Decision Problems*, 1985], or as doubly-stochastic models in stochastic processes in particular [Pinsky and Karlin, *An Introduction to Stochastic Modeling*, 2010; Schnoerr et al., *Cox process representation and inference for stochastic reaction-diffusion processes*, 2016]. Such superpositions of distributions are also readily used in direct applications to describe non-Gaussian data in quantitative finance [Delpini and Bormetti, *Minimal model of financial stylized facts* 2011; Langrene et al., *Switching to non-affine stochastic volatility: A closed-form expansion for the Inverse Gamma model*, 2015] or econometrics models [Nelson, *ARCH models as diffusion approximations*, 1990].\n* We have carefully re-read and revised the manuscript, and restructured the presentation of our results. The quantity $\\boldsymbol g$ is defined in Eq. 6, right after its appearance in Eq. 4, whilst $\\boldsymbol h$ is defined in line 141 via the quantities $h_\\pm$ and $\\omega_\\pm$ introduced by Eq. 5 (in the new version of our manuscript we have now used the notation $\\coloneqq$ to clarify that we are indeed defining the quantities therein). Similarly, all order parameters are defined as the solutions of Eq. 7. We hope that the new version of the manuscript will be more clear, and we thank the referee for pointing out places where the clarity of the text could be improved. "}, {"review_id": "SpjVvYIRjM", "reviewer": "Reviewer_7S1w", "summary": "The paper is concerned with binary classification, when the data comes from two point clouds that are superposition of Gaussian distribution. This model allows for data distribution with fat tails. The authors analyse the performance of empirical risk minimization in the high-dimensional regime where the number of training samples and the dimension jointly diverge. Using the replica method from statistical physics, they reduce the computation of e.g the training loss / generalization error to the resolution of self-consistent equations.\nIn the third section of the paper, the authors apply their main result on experiments with synthetic data.", "questions": "Could the computations be easily extended to  : \n1) other types of estimators, e.g. Bayesian estimators that sample from a posterior distribution instead of ERM \n2) more generic covariance for the Gaussian distributions. For instance instead of $N(\\mu, \\Delta I_d)$, use $N(\\mu, \\Delta \\Sigma)$ where $\\Sigma is a generic covariance matrix.", "limitations": "n/a", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is clearly written, and the mathematical derivations are easy to follow. The application of the replica method on a mixture of superposition of Gaussian distributions is new, to the best of my knowledge.", "weaknesses": "The paper lacks an experiment on real data to showcase situations in which the data model used in this paper (superposition of Gaussian) is more realistic / useful than simply using a mixture of two Gaussians.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the referee for her/his positive comments on our paper. We summarise our answers to her/his questions below.\n* As foreseen by the referee, the method can indeed be applied to the study of an estimator obtained by minimisation of a proper convex function, and in particular, to the study of the Bayes optimal estimator itself (in this case, the convexity of the problem in the Bayes optimal setting is guaranteed by the so-called \"Nishimori conditions\" in the statistical physics jargon). As the Bayes optimal estimator is indeed an important reference for our results, stimulated by the question of the referee we have derived the expression for the Bayes optimal error in a new Appendix A.4 and added, as guide for comparison, the corresponding error curves in our plots.\n* The case suggested by the referee is indeed within the reach of our theory, which can be easily generalised to consider clouds of the type $P(\\boldsymbol{x}|\\boldsymbol \\mu)=\\mathbb E[\\mathcal N(\\boldsymbol x;\\boldsymbol\\mu,\\Delta\\boldsymbol\\Sigma)]$, possibly with a different covariance matrix for each cloud (the case of identical $\\boldsymbol\\Sigma$ can be reduced to the current setting by a simple change of variables). For example, by supposing that the two clouds are described by a density $P_\\pm(\\boldsymbol x|\\boldsymbol\\mu_\\pm)=\\mathbb E[\\mathcal N(\\boldsymbol x;\\boldsymbol\\mu_\\pm,\\Delta\\boldsymbol\\Sigma)]$, fixed-point equations similar to the ones presented in the paper would be obtained with respect to scalar order parameters, although two order parameters $\\boldsymbol q_\\pm$ would be needed in this case, each taking into account the correlation between the estimator $\\boldsymbol w^\\star$ and each one of the covariances. To avoid this complication, we restricted our analysis to the simplest case for illustrative purposes: we however added a line in the main text commenting about this more general setting pointed out by the referee. We mention here that the case in which the matrix $\\boldsymbol\\Sigma$ is assumed to be random is, instead, much more challenging to analyse: the dependence of the order parameters on the additional stochasticity is, in this case, much more complicated than the one appearing in Eq. 47 and does not allow the simple factorisation given therein. "}]

Author rebuttal:
*General remarks* We thank the reviewers for their helpful feedback which helped us improve the readability and clarity of our work, and better express the importance of our contribution. To take into account their comments, we have prepared a new version of our manuscript, in which, beyond addressing the referees’ comments, we provide additional, important results, namely 
* a generalisation of our result to $K$ classes (to appear in the Appendix);
* a derivation of an analytical formula for the data separability threshold;
* a derivation of an analytical formula for the Bayes-optimal error, and the comparison of the bound with the ERM results (we provide, as an example, the new Fig. 1 where the Bayes-optimal error bound is given as a dashed line for each value $a$). 

We have answered the questions and addressed specific concerns of each reviewer in the individual answers below.

Meta-review: This paper provides statistical rates for classification in high dimensions, with the main contribution being that the data distribution is a superposition of Gaussians, and thus flexible enough to capture many phenomenon heretofore avoided in this literature, such as heavy tails. I am happy to recommend acceptance, though I urge the authors to adjust the manuscript during the camera ready phase to reflect the fruitful reviewer discussions.
Decision: Accept (poster)