# Are Gaussian data all you need? Extents and limits of universality
in high-dimensional generalized linear estimation

Luca Pesce  Ecole Polytechnique Fédérale de Lausanne (EPFL). Information,
Learning and Physics (IdePHICS) lab.
CH-1015 Lausanne, Switzerland. Florent Krzakala  Ecole Polytechnique Fédérale
de Lausanne (EPFL). Information, Learning and Physics (IdePHICS) lab.
CH-1015 Lausanne, Switzerland. Bruno Loureiro  Département d’Informatique,
École Normale Supérieure - PSL & CNRS, 45 rue d’Ulm,
F-75230 Paris cedex 05, France.  Ludovic Stephan  Ecole Polytechnique Fédérale
de Lausanne (EPFL). Information, Learning and Physics (IdePHICS) lab.
CH-1015 Lausanne, Switzerland.

###### Abstract

In this manuscript we consider the problem of generalized linear estimation on
Gaussian mixture data with labels given by a single-index model. Our first
result is a sharp asymptotic expression for the test and training errors in
the high-dimensional regime. Motivated by the recent stream of results on the
Gaussian universality of the test and training errors in generalized linear
estimation, we ask ourselves the question: "_when is a single Gaussian enough
to characterize the error?_ ". Our formula allow us to give sharp answers to
this question, both in the positive and negative directions. More precisely,
we show that the sufficient conditions for Gaussian universality (or lack of
thereof) crucially depend on the alignment between the target weights and the
means and covariances of the mixture clusters, which we precisely quantify. In
the particular case of least-squares interpolation, we prove a strong
universality property of the training error, and show it follows a simple,
closed-form expression. Finally, we apply our results to real datasets,
clarifying some recent discussion in the literature about Gaussian
universality of the errors in this context.

##  1 Introduction

It is commonsense in machine learning that structure in the data is an
important ingredient for successful learning. Quantifying this statement, and
in particular how structure in the features impact the training and
generalization errors the most, is an important endeavor in the broad program
of "seeing through" the modern machine learning black box. On the theoretical
side, there has been some important recent progress in this direction in the
context of generalized linear estimation. For instance, a recent line of work
on linear regression trained on Gaussian data has shown that good
generalization can arise even in the "overparametrized regime" where the
training error is exactly zero [1, 2, 3]. This _benign overfitting_ property
crucially depends on the covariance structure, occurring when the signal
components of the target align with a lower-dimensional of the data, leaving
space for the noise to spread along the higher-dimensional orthogonal subspace
[1]. Analogous conclusions hold, under similar conditions, to generalized
linear tasks [4, 5, 6]. Indeed, this is only one example of many surprising
insights learned from the study of generalized linear models on Gaussian data
over the past few years [7, 8, 9, 10]. Despite the seemingly constraining
assumption on the distribution of the features, a recent line of work provides
strong evidence for the _Gaussian universality_ of the training and
generalization errors in generalized estimation in different settings. These
includes rigorous results for non-Gaussian designs [11, 12], random feature
maps [13, 14, 15, 16], neural tangent features [17], Gaussian mixtures with
random labels [18], as well as extensive numerical evidence for other feature
maps [15, 19] and even real datasets [19, 20, 21]. These works beg the
question "when are Gaussian features a good model for learning?". Our aim is
to give precise answers to this question in the context of generalized linear
estimation on a popular model for multi-modal data, known to be able to
approximate any distribution: the Gaussian mixture model.

###  1.1 Main results

Our main contributions in this work are as follows:

  * •

Exact asymptotics of GLMs: We provide the exact asymptotic limit for the
training and test errors of a generalized linear model with convex loss in
high dimensions, when the data is drawn from a Gaussian mixture model with a
single-index target. These asymptotics are based on the so-called _replica
method_ from statistical physics, and follow the same line as [22], which
considered instead the task of learning the mixture labels.

  * •

Universality of training and test errors: We provide a set of sufficient
conditions on the target weights $\theta_{0}$ such that both the asymptotic
training and test errors for a Gaussian mixture model are independent from the
cluster means. In particular, these conditions are satisfied by a target whose
direction is uniform on $\mathbb{S}^{d-1}$. In the case of ridge regression,
we show an even stronger result: namely, the training loss is also independent
from the cluster covariances, and reduces to that of a single Gaussian with
identity covariance.

  * •

The importance of Homoscedasticity: In the particular case of a homoscedastic
Gaussian mixture (a mixture of Gaussians that share the same covariance
matrix), we further demonstrate universality results that can actually be
observed on real data after a random feature map (see e.g. Fig. 1 and 2). We
also unveil the universal behavior of the linear separability transition, a
phenomenon studied in detail for pure Gaussian data in [10] and that appears
to be universal for a homoscedastic mixture.

  * •

Breaking universality: In contrast to the results of the previous paragraph,
we show that there are two ways to break Gaussian universality. First,
strongly heteroscedasticity can break the universal behavior. Second, in the
homoscedastic case we show that the correlation between the data and the task
matters: even a small correlation between the target weights and the cluster
means suffice to break universality, in the sense that the asymptotic errors
of a model trained in a Gaussian mixture differs from the one of a model
trained on Gaussian data. Rather than the structure of the data itself, what
appears to matter is thus the correlation between this structure and the task
to be learned.

###  1.2 Related works

#### Exact asymptotics:

An appealing feature of Gaussian data is that the asymptotic performance of
different models can be sharply characterized in the proportional high-
dimensional limit where $n,d\to\infty$ at fixed _sample complexity_
$\alpha\coloneqq\nicefrac{{n}}{{d}}$. This is particularly the case for ridge
regression, because of the close connection to a random matrix theory problem;
see e.g. [23, 2, 3]. Beyond the quadratic case, there exists many asymptotic
rigorous studies, for instance [24] studied M-estimators, [25] the performance
of the LASSO estimator, while [19] provided a general result for convex losses
and penalties with arbitrary covariances.

In the case of Gaussian mixtures, most of the effort has been geared towards
classification, i.e. recovering the cluster label instead of teacher-generated
ones. For binary classification, examples include [26, 27, 28], the latter of
which also shows an equivalence between classification and a single-index
model. In the multi-class setting, [29] studied the performance of ridge
regression classifiers; the most general result in this line is [22], which
considers any convex (not necessarily separable) loss.

#### Gaussian universality:

Remarkably, this model is also able to capture the errors of particular
classes of non-Gaussian features. This _Gaussian universality property_ (GEP)
[30] was proven to hold for generalized linear estimation with random features
[13, 14, 15, 16, 31], neural tangent features [17] and kernel features [32,
21, 33, 34, 35]. [36] showed that Gaussian universality is preserved on the
random features model when the weights are trained at order one steps, but
break if an extensive number of steps are taken. Beyond the realm of theorems,
[19] provided numerical evidence of Gaussian universality for a broader class
of realistic features from trained neural networks. On a close line, [18]
studied the Gaussian universality for pure random binary labels, which was an
important source of inspiration for the present work, while [21, 37, 20, 19,
34, 38] has numerically shown that for ridge regression in particular, the
Gaussian formula captured the learning curves of some simple real datasets.

##  2 Setting & motivation

Let $(\bm{x}^{\nu},y^{\nu})\in\mathbb{R}^{d}\times\mathcal{Y}$ denote
$\nu\\!=\\!1,\cdots,n$ pairs of independently sampled training points. We
shall be interested in studying the properties of generalized linear
estimation $\hat{y}(\bm{x})=\hat{f}(\bm{\theta}^{\top}\bm{x})$ with weights
$\bm{\theta}\in\mathbb{R}^{d}$ learned from the training data by minimizing
the following empirical risk:

$\displaystyle\hat{\mathcal{R}}^{\lambda}_{n}(\bm{\theta})=\frac{1}{n}\sum\limits_{\nu=1}^{n}\ell\left(y^{\nu},\frac{\bm{\theta}^{\top}\bm{x}}{\sqrt{d}}\right)+\lambda
r(\bm{\theta})$ (1)

where $\ell:\mathcal{Y}\times\mathbb{R}\to\mathbb{R}_{+}$ is a convex loss
function and $r:\mathbb{R}^{d}\to\mathbb{R}$ is a convex penalty. For example,
this includes the particular case of ridge regression where
$\ell(y,\hat{y})=(y-\hat{y})^{2}$ and
$r(\vec{\theta})=||\vec{\theta}||_{2}^{2}$. The key quantities of interested
in the following will be the training and generalization error, defined as:

$\displaystyle\varepsilon_{\text{tr}}(\hat{\vec{\theta}})$
$\displaystyle=\frac{1}{n}\sum\limits_{\nu=1}^{n}\ell\left(y^{\nu},\frac{\hat{\bm{\theta}}^{\top}\bm{x}^{\nu}}{\sqrt{d}}\right)$
(2) $\displaystyle\varepsilon_{\text{gen}}(\hat{\vec{\theta}})$
$\displaystyle=\mathbb{E}\left[g\left(y_{\text{new}},\hat{f}\left(\frac{\hat{\vec{\theta}}^{\top}\vec{x}_{\text{new}}}{\sqrt{d}}\right)\right)\right]$
(3)

where $g:\mathcal{Y}\times\mathbb{R}\to\mathbb{R}_{+}$ is a performance metric
of the choice of the statistician, not necessarily equal to the loss function
$\ell(y,\hat{y})$. For example, in the case of binary classification with
$\mathcal{Y}=\\{-1,+1\\}$, we can take $\ell(y,\hat{y})=\log(1+e^{-y\hat{y}})$
to be the logistic loss, while taking
$g(y,\hat{y})=\mathbb{I}\left[y\neq\hat{y}\right]$ to be the classification
error. Note that in eq. (3) the expectation is taken over a new data pair
$(\bm{x}_{\text{new}},y_{\text{new}})$ which we assume is independently drawn
from the same distribution of the training data.

In particular, we will be interested in characterizing these errors under the
assumption that the labels have been generated by the following target
distribution:

$\displaystyle y^{\nu}\sim
P_{0}\left(\cdot\,\big{|}\,\frac{\bm{\theta_{0}}^{\top}\bm{x}^{\nu}}{\sqrt{d}}\right)$
(4)

for some fixed vector $\bm{\theta}_{0}\in\mathbb{R}^{d}$ and distribution
$P_{0}$. A common choice for $P_{0}$ is
$y^{\nu}=f_{0}\left(\frac{\bm{\theta_{0}}^{\top}\bm{x}^{\nu}}{\sqrt{d}}+\xi\right)$
for additive Gaussian noise $\xi\sim\mathcal{N}(0,\Delta)$. This setting is
sometimes refereed to as a _teacher-student setting_. We will sometimes adopt
this convenient terminology, refereeing to the target distribution $P_{0}$ as
the _teacher_ and $\bm{\theta}_{0}\in\mathbb{R}^{d}$ as the _teacher weights_.
Similarly, we will sometimes refer to the model $\hat{y}$ as the _student_ and
$\hat{\bm{\theta}}\in\mathbb{R}^{d}$ the _student weights_. As previously
mentioned, we shall be considering both regression $\mathcal{Y}=\mathbb{R}$
and binary classification $\mathcal{Y}=\\{-1,+1\\}$.

This model has been the subject of a plethora of works in the high-dimensional
statistics literature over the past few years, in particular under the
_Gaussian design_ assumption:

######  Model 2.1 (Gaussian covariate model).

In the Gaussian covariate model (GCM), we assume the inputs are independently
drawn from a Gaussian distribution:

$\displaystyle\bm{x}^{\nu}\sim\mathcal{N}\left(\frac{\vec{\mu}}{\sqrt{d}},\Sigma\right),\qquad\nu=1,\cdots,n.$
(5)

We denote $\mathcal{G}_{\vec{\theta}_{0},\Sigma,\vec{\mu}}$ the teacher-
student problem under this data assumption.

A key motivation for this work is the common intuition that data from standard
classification tasks such as MNIST are closer to multi-modal distributions
than to a single-mode Gaussian. Therefore, in this manuscript we ask ourselves
the question: "when is Gaussian data all you need?". In particular, we focus
our attention to a prototypical distribution to model multi-modal data (and a
universal approximator of densities): the $K$ cluster Gaussian mixture:

######  Model 2.2 (Gaussian mixture model).

In the Gaussian mixture model (GMM), we assume the inputs are independently
drawn from a mixture of $K$ Gaussians:

$\displaystyle\bm{x}^{\nu}\sim\sum\limits_{c\in\mathcal{C}}p_{c}\,\mathcal{N}\left(\frac{\bm{\mu}_{c}}{\sqrt{d}},\Sigma_{c}\right),$
$\displaystyle\nu=1,\cdots,n.$ (6)

where $\mathcal{C}\coloneqq\\{1,\cdots,K\\}$ is the set of possible clusters,
and $p_{c}\in[0,1]$ is the probability of belonging to cluster
$c\in\mathcal{C}$, whose means and covariance are given by
$(\vec{\mu}_{c},\Sigma_{c})$.

We note that despite Model 2.1 being a special case $K=1$ of Model 2.2 when
the labels $y^{\nu}$ do not depend on the input cluster, it will be
instructive to treat the $K=1$ and $K>1$ case as two different models.

##  3 Main theoretical results

Figure 1: An illustration of Gaussian universality with vanishing
regularization $\lambda=0^{+}$ for a selection of datasets (MNIST, Fashion-
MNIST, Cifar10) with a random teacher function, after a random feature map:
Generalization (left) and training (right) errors as a function of the number
of samples per dimension $\alpha=n/d$ for ridge regression (top pannels) and
logistic classification (bottom pannels). The solid purple line is the exact
asymptotics formula for Gaussian covariates with identity covariance, while
the vertical green lines are the threshold values $\alpha^{\star}(\Delta)$ at
which an unique minimizer of the loss starts to exist. Dots show numerical
simulations for different real datasets with random features maps. All data
follows the Gaussian predictions for the training loss, illustrating Theorems
3.4 and 3.6. In particular the separability, or interpolation, threshold is
the one of the Gaussian model (see corollary 3.5). Each learning task was ran
for two different value for the noise variance corrupting the labels - crosses
are associated to $\Delta=3$ while pluses to $\Delta=0.3$. Error bars are
built using standard deviation over 30 runs.

In this section, we introduce our main theoretical results concerning
universality of high-dimensional generalized linear estimation of GMMs. Our
result builds on a long line of works providing an exact asymptotic
characterization of empirical risk minimizers (2) on the proportional high-
dimensional limit for Model 2.1 [23, 2, 24, 19]. In particular, closer to our
derivation are the rigorous results in [22, 19]. We prove that the training
and generalization error concentrate in high-dimensions in a deterministic
expression given by the solution of a set of self-consistent equations. Then
we analyze Gaussian universality, provably characterizing a set of sufficient
conditions the learning task must respect such that the test and training
errors of Model 2.2 asymptotically agrees with the ones from Model 2.1.

Since we deal with sequences of random variables, we will need a rigorous
definition of convergence. For two sequences of numbers $(a_{n}),(b_{n})$, we
write

$a_{n}\simeq b_{n}\quad\text{iff}\quad\lim_{n\to\infty}a_{n}-b_{n}=0.$ (7)

Accordingly, for two sequences of random variables $(X_{n}),(Y_{n})$, we
define closeness in probability by

$X_{n}\overset{P}{\simeq}Y_{n}\quad\text{iff}\quad
X_{n}-Y_{n}\overset{P}{\to}0,$ (8)

where $\overset{P}{\to}$ denotes convergence in probability.

###  3.1 Exact asymptotics

Our first result is to give closed-form asymptotic characterization of the
performance of the minimizer of 2 for the Gaussian Mixture model 2.2,
generalizing the rigorous results of [22]:

######  Proposition 3.1.

(Exact asymptotics, informal statement) Consider the empirical risk
minimization problem introduced in eq. (1) under Gaussian mixture data given
by Model 2.2. For any pseudo-Lispchitz performance metric $g$, the training
and generalization errors (2) converge in the high-dimensional limit of
$n,d\to\infty$ with fixed ratio $\alpha=\nicefrac{{n}}{{d}}$ to deterministic
expressions which are entirely determined by the solution of a set of self-
consistent replica saddle-point equations (33).

$\displaystyle\varepsilon_{\text{tr}}(\bm{\hat{\theta}})$
$\displaystyle\overset{P}{\simeq}\;\varepsilon_{\text{tr}}(\bm{\theta}_{0},\\{\bm{\mu}_{c}\\}_{c\in\mathcal{C}},\\{\Sigma_{c}\\}_{c\in\mathcal{C}})$
(9) $\displaystyle\varepsilon_{\text{gen}}(\bm{\hat{\theta}})$
$\displaystyle\overset{P}{\simeq}\;\varepsilon_{\text{gen}}(\bm{\theta}_{0},\\{\bm{\mu}_{c}\\}_{c\in\mathcal{C}},\\{\Sigma_{c}\\}_{c\in\mathcal{C}})$

The key difference between Prop. 3.1 and Thm. 1 from [22] is the distribution
of the labels. While the labels in [22] are given by the GMM cluster index, in
Prop. 3.1 we consider labels generated by the target function (4). While we do
not provide a formal proof of these results, note that the generic proof
scheme from [22], that maps the solution to the study of a so-called
approximate message passing algorithm [39], can be readily adapted to our
setting. Indeed, the approximate message passing scheme in our scenario is the
same, and the only difference in the proof is to include a teacher in the its
asymptotic analysis, similarly as was done in [40].

Once we specified the GMM, and properly defined the ERM in eq. (1) we can
evaluate the expressions in eqs. (2),(3) as a function of low-dimensional
quantities which define the sufficient statistics of the asymptotic errors,
and are also known as order parameters. We give the detail of the computation
in Appendix A.

###  3.2 Uncorrelated teachers

An in-depth comparison of the asymptotic expressions for the errors of Models
2.1 & 2.2 reveals that their key difference lies on the way the leading target
direction $\bm{\theta}_{0}$ correlate with the cluster means and covariances.
A first step towards universality is therefore to characterize under which
conditions the asymptotic errors are independent of the means. We make the
following assumptions:

######  Assumption 1.

The teacher $\vec{\theta}_{0}$ respects,
$\forall(c,c^{\prime})\in\mathcal{C}\times\mathcal{C}$:

$\displaystyle\lim_{n,d\to\infty}\frac{\vec{\theta}_{0}^{\top}\vec{\mu}_{c}}{d}=0$
(10)
$\displaystyle\lim_{n,d\to\infty}\frac{1}{d}\vec{\theta}_{0}^{\top}\Sigma_{c^{\prime}}\left(\lambda+\sum_{c\in\mathcal{C}}\hat{V}^{\star}_{c}\Sigma_{c}\right)^{-1}\bm{\mu}_{c}\to
0$ (11)

where $\\{\hat{V}^{\star}_{c}\\}_{c=1}^{K}$ are the fixed points of the
(replica) saddle point equations describing the centered GMM problem.

######  Assumption 2.

The loss function, and the teacher distribution are both symmetric:

$\displaystyle\ell(x,y)$ $\displaystyle=\ell(-x,-y)$ (12) $\displaystyle
P_{0}\left(y|\tau\right)$ $\displaystyle=P_{0}\left(-y|-\tau\right),$ (13)

and the regularization is an $\ell^{2}$ penalty
$\nicefrac{{\lambda}}{{2}}\lVert\cdot\rVert_{2}^{2}$.

######  Proposition 3.2.

(Mean Universality)

Under Assumptions 1 and 2,the cluster means
$\\{\vec{\mu}_{c}\\}_{c\in\mathcal{C}}$ are not relevant in high-dimensional
ERM estimation:

$\displaystyle\varepsilon^{\rm
GMM}_{gen}\left(\\{\vec{\mu}_{c}\\}_{c=1}^{K},\\{\Sigma_{c}\\}_{c=1}^{K}\right)$
$\displaystyle\simeq\varepsilon^{\rm
GMM}_{gen}\left(\vec{0},\\{\Sigma_{c}\\}_{c=1}^{K}\right)$ (14)
$\displaystyle\varepsilon^{\rm
GMM}_{tr}\left(\\{\vec{\mu}_{c}\\}_{c=1}^{K},\\{\Sigma_{c}\\}_{c=1}^{K}\right)$
$\displaystyle\simeq\varepsilon^{\rm
GMM}_{tr}\left(\vec{0},\\{\Sigma_{c}\\}_{c=1}^{K}\right)$ (15)

Therefore, intuitively mean universality can be achieved when the label
generation process is uncorrelated with the data structure, see Appendix B for
a detailed discussion. The assumptions on the target and loss function are not
restrictive, and are easily satisfied by odd target activation $f_{0}$ and
margin-based losses of the form $l(y,z)=\ell(yz)$. Stronger universality can
be shown by doing simplifying assumption on the data structure. Indeed, if the
mixture is homogeneous (in the sense that the all the covariance are
identical, a condition often called homoscedasticity in statistics) we have
Gaussian universality:

######  Theorem 3.3.

(Gaussian Universality of homoscedastic GMMs) Under the assumption of 3.2,
consider an homoscedastic GMM:

$\Sigma_{c}=\Sigma\qquad\forall c\in\mathcal{C}$

Then for all $\alpha$ and $\Delta$, the errors of the GMM are asymptotically
equal to those of a GCM:

$\displaystyle\varepsilon^{\rm
GMM}_{\text{gen}}\left(\\{\vec{\mu}_{c}\\}_{c=1}^{K},\Sigma\right)$
$\displaystyle\simeq\varepsilon^{\rm
GCM}_{\text{gen}}(\vec{\theta}_{0},\Sigma,\vec{0})$ (16)
$\displaystyle\varepsilon^{\rm
GMM}_{\text{tr}}\left(\\{\vec{\mu}_{c}\\}_{c=1}^{K},\Sigma\right)$
$\displaystyle\simeq\varepsilon^{\rm
GCM}_{\text{tr}}(\vec{\theta}_{0},\Sigma,\vec{0})$ (17)

This results follow in a straightforward way from Theorem. 3.2. We have mapped
under controllable assumptions a GMM problem to a simpler Gaussian one.

Additionally, if we consider vanishing regularization, we can prove that
errors are independent from the shape of the covariance, and therefore we can
take an isotropic mixture:

######  Theorem 3.4.

(Covariance universality for $\lambda=0^{+}$) Under the assumption of Thm.
3.3, assume that it exists an unique minimizer of the empirical risk (1) with
zero regularization. Then, the test & training errors a homoscedastic GMM
(2.2) estimation problem asymptotically coincide with those of a centered,
isotropic Gaussian model $\mathcal{G}_{(\vec{\theta}_{0},\mathbf{I},\vec{0})}$
(2.1).

The proof of this general result is given in Appendix B.

Since the training error for losses such as logistic or hinge characterize the
separability transition at $\lambda=0^{+}$, an interesting consequence of Thm.
3.4 is the following:

######  Corollary 3.5.

(Universality of linear separability for homoscedastic GMMs) The location of
the separability / interpolation transition $\alpha_{c}$ above which the data
stop to be linearly separable is the same for homoscedastic GMMs and the
Gaussian model.

This universality is particularly interesting in light of the detailed study
of the separability transition for random teacher weights and Gaussian data in
[10]. Similar universality phenomena was observed for the reconstruction
transition in linear estimation in [41].

Surprisingly, if we further consider a square loss minimization, the
estimation of any GMM (homoscedastic or not!) under mean universality
condition can be mapped to those of a trivial Gaussian problem:

######  Theorem 3.6.

(Strong universality of the square loss for $\lambda=0^{+}$) Consider
underparametrized learning of GMMs with general means and covariance
respecting the assumptions of Theorem. 3.2. Set in eq. (1):
$l(y,\hat{y})=(y-\hat{y})^{2}$ and consider a generic noisy linear teacher
activation $P_{0}(\tau;\Delta)=\mathcal{N}(\tau,\Delta)$, with
$\xi\sim\mathcal{N}(0,1)$. Then, the training error for a GMM estimation
problem is given by

$\displaystyle\varepsilon^{\rm
GMM}_{tr}\left(\\{\vec{\mu}_{c}\\}_{c=1}^{K},\\{\Sigma_{c}\\}_{c=1}^{K}\right)$
$\displaystyle\simeq\frac{(\alpha-1)\Delta}{\alpha}$ (18)

Note that the strong universality statement does not hold for the test error
(a counterexample is discussed in Appendix C using a strongly heteroscedastic
case). However, as we see from Fig. 1, it seems surprisingly true that random
teacher regression on real data follows the Gaussian asymptotic prediction in
the underparametrized region. Moreover, we observe in the lower panel of Fig.
1 that even for non-quadratic losses we can draw a similar conclusion, and we
investigate this further in next section.

Finally, we note that it the limit of infinite noise and binary labels, our
results give back the ones observed for purely random Rachemacher labels
proven in [18].

Figure 2: An illustration of Gaussian universality with finite regularization
$\lambda$ for a selection of datasets (MNIST, Fashion-MNIST, Cifar10) with a
random teacher function, after a random feature map: Generalization (left) and
training (right) errors as a function of the number of samples per dimension
$\alpha=n/d$. In the upper panel we show ridge regression on Cifar10
preprocessed with RF and $\tanh$ activation, while the lower one is logistic
regression on MNIST with RF and $\operatorname{erf}$ activation. The dashed
curves are the exact asymtotics prediction of the Gaussian theory by matching
covariance. The good agreement illustrates the property of theorem 3.3. The
different colours represents different regularization strength
$\lambda\in\\{0.1,1,10\\}$ respectively in yellow, blue, and red. Error bars
are built over $30$ runs.

###  3.3 Correlated teachers

We now consider the general case where the target weights correlate with the
structure in the data, relaxing the assumptions in Theorem. 3.2. Can we still
say something? Our first result shows that the answer is yes. We focus on a
simple controlled setting in which we can express the ERM performance in a
closed form for any teacher vector:

######  Theorem 3.7.

(Exact asymptotics for isotropic covariance) Consider a ridge regression task
with a 2-clusters GMM (2.2). Note that, without loss of generality we can take
$\vec{\mu}_{+}=-\vec{\mu}_{-}=\vec{\mu}$. Assume isotropic covariances:

$\displaystyle\qquad\Sigma_{+}=\Sigma_{-}=I_{d},$ (19)

and denote:

$\displaystyle\rho$
$\displaystyle=\lim\limits_{d\to\infty}\frac{1}{d}\vec{\theta}_{0}^{\top}\Sigma\vec{\theta}_{0},\qquad\gamma=\lim\limits_{d\to\infty}\frac{1}{d}||\vec{\mu}||_{2}^{2}$
(20) $\displaystyle\pi$
$\displaystyle=\lim\limits_{d\to\infty}\frac{1}{d}\vec{\mu}^{\top}\vec{\theta}_{0}\qquad
P_{0}(\tau;\Delta)=\mathcal{N}(\tau,\Delta).$ (21)

Defining:

$\displaystyle A(\eta)$
$\displaystyle=\frac{(\eta-1)^{2}\left(\gamma\eta^{2}+2\eta-1\right)}{\left(\Delta+(\eta-1)^{2}\rho\right)(1+\gamma\eta)^{2}}\,\,\,\rm{with}$
(22) $\displaystyle\eta$
$\displaystyle=1+\frac{1}{2}\left(\alpha-1+\lambda-\sqrt{4\lambda+\left(\alpha-1+\lambda\right)^{2}}\right)$

The asymptotic errors admit a closed form expression in terms of $\eta$:

$\displaystyle\varepsilon^{\rm GMM}_{\text{gen}}$ $\displaystyle\simeq
g(\alpha,\Delta,\rho,\eta)(1-\pi^{2}A(\eta))$ (23)
$\displaystyle\varepsilon^{\rm GMM}_{\text{tr}}$ $\displaystyle\simeq
t(\alpha,\Delta,\rho,\eta)(1-\pi^{2}A(\eta))$ (24)

where $g$ and $t$ are the asymptotic limit of the generalisation
$\varepsilon^{\rm GCM}_{\text{gen}}(\vec{\theta}_{0},I,\vec{0})$ and training
$\varepsilon^{\rm GCM}_{\text{tr}}(\vec{\theta}_{0},I,\vec{0})$ error for a
single Gaussian model with unit covariance and uncorrelated teacher:

$\displaystyle g(\alpha,\Delta,\rho,\eta)$
$\displaystyle=\frac{\alpha\left(\Delta+(\eta-1)^{2}\rho\right)}{\left(\alpha-\eta^{2}\right)}$
(25) $\displaystyle t(\alpha,\Delta,\rho,\eta)$
$\displaystyle=\frac{(\alpha-\eta)^{2}\left(\Delta+(\eta-1)^{2}\rho\right)}{\alpha\left(\alpha-\eta^{2}\right)}$
(26)

Figure 3: An illustration of Gaussian universality for Training, and lack
thereof for Generalization at vanishing regularization. Instead of a random
feature maps as in Fig. 1, we pre-processed grayscale Cifar10 with wavelet
scattering (blue dots) or orthogonal hadamard projections (yellow dots), or
used directly the raw data (red dots). Using again a random teacher function,
the training error is found to agree perfectly with the universal Gaussian
asymptotics predictions with identity covariance and $\lambda=0^{+}$, as for
theorem 3.6, but the generalization is found to show clear deviation. This
highlight the roles of heteroscedasticity in the data. Error bars are built
over $30$ runs.

The full closed form expressions and the extension to covariances of the form
$\Sigma=\sigma I_{d}$ are derived in Appendix C. Note that the correction
factor to the Gaussian performance scales with $\pi^{2}$, a measure of
correlation between teacher vector and data structure. Hence, we would be
tempted to state that targets that correlate with the data structure, i.e.
$\pi\neq 0$, always break Gaussian universality. Although this is usually the
case, in the limit of vanishing regularization we are in the position to
present an interesting corollary of Theorem. 3.7:

######  Corollary 3.8.

(Restoration of universality for $\lambda=0^{+}$) Consider the same setting of
Theorem. 3.7, further consider underparametrized learning in the limit of
vanishing regularization. Then, the test and training errors for a GMM
estimation are equal to a Gaussian one for any teacher vector (that is eqs
(23) and (24) at $\eta=0$):

$\displaystyle\varepsilon^{\rm GMM}_{\text{gen}}$
$\displaystyle=\Delta\frac{\alpha}{\alpha-1}$ (27)
$\displaystyle\varepsilon^{\rm GMM}_{\text{tr}}$
$\displaystyle=\Delta\frac{(\alpha-1)}{\alpha}$ (28)

We thus restore Gaussian universality for correlated teachers in the
underparametrized regime, in a similar fashion to what Theorem. 3.4 is stating
for general convex losses and covariances under the mean universality
condition. The universality property for correlated teachers is valid also for
more general homoscedastic mixtures with identity covariance. For the sake of
brevity we refer to Appendix C where we discuss in detail this interesting
extension.

##  4 Illustration on synthetic and real datasets

In this section we investigate the consequences of our main theoretical
results. First, we consider the case of real data, illustrating the
applicability of our universality theorems in cases in which the target is not
correlated with the data structure. Based on these observations we move on
studying simple synthetic settings described in Theorem 3.7 and for which we
can derive analytical results and systematically probe Gaussian universality.
All the code used in our experiments are available in a GitHub repository.

###  4.1 Random teacher universality

We analyze Gaussian universality of real data under the random teacher
function assumptions of Prop. 3.2. We consider three standard datasets: MNIST,
Fashion-MNIST, and grayscale CIFAR-10, as well as synthetic data.

Figure 4: On the importance of the correlation of the task and the structure
of the data: Left & Center : Training errors achieved in the ERM for
estimation of 2-cluster GMM with $\vec{\mu}_{+}=-\vec{\mu}_{-}$ and
$\Sigma_{+}=\Sigma_{-}=I_{d}$, plotted versus a correlation parameter
$\Omega$: we build a series of learning tasks with teacher weights dependent
on $\Omega$, namely we take
$\vec{\theta}_{0}(\Omega)=\Omega\vec{\mu}_{\perp}+\sqrt{1-\Omega^{2}}\vec{\mu}_{+}$,
with $\vec{\mu}_{+}^{\top}\vec{\mu}_{\perp}=0$. We fix the dimension to be
$d=500$, and $(\alpha,\lambda)=(1.2,0.7)$. The solid black line is the
theoretical prediction coming from Gaussian asymptotics. The solid red line is
the theoretical prediction for the GMM performance, while the orange dots are
the numerical simulations which agree as expected with the theoretical
prediction. In the left panel we have real labels and perform ridge
regression, while on the central one we consider binary labels and perform
logistic regression. The error bars are built using standard deviation over
$30$ runs. Right: Geometrical intuition plotted for $d=2$. Three hyperplanes
(lines in 2D) are displayed correspondent to $\Omega=\\{0,0.5,1\\}$,
respectively in green, blue, and orange. As $\Omega$ increase the labels
become more uncorrelated with the data structure.

#### Universality with random feature maps —

To go beyond simple linear fit, we pre-process the data with a random feature
map [42] as follows: take an image
$\vec{\omega}_{\nu}\in\mathbb{R}^{d^{\prime}}$ and map it to
$\vec{x}_{\nu}=\sigma(F\vec{\omega}_{\nu})$, where the elements
$F\in\mathbb{R}^{d\times d^{\prime}}$ are i.i.d drawn from $\mathcal{N}(0,1)$
and $\sigma(\cdot)$ is an activation function. As shown in [42], for
$d\to\infty$ this converge to a kernel method. For each dataset we take a
different non-linearity ($\operatorname{erf}$,$\tanh$, and $\rm{sign}$ as
described in Fig. 1). We then build new labels by plugging in eq. (4)
$\vec{\theta}_{0}\sim\mathcal{N}(\vec{0},I_{d})$, and we set for regression
tasks $f_{0}(\tau)=\tau+\Delta$, while for classification tasks
$f_{0}(\tau)=\rm{sign}(\tau+\Delta)$.

First, we analyze the vanishing regularization case in Fig. 1: random teacher
regression on pre-processed real datasets respects the strong universality of
the training loss as stated in Theorem. 3.6. Interestingly, we also observe a
perfect match between the Gaussian asymptotics prediction and the simulations
in Fig. 1 for the generalization error. This suggests that after the random
features maps the data is sufficiently close to a homoscedastic mixture.
Further, the lower panel of Fig. 1 shows that the strong universality property
seems to hold beyond square loss minimization.

We analyze as well the finite regularization setting in Fig. 2: we compare the
simulations on real data with the prediction of the exact Gaussian
asymptotics: we match the covariance of each pre-processed real dataset and
compute the performance of the ERM estimator thanks to the deterministic
replica formula. 3.1. The predictions of the Gaussian theory matches the
numerical simulations. As discussed for Fig. 1, it seems that the
homoscedasticity assumption in Theorem. 3.3 can be sometimes relaxed.

#### Universality of the double descent phenomena —

One finding in modern machine learning that goes against the classical
statistical theory wisdom is the double descent behaviour of learning curves
[43, 44, 45, 13], see upper panel of Fig. 1: the test error does not
deteriorate as the number of parameters is increased with a characteristic
divergence at $\alpha=1$, known as the interpolation peak. As shown in Fig. 1
and Fig. 2, the generalization error and its characteristic "double descent"
behavior for $\alpha>1$ is universal for homoscedastic data, while the
training error appears universal even for heteroscedastic data (see Fig. 3).

#### Universality of linear separability —

Recently, [10] investigated the linear separability of Gaussian data with a
random teacher and noise $\Delta$, generalizing the classical results by [46]
to a single-index target. [10] has shown the existence of a critical phase
transition $\alpha_{c}(\Delta)$ that goes continuously from
$\alpha_{c}(\infty)=2$ (for infinite noise) to $\alpha_{c}>2$ for finite
$\Delta$ (with $\alpha_{c}\to\infty$ as $\Delta\to 0$). As discussed in
corollary 3.5, this transition is universal for homoscedastic mixtures. The
interest of this transition thus extends way beyond Gaussian data and in fact
agrees with the real data experiment in Figs. 1 and 2. As for the double
descent phenomena, we believe it is a very interesting consequence of our
theorems that the theoretical works mentioned above are valid way beyond the
simple Gaussian assumption.

#### Non-universal behavior for generalization —

Finally, we implemented different transforms beyond random features for the
pre-processing step: we considered the wavelet scattering transform [47],
orthogonal random projections [48], and even no transform at all. Without the
shuffling of the random projection, the fact that the data are complex and
probably rather heteroscedastic than homoscedastic, we expected a weaker form
of universality. We present the results of ridge regression with labels
generated by a random teacher vector in Fig. 3. The strong universality
statement in Theorem. 3.6, that did not required any assumption on the data,
remains valid and we observed a perfect collapse of the training data.
However, as expected, the generalization error shows clear deviations with
respect to the Gaussian behavior. Presumably these transforms do not
homogenize the data enough such that Gaussian universality for the test error
to hold.

###  4.2 Correlated teacher

Previously, we showed that ridge regression with a random teacher model, even
when the data is structured, can lead to universal behaviour. We now consider
the dual task: take a simple homogeneous model for the data and study
correlated target weights. Indeed, some recent works have studied examples of
the lack of universality between Gaussian mixtures and Gaussian models [49,
50]; we want to follow this direction and use the setting described in
Theorem. 3.7: we consider a 2-cluster GMM with opposite means of norm
$\vec{\mu}$ and same covariance $\Sigma$. We build a series of learning tasks
at fixed $(\alpha,\lambda)$ varying the overlap between the teacher vector and
the cluster means as follows:

$\displaystyle\vec{\theta}_{0}(\Omega)=\Omega\vec{\mu}_{\perp}+\sqrt{1-\Omega^{2}}\vec{\mu}$
(29) $\displaystyle\text{such that}\qquad\vec{\mu}^{\top}\vec{\mu}_{\perp}=0$
(30)

The results are presented in Fig. 4 for both ridge and logistic regression.
They clearly show that a small correlation between the target weights and the
mixture means breaks Gaussian Universality. A neat geometrical intuition of
the result is given in Fig. 4: as we decrease $\Omega$, we generate labels
which are more and more correlated with the data distribution, and
consequently there is a correction factor to the Gaussian prediction as
Theorem. 3.7 predicts. We refer to Appendix. C for a more detailed discussion
on correlated teachers. We conclude that universality is broken in tasks where
the labels are correlated with the structure. Note, however, that as proven in
Corollary. 3.8: the discrepancy between the GMM prediction and the Gaussian
one goes to zero for any correlation measure $\Omega$ as $\lambda\to 0^{+}$,
where the universality is restored as suggested by Thm. 3.6.

##  5 Acknowledgments

We thank Yatin Dandi, Federica Gerace, Sebastian Goldt, Denny Wu & Lenka
Zdeborová for useful discussions. We acknowledge funding from the Swiss
National Science Foundation grant SNFS OperaGOST, $200021\\_200390$ and the
Choose France - CNRS AI Rising Talents program.

## References

  * [1] Peter L. Bartlett, Philip M. Long, Gá bor Lugosi, and Alexander Tsigler.  Benign overfitting in linear regression.  Proceedings of the National Academy of Sciences, 117(48):30063–30070, apr 2020.
  * [2] Trevor Hastie, Andrea Montanari, Saharon Rosset, and Ryan J. Tibshirani.  Surprises in high-dimensional ridgeless least squares interpolation.  Annals of Statistics, 50(2):949–986, April 2022.
  * [3] Denny Wu and Ji Xu.  On the optimal weighted \ell_2 regularization in overparameterized linear regression.  In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 10112–10123. Curran Associates, Inc., 2020.
  * [4] Vidya Muthukumar, Adhyyan Narang, Vignesh Subramanian, Mikhail Belkin, Daniel Hsu, and Anant Sahai.  Classification vs regression in overparameterized regimes: Does the loss function matter?  Journal of Machine Learning Research, 22(222):1–69, 2021.
  * [5] Ke Wang and Christos Thrampoulidis.  Benign overfitting in binary classification of gaussian mixtures.  In ICASSP 2021 - 2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 4030–4034, 2021.
  * [6] Ohad Shamir.  The implicit bias of benign overfitting.  In Po-Ling Loh and Maxim Raginsky, editors, Proceedings of Thirty Fifth Conference on Learning Theory, volume 178 of Proceedings of Machine Learning Research, pages 448–478. PMLR, 02–05 Jul 2022.
  * [7] Elizabeth Gardner and Bernard Derrida.  Optimal storage properties of neural network models.  Journal of Physics A: Mathematical and general, 21(1):271, 1988\.
  * [8] Anders Krogh and John A Hertz.  A simple weight decay can improve generalization.  In Advances in neural information processing systems, pages 950–957, 1992.
  * [9] David Donoho and Jared Tanner.  Observed universality of phase transitions in high-dimensional geometry, with implications for modern data analysis and signal processing.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 367(1906):4273–4293, 2009.
  * [10] Emmanuel J Candès, Pragya Sur, et al.  The phase transition for the existence of the maximum likelihood estimate in high-dimensional logistic regression.  The Annals of Statistics, 48(1):27–42, 2020.
  * [11] Andrea Montanari and Phan-Minh Nguyen.  Universality of the elastic net error.  In 2017 IEEE International Symposium on Information Theory (ISIT), pages 2338–2342, 2017.
  * [12] Ashkan Panahi and Babak Hassibi.  A universal analysis of large-scale regularized least squares solutions.  In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
  * [13] Song Mei and Andrea Montanari.  The generalization error of random features regression: Precise asymptotics and double descent curve, 2019.
  * [14] Federica Gerace, Bruno Loureiro, Florent Krzakala, Marc Mezard, and Lenka Zdeborova.  Generalisation error in learning with random features and the hidden manifold model.  In Hal Daumé III and Aarti Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 3452–3462. PMLR, 13–18 Jul 2020\.
  * [15] Sebastian Goldt, Bruno Loureiro, Galen Reeves, Florent Krzakala, Marc Mézard, and Lenka Zdeborová.  The gaussian equivalence of generative models for learning with shallow neural networks.  In MSML, 2021.
  * [16] Hong Hu and Yue M. Lu.  Universality laws for high-dimensional learning with random features.  IEEE Transactions on Information Theory, pages 1–1, 2022.
  * [17] Andrea Montanari and Basil N. Saeed.  Universality of empirical risk minimization.  In Po-Ling Loh and Maxim Raginsky, editors, Proceedings of Thirty Fifth Conference on Learning Theory, volume 178 of Proceedings of Machine Learning Research, pages 4310–4312. PMLR, 02–05 Jul 2022.
  * [18] Federica Gerace, Florent Krzakala, Bruno Loureiro, Ludovic Stephan, and Lenka Zdeborová.  Gaussian universality of linear classifiers with random labels in high-dimension, 2022.
  * [19] Bruno Loureiro, Cedric Gerbelot, Hugo Cui, Sebastian Goldt, Florent Krzakala, Marc Mezard, and Lenka Zdeborová.  Learning curves of generic features maps for realistic datasets with a teacher-student model.  In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan, editors, Advances in Neural Information Processing Systems, volume 34, pages 18137–18151. Curran Associates, Inc., 2021.
  * [20] Arthur Jacot, Berfin Simsek, Francesco Spadaro, Clement Hongler, and Franck Gabriel.  Kernel alignment risk estimator: Risk prediction from training data.  In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 15568–15578. Curran Associates, Inc., 2020.
  * [21] Blake Bordelon, Abdulkadir Canatar, and Cengiz Pehlevan.  Spectrum dependent learning curves in kernel regression and wide neural networks.  2020\.
  * [22] Bruno Loureiro, Gabriele Sicuro, Cédric Gerbelot, Alessandro Pacco, Florent Krzakala, and Lenka Zdeborová.  Learning gaussian mixtures with generalised linear models: Precise asymptotics in high-dimensions.  2021\.
  * [23] Edgar Dobriban and Stefan Wager.  High-Dimensional Asymptotics of Prediction: Ridge Regression and Classification.  The Annals of Statistics, 46(1):247–279, 2018.
  * [24] Christos Thrampoulidis, Ehsan Abbasi, and Babak Hassibi.  Precise Error Analysis of Regularized $M$ -Estimators in High Dimensions.  IEEE Transactions on Information Theory, 64(8):5592–5628, August 2018.
  * [25] Michael Celentano, Andrea Montanari, and Yuting Wei.  The lasso with general gaussian designs with applications to hypothesis testing, 2020.
  * [26] Xiaoyi Mai, Zhenyu Liao, and Romain Couillet.  A Large Scale Analysis of Logistic Regression: Asymptotic Performance and New Insights.  In ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 3357–3361, May 2019.  ISSN: 2379-190X.
  * [27] Francesca Mignacco, Florent Krzakala, Yue M. Lu, and Lenka Zdeborová.  The role of regularization in classification of high-dimensional noisy gaussian mixture.  2020\.
  * [28] Zeyu Deng, Abla Kammoun, and Christos Thrampoulidis.  A model of double descent for high-dimensional binary linear classification.  Information and Inference: A Journal of the IMA, 11(2):435–495, June 2022.
  * [29] Christos Thrampoulidis, Samet Oymak, and Mahdi Soltanolkotabi.  Theoretical insights into multiclass classification: a high-dimensional asymptotic view.  In Proceedings of the 34th International Conference on Neural Information Processing Systems, NIPS’20, pages 8907–8920, Red Hook, NY, USA, December 2020. Curran Associates Inc.
  * [30] Sebastian Goldt, Marc Mézard, Florent Krzakala, and Lenka Zdeborová.  Modeling the influence of data structure on learning in neural networks: The hidden manifold model.  Phys. Rev. X, 10:041044, Dec 2020.
  * [31] Dominik Schröder, Hugo Cui, Daniil Dmitriev, and Bruno Loureiro.  Deterministic equivalent and error universality of deep random features learning, 2023.
  * [32] Noureddine El Karoui.  Spectrum estimation for large dimensional covariance matrices using random matrix theory.  The Annals of Statistics, pages 2757–2790, 2008.
  * [33] Song Mei, Theodor Misiakiewicz, and Andrea Montanari.  Generalization error of random feature and kernel methods: Hypercontractivity and kernel matrix concentration.  Applied and Computational Harmonic Analysis, 59:3–84, 2022.  Special Issue on Harmonic Analysis and Machine Learning.
  * [34] Hugo Cui, Bruno Loureiro, Florent Krzakala, and Lenka Zdeborová.  Generalization error rates in kernel regression: The crossover from the noiseless to noisy regime.  In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan, editors, Advances in Neural Information Processing Systems, volume 34, pages 10131–10143. Curran Associates, Inc., 2021.
  * [35] Hugo Cui, Bruno Loureiro, Florent Krzakala, and Lenka Zdeborová.  Error rates for kernel classification under source and capacity conditions, 2022.
  * [36] Jimmy Ba, Murat A. Erdogdu, Taiji Suzuki, Zhichao Wang, Denny Wu, and Greg Yang.  High-dimensional asymptotics of feature learning: How one gradient step improves the representation, 2022.
  * [37] Stefano Spigler, Mario Geiger, and Matthieu Wyart.  Asymptotic learning curves of kernel methods: empirical data versus teacher–student paradigm.  Journal of Statistical Mechanics: Theory and Experiment, 2020(12):124001, 2020.
  * [38] Alexander Wei, Wei Hu, and Jacob Steinhardt.  More than a toy: Random matrix models predict how real-world neural representations generalize.  In Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato, editors, Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 23549–23588. PMLR, 17–23 Jul 2022.
  * [39] David L Donoho, Arian Maleki, and Andrea Montanari.  Message-passing algorithms for compressed sensing.  Proceedings of the National Academy of Sciences, 106(45):18914–18919, 2009.
  * [40] Elisabetta Cornacchia, Francesca Mignacco, Rodrigo Veiga, Cédric Gerbelot, Bruno Loureiro, and Lenka Zdeborová.  Learning curves for the multi-class teacher-student perceptron.  Machine Learning: Science and Technology, 2022.
  * [41] Ehsan Abbasi, Fariborz Salehi, and Babak Hassibi.  Universality in learning from linear measurements.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
  * [42] Ali Rahimi and Benjamin Recht.  Random Features for Large-Scale Kernel Machines.  In Advances in Neural Information Processing Systems, volume 20. Curran Associates, Inc., 2007.
  * [43] Manfred Opper and Wolfgang Kinzel.  Statistical mechanics of generalization.  In Models of neural networks III, pages 151–209. Springer, 1996\.
  * [44] Mikhail Belkin, Daniel Hsu, Siyuan Ma, and Soumik Mandal.  Reconciling modern machine-learning practice and the classical bias–variance trade-off.  Proceedings of the National Academy of Sciences, 116(32):15849–15854, jul 2019.
  * [45] S Spigler, M Geiger, S d’Ascoli, L Sagun, G Biroli, and M Wyart.  A jamming transition from under- to over-parametrization affects generalization in deep learning.  Journal of Physics A: Mathematical and Theoretical, 52(47):474001, oct 2019.
  * [46] Thomas M Cover.  Geometrical and statistical properties of systems of linear inequalities with applications in pattern recognition.  IEEE transactions on electronic computers, (3):326–334, 1965.
  * [47] Mathieu Andreux, Tomás Angles, Georgios Exarchakis, Roberto Leonarduzzi, Gaspar Rochette, Louis Thiry, John Zarka, Stéphane Mallat, Joakim Andén, Eugene Belilovsky, Joan Bruna, Vincent Lostanlen, Muawiz Chaudhary, Matthew J. Hirn, Edouard Oyallon, Sixin Zhang, Carmine Cella, and Michael Eickenberg.  Kymatio: Scattering transforms in python.  Journal of Machine Learning Research, 21(60):1–6, 2020.
  * [48] Krzysztof M Choromanski, Mark Rowland, and Adrian Weller.  The unreasonable effectiveness of structured random orthogonal embeddings.  In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
  * [49] Umberto M. Tomasini, Antonio Sclocchi, and Matthieu Wyart.  Failure and success of the spectral bias prediction for Laplace Kernel Ridge Regression: the case of low-dimensional data.  pages 21548–21583. PMLR, June 2022.
  * [50] Alessandro Ingrosso and Sebastian Goldt.  Data-driven emergence of convolutional structure in neural networks.  Proceedings of the National Academy of Science, 119(40):e2201854119, October 2022.
  * [51] Marc Mezard, Giorgio Parisi, and Miguel Angel Virasoro.  Spin Glass Theory And Beyond: An Introduction To The Replica Method And Its Applications.  World Scientific Publishing Company, November 1987.  Google-Books-ID: DwY8DQAAQBAJ.
  * [52] Michel Talagrand.  What Is a Quantum Field Theory? Cambridge University Press, Cambridge, 2022.
  * [53] 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.  Journal of Machine Learning Research, 12:2825–2830, 2011.

##  Appendix A Replica computation

###  A.1 Formal statement of the theorem

We first provide the full statement of Proposition 3.1. Consider a
minimization problem of the form

$\bm{\hat{\theta}}=\operatorname*{arg\,min}_{\bm{\theta\in\mathbb{R}^{d}}}\frac{1}{n}\sum\limits_{\nu=1}^{n}\ell\left(y^{\nu},\bm{\theta}^{\top}\bm{x}^{\nu}\right)+r(\bm{\theta}),$
(31)

where the data $(\bm{x}^{\nu},y^{\nu})$ is generated according to the
following Gaussian mixture model:

$\bm{x}^{\nu}\sim\sum_{c\in\mathcal{C}}p_{c}\,\mathcal{N}\left(\frac{\bm{\mu}_{c}}{\sqrt{d}},\Sigma_{c}\right)\quad\text{and}\quad
y^{\nu}\sim P_{0}(\cdot\,|\,\bm{\theta}_{0}^{\top}\bm{x}^{\nu})$ (32)

and assume the following:

  1. 1.

The functions $\ell(y,\cdot)$ and $r$ are continuous and coercive, and the
function $\ell(y,\cdot)+r(\cdot)$ is strongly convex,

  2. 2.

The covariance matrices $\Sigma_{c}$ are positive definite, and their spectral
norms are uniformly bounded,

  3. 3.

The means $\bm{\mu_{c}}/\sqrt{d}$ and the teacher vector $\bm{\theta}_{0}$ are
uniformly bounded,

  4. 4.

the number of clusters $|\mathcal{C}|$ is finite,

  5. 5.

the distribution $P_{0}$ is sub-gaussian with uniformly bounded norm.

Then, as $n,d\to\infty$ with $n/d\to\alpha>0$, we have

$\displaystyle\varepsilon_{\text{tr}}(\bm{\hat{\theta}})$
$\displaystyle\overset{P}{\simeq}\;\sum_{c\in\mathcal{C}}p_{c}\,\mathbb{E}_{\omega^{(s)}_{c},\omega^{(t)}_{c},y}\left[\ell\left(y,\operatorname{prox}_{V_{c}^{\star}\ell(y,\cdot)}(\omega^{(s)}_{c})\right)\right]$
$\displaystyle=:\varepsilon_{\text{tr}}(\bm{\theta}_{0},\\{\bm{\mu}_{c}\\}_{c\in\mathcal{C}},\\{\Sigma_{c}\\}_{c\in\mathcal{C}})$
(33) $\displaystyle\varepsilon_{\text{tr}}(\bm{\hat{\theta}})$
$\displaystyle\overset{P}{\simeq}\;\sum_{c\in\mathcal{C}}p_{c}\,\mathbb{E}_{\omega^{(s)}_{c},\omega^{(t)}_{c},y}\left[\ell\left(y,\omega^{(s)}_{c}\right)\right]$
$\displaystyle=:\varepsilon_{\text{tr}}(\bm{\theta}_{0},\\{\bm{\mu}_{c}\\}_{c\in\mathcal{C}},\\{\Sigma_{c}\\}_{c\in\mathcal{C}})$

where

$\begin{pmatrix}\omega_{c}^{(s)}\\\
\omega_{c}^{(t)}\end{pmatrix}\sim\mathcal{N}\left(\begin{bmatrix}\pi_{c}\\\
h^{\star}_{c}\end{bmatrix},\begin{bmatrix}\rho&m^{\star}_{c}\\\
m^{\star}_{c}&q^{\star}_{c}\end{bmatrix}\right)\quad\text{and}\quad y\sim
P_{0}\left(\cdot\,\big{|}\,\omega_{c}^{(t)}\right).$ (34)

and the $\operatorname{prox}$ function is

$\operatorname{prox}_{f}(x)=\min_{z\in\mathbb{R}}\left[\frac{1}{2}(z-x)^{2}+f(z)\right]$
(35)

The overlaps used in the equation are defined as follows:

$\rho_{c}=\frac{1}{d}\bm{\theta}_{0}^{\top}\Sigma_{c}\bm{\theta}_{0},\quad\pi_{c}=\frac{1}{d}\bm{\theta}_{0}^{\top}\bm{\mu}_{c},$
(36)

and
$(V^{\star}_{c},q^{\star}_{c},m^{\star}_{c},h^{\star}_{c})_{c\in\mathcal{C}}$
are the unique fixed point of the following set of self-consistent _replica
saddle-point equations_ :

$\begin{split}&\begin{cases}\hat{V}_{c}&=\alpha
p_{c}\mathbb{E}_{\xi_{c}\sim\mathcal{N}(0,1)}\left[\int\mathrm{d}y~{}\mathcal{Z}_{0}\left(y,\pi_{c}+\frac{m_{c}}{\sqrt{q}_{c}}\xi_{c},\rho-\frac{m_{c}^{2}}{q_{c}}\right)\partial_{\omega}f_{\ell}(y,h_{c}+\sqrt{q_{c}}\xi_{c},V_{c})\right]\\\
\hat{q}_{c}&=\alpha
p_{c}\mathbb{E}_{\xi_{c}\sim\mathcal{N}(0,1)}\left[\int\mathrm{d}y~{}\mathcal{Z}_{0}\left(y,\pi_{c}+\frac{m_{c}}{\sqrt{q}_{c}}\xi_{c},\rho-\frac{m_{c}^{2}}{q_{c}}\right)f_{\ell}(y,h_{c}+\sqrt{q_{c}}\xi_{c},V_{c})^{2}\right]\\\
\hat{m}_{c}&=\alpha
p_{c}\mathbb{E}_{\xi_{c}\sim\mathcal{N}(0,1)}\left[\int\mathrm{d}y~{}\partial_{\omega}\mathcal{Z}_{0}\left(y,\pi_{c}+\frac{m_{c}}{\sqrt{q}_{c}}\xi_{c},\rho-\frac{m_{c}^{2}}{q_{c}}\right)f_{\ell}(y,h_{c}+\sqrt{q_{c}}\xi_{c},V_{c})\right]\\\
\hat{h}_{c}&=\alpha
p_{c}\mathbb{E}_{\xi_{c}\sim\mathcal{N}(0,1)}\left[\int\mathrm{d}y~{}\mathcal{Z}_{0}\left(y,\pi_{c}+\frac{m_{c}}{\sqrt{q}_{c}}\xi_{c},\rho-\frac{m_{c}^{2}}{q_{c}}\right)f_{\ell}(y,h_{c}+\sqrt{q_{c}}\xi_{c},V_{c})\right]\end{cases}\\\\[9.00002pt]
&\begin{cases}V_{c}&=\mathbb{E}_{\\{\bm{\xi}_{c}\\}~{}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}~{}\mathcal{N}(\bm{0},I_{d})}\left[\bm{\hat{\eta}}^{\top}\hat{q}_{c}^{-1/2}\Sigma_{c}^{1/2}\bm{\xi}_{c}\right]\\\
q_{c}&=\mathbb{E}_{\\{\bm{\xi}_{c}\\}~{}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}~{}\mathcal{N}(\bm{0},I_{d})}\left[\bm{\hat{\eta}}^{\top}\Sigma_{c}\bm{\hat{\eta}}\right]\\\
m_{c}&=\mathbb{E}_{\\{\bm{\xi}_{c}\\}~{}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}~{}\mathcal{N}(\bm{0},I_{d})}\left[\bm{\theta}_{0}^{\top}\Sigma_{c}\bm{\hat{\eta}}\right]\\\
h_{c}&=\mathbb{E}_{\\{\bm{\xi}_{c}\\}~{}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}~{}\mathcal{N}(\bm{0},I_{d})}\left[\bm{\mu}_{c}^{\top}\bm{\hat{\eta}}\right]\end{cases}\end{split}$
(37)

and we have defined the following auxiliary functions:

$\displaystyle\mathcal{Z}_{0}(y,\omega,V)$
$\displaystyle=\int\frac{\mathrm{d}\lambda}{\sqrt{2\pi
V}}P_{0}(y|\lambda)~{}e^{-\frac{(\lambda-\omega)^{2}}{2V}}$ (38)
$\displaystyle f_{\ell}(y,\omega,V)$
$\displaystyle=-V^{-1}\partial_{\omega}\mathcal{M}_{V\ell(y,\cdot)}(\omega)=V^{-1}\left(\operatorname{prox}_{V\ell(y,\cdot)}(\omega)-\omega\right).$
(39)

Finally, the auxiliary variable $\bm{\hat{\eta}}$ is a function of
$\\{\bm{\xi}_{c},\hat{V}_{c},\hat{m}_{c},\hat{q}_{c},\hat{h}_{c}\\}$:

$\bm{\hat{\eta}}=\operatorname{prox}_{r(\hat{\Sigma}^{-1/2}\cdot)}\left(\hat{\Sigma}^{-1/2}\left(\sum_{c\in\mathcal{C}}\bm{\hat{\mu}}_{c}+\sqrt{\hat{q}_{c}}\Sigma_{c}^{1/2}\bm{\xi}_{c}\right)\right)$
(40)

with

$\bm{\hat{\mu}_{c}}=\hat{h}_{c}\bm{\mu_{c}}+\hat{m}_{c}\Sigma_{c}\bm{\theta}_{0}\quad\text{and}\quad\hat{\Sigma}=\sum_{c\in\mathcal{C}}\hat{V}_{c}\Sigma_{c}$
(41)

We will not prove this theorem, since the proof is virtually equivalent to the
one in [22, 40]. Instead, the next sections are dedicated to the derivation of
these equations, using the so-called _replica method_ from statistical physics
[51]. Simpler particular cases of these equations can be found in Appendix
A.9.

###  A.2 Gibbs measure and free energy

The starting point for our replica computation is to define the Gibbs measure
over the weights $W\in\mathbb{R}^{K\times d}$:

$\displaystyle\mu_{\beta}(\bm{\theta})=\frac{1}{\mathcal{Z}_{\beta}}e^{-\beta\left[\sum\limits_{\nu=1}^{n}\ell\left(y^{\nu},\bm{\theta}^{\top}\bm{x}^{\nu}\right)+\lambda
r(\bm{\theta})\right]}=\frac{1}{\mathcal{Z}_{\beta}}e^{-\beta
r(\bm{\theta})}\prod\limits_{\nu=1}^{n}e^{-\beta\ell(y^{\nu},\bm{\theta}^{\top}\bm{x}^{\nu})}$
(42)

where $\beta>0$ is a parameter we eventually want to send to $\beta\to 0^{+}$,
and $\mathcal{Z}_{\beta}\in\mathbb{R}$ is the partition function (the
normalisation of the Gibbs measure). For convenience, we will define the
following useful notation:

$\displaystyle P_{\bm{\theta}}(\bm{\theta})\equiv e^{-\beta r(\bm{\theta})},$
$\displaystyle P_{\ell}(y^{\nu}|\bm{\theta}^{\top}\bm{x}^{\nu})\equiv
e^{-\beta\ell(y^{\nu},\bm{\theta}^{\top}\bm{x}^{\nu})}$ (43)

which allow us to write the Gibbs measure as:

$\displaystyle\mu_{\beta}(\bm{\theta})=\frac{1}{\mathcal{Z}_{\beta}}P_{\bm{\theta}}(\bm{\theta})\prod\limits_{\nu=1}^{n}P_{\ell}(y^{\nu}|\bm{\theta}^{\top}\bm{x}^{\nu})$
(44)

When $\beta\to\infty$, the Gibbs measure $\mu_{\beta}$ will concentrate around
the $\bm{\theta}$ that minimize the empirical risk (1). As a result, the free
energy density defined as

$\displaystyle-\beta
f_{\beta}=\lim\limits_{d\to\infty}\frac{1}{d}\mathbb{E}\log\mathcal{Z}_{\beta}$
(45)

where the limit is taken with $n/d=\alpha$ fixed, and the expectation is over
the distribution of the data, will concentrate around the minimum risk. In
order to take the expectation explicitly, we use the replica trick:

$\displaystyle\log{\mathcal{Z}_{\beta}}=\lim\limits_{s\to
0^{+}}\frac{1}{s}\partial_{s}\mathcal{Z}_{\beta}^{s}|_{s=0}$ (46)

Therefore, the replica computation boils down to the computation of the
averaged replicated partition function:

$\displaystyle\mathbb{E}\mathcal{Z}_{\beta}^{s}$
$\displaystyle=\mathbb{E}\left[\int_{\mathbb{R}^{d}}\mathrm{d}\bm{\theta}~{}P_{\bm{\theta}}(\bm{\theta})\prod_{\nu=1}^{n}P_{\ell}\left(y^{\nu}|\bm{\theta}^{\top}\bm{x}^{\nu}\right)\right]^{s}$
(47)
$\displaystyle=\prod\limits_{\nu=1}^{n}\mathbb{E}_{(\bm{x}^{\nu},y^{\nu})}\int\prod\limits_{a=1}^{s}\mathrm{d}\bm{\theta}^{a}~{}P_{\bm{\theta}}(\bm{\theta}^{a})P_{\ell}\left(y^{\nu}|{\bm{\theta}^{a}}^{\top}\bm{x}^{\nu}\right)$
(48)

###  A.3 Taking the average over the data

The main difference in this computation with respect to the usual binary
Gaussian mixture is that the labels are generated by a teacher target
function. In other words, the joint distribution of the data is given by:

$\displaystyle
P_{\bm{\theta}_{0}}(\bm{x},y)=P_{0}(y|\tau)\sum\limits_{c\in\mathcal{C}}\rho_{c}\,\mathcal{N}\left(\bm{x};\frac{\bm{\mu}_{c}}{\sqrt{d}},I_{d}\right)$
(49)

where $P_{0}(y|\tau)$ is the probability induced by the teacher activation
$f_{0}(\tau)$. The expression in (48) can then be written as:

$\displaystyle\mathbb{E}\mathcal{Z}_{\beta}^{s}$
$\displaystyle=\int\prod\limits_{a=1}^{s}\mathrm{d}\bm{\theta}^{a}~{}P_{\bm{\theta}}(\bm{\theta}^{a})\left(\int_{\mathbb{R}}\mathrm{d}{y}\int_{\mathbb{R}^{d}}\mathrm{d}{\bm{x}}P_{0}(\bm{y}|\bm{\theta}_{0}^{\top}\bm{x})P_{\bm{x}}(\bm{x})\prod\limits_{a=1}^{s}P_{\ell}\left(y^{\nu}|{\bm{\theta}^{a}}^{\top}\bm{x}^{\nu}\right)\right)^{n}$
(50)

We make a change of variables by introducing the _local fields_
$\tau=\bm{\theta}_{0}^{\top}\bm{x}$ and
$\lambda^{a}=\bm{\theta}^{a\top}\bm{x}$. The distribution of the local fields
will be a Gaussian mixture itself, and we can compute the moments for each
mode:

$\begin{split}\pi_{c}&\equiv\mathbb{E}_{\bm{x}^{\nu}}\left[\tau_{c}\right]=\frac{1}{d}\bm{\theta}_{0}^{\top}\bm{\mu}_{c}\\\
\rho_{c}&\equiv\text{Var}_{\bm{x}^{\nu}}\left[\tau_{c}\right]=\frac{1}{d}\bm{\theta}_{0}^{\top}\Sigma_{c}\bm{\theta}_{0}\\\
h^{a}_{c}&\equiv\mathbb{E}_{\bm{x}^{\nu}}\left[\lambda^{a}_{c}\right]=\frac{1}{d}{\bm{\theta}^{a}}^{\top}\bm{\mu}_{c}\\\
q_{c}^{ab}&\equiv\text{Cov}_{\bm{x}^{\nu}}\left[\lambda_{c}^{a},\lambda_{c}^{b}\right]=\frac{1}{d}{\bm{\theta}^{a}}^{\top}\Sigma_{c}\bm{\theta}^{b}\\\
m_{c}^{a}&\equiv\text{Cov}_{\bm{x}^{\nu}}\left[\tau_{c},\lambda_{c}^{a}\right]=\frac{1}{d}\bm{\theta}_{0}^{\top}\Sigma_{c}\bm{\theta}^{a}\end{split}$
(51)

Equivalently, the local fields distribution is the following low-dimensional
Gaussian mixture distribution:

$\displaystyle
P(\tau,\bm{\lambda})=\sum\limits_{c\in\mathcal{C}}p_{c}~{}\mathcal{N}\left(\begin{bmatrix}\pi_{c}\\\
\bm{h}_{c}\end{bmatrix},\begin{bmatrix}\rho_{c}&\bm{m}_{c}^{\top}\\\
\bm{m}_{c}&Q_{c}\end{bmatrix}\right),$ (52)

where we have defined the vectors $\bm{m},\bm{h}\in\mathbb{R}^{s}$ with
entries $m^{a}$, $h^{a}$ and the matrix $Q\in\mathbb{R}^{s\times s}$ with
entries $q^{ab}$. Therefore, we can factorize the partition function by only
integrating over the local field distribution:

$\displaystyle\mathbb{E}\mathcal{Z}_{\beta}^{s}$
$\displaystyle=\int\prod\limits_{a=1}^{s}\mathrm{d}\bm{\theta}^{a}~{}P_{\bm{\theta}}(\bm{\theta}^{a})\left[\int\mathrm{d}y\int\mathrm{d}\tau~{}P_{0}(y|\tau)\int\left(\prod\limits_{a=1}^{s}\mathrm{d}\lambda^{a}P_{\ell}(y|\lambda^{a})\right)P(\tau,\bm{\lambda}|\bm{m},\bm{h},q)\right]^{n}$
(53)

###  A.4 Writing as a saddle-point problem

Note that $(\pi_{c},\rho_{c})$ are fixed inputs in the problem. By Fourier
transform arguments, we can write

$\delta\left(m_{c}^{a}-\frac{1}{d}\bm{\theta}_{0}^{\top}\Sigma_{c}\bm{\theta}^{a}\right)=\int_{i\mathbb{R}}\frac{\mathrm{d}{\hat{m}_{c}^{a}}}{2\pi}\,e^{\hat{m}_{c}^{a}\left(dm_{c}^{a}-\bm{\theta}_{0}^{\top}\Sigma_{c}\bm{\theta}^{a}\right)},$
(54)

where the integral is on the imaginary line. By doing the same arguments for
the $q_{c}^{ab}$ and $h_{c}^{a}$, it ensues that

$\displaystyle\mathbb{E}\mathcal{Z}_{\beta}^{s}$
$\displaystyle=\int\prod_{c\in\mathcal{C}}\prod\limits_{a=1}^{s}\frac{\mathrm{d}m_{c}^{a}\hat{m}_{c}^{a}}{2\pi}\frac{\mathrm{d}h_{c}^{a}\hat{h}_{c}^{a}}{2\pi}\prod\limits_{1\leq
a\leq b\leq
s}\frac{\mathrm{d}q_{c}^{ab}\mathrm{d}\hat{q}_{c}^{ab}}{2\pi}e^{d\Phi^{(s)}_{\rm{rs}}(\bm{m},\hat{\bm{m}},q,\hat{q})}$
(55)

where we have defined the free energy potential:

$\displaystyle\Phi^{(s)}_{\rm{rs}}(\bm{m},\hat{\bm{m}},q,\hat{q})$
$\displaystyle=\sum\limits_{c\in\mathcal{C}}\left[\sum\limits_{a=1}^{s}\hat{m}_{c}^{a}m_{c}^{a}+\sum\limits_{a=1}^{s}\hat{h}_{c}^{a}h_{c}^{a}+\sum\limits_{1\leq
a\leq b\leq
s}\hat{q}_{c}^{ab}q_{c}^{ab}\right]-\alpha\Psi^{(s)}_{g}(\bm{m},\bm{h},q)-\Psi^{(s)}_{\theta}(\hat{m},\hat{h},\hat{q})$
$\displaystyle\Psi^{(s)}_{\ell}(\bm{m},\bm{h},q)$
$\displaystyle\equiv\log\int\mathrm{d}y\int\mathrm{d}\tau~{}P_{0}(y|\tau)\int\prod\limits_{a=1}^{s}\mathrm{d}\lambda^{a}P_{\ell}(y|\lambda^{a})P(\tau,\bm{\lambda}|\bm{m},\bm{h},q)$
$\displaystyle\Psi^{(s)}_{\theta}(\hat{\bm{m}},\hat{\bm{h}},\hat{q})$
$\displaystyle\equiv\frac{1}{d}\log\int\prod\limits_{a=1}^{s}\mathrm{d}\bm{\theta}^{a}~{}P_{\bm{\theta}}(\bm{\theta})\prod\limits_{c\in\mathcal{C}}e^{\sum\limits_{a=1}^{s}\left[\hat{h}_{c}^{a}{\bm{\theta}^{a}}^{\top}\bm{\mu}_{c}+\hat{m}^{a}\bm{\theta}_{0}^{\top}\Sigma_{c}\bm{\theta}^{a}\right]+\sum\limits_{1\leq
a\leq b\leq
s}\hat{q}_{c}^{ab}{\bm{\theta}^{a}}^{\top}\Sigma_{c}\bm{\theta}^{b}}$ (56)

Therefore, as we take the $d\to\infty$ limit, we can apply the saddle-point
method [51] to compute $f_{\beta}$

$\displaystyle-\beta
f_{\beta}=\underset{\bm{m},\hat{\bm{m}},\bm{h},\hat{\bm{h}},q,\hat{q}}{\operatorname{extr}}~{}\lim\limits_{s\to
0^{+}}\frac{\Phi^{(s)}(\bm{m},\hat{\bm{m}},\bm{h},\hat{\bm{h}},q,\hat{q})}{s}$
(57)

######  Remark A.1.

We have introduced the parameters $\hat{q},\hat{m},\hat{h}$ as pure imaginary
numbers, but the optimization problem considers them as real numbers. This
stems from the fact that the function $\Phi$ is holomorphic, so the integral
is independent from the contour of integration. More details can be found in
[52].

###  A.5 Replica symmetric ansatz

In order to make progress with the $s\to 0^{+}$ limit, we make the following
replica symmetric ansatz:

$\displaystyle m_{c}^{a}=m_{c}$ $\displaystyle\hat{m}^{a}_{c}=\hat{m}_{c},$
$\displaystyle a=1,\cdots,s$ $\displaystyle h_{c}^{a}=h_{c}$
$\displaystyle\hat{h}^{a}_{c}=\hat{h}_{c},$ $\displaystyle a=1,\cdots,s$
$\displaystyle q_{c}^{aa}=r_{c},$
$\displaystyle\hat{q}^{aa}_{c}=-\frac{1}{2}\hat{r}_{c},$ $\displaystyle
a=1,\cdots,s$ $\displaystyle q_{c}^{ab}=q_{c},$
$\displaystyle\hat{q}^{aa}_{c}=\hat{q}_{c},$ $\displaystyle 1\leq a<b\leq s$
(58)

for all $c\in\mathcal{C}$. Since $\ell$ and $r$ are convex, the function
$\Phi^{(s)}$ has a unique saddle point, which must therefore coincide with the
replica symmetric ansatz. We also define

$\displaystyle
V_{c}=r_{c}-q_{c},\qquad\text{and}\qquad\hat{V}_{c}=\hat{r}_{c}-\hat{q}_{c}$
(59)

By inserting this ansatz above, we can take the $s\to 0^{+}$ limit. This is
done through a classical but computationally heavy method known as the
_Hubbard-Stratonovich transform_ ; details can be found in [14], Appendix C.
We simply reproduce the final results here:

$\displaystyle-\beta
f_{\beta}=\underset{\\{m_{c},\hat{m_{c}},h_{c},\hat{h}_{c},q_{c},\hat{q_{c}},V_{c},\hat{V}_{c}\\}}{\text{extr}}~{}\Phi_{\text{rs}}(\\{m_{c},\hat{m_{c}},h_{c},\hat{h}_{c},q_{c},\hat{q_{c}},V_{c},\hat{V}_{c}\\})$
(60)

The replicated free energy $\Phi_{\text{rs}}$ is given by the following
formula:

$\displaystyle\Phi_{\text{rs}}(\\{m_{c},\hat{m_{c}},h_{c},\hat{h}_{c},q_{c},\hat{q_{c}},V_{c},\hat{V}_{c}\\})=$
$\displaystyle\sum\limits_{c\in\mathcal{C}}\left[\frac{1}{2}\left(\hat{V}_{c}q_{c}-\hat{q}_{c}V_{c}\right)-\frac{1}{2}\hat{V}_{c}V_{c}+\hat{m}_{c}m_{c}+\hat{h}_{c}h_{c}\right]$
(61)
$\displaystyle-\alpha\Psi_{\ell}(\\{m_{c},h_{c},q_{c},V_{c}\\})-\Psi_{\bm{\theta}}(\\{\hat{m_{c}},\hat{h}_{c},\hat{q_{c}},\hat{V}_{c}\\})$

where we have decomposed the contributions coming from the loss
($\Psi_{\ell}$) and the regularization ($\Psi_{\bm{\theta}}$):

$\displaystyle\Psi_{\ell}$
$\displaystyle=\sum\limits_{c\in\mathcal{C}}\mathbb{E}_{\xi_{c}~{}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}~{}\mathcal{N}(0,1)}\int\mathrm{d}y\mathcal{Z}_{0}\left(y,\frac{m_{c}}{\sqrt{q_{c}}}\xi+\pi_{c},\rho-\frac{m_{c}^{2}}{q_{c}}\right)\log\mathcal{Z}_{\ell}(y,h_{c}+\sqrt{q_{c}}\,\xi_{c},V_{c})$
(63) $\displaystyle\Psi_{\bm{\theta}}$
$\displaystyle=\frac{1}{d}\,\mathbb{E}_{\bm{\xi}_{c}~{}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}~{}\mathcal{N}(\bm{0},I_{d})}\log\mathcal{Z}_{\bm{\theta}}\left(\sum\limits_{c\in\mathcal{C}}\hat{V}_{c}\Sigma_{c},\sum\limits_{c\in\mathcal{C}}\hat{h}_{c}\bm{\mu}_{c}+\hat{m}_{c}\Sigma_{c}\bm{\theta}_{0}+\sqrt{\hat{q}_{c}}\Sigma_{c}^{1/2}\bm{\xi}_{c}\right),$
(64)

and defined the following auxiliary free energies:

$\displaystyle\mathcal{Z}_{\ell/0}(y,\omega,V)$
$\displaystyle=\int\frac{\mathrm{d}\lambda}{\sqrt{2\pi
V}}P_{\ell/0}(y|\lambda)~{}e^{-\frac{(\lambda-\omega)^{2}}{2V}}$ (65)
$\displaystyle\mathcal{Z}_{\bm{\theta}}(A,\bm{b})$
$\displaystyle=\int\mathrm{d}\bm{\theta}e^{-\beta
r(\bm{\theta})}e^{-\frac{1}{2}\bm{\theta}^{\top}A\bm{\theta}+\bm{b}^{\top}\bm{\theta}}$
(66)

###  A.6 Taking the zero temperature limit

In order to take the $\beta\to\infty$ limit, we make the following rescalings:

$\displaystyle V_{c}$ $\displaystyle\to\beta^{-1}V_{c}$
$\displaystyle\hat{V}_{c}$ $\displaystyle\to\beta\hat{V}_{c}$
$\displaystyle\hat{q}_{c}$ $\displaystyle\to\beta^{2}\hat{q}_{c}$
$\displaystyle\hat{m}_{c}$ $\displaystyle\to\beta\hat{m}_{c}.$ (67)

It is easy to check how this rescaling affects $\beta^{-1}\Phi^{(\text{rs})}$,
so we only need to consider $\Psi_{\ell}$ and $\Psi_{\bm{\theta}}$.

We start with the latter: letting

$\mathcal{L}_{\bm{\theta}}(\bm{\theta})=\frac{1}{2}\,\bm{\theta}^{\top}A\bm{\theta}-\bm{b}^{\top}\bm{\theta}+r(\bm{\theta})$

by the Laplace method for any $A,\bm{b}$,

$\displaystyle\lim\limits_{\beta\to\infty}\frac{1}{\beta}\log
Z_{\bm{\theta}}(\beta
A,\beta\bm{b})=-\inf_{\bm{\theta}}\mathcal{L}_{\bm{\theta}}(\bm{\theta})=-\mathcal{L}_{\bm{\theta}}(\bm{\eta})$
(68)

where

$\bm{\eta}=\operatorname{prox}_{r(A^{1/2}\cdot)}(A^{1/2}\bm{b}).$ (69)

As a result, every integral involved in $\Psi_{\bm{\theta}}$ (and, later, its
partial derivatives) concentrates around its value at $\bm{\hat{\eta}}$
defined in (40). For the term $\Psi_{\ell}$, the term $\mathcal{Z}_{0}$ is
left unchanged, and by the same reasoning as above

$\lim\limits_{\beta\to\infty}\frac{1}{\beta}\log\mathcal{Z}_{\ell}(y,\omega,\beta
V)=-V^{-1}\mathcal{M}_{V\ell(y,\cdot)}(\omega).$ (70)

###  A.7 Saddle-point equations

The saddle-point equations are obtained by taking the derivatives of the free
energy potential with respect to the overlap parameters. We obtain a set of
self-consistent equations which we should solve in order to find a fixed
point:

$\displaystyle\begin{cases}\hat{V}_{c}&=\alpha
p_{c}\mathbb{E}_{\xi_{c}\sim\mathcal{N}(0,1)}\left[\int\mathrm{d}y~{}\mathcal{Z}_{0}\left(y,\pi_{c}+\frac{m_{c}}{\sqrt{q}_{c}}\xi_{c},\rho-\frac{m_{c}^{2}}{q_{c}}\right)\partial_{\omega}f_{\ell}(y,h_{c}+\sqrt{q_{c}}\xi_{c},V_{c})\right]\\\
\hat{q}_{c}&=\alpha
p_{c}\mathbb{E}_{\xi_{c}\sim\mathcal{N}(0,1)}\left[\int\mathrm{d}y~{}\mathcal{Z}_{0}\left(y,\pi_{c}+\frac{m_{c}}{\sqrt{q}_{c}}\xi_{c},\rho-\frac{m_{c}^{2}}{q_{c}}\right)f_{\ell}(y,h_{c}+\sqrt{q_{c}}\xi_{c},V_{c})^{2}\right]\\\
\hat{m}_{c}&=\alpha
p_{c}\mathbb{E}_{\xi_{c}\sim\mathcal{N}(0,1)}\left[\int\mathrm{d}y~{}\partial_{\omega}\mathcal{Z}_{0}\left(y,\pi_{c}+\frac{m_{c}}{\sqrt{q}_{c}}\xi_{c},\rho-\frac{m_{c}^{2}}{q_{c}}\right)f_{\ell}(y,h_{c}+\sqrt{q_{c}}\xi_{c},V_{c})\right]\\\
\hat{h}_{c}&=\alpha
p_{c}\mathbb{E}_{\xi_{c}\sim\mathcal{N}(0,1)}\left[\int\mathrm{d}y~{}\mathcal{Z}_{0}\left(y,\pi_{c}+\frac{m_{c}}{\sqrt{q}_{c}}\xi_{c},\rho-\frac{m_{c}^{2}}{q_{c}}\right)f_{\ell}(y,h_{c}+\sqrt{q_{c}}\xi_{c},V_{c})\right]\end{cases}$
(71)
$\displaystyle\begin{cases}V_{c}&=\mathbb{E}_{\\{\bm{\xi}_{c}\\}~{}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}~{}\mathcal{N}(\bm{0},I_{d})}\left[\bm{\hat{\eta}}^{\top}\hat{q}_{c}^{-1/2}\Sigma_{c}^{1/2}\bm{\xi}_{c}\right]\\\
q_{c}&=\mathbb{E}_{\\{\bm{\xi}_{c}\\}~{}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}~{}\mathcal{N}(\bm{0},I_{d})}\left[\bm{\hat{\eta}}^{\top}\Sigma_{c}\bm{\hat{\eta}}\right]\\\
m_{c}&=\mathbb{E}_{\\{\bm{\xi}_{c}\\}~{}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}~{}\mathcal{N}(\bm{0},I_{d})}\left[\bm{\theta}_{0}^{\top}\Sigma_{c}\bm{\hat{\eta}}\right]\\\
h_{c}&=\mathbb{E}_{\\{\bm{\xi}_{c}\\}~{}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}~{}\mathcal{N}(\bm{0},I_{d})}\left[\bm{\mu}_{c}^{\top}\bm{\hat{\eta}}\right]\end{cases}$
(72)

where all the relevant quantities have been defined in Appendix A.1.

###  A.8 Training and generalization errors

It now remains to compute the training and generalization errors from the free
energy. Recalling that

$\varepsilon_{\text{tr}}(\hat{\vec{\theta}})=\frac{1}{n}\sum\limits_{\nu=1}^{n}\ell\left(y^{\nu},\hat{\bm{\theta}}^{\top}\bm{x}^{\nu}\right),$

we can write using the definition of the free energy (45):

$\lim_{n\to\infty}\varepsilon_{\text{tr}}=\lim_{\beta\to\infty}\partial_{\beta}f_{\beta}-r(\bm{\hat{\theta}}).$

Computing explicitly the derivative and averaging again over the data yields

$\lim_{n\to\infty}\varepsilon_{\text{tr}}=-\lim_{\beta\to\infty}\partial_{\beta}\Psi_{\ell},$

where $\Psi_{\ell}$ is the free energy contribution of the loss defined in
(63). Writing explicitly this derivative,

$-\partial_{\beta}\Psi_{\ell}=\sum_{c\in\mathcal{C}}\mathbb{E}_{\xi_{c}\sim\mathcal{N}(0,1)}\int_{\mathbb{R}}\mathrm{d}{y}\frac{\mathcal{Z}_{0}\left(y,\omega_{c}^{(0)},\rho-\frac{m_{c}^{2}}{q_{c}}\right)}{\mathcal{Z}_{\ell}(y,\omega_{c}^{(\ell)},V_{c})}\underbrace{\int_{\mathbb{R}}\frac{\mathrm{d}{\lambda}}{\sqrt{2\pi
V_{c}}}e^{-(\lambda-\omega_{c}^{(\ell)})V_{c}^{-1}(\lambda-\omega_{c}^{(\ell)},V_{c}))-\beta\ell(y,\lambda)}\ell(y,\lambda)}_{\tilde{\mathcal{Z}}_{\ell}(y,\omega_{c}^{(\ell)},V_{c})}$
(73)

where

$\omega_{c}^{(0)}=\pi_{c}+\frac{m_{c}}{\sqrt{q_{c}}}\xi_{c}\quad\text{and}\quad\omega_{c}^{(\ell)}=h_{c}+\sqrt{q_{c}}\xi_{c}.$

We can now make the change of variables in (67), and use again Laplace’s
approximation: if

$\eta_{c}=\operatorname{prox}_{V\ell(y,\cdot)}(\omega_{c}^{(\ell)})\quad\text{and}\quad\mathcal{L}(\lambda)=-(\lambda-\omega_{c}^{(\ell)})V_{c}^{-1}(\lambda-\omega_{c}^{(\ell)})-\beta\ell(y,\lambda),$

then

$\mathcal{Z}_{\ell}(y,\omega_{c}^{(\ell)},V_{c})\sim
e^{-\beta\mathcal{L}(\eta_{c})}\quad\text{and}\quad\tilde{\mathcal{Z}}_{\ell}(y,\omega_{c}^{(\ell)},V_{c})\sim
e^{-\beta\mathcal{L}(\eta_{c})}\ell(y,\eta_{c}),$

and all of the overlaps will concentrate around the solutions of (71), (72).
Finally, the term containing $\mathcal{Z}_{0}$ is an expectation of $y$
according to $P_{0}(y\,\mid\,\tau_{c})$, where

$\tau_{c}=\omega_{c}^{(0)}+\sqrt{\rho-\frac{m_{c}^{2}}{q_{c}}}\xi^{(0)}_{c}$

and $\xi_{c}^{(0)}$ is a Gaussian variable independent from everything else.
Putting all together, we can write

$\lim_{n\to\infty}\varepsilon_{\text{tr}}(\bm{\hat{\theta}})=\sum_{c\in\mathcal{C}}p_{c}\,\mathbb{E}_{\nu_{c},\tau_{c},y}\left[\ell\left(y,\operatorname{prox}_{V_{c}^{\star}\ell(y,\cdot)}(\nu_{c})\right)\right]$
(74)

with

$\begin{pmatrix}\nu_{c}\\\
\tau_{c}\end{pmatrix}\sim\mathcal{N}\left(\begin{bmatrix}\pi_{c}\\\
h_{c}\end{bmatrix},\begin{bmatrix}\rho_{c}&m_{c}^{\star}\\\
m_{c}^{\star}&q_{c}^{\star}\end{bmatrix}\right)\quad\text{and}\quad y\sim
P_{0}(\cdot\,|\,\tau_{c}).$ (75)

The generalization error is much simpler to obtain: since
$x_{\text{new}},y_{\text{new}}$ are independent from the estimator
$\bm{\hat{\theta}}$, we simply have

$\lim_{n\to\infty}\varepsilon_{\text{gen}}(\bm{\hat{\theta}})=\sum_{c\in\mathcal{C}}p_{c}\,\mathbb{E}_{\nu_{c},\tau_{c},y}\left[\ell\left(y,\nu_{c}\right)\right]$
(76)

where $\nu_{c},\tau_{c},y$ follow the same distribution as above.

###  A.9 Examples

#### Ridge penalty

Consider a particular case of the general equations reported above: the case
of a ridge penalty

$r(\bm{\theta})=\frac{\lambda}{2}||\bm{\theta}||^{2}_{2}.$

We then have:

$\displaystyle\operatorname{prox}_{r}(\bm{x})=(1+\lambda)^{-1}\bm{x}$ (77)

which simplify the prior equations considerably. Indeed, we can now compute
every expectation in (72), which yields the following fixed-point equations:

$\displaystyle\begin{cases}\hat{V}_{c}&=-\alpha
p_{c}\mathbb{E}_{\xi_{c}\sim\mathcal{N}(0,1)}\left[\int\mathrm{d}y~{}\mathcal{Z}_{0}\left(y,\pi_{c}+\frac{m_{c}}{\sqrt{q}_{c}}\xi_{c},\rho-\frac{m_{c}^{2}}{q_{c}}\right)\partial_{\omega}f_{\ell}(y,h_{c}+\sqrt{q_{c}}\xi_{c},V_{c})\right]\\\
\hat{q}_{c}&=\alpha
p_{c}\mathbb{E}_{\xi_{c}\sim\mathcal{N}(0,1)}\left[\int\mathrm{d}y~{}\mathcal{Z}_{0}\left(y,\pi_{c}+\frac{m_{c}}{\sqrt{q}_{c}}\xi_{c},\rho-\frac{m_{c}^{2}}{q_{c}}\right)f_{\ell}(y,h_{c}+\sqrt{q_{c}}\xi_{c},V_{c})^{2}\right]\\\
\hat{m}_{c}&=\alpha
p_{c}\mathbb{E}_{\xi_{c}\sim\mathcal{N}(0,1)}\left[\int\mathrm{d}y~{}\partial_{\omega}\mathcal{Z}_{0}\left(y,\pi_{c}+\frac{m_{c}}{\sqrt{q}_{c}}\xi_{c},\rho-\frac{m_{c}^{2}}{q_{c}}\right)f_{\ell}(y,h_{c}+\sqrt{q_{c}}\xi_{c},V_{c})\right]\\\
\hat{h}_{c}&=\alpha
p_{c}\mathbb{E}_{\xi_{c}\sim\mathcal{N}(0,1)}\left[\int\mathrm{d}y~{}\mathcal{Z}_{0}\left(y,\pi_{c}+\frac{m_{c}}{\sqrt{q}_{c}}\xi_{c},\rho-\frac{m_{c}^{2}}{q_{c}}\right)f_{\ell}(y,h_{c}+\sqrt{q_{c}}\xi_{c},V_{c})\right]\end{cases}$
(78)
$\displaystyle\begin{cases}V_{c}&=\frac{1}{d}\operatorname{tr}\left[\Sigma_{c}\left(\lambda
I_{d}+\hat{\Sigma}\right)^{-1}\right]\\\ \vskip 1.79997pt\cr
q_{c}&=\frac{1}{d}\operatorname{tr}\left[\left(\sum\limits_{c^{\prime}\in\mathcal{C}}\hat{q}_{c^{\prime}}\Sigma_{c^{\prime}}+\sum\limits_{c^{\prime},c^{\prime\prime}\in\mathcal{C}}\bm{\hat{\mu}}_{c^{\prime}}\bm{\hat{\mu}}_{c^{\prime\prime}}^{\top}\right)\Sigma_{c}\left(\lambda
I_{d}+\hat{\Sigma}\right)^{-2}\right]\\\ \vskip 1.79997pt\cr
m_{c}&=\frac{1}{d}\operatorname{tr}\left[\left(\sum\limits_{c^{\prime}\in\mathcal{C}}\bm{\hat{\mu}}_{c^{\prime}}\bm{\theta}_{0}^{\top}\right)\Sigma_{c}\left(\lambda
I_{d}+\hat{\Sigma}\right)^{-1}\right]\\\ \vskip 1.79997pt\cr
h_{c}&=\frac{1}{d}\operatorname{tr}\left[\left(\sum\limits_{c^{\prime}\in\mathcal{C}}\bm{\hat{\mu}}_{c^{\prime}}\bm{\mu}_{c}^{\top}\right)\left(\lambda
I_{d}+\hat{\Sigma}\right)^{-1}\right]\end{cases}$ (79)

#### Gaussian covariate model

The equations for the Gaussian covariate model can be found in [19]; they also
correspond to taking $|\mathcal{C}|=1$ in the ones above. We reproduce them
here for completeness:

$\displaystyle\begin{cases}\hat{V}&=-\alpha\mathbb{E}_{\xi\sim\mathcal{N}(0,1)}\left[\int\mathrm{d}y~{}\mathcal{Z}_{0}\left(y,\pi+\frac{m}{\sqrt{q}}\xi,\rho-\frac{m^{2}}{q}\right)\partial_{\omega}f_{\ell}(y,h+\sqrt{q}\xi,V)\right]\\\
\hat{q}&=\alpha\mathbb{E}_{\xi\sim\mathcal{N}(0,1)}\left[\int\mathrm{d}y~{}\mathcal{Z}_{0}\left(y,\pi+\frac{m}{\sqrt{q}}\xi,\rho-\frac{m^{2}}{q}\right)f_{\ell}(y,h+\sqrt{q}\xi,V)^{2}\right]\\\
\hat{m}&=\alpha\mathbb{E}_{\xi\sim\mathcal{N}(0,1)}\left[\int\mathrm{d}y~{}\partial_{\omega}\mathcal{Z}_{0}\left(y,\pi+\frac{m}{\sqrt{q}}\xi,\rho-\frac{m^{2}}{q}\right)f_{\ell}(y,h+\sqrt{q}\xi,V)\right]\\\
\hat{h}&=\alpha\mathbb{E}_{\xi\sim\mathcal{N}(0,1)}\left[\int\mathrm{d}y~{}\mathcal{Z}_{0}\left(y,\pi+\frac{m}{\sqrt{q}}\xi,\rho-\frac{m^{2}}{q}\right)f_{\ell}(y,h+\sqrt{q}\xi,V)\right]\end{cases}$
(80)
$\displaystyle\begin{cases}V&=\frac{1}{d}\operatorname{tr}\left[\Sigma\left(\lambda
I_{d}+\hat{\Sigma}\right)^{-1}\right]\\\ \vskip 1.79997pt\cr
q&=\frac{1}{d}\operatorname{tr}\left[\left(\hat{q}\Sigma+\bm{\hat{\mu}}\bm{\hat{\mu}}^{\top}\right)\Sigma\left(\lambda
I_{d}+\hat{\Sigma}\right)^{-2}\right]\\\ \vskip 1.79997pt\cr
m&=\frac{1}{d}\operatorname{tr}\left[\bm{\hat{\mu}}\bm{\theta}_{0}^{\top}\Sigma\left(\lambda
I_{d}+\hat{\Sigma}\right)^{-1}\right]\\\ \vskip 1.79997pt\cr
h&=\frac{1}{d}\operatorname{tr}\left[\bm{\hat{\mu}}\bm{\mu}^{\top}\left(\lambda
I_{d}+\hat{\Sigma}\right)^{-1}\right]\end{cases}$ (81)

where this time

$\bm{\hat{\mu}}=\hat{h}\bm{\mu}+\hat{m}\Sigma\bm{\theta}_{0}\qquad\text{and}\qquad\hat{\Sigma}=\hat{V}\Sigma.$
(82)

The errors are given by

$\displaystyle\lim_{n\to\infty}\varepsilon_{\text{tr}}(\bm{\hat{\theta}})$
$\displaystyle=\mathbb{E}_{\nu,\tau,y}\left[\ell\left(y,\operatorname{prox}_{V^{\star}\ell(y,\cdot)}(\nu)\right)\right]$
(83)
$\displaystyle\lim_{n\to\infty}\varepsilon_{\text{gen}}(\bm{\hat{\theta}})$
$\displaystyle=\mathbb{E}_{\nu,\tau,y}\left[\ell\left(y,\nu\right)\right]$
(84)

with

$\begin{pmatrix}\nu\\\
\tau\end{pmatrix}\sim\mathcal{N}\left(\begin{bmatrix}\pi\\\
h\end{bmatrix},\begin{bmatrix}\rho&m^{\star}\\\
m^{\star}&q^{\star}\end{bmatrix}\right)\quad\text{and}\quad y\sim
P_{0}(\cdot\,|\,\tau).$ (85)

#### Ridge regression

We place ourselves in the ridge regression case, where

$\displaystyle P_{0}(\tau|\Delta)$ $\displaystyle=\mathcal{N}(\tau,\Delta)$
$\displaystyle l(y,\hat{y})$ $\displaystyle=(y-\hat{y})^{2}.$ (86)

In this case, the equations in (79) simplify even further, yielding

$\displaystyle\begin{cases}\hat{V}_{c}&=\frac{\alpha p_{c}}{1+V_{c}}\\\
\hat{q}_{c}&=\frac{\alpha
p_{c}}{(1+V_{c})^{2}}(\rho_{c}+\Delta+q_{c}-2m_{c}+(h_{c}-\pi_{c})^{2})\\\
\hat{m}_{c}&=\frac{\alpha p_{c}}{1+V_{c}}\\\ \hat{h}_{c}&=\frac{\alpha
p_{c}(\pi_{c}-h_{c})}{1+V_{c}}\end{cases}$ (87)
$\displaystyle\begin{cases}V&=\frac{1}{d}\operatorname{tr}\left[\Sigma\left(\lambda
I_{d}+\hat{\Sigma}\right)^{-1}\right]\\\ \vskip 1.79997pt\cr
q&=\frac{1}{d}\operatorname{tr}\left[\left(\hat{q}\Sigma+\bm{\hat{\mu}}\bm{\hat{\mu}}^{\top}\right)\Sigma\left(\lambda
I_{d}+\hat{\Sigma}\right)^{-2}\right]\\\ \vskip 1.79997pt\cr
m&=\frac{1}{d}\operatorname{tr}\left[\bm{\hat{\mu}}\bm{\theta}_{0}^{\top}\Sigma\left(\lambda
I_{d}+\hat{\Sigma}\right)^{-1}\right]\\\ \vskip 1.79997pt\cr
h&=\frac{1}{d}\operatorname{tr}\left[\bm{\hat{\mu}}\bm{\mu}^{\top}\left(\lambda
I_{d}+\hat{\Sigma}\right)^{-1}\right]\end{cases}$ (88)

The training and generalization error also benefit from a very simple
expression

$\displaystyle\varepsilon_{\text{tr}}=\sum_{c}\frac{\hat{q}_{c}^{\star}}{\alpha}$
$\displaystyle\varepsilon_{\text{tr}}=\sum_{c}\frac{\hat{q}^{\star}_{c}(1+V_{c}^{\star})^{2}}{\alpha}.$
(89)

##  Appendix B Universality of Gaussian Mixture Models

In this section we prove the main results shown in in Sec. 3 regarding
universality properties of Gaussian Mixture Models.

###  B.1 Mean Universality : proof of Proposition 3.2

Let us rewrite the assumptions here. Let $\\{\hat{V}^{\star}_{c}\\}_{c=1}^{K}$
be the fixed points of the (replica) saddle point equations describing the
centered Gaussian mixture problem, i.e. the solutions of of eqs. (78), (79).
We assume that the teacher vector and the data structure respect the
following:

$\displaystyle\lim_{n,d\to\infty}\frac{\vec{\theta}_{0}^{\top}\vec{\mu}_{c}}{d}=0\qquad\forall
c\in\mathcal{C}$ (90)
$\displaystyle\lim_{n,d\to\infty}\frac{1}{d}\vec{\theta}_{0}^{\top}\Sigma_{c^{\prime}}\left(\lambda+\sum_{c\in\mathcal{C}}\hat{V}^{\star}_{c}\Sigma_{c}\right)^{-1}\bm{\mu}_{c}\to
0\qquad\forall(c,c^{\prime})\in\mathcal{C}\times\mathcal{C}$ (91)

and the loss and teacher are both symmetric:

$\displaystyle\ell(x,y)$ $\displaystyle=\ell(-x,-y)$ (92) $\displaystyle
P_{0}\left(y|\tau\right)$ $\displaystyle=P_{0}\left(-y|-\tau\right).$ (93)

In the saddle-point equations (78), (79), the mean vectors appear always
coupled with the overlaps $\\{h_{c},\hat{h}_{c}\\}_{c\in\mathcal{C}}$. Hence,
to prove Prop. 3.2 it suffices to prove the following result:

######  Lemma B.1.

If $(a+b+c)$ hold, then $\\{h_{c}=\hat{h}_{c}=0\\}_{c\in\mathcal{C}}$ is a
fixed point for the problem.

First consider the update equations for $\\{h_{c}\\}_{c\in\mathcal{C}}$:

$\displaystyle h_{c}$
$\displaystyle=\frac{1}{d}\sum_{l}\text{tr}\left[\left(\hat{h}_{l}\mu_{l}\mu_{c}^{\top}+\hat{m}_{l}\Sigma_{l}\bm{\theta}_{0}\mu_{c}^{\top}\right)\left(\lambda
I_{d}+\sum_{l^{\prime}}\hat{V}_{l^{\prime}}\Sigma_{l^{\prime}}\right)^{-1}\right]$
(94)
$\displaystyle\underset{a+c}{\simeq}\quad\frac{1}{d}\sum_{l}\text{tr}\left[\left(\hat{h}_{l}\mu_{l}\mu_{c}^{\top}\right)\left(\lambda
I_{d}+\sum_{l^{\prime}}\hat{V}_{l^{\prime}}\Sigma_{l^{\prime}}\right)^{-1}\right]$
(95)

By continuity of the saddle-point equations, if at the fixed point
$\\{\hat{h}_{c}^{\star}\to 0\\}_{c\in\mathcal{C}}$ holds, we also easily have
that $\\{h_{c}^{\star}\to 0\\}_{c\in\mathcal{C}}$. Now assume
$\\{h_{c}^{\star}=0\\}_{c\in\mathcal{C}}$ at the fixed point. We exploit
symmetry argument inherent to the update functions to show that
$\\{\hat{h}_{c}^{\star}=0\\}_{c\in\mathcal{C}}$ under weak assumptions on the
teacher and loss functions. We write the updates using assumption c):

$\displaystyle\hat{h}_{c}$
$\displaystyle=\alpha\mathbb{E}_{\xi\sim\mathcal{N}(0,1)}\left[\int\mathrm{d}y~{}\mathcal{Z}_{0}\left(y,\frac{m_{c}}{\sqrt{q_{c}}}\xi+\pi_{c},\rho_{c}-\frac{m_{c}^{2}}{q_{c}}\right)f_{\ell}(y,h_{c}+\sqrt{q_{c}}\xi,V)\right]$
(96) $\displaystyle\underset{h_{c}\to
0}{\simeq}\quad\alpha\mathbb{E}_{\xi\sim\mathcal{N}(0,1)}\left[\int\mathrm{d}y~{}\mathcal{Z}_{0}\left(y,\frac{m_{c}}{\sqrt{q_{c}}}\xi,\rho_{c}-\frac{m_{c}^{2}}{q_{c}}\right)f_{\ell}(y,\sqrt{q_{c}}\xi,V_{c})\right]$
(97)

Reminding the definition of the teacher measure term:

$\displaystyle
Z_{0}(y,\omega,V)=\mathbb{E}_{\lambda\sim\mathcal{N}(\omega,V)}\left[P_{0}(y|\lambda)\right]$
$\displaystyle\lambda\equiv\frac{\vec{x}^{\top}\vec{\theta}_{0}}{\sqrt{d}}$
(98)

we see that the symmetry conditions impose basic restrictions on the label
generation:

$\displaystyle
f)\quad\mathbb{E}_{\lambda\sim\mathcal{N}(\omega,V)}\left[P_{0}(y|\lambda)\right]=\mathbb{E}_{\lambda\sim\mathcal{N}(-\omega,V)}\left[P_{0}(-y|\lambda)\right]$
(99)

and hence $\mathcal{Z}_{0}$ is even in its second argument. Additionally, we
have

$\displaystyle
f_{\ell}(y,\omega,V)=\frac{1}{V}(\operatorname{prox}\,_{Vl(y,\cdot)}(\omega)-\omega)$
$\displaystyle\operatorname{prox}\,_{V\ell(y,\cdot)}(w)=\operatorname*{arg\,min}_{z}\left[\frac{1}{2V}(z-w)^{2}+\ell(y,z)\right]$
(100)

The symmetry condition on $\ell$ then implies that

$\displaystyle\operatorname{prox}\,_{V\ell(-y,\cdot)}(-\omega)=-\operatorname{prox}\,_{V\ell(y,\cdot)}(\omega),$
(101)

so $f_{\ell}$ is odd in its second argument. All that’s left to notice is that
(97) is an Gaussian integral of an odd function, hence it is equal to 0.

###  B.2 Gaussian Universality : proof of Proposition 3.3

We proved that the fixed point respects
$\\{h_{c}^{\star},\hat{h}_{c}^{\star}=0,0\\}_{c\in\mathcal{C}}$ under the
hypothesis of Prop. 3.2. Assume now that the mixture is homogeneous:

$\displaystyle\Sigma_{c}=\Sigma\qquad\forall c\in\mathcal{C}$ (102)

Then we have

$\displaystyle\varepsilon^{\rm
GMM}_{gen}\left(\\{\vec{\mu}_{c}\\}_{c=1}^{K},\\{\Sigma\\}_{c=1}^{K}\right)$
$\displaystyle\simeq\varepsilon^{\rm
GMM}_{gen}\left(\vec{0},\\{\Sigma\\}_{c=1}^{K}\right)$ (103)
$\displaystyle\varepsilon^{\rm
GMM}_{tr}\left(\\{\vec{\mu}_{c}\\}_{c=1}^{K},\\{\Sigma\\}_{c=1}^{K}\right)$
$\displaystyle\simeq\varepsilon^{\rm
GMM}_{tr}\left(\vec{0},\\{\Sigma\\}_{c=1}^{K}\right)$ (104)

But the right-hand side of those equations corresponds to a distribution of
the form

$\bm{x}\sim\sum_{c\in\mathcal{C}}\mathcal{N}(\bm{0},\Sigma)=\mathcal{N}(\bm{0},\Sigma)!$

This proves Proposition 3.3

###  B.3 Covariance universality

Surprisingly in the limit of vanishing regularization, $\lambda\to 0^{+}$, the
covariance is not relevant for the high dimensional learning problem. We show
that when a unique minimizer of the loss exists, and we can take safely the
limit $\lambda\to 0^{+}$ in the saddle point equations, the covariance
$\Sigma$ disappears completely from the overlap expression. Indeed assuming
the minimizer $\bm{\hat{\theta}}$ is unique, we can safely simplify expression
of the type:

$\displaystyle\lim_{\lambda\to 0^{+}}\Sigma^{n}(\lambda
I_{d}+\hat{V}\Sigma)^{n^{\prime}}=\frac{1}{\hat{V}}\Sigma^{n-n^{\prime}}$
(105)

Then by plugging in the $\lambda\to 0^{+}$ simplification in eqs. (80), (81)
we have complete independence from the covariance matrix $\Sigma$:

$\displaystyle\begin{cases}\hat{V}&=-\alpha\mathbb{E}_{\xi\sim\mathcal{N}(0,1)}\left[\int\mathrm{d}y~{}\mathcal{Z}_{0}\left(y,\frac{m_{c}}{\sqrt{q_{c}}}\xi,\rho_{c}-\frac{m_{c}^{2}}{q_{c}}\right)\partial_{\omega}f_{\ell}(y,\sqrt{q_{c}}\xi,V_{c})\right]\\\
\hat{q}&=\alpha\mathbb{E}_{\xi\sim\mathcal{N}(0,1)}\left[\int\mathrm{d}y~{}\mathcal{Z}_{0}\left(y,\frac{m_{c}}{\sqrt{q}}\xi,\rho-\frac{m^{2}}{q}\right)f_{\ell}(y,\sqrt{q}\xi,V)^{2}\right]\\\
\hat{m}&=\alpha\mathbb{E}_{\xi\sim\mathcal{N}(0,1)}\left[\int\mathrm{d}y~{}\partial_{\omega}\mathcal{Z}_{0}\left(y,\frac{m}{\sqrt{q}}\xi,\rho-\frac{m^{2}}{q}\right)f_{\ell}(y,\sqrt{q}\xi,V)\right]\end{cases}$
(106) $\displaystyle\begin{cases}V&=1\\\ q&=\hat{q}+\rho\hat{m}^{2}\\\
m&=\hat{m}\rho\end{cases}$ (107)

which concludes the proof.

###  B.4 Strong universality

We never made any specific assumption on the loss up to now, apart the very
general symmetry condition in Assumption 2. In this section we show strong
universality for square loss regression for any GMM estimation problem
respecting the mean universality property (Prop. 3.2). We assume to consider a
ridge regression problem in the underparametrized regime $\alpha>1$:

$\displaystyle P_{0}(\tau|\Delta)=\mathcal{N}(\tau,\Delta)$
$\displaystyle\ell(y,\hat{y})=(y-\hat{y})^{2}$ (108)

the replicas then correspond to equations (87), (88). In the $\lambda\to
0^{+}$ limit, the equation simplify greatly:

$\displaystyle\begin{cases}\hat{V}_{c}&=\frac{\alpha p_{c}}{1+V_{c}}\\\
\hat{q}_{c}&=\frac{\alpha
p_{c}}{(1+V_{c})^{2}}(\rho_{c}+\Delta+q_{c}-2m_{c})\\\
\hat{m}_{c}&=\frac{\alpha p_{c}}{1+V_{c}}\end{cases}$ (109)
$\displaystyle\begin{cases}V_{c}&=\frac{1}{d}\operatorname{tr}\left[\Sigma_{c}(\sum_{l\in\mathcal{C}}\hat{V}_{l}\Sigma_{l})^{-1}\right]\\\
q_{c}&=\frac{1}{d}\sum_{l\in\mathcal{C}}^{K}\operatorname{tr}\left[\hat{q}_{l}\Sigma_{l}\Sigma_{c}(\sum_{l^{\prime}\in\mathcal{C}}\hat{V}_{l^{\prime}}\Sigma_{l^{\prime}})^{-2}\right]+\frac{1}{d}\sum_{(l,l)^{\prime}\in\mathcal{C}\times\mathcal{C}}\operatorname{tr}\left[\hat{m}_{l}\hat{m}_{l^{\prime}}\Sigma_{l}\bm{\theta}_{0}\vec{\theta}_{0}^{\top}\Sigma_{l^{\prime}}\Sigma_{c}(\sum_{l^{\prime\prime}\in\mathcal{C}}\hat{V}_{l^{\prime\prime}}\Sigma_{l^{\prime\prime}})^{-2}\right]\\\
m_{c}&=\frac{1}{d}\sum_{l=1}^{K}\operatorname{tr}\left[\hat{m}_{l}\Sigma_{l}\bm{\theta}_{0}\bm{\theta}_{0}^{T}\Sigma_{c}(\sum_{l^{\prime}}\hat{V}_{l^{\prime}}\Sigma_{l^{\prime}})^{-1}\right]\end{cases}$
(110)

Consider the equation for the overlaps $\\{q_{c}\\}_{c\in\mathcal{C}}$:

$\displaystyle q_{c}$
$\displaystyle=\frac{1}{d}\sum_{l=1}^{K}\operatorname{tr}\left[\hat{q}_{l}\Sigma_{l}\Sigma_{c}(\sum_{l^{\prime}}\hat{V}_{l^{\prime}}\Sigma_{l^{\prime}})^{-2}\right]+\frac{1}{d}\sum_{l,k=1}^{K}\operatorname{tr}\left[\hat{m}_{l}\hat{m}_{k}\Sigma_{l}\vec{\theta}_{0}\bm{\theta}_{0}^{T}\Sigma_{k}\Sigma_{c}(\sum_{l^{\prime}}\hat{V}_{l^{\prime}}\Sigma_{l^{\prime}})^{-2}\right]$
(111) $\displaystyle\equiv R_{c}+G_{c}$ (112)

The error metrics can be decomposed as:

$\displaystyle\varepsilon_{\text{tr}}$
$\displaystyle=\sum_{c\in\mathcal{C}}\frac{p_{c}}{(1+V_{c})^{2}}(R_{c}+\Delta)+\sum_{c\in\mathcal{C}}\frac{p_{c}}{(1+V_{c})^{2}}(\rho_{c}+G_{c}-2m_{c})$
(113) $\displaystyle\varepsilon_{\text{gen}}$
$\displaystyle=\sum_{c\in\mathcal{C}}p_{c}(R_{c}+\Delta)+\sum_{c\in\mathcal{C}}p_{c}(\rho_{c}+G_{c}-2m_{c})$
(114)

We focus on the second term and show that it is equal to zero. Looking at eqs.
(110),(109), we note that:

$\displaystyle\hat{m}_{c}=\hat{V}_{c}\qquad\forall c\in\mathcal{C}$ (115)

By using eq. (115) we can simplify the equations for
$\\{\hat{m}_{c},G_{c}\\}$:

$\displaystyle G_{c}=\rho_{c},\qquad m_{c}=\rho_{c}\qquad\forall
c\in\mathcal{C}$ (116)

Plugging this relations in the saddle point equations we obtain:

$\displaystyle\begin{cases}\hat{V}_{c}&=\frac{\alpha
p_{c}}{1+V_{c}}=\hat{m}_{c}\\\ \hat{q}_{c}&=\frac{\alpha
p_{c}}{(1+V_{c})^{2}}(\Delta+R_{c})\end{cases}$ (117)
$\displaystyle\begin{cases}V_{c}&=\frac{1}{d}\operatorname{tr}\left[\Sigma_{c}(\sum_{l\in\mathcal{C}}\hat{V}_{l}\Sigma_{l})^{-1}\right]\\\
q_{c}&=\frac{1}{d}\sum_{l=1}^{K}\operatorname{tr}\left[\hat{q}_{l}\Sigma_{l}\Sigma_{c}(\sum_{l^{\prime}}\hat{V}_{l^{\prime}}\Sigma_{l^{\prime}})^{-2}\right]+\rho_{c}\equiv
R_{c}+\rho_{c}\\\ m_{c}&=\rho_{c}\end{cases}$ (118)

The fixed point of the equations does not depend on
$\\{\rho_{c}\\}_{\in\mathcal{C}}$ anymore and we are left with equations only
for $\\{\hat{V}_{c},\hat{q}_{c},V_{c},R_{c}\\}_{c\in\mathcal{C}}$:

$\displaystyle\begin{cases}\hat{V}_{c}&=\frac{\alpha p_{c}}{1+V_{c}}\\\
\hat{q}_{c}&=\frac{\alpha p_{c}}{(1+V_{c})^{2}}(\Delta+R_{c})\end{cases}$
(119)
$\displaystyle\begin{cases}V_{c}&=\frac{1}{d}\operatorname{tr}\left[\Sigma_{c}(\sum_{l\in\mathcal{C}}\hat{V}_{l}\Sigma_{l})^{-1}\right]\\\
R_{c}&=\frac{1}{d}\sum_{l=1}^{K}\operatorname{tr}\left[\hat{q}_{l}\Sigma_{l}\Sigma_{c}(\sum_{l^{\prime}}\hat{V}_{l^{\prime}}\Sigma_{l^{\prime}})^{-2}\right]\end{cases}$
(120)

The errors can be computed in a closed form with a series of algebraic
manipulation. One can verify by algebraic manipulation the following
relations:

$\displaystyle\sum_{c}V_{c}\hat{V}_{c}$
$\displaystyle=\sum_{c}\frac{1}{d}\operatorname{tr}\left[\hat{V}_{c}\Sigma_{c}(\sum_{l\in\mathcal{C}}\hat{V}_{l}\Sigma_{l})^{-1}\right]=1$
(121) $\displaystyle\sum_{c}R_{c}\hat{V}_{c}-V_{c}\hat{q}_{c}$
$\displaystyle=\frac{1}{d}\sum_{l=1}^{K}\operatorname{tr}\left[\hat{q}_{l}\Sigma_{l}(\sum_{c}\hat{V}_{c}\Sigma_{c})(\sum_{l^{\prime}}\hat{V}_{l^{\prime}}\Sigma_{l^{\prime}})^{-2}\right]-\frac{1}{d}\operatorname{tr}\left[\sum_{c}\hat{q}_{c}\Sigma_{c}(\sum_{l\in\mathcal{C}}\hat{V}_{l}\Sigma_{l})^{-1}\right]=0$
(122)

Now we express everything in eq. (121) in terms of the non-hatted overlaps to
obtain:

$\displaystyle\star)\,\sum_{c}\frac{\alpha p_{c}}{1+V_{c}}V_{c}=1$ (123)
$\displaystyle\\#)\,\sum_{c}\frac{\alpha
p_{c}}{(1+V_{c})^{2}}R_{c}=\sum_{c}\frac{\alpha
p_{c}}{(1+V_{c})^{2}}V_{c}\Delta$ (124)

We remark that we can write the training and generalization error can be
written as:

$\displaystyle\varepsilon_{\text{tr}}=\sum_{c}\frac{\hat{q}_{c}}{\alpha}$
$\displaystyle\varepsilon_{\text{tr}}=\sum_{c}\frac{\hat{q}^{\star}_{c}(1+V_{c}^{\star})^{2}}{\alpha}$
(125)

So if we compute the quantity $\sum_{c}\hat{q}_{c}$ we conclude. By plugging
in $(\\#)$ into the expression above we obtain:

$\displaystyle\sum_{c}\hat{q}_{c}$ $\displaystyle=\sum_{c}\frac{\alpha
p_{c}}{(1+V_{c})^{2}}(R_{c}+\Delta)$ (126)
$\displaystyle\underset{(\\#)}{=}\Delta\sum_{c}\frac{\alpha
p_{c}}{(1+V_{c})^{2}}+\Delta\sum_{c}\frac{\alpha p_{c}}{(1+V_{c})^{2}}V_{c}$
(127) $\displaystyle=\Delta\sum_{c}\frac{\alpha p_{c}}{1+V_{c}}$ (128)

and finally using relation $(\star)$ we prove the theorem:

$\displaystyle\varepsilon_{\text{tr}}$
$\displaystyle=\sum_{c}\frac{\hat{q}_{c}^{\star}}{\alpha}$ (129)
$\displaystyle=\sum_{c}\frac{\hat{q}_{c}^{\star}}{\alpha}\pm\frac{\Delta}{\alpha}$
(130) $\displaystyle\underset{(\star)}{=}\Delta(1-\frac{1}{\alpha})$ (131)

On the other hand, the generalization error does not respect the strong
universlaity statement as we analyze in the next section (See Fig. 5).

##  Appendix C Non-universality of Gaussian Mixture Models

In the previous section we enumerated a series of results unveiling
universality of GMM. Now we want to study the dual task: when GMM model are
not universal? We have two ways to break universality: a) allow strong
heterogeneity in the data structure; b) consider labels which are strongly
correlated with the data structure. The plan for this section is to review
more in detail these two processes for universality breaking.

###  C.1 Strongly heterogeneous mixtures

We proved in Theorem. 3.6 a strong universality statement for the training
loss of ridge regression. However in this section we want to clarify that the
theorem is not valid beyond its assumption for general mixtures. We prove this
by analyzing a counterexample: a strongly heterogeneous 2-clusters Gaussian
Mixture:

$\displaystyle\Sigma_{+}$
$\displaystyle=\text{diag}(0.1,\dots,0.1,1.9,\dots,1.9)\qquad p_{+}=0.8$ (132)
$\displaystyle\Sigma_{-}$
$\displaystyle=\text{diag}(1.9,\dots,1.9,0.1,\dots,0.1)\qquad p_{-}=0.2$ (133)

We present in Fig. 5 the comparison of the heterogeneous GMM performance with
the Gaussian theory: although the training errors coincide as predicted by
Theorem. 3.6, the generalization errors are different. However, we remark that
real data after preprocessing seem homogeneous enough to obtain a good
agreement with the exact Gaussian asymtptotics as we see in Fig. 1.

Figure 5: Training and generalization error for ridge regression on the GMM
defined in eq. (132). Left: We compare the performance of the Gaussian
asymtptotics (solid black line) with the GMM one (solid red line), while the
orange dots represents simulations which agrees with the theoretical
predictionas predicted by Theorem. 3.1. Right: Two dimensional toy plot of
$10000$ samples coming from the heterogeneous Gaussian mixture.

###  C.2 Correlated teachers

In this section we investigate in deeper detail the controlled setting of
Theorem. 3.7. We do not necessarily consider estimation problems which
respects the assumption of mean universality property (Prop. 3.2), and in fact
we precisely analyze the consequences if we relax these conditions. We perform
ridge regression on a two-mixtures GMM with opposite means and same covariance
matrix proportional to $I_{d}$. In this scenario the exact asymptotics for the
performance of the ERM estimator admits a closed form expression. We remind
here the main parameters for the theoretical analysis are:

$\displaystyle\rho=\frac{1}{d}\vec{\theta}_{0}^{\top}\Sigma\vec{\theta}_{0}=\frac{1}{d}\lVert\vec{\theta}_{0}\rVert_{2}^{2}\qquad\gamma=\frac{1}{d}||\vec{\mu}||_{2}^{2}\qquad\Delta=\mathbb{E}_{\xi}\xi^{2}\qquad\bm{\hat{\mu}_{\pm}}=\hat{h}_{\pm}\bm{\mu}_{\pm}+\hat{m}_{\pm}\Sigma\bm{\theta}_{0}$
(134)

In this simple setting we can simplify the general mixture equations in eqs.
(80),(81) without doing any assumption on the teacher. We first map the
equation to a single Gaussian problem
$\mathcal{G}_{\vec{\theta},I_{d},\vec{\mu}}$. We have to slightly modify the
proof of Prop. 3.3 presented in Sec. B.2, indeed we cannot blindly assume that
the mean overlaps $\\{h_{+},h_{-}\\}$ are zero but one can show that at the
fixed point the overlaps respect:

$\displaystyle\hat{h}_{+}^{\star}=-\hat{h}_{-}^{\star},\qquad\hat{m}_{+}^{\star}=\hat{m}_{-}^{\star}\qquad\to\qquad\vec{\hat{\mu}}_{+}^{\star}=\vec{\hat{\mu}}_{-}^{\star}$
(135) $\displaystyle\hat{V}=\hat{V}_{+}+\hat{V}_{-}$
$\displaystyle\hat{q}=\hat{q}_{+}+\hat{q}_{-}$
$\displaystyle\hat{m}=\hat{m}_{+}+\hat{m}_{-}$
$\displaystyle\hat{h}=\hat{h}_{+}-\hat{h}_{-}$ (136) $\displaystyle
V=V_{+}=V_{-}$ $\displaystyle q=q_{+}=q_{-}$ $\displaystyle m=m_{+}=m_{-}$
$\displaystyle h=h_{+}=-h_{-}$ (137)

the fixed point of the replica equations for the mixture defined above are
mapped to the one of a single Gaussian problem for the overlaps
$\\{V,m,q,h,\hat{V},\hat{m},\hat{q},\hat{h}\\}$. Now we can simplify them even
further by plugging in the assumption on the covariance, all the traces
simplify and we get:

$\begin{cases}V&=\frac{1}{\lambda+\hat{V}}\\\
q&=V^{2}\left(\hat{q}+\rho\hat{m}^{2}+2\pi\hat{m}\hat{h}+\gamma\hat{h}^{2}\right)\\\
m&=V\left(\rho\hat{m}+\pi\hat{h}\right)\\\
h&=V\left(\gamma\hat{h}+\pi\hat{m}\right)\end{cases}$ (138)

These equations are actually solvable! Define

$\eta:=V\hat{V}=\frac{\alpha V}{1+V};$

Now, we know that the generalization error satisfies:

$V^{2}\hat{q}=\frac{V^{2}\alpha}{(1+V)^{2}}\varepsilon_{\text{gen}}=\frac{\eta^{2}}{\alpha}(\varepsilon_{\text{gen}}).$

We can plug this into the equation for $q$ to get

$\displaystyle\varepsilon_{\text{gen}}$
$\displaystyle=\rho+\underbrace{\left(\frac{\eta^{2}}{\alpha}(\varepsilon_{\text{gen}})+\rho(V\hat{m})^{2}+2\pi(V\hat{m})(V\hat{h})+\gamma(V\hat{h})^{2}\right)}_{q}-2\underbrace{\left(\rho
V\hat{m}+\pi V\hat{h}\right)}_{m}+(\pi-h)^{2}$ (139)
$\displaystyle\varepsilon_{\text{tr}}$
$\displaystyle=\frac{\varepsilon_{\text{gen}}}{(1+V)^{2}}$ (140)

with the relations:

$V=\frac{\eta}{\alpha-\eta},\qquad V\hat{m}=\eta,\quad V\hat{h}=\eta(\pi-h)$

and the expression of $\pi-h$ as a function of $\eta$:

$\displaystyle\pi-h=\pi-V(\gamma\hat{h}+\pi\hat{m})=\pi-\gamma\eta(\pi-h)-\pi\eta\iff\pi-h=\pi\frac{1-\eta}{1+\gamma\eta}$
(141)

we simplify everything in terms only of $\eta$ to get:

$\displaystyle\varepsilon_{\text{gen}}$
$\displaystyle=\frac{\alpha\left(\Delta+(\eta-1)^{2}\rho\right)}{\left(\alpha-\eta^{2}\right)}\left(1-\pi^{2}\frac{(\eta-1)^{2}\left(\gamma\eta^{2}+2\eta-1\right)}{\left(\Delta+(\eta-1)^{2}\rho\right)(1+\gamma\eta)^{2}}\right)$
(142) $\displaystyle\varepsilon_{\text{tr}}$
$\displaystyle=\frac{(\alpha-\eta)^{2}\left(\Delta+(\eta-1)^{2}\rho\right)}{\alpha\left(\alpha-\eta^{2}\right)}\left(1-\pi^{2}\frac{(\eta-1)^{2}\left(\gamma\eta^{2}+2\eta-1\right)}{\left(\Delta+(\eta-1)^{2}\rho\right)(1+\gamma\eta)^{2}}\right)$
(143)

All that remains is to solve for $\eta$ using the equations for $V$ and
$\hat{V}$, which yields

$\eta=1+\frac{1}{2}\left(\alpha-1+\lambda-\sqrt{4\lambda+(\alpha-1+\lambda)^{2}}\right).$
(144)

#### Vanishing regularization

We can take $\lambda=0$ inside (144) even when $\alpha<1$, and we find

$\displaystyle\eta(\lambda=0)=\min(\alpha,1)$ (145)

Finally, this yields

$\displaystyle\varepsilon_{gen}$
$\displaystyle=\begin{cases}\frac{\alpha\Delta}{\alpha-1}\qquad&\alpha\geq
1\\\
\frac{\pi^{2}(\alpha-1)\left(\alpha^{2}\gamma+2\alpha-1\right)}{(\alpha\gamma+1)^{2}}-\frac{(\alpha-1)^{2}\rho+\Delta}{\alpha-1}\qquad&\text{else}\end{cases}$
(146) $\displaystyle\varepsilon_{\text{train}}$
$\displaystyle=\begin{cases}\frac{(\alpha-1)\Delta}{\alpha}\qquad&\alpha\geq
1\\\ 0\qquad&\text{else}\\\ \end{cases}$ (147)

We retrieve the results of Corollary. 3.8: even with correlated teachers we
show that in the underparametrized regime we have Gaussian universality.

#### Extension of Corollary. 3.8:

The Gaussian universality result for correlated teacher at vanishing
regularization can be extended as well to general balanced mixtures with
homoscedastic covariance, i.e. the case where

$p_{c}=\frac{1}{|\mathcal{C}|},\quad\Sigma_{c}=\Sigma,\quad\lambda=0.$

In order to do so, we consider the general saddle point equation for ridge
regression in eq. (87),(88), and plug in the assumptions:

$\displaystyle\begin{cases}\hat{V}&=\frac{\alpha}{|\mathcal{C}|(1+V)}\\\
\hat{q}_{c}&=\frac{\alpha}{|\mathcal{C}|(1+V)^{2}}(\rho+\Delta+q_{c}-2m+(h_{c}-\pi_{c})^{2})\\\
\hat{m}&=\frac{\alpha}{|\mathcal{C}|(1+V)}\\\
\hat{h}_{c}&=\frac{\alpha(\pi_{c}-h_{c})}{|\mathcal{C}|(1+V)}\end{cases}$
(148)
$\displaystyle\begin{cases}V&=\frac{1}{d}\operatorname{tr}\left[\Sigma\hat{\Sigma}^{-1}\right]\\\
\vskip 1.79997pt\cr
q_{c}&=\frac{1}{d}\operatorname{tr}\left[\left(\hat{q}_{c}\Sigma+\bm{\hat{\mu}}\bm{\hat{\mu}}^{\top}\right)\Sigma\hat{\Sigma}^{-2}\right]\\\
\vskip 1.79997pt\cr
m&=\frac{1}{d}\operatorname{tr}\left[\bm{\hat{\mu}}\bm{\theta}_{0}^{\top}\Sigma\hat{\Sigma}^{-1}\right]\\\
\vskip 1.79997pt\cr
h_{c}&=\frac{1}{d}\operatorname{tr}\left[\bm{\hat{\mu}}\bm{\mu}_{c}^{\top}\hat{\Sigma}^{-1}\right]\end{cases}$
(149)

In particular, these assumptions imply that
$V_{c},\hat{V}_{c},m_{c},\hat{m}_{c}$ do not depend on the class labels $c$,
and hence

$\hat{\Sigma}=\hat{V}\Sigma.$

Now, consider the assignment

$\hat{h}_{c}=0\quad\text{and}\quad h_{c}=\pi_{c};$ (150)

is is easy to check that they are a fixed point of the above saddle-point
equations, since then

$h_{c}=\frac{\hat{m}\pi_{c}}{\hat{V}}=\pi_{c}$

Plugging in this relation in the expression of the update for
$\\{\hat{q}_{c}\\}_{c\in\mathcal{C}}$ we obtain:

$\displaystyle\hat{q}_{c}^{\star}=\frac{\alpha
p_{c}}{(1+V_{c}^{\star})^{2}}(\rho_{c}+\Delta+q_{c}^{\star}-2m_{c}^{\star}+(h_{c}^{\star}-\pi_{c})^{2})=\frac{\alpha
p_{c}}{(1+V_{c}^{\star})^{2}}(\rho_{c}+\Delta+q_{c}^{\star}-2m_{c}^{\star})$
(151)

Hence by looking at the expression of generalization and training error for
ridge regression in eq. (125) we conclude that they will not depend on the
values of the mean overlaps $\\{h_{c}^{\star}\\}_{c\in\mathcal{C}}$ and we can
prove Gaussian universality in the same way as we did in the proof of Theorem.
3.4.

Figure 6: Theoretical prediction from eq. (158) for the difference between the
generalization error computed from the Gaussian theory and GMM one plotted as
a function of the number of samples used in the fit in Algorithm. 1. We use
different regularization strength $\Lambda\in\\{1e-10,0.1,1,10\\}$, increasing
clockwise in the figure.

###  C.3 Interpolating teachers

Different works observed over different datasets that if we find an
interpolating teacher vector, in such a way that we can keep the real labels
to study ERM performance, the theoretical predictions coming from Gaussian
asymptotics would agree with simulations, [19] among others. This observation
apparently contradicts the result of Theorem. 3.7 in which we show that
correlated labels (as we expect real labels to be) break Gaussian
universality. We try to motivate this in the same setting as Theorem. 3.7:

######  Theorem C.1.

(Real labels universality) Consider the same setting of Theorem. 3.7. Fix now
the teacher vector to be the maximally correlated one $\vec{\theta}=\vec{\mu}$
and find an interpolating teacher $\tilde{\vec{\theta}}$ using $\tilde{n}$
samples coming from the GMM. If we compare the exact asymptotics of Gaussian
and GMM for regularization $\Lambda$ at fixed $\alpha=\frac{\tilde{n}}{d}$ we
obtain a discrepancy:

$\displaystyle\Delta\varepsilon_{\text{gen}}=\frac{4(E-1)\alpha^{2}\left(E^{2}\gamma+2e-1\right)}{(\alpha+1)^{2}(E\gamma+1)^{2}}$
(152)

with E($\alpha$,$\Lambda$) solution of:

$\displaystyle
E=1+\frac{1}{2}\left(\alpha-1+\Lambda-\sqrt{4\Lambda+(\alpha-1+\Lambda)^{2}}\right).$
(153)

The result above solve the apparent contradiction: if we strongly
overparametrize the fitting model, interpolating teachers can be uncorrelated
with the data structure.

Algorithm 1 Get theoretical learning curves keeping real labels

Data generation: Generate the data from the 2-clusters GMM described in C.1
and set:
$\displaystyle\vec{\mu}_{+}=-\vec{\mu}_{-}=\sqrt{\gamma}(1,\dots,1)^{\top},\qquad\Sigma_{+}=\Sigma_{-}=I_{d},\qquad\vec{\theta}_{0}=\vec{\mu}_{+}$
(154)

Fit interpolating teacher Perform ridge regression for $\lambda\to 0$. Take
$\alpha=\frac{n}{d}<1$ so that the ERM estimator will interpolate the data.

Run exact asymptotics Compute the theoretical performance of the associated
2-GMM and equivalent GCM from the estimated $\hat{\theta}$ for the same
$\alpha$ with a fixed regularization parameter $\Lambda$.

The proof of the results relies strongly on the fact that we can analyze
analytically the performance in this scenario being it the same as the
previous section. Reminding from eq. (145) that in the limit $\lambda\to 0$ we
have in the overparametrized setting $\eta(\alpha,\lambda=0)=\alpha$, we can
write:

$\displaystyle h_{\textit{erm}}$
$\displaystyle=\frac{\vec{\mu}^{\top}\hat{\vec{\theta}}}{d}=\gamma\alpha\frac{\gamma+1}{1+\gamma\alpha}$
(155) $\displaystyle q_{\textit{erm}}$
$\displaystyle=\frac{\hat{\vec{\theta}}^{\top}\Sigma\hat{\vec{\theta}}}{d}=\alpha\left(-\frac{\Delta}{\alpha-1}-\frac{(\alpha-1)\pi^{2}}{\alpha\gamma+1}+\rho\right)$
(156)

We summarize the algorithmic procedure in in Algorithm. 1. In order to compute
the generalization error from eq. (142) we need the quantities
$\\{\rho,\pi\\}$. Note that in Algorithm. 1 we compute the theoretical curve
for the GMM and the GCM from the estimated $\hat{\vec{\theta}}$. This
translates into the fact that the estimated overlaps from numerical
simulations in the previous step become the equivalent to:

$\displaystyle\pi=h_{erm}$ $\displaystyle\rho=q_{erm}$ (157)

We fix the regularization for the theoretical prediction to be $\Lambda$, and
for simplicity we use the same sample complexity parameter $\alpha$. We are
now in a position to compare Gaussian and GMM theoretical prediction, more
precisely we analyze the difference of the generalization errors for the two
data model. Recalling the expression for the generalization error in eq.
(142), we just need to plug the estimated overalps in eq. (157) to obtain:

$\displaystyle\Delta\varepsilon_{\text{gen}}=\frac{4\alpha^{2}(E-1)\left(E^{2}\gamma+2e-1\right)}{(\alpha+1)^{2}(E\gamma+1)^{2}}$
(158)

where $E(\alpha,\Lambda)$ is the solution of eq. (144), retrieving the result
in Theorem. C.1. For the sake of clarity we stated the theorem in a simple
setting, indeed once we estimated $\hat{\vec{\theta}}$ we could have decided
to change the value of the sample complexity, as we are now interested in
running theory curves. However, to simplify the equations we assumed to run
the theoretical prediction for the same value of $\alpha$ as in the ERM fit.
We present the results in Fig. 6 for $\gamma=1$: although the true labels are
maximally correlated with the data structure, we see that the Gaussian
theoretical predictions with the interpolating teacher $\vec{\hat{\theta}}$
are very close to the GMM one. We remark that the upper-left plot in Fig. 1 is
a nice characterization of Theorem. (3.8): as we reach the underparametrized
regime we restore Gaussian universality for $\Lambda\to 0^{+}$ even for
correlated teachers.

##  Appendix D Details on real dataset simulations

In this section we report the procedure to create the random regression task
on real data, see Algorithm. 2. The implementation of the different numerical
simulations described in this work are available in a GitHub respository.

Algorithm 2 Random teacher regression on real data

Data processing: Load data and perform the preprocessing step with a transform
matrix $F\in\mathbb{R}^{d\times d^{\prime}}$, with $d^{\prime}$ dimension of
the images in the chosen dataset. We choose for all the figures in this
manuscript $d=2000$.

Match covariance Compute $\hat{\Sigma}=\frac{1}{n}XX^{\top}$

Create new labels Forget the real labels associated with the dataset and
create new label according to:  $\displaystyle y_{\nu}$
$\displaystyle=\begin{cases}&\frac{\vec{\theta}_{0}^{\top}\vec{x}_{\nu}}{\sqrt{d}}+\sqrt{\Delta}\xi\qquad\text{Ridge
regression}\\\
&\rm{sign}\left(\frac{\vec{\theta}_{0}^{\top}\vec{x}_{\nu}}{\sqrt{d}}+\sqrt{\Delta}\xi\right)\qquad\text{Logistic
regression}\end{cases}$ (159) $\displaystyle\vec{\theta}_{0}$
$\displaystyle\sim\mathcal{N}(\vec{0},I_{d})\qquad\xi\sim\mathcal{N}(0,1)$
(160)

Run learning curves Fix the regularization parameter $\lambda$.

for $\alpha$ in a given range do

Simulation Solve the ERM in eq. (1) using sklearn package LogisticRegression
[53] or the Moore-Penrose pseudo-inverse for the ridge estimator:
$\displaystyle\hat{\vec{\theta}}$
$\displaystyle=\begin{cases}&X^{\top}(XX^{\top}+\lambda I_{d})^{-1}Y\qquad
n<d\\\ &(X^{\top}X+\lambda I_{d})^{-1}X^{\top}y\qquad\qquad n>d\end{cases}$
(161)

Gaussian theory Solve saddle point equations in eqs. (80),(81) defining
Gaussian asymptotics for the model. 2.1.

Compute erros Compute the training loss and generalization error using the
metrics:  $\displaystyle g(y,\hat{y})$
$\displaystyle=\begin{cases}(y-\hat{y})^{2}\qquad\text{Ridge regression}\\\
\mathbb{P}(y\neq\hat{y})\qquad\text{Logistic regression}\end{cases}$ (162)

end for