The role of regularization in classification of high-dimensional noisy Gaussian mixture

I Introduction

where both ${\bf z}_{i}$ and ${\bf v^{*}}$ have components taken in ${\cal N}(0,1)$. The labels $y_{i}\in{\pm 1}$ are generated randomly with a fraction $\rho$ of $+1$ (and $1-\rho$ of $-1$). We focus on the high-dimensional limit where $n,d\!\to\!\infty$ while $\alpha=n/d$, $\rho$ and $\Delta$ are fixed. The factor $\sqrt{d}$ in (1) is such that a classification better than random is possible, yet even the oracle-classifier that knows exactly the centroid $\frac{{\bf v}^{*}}{\sqrt{d}}$ only achieves a classification error bounded away from zero. We focus on ridge regularized learning performed by the empirical risk minimization of the loss:

The unsupervised version of the problem is the standard Gaussian mixture modeling problem in statistics . For the supervised model considered here, recently computed rigorously the performance of the Bayes-optimal estimator (that knows the generative model of the data, but does not have access to the vector ${\bf v}^{*}$) for the case of equally sized clusters. We generalize these results for arbitrary cluster sizes to provide a baseline for the estimators obtained by empirical risk minimization.

The model was recently under investigation in a number of papers. In , the authors study the same data generative model in the particular case of equally sized clusters, and analyze non-regularized losses under the assumption that the data are not linearly separable. They conclude that in that case the square loss is a universally optimal loss function. Our study of the regularized losses shows that the performance of the non-regularized square loss can be easily, and drastically improved. studied the logistic loss, again without regularization and for two clusters of equal size, and derive the linear separability condition in this case.

Secondly, we present a systematic investigation of the effects of regularization and of the cluster size, discussing in particular how far estimators obtained by empirical risk minimization fall short of Bayes-optimal one, with surprising conclusions where we illustrate the effect of strong and weak regularizations. In particular, when data are linearly separable, Rosset et al. proves that all monotone non-increasing loss functions depending on the margin find a solution maximizing the margin. This is indeed exemplified in our model by the fact that for $\alpha<\alpha^{*}(\Delta,\rho)$ (the location of transition for linear separability) the hinge, and logistic losses converge to the same test error as the regularization tends to zero. This is related to the implicit regularization of gradient descent for the non-regularized minimization , and we discuss this in connection with the “double-descent” phenomenon that is currently the subject of intense studies .

The existence of a sharp transition for perfect separability in the model, with and without bias, is interesting in itself. Recently analyzed the maximum likelihood estimate (MLE) in high-dimensional logistic regression. While they analyzed Gaussian data (whereas we study Gaussian mixture) their results on the existence of the MLE being related to the separability of the data and having a sharp phase transition are of the same nature as ours, and similar to earlier works in statistical physics .

Finally, we note that the formulas proven here can also be obtained from the heuristic replica theory from statistical physics. Indeed, a model closely related to ours was studied in this literature and our rigorous solution thus provides a further example of a rigorous proof of a result obtained by this technique.

All these results show that the Gaussian mixtures model studied here allows to discuss, illustrate, and clarify in a unified fashion many phenomena that are currently the subject of intense scrutiny in high-dimensional statistics and machine learning.

II Main theoretical results

Our first result is a rigorous analytical formula for the generalization classification error obtained by the empirical risk minimization of (2). Define $q$ as the length of the vector $\bf w$ and $m$ as its overlap with ${\bf v^{*}}$, both rescaled by the dimensionality $d$

In the high dimensional limit when $n,d\to\infty$ with a fixed ratio $\alpha=n/d$, the length $q$ and overlap $m$ of the vector ${\bf w}$ obtained by the empirical risk minimization of (2) with a convex loss converge to deterministic quantities given by the unique fixed point of the system:

where $h\sim{\cal N}(m+yb,\Delta q)$, $\rho\in(0,1)$ is the probability with which $y_{i}=1$, and $v$ is the solution of

and the bias $b$, defined in (2), is the solution of the equation

This is proven in the next section using Gordon’s minimax approach. Once the fixed point values of the overlap $m$ and length $q$ are known, then we can express the asymptotic values for the generalization error and the training loss:

In the same limit as in theorem 11, the generalization error expressed as fraction of wrong labeled instances is given by

where $Q(x)=\frac{1}{\sqrt{2\pi}}\int_{x}^{\infty}e^{-t^{2}/2}dt$ is the Gaussian tail function. The value of the training loss rescaled by the data dimension reads

The details on (12) and (13) are provided in Appendices A and C.

II.2 MLE and Bayes-optimal estimator

The maximum likelihood estimation (MLE) for the considered model corresponds to the optimization of the non-regularized logistic loss. This follows directly from the Bayes formula:

where $c=\frac{2}{\Delta}y\left(\tfrac{1}{\sqrt{d}}{\rm v}^{\top}{\rm x}+\frac{\Delta}{2}\log\frac{\rho}{1-\rho}\right)$, therefore a simple redefinition of the variables leads to the logistic cost function that turns out to be the MLE (or rather the maximum a posteriori estimator if one allows the learning of a bias to account for the possibility of different cluster sizes).

The Bayes-optimal estimator is the “best” possible one in the sense that it minimizes the number of errors for new labels. It can be computed as

where $\{\bf X,\bf y\}$ is the training set and $\bf x_{\text{new}}$ is a previously unseen data point. In the Bayes-optimal setting, the model generating the data (1) and the prior distributions ${\rm p}_{y}$, ${\rm p}_{\bf z}$, ${\rm p}_{\bf v^{*}}$ are known. Therefore, we can compute the posterior distribution in (15):

Hence, we can compute the Bayes-optimal generalization error using

where $m_{BO}=q_{\rm BO}=\tfrac{\alpha}{\Delta+\alpha}$ and $b_{BO}=\tfrac{\Delta}{2}\log\tfrac{\rho}{1-\rho}$. This formula is derived in the Appendix B. The case $\rho=1/2$ was also discussed in .

Finally, it turns out that in this problem, one can reach the performances of the Bayes-optimal estimator, usually difficult to compute, efficiently using a simple plug-in estimator akin to applying the Hebb’s rule . Consider indeed the weight vector averaged over the training samples, each multiplied by its label and rescaled by $\sqrt{d}$

It is straightforward to check that, for ${\bf\hat{w}}_{\rm Hebb}$, one has in large dimension $m=1$ and $q=(1+\tfrac{\Delta}{\alpha})$. If one further optimizes the bias (for instance by cross validation) and uses its optimal value $b=\tfrac{\Delta q}{2m}\log\tfrac{\rho}{1-\rho}$, plugging these in eq. (12) one reaches Bayes-optimal performance $\varepsilon_{\rm{gen}}^{\rm Hebb}=\varepsilon_{\rm{gen}}^{\rm BO}$. Since there exists a plug-in estimator that reaches the Bayes-optimal performance, it is particularly interesting to see how the ones obtained by empirical risk minimization compare with the optimal results.

II.3 High-Dimensional Landscapes of Training Loss

Our analysis also leads to an analytical characterization of the high-dimensional landscapes of the training loss. First, we let

to denote the normalized training loss when we restrict the weight vector to satisfy the two conditions in (21). In what follows, we refer to ${\cal L}_{\lambda}(q,m,b)$ as the “local training loss” at fixed values of $q,m$ and $b$. The “global training loss” can then be obtained as

where the constraint that $m^{2}\leq q$ is due to the Cauchy-Schwartz inequality: $\lvert m\rvert=\frac{\lvert{\bf w}^{\top}{\bf v}^{\ast}\rvert}{d}\leq\frac{\lVert{\bf w}\rVert}{\sqrt{d}}\frac{\lVert{\bf v}^{\ast}\rVert}{\sqrt{d}}=\sqrt{q}$.

In the high-dimensional limit when $n,d\to\infty$ with a fixed ratio $\alpha=n/d$, many properties of the local training loss can be characterized by a deterministic function, defined as

Here, for any $\gamma\geq 0$, $v_{\gamma}$ denotes a random variable whose cumulative distribution function is given by

Moreover, $\gamma^{\ast}$ in (23) is the unique solution to the equation

Let $\Omega$ be an arbitrary compact subset of $\left\{(q,m,b):m^{2}\leq q\right\}$. We define

For any constant $\delta>0$ and as $n,d\to\infty$ with $\alpha=n/d$ fixed, it holds that

where ${\cal L}_{\lambda}^{\ast}$ is the global training loss defined in (22).

The characterization in (27) shows that the global training loss will concentrate around the fixed value ${\cal E}_{\lambda}^{\ast}$. Meanwhile, (26) implies that the deterministic function ${\cal E}_{\lambda}(q,m,b)$ serves as a high-probability lower bound of the local training loss ${\cal L}_{\lambda}(\Omega)$ over any given compact subset $\Omega$. This latter property allows us to study the high-dimensional landscapes of the training loss as we move along the 3-dimensional space of the parameters $q,m$ and $b$.

In particular, by studying ${\cal E}_{\lambda}(q,m,b)$, we can obtain the phase transition boundary characterizing the critical value of $\alpha$ below which the training data become perfectly separable.

and $f(x)$ is the probability density function of ${\cal N}(0,1)$.

III Proof Sketches

In this section, we sketch the proof steps behind our main results presented in Section II. The full technical details are given in the Appendix C.

Roughly speaking, our proof strategy consists of three main ingredients: (1) Using Gordon’s minimax inequalities , we can show that the random optimization problem associated with the local training loss in (21) can be compared against a much simpler optimization problem (see (32) in Section III.1) that is essentially decoupled over its coordinates; (2) we show in Section III.2 that the aforementioned simpler problem concentrates around a well-defined deterministic limit as $n,d\to\infty$; and (3) by studying properties of the deterministic function, we reach the various characterizations given in Theorem 11, Proposition 1 and Proposition 2.

For example, for the square, logistic, and hinge losses defined in Section I, their corresponding convex conjugates are given by

where $H(u)\overset{\text{def}}{=}-u\log u-(1-u)\log(1-u)$ is the binary entropy function, and

Substituting (29) into (21) and recalling the data model (1), we can rewrite (21) as the following minimax problem

where ${\cal S}_{q,m}\overset{\text{def}}{=}\left\{{\bf w}:q=\frac{1}{d}\lVert{\bf w}\rVert^{2}\text{ and }m=\frac{1}{d}{\bf w}^{\top}{\bf v}^{\ast}\right\}$.

For every $(q,m,b)$ satisfying $q>m^{2}$, let

where $\Delta_{d}\overset{\text{def}}{=}(Q_{d}/d)\Delta$ with $Q_{d}\sim\chi^{2}_{d}$,

and $\bm{s}\sim\mathcal{N}(0,\bm{I}_{n})$ is an i.i.d. Gaussian random vector. Then for any constant $c$ and $\delta>0$, we have

The proof of Proposition 3, which can be found in the Appendix C.1, is based on an application of Gordon’s comparison inequalities for Gaussian processes . Similar techniques have been used by the authors of to study the Gaussian mixture model for the non-regularized logistic loss for two clusters of the same size.

III.2 Asymptotic Characterizations

The definition of ${\cal E}^{(d)}_{\lambda}(q,m,b)$ in (32) still involves an optimization with an $n$-dimensional vector $\bm{u}$, but it can be simplified to a one-dimensional optimization problem with respect to a Lagrange multiplier $\gamma$:

One can show that the problem in (36) reaches its maximum at a point $\gamma^{\ast}$ that is the unique solution to

In the asymptotic limit, as $n,d\to\infty$, both (38) and (39) converge towards their deterministic limits:

Substituting these identities into (40) and (41) then gives us the characterizations (25) and (23) as stated in Section II.

Finally, the fixed point characterizations given in Theorem 11 can be obtained by taking derivatives of ${\cal E}_{\lambda}(q,m,b)$ with respect to $q,m,b$ and setting them to . Similarly, the phase transition curve given in Proposition 2 can be obtained by quantifying the conditions under which the deterministic function ${\cal E}_{\lambda}(q,m,b)$ reaches its minimum at a finite point. We give more details in Appendix C.4 - C.5.

III.3 Interpretation from the replica method

These same equations can be independently derived from the non-rigorous replica methods from statistical physics , a technique that has proven useful in the study of high-dimensional statistical models, for instance following . Alternatively, these equations can also be seen as a special case of the State Evolution equation of the Approximate Message Passing algorithm . Both interpretations can be useful, since the various quantities enjoy additional heuristic interpretations that allow us to obtain further insight. For instance, the parameter $\gamma$ in (6) is connected to the rescaled variance of the estimator ${\bf w}$:

The zero temperature limit of the fixed point equations obtained with the replica method corresponds to the loss minimization . In this limit, the behaviour of the rescaled variance $V$ at zero penalty ($\lambda=0$) is an indicator of data separability. In the non-separable regime, the minimizer of the loss is unique and $V\rightarrow 0$ at temperature $T=0$. The parameter $\gamma$ turns out to be simply $\gamma=\tfrac{V}{T}$. However, in the regime where data are separable there is a degeneracy of solutions at $\lambda=0$, and the variance is finite: $V>0$. Hence the parameter $\gamma$ has a divergence at the transition, and this provides a very easy way to compute the location of the phase transition.

IV Consequences of the formulas

In this section we evaluate the above formulas and investigate how does the test error depend on the regularization parameter $\lambda$, the fraction taken by the smaller cluster $\rho$, the ratio between the number of samples and the dimension $\alpha$ and the cluster variance $\Delta$. The details on the evaluation and iteration of the fixed point equations in Theorem 11 are provided in Appendices D and F respectively. Keeping in mind that minimization of the non-regularized logistic loss corresponds in the considered model to the maximum likelihood estimation (MLE), we thus pay a particular attention to it as a benchmark of what the most commonly used method in statistics would achieve in this problem. Another important benchmark is the Bayes-optimal performance that provides a threshold that no algorithm can improve.

Fig. 1 summarizes how the regularization parameter $\lambda$ and the cluster size $\rho$ influence the generalization performances. The left panel of Fig. 1 is for the symmetric case $\rho\!=\!0.5$, the right panel for the non-symmetric case $\rho=0.2$. Let us define as $\alpha^{*}$ the value of $\alpha$ such that for $\alpha<\alpha^{*}$ the training loss for hinge and logistic goes to zero (in other words, the data are linearly separable . In the left part of Fig. 1 we depict (in green) the performance of the non-regularized logistic loss a.k.a. the maximum likelihood. For $\alpha>\alpha^{*}(\rho,\Delta)$ the training data are not linearly separable and the minimum training loss is bounded away from zero. For $\alpha<\alpha^{*}(\rho,\Delta)$ the data are linearly separable, in which case properly speaking the maximum likelihood is ill-defined , the curve that we depict is the limiting value reached as $\lambda\to 0^{+}$. The points are results of simulations with a standard scikitlearn package. As shown in , even though the logistic estimator does not exist, gradient descent actually converges to the max-margin solution in this case, or equivalently to the least norm solution corresponding to $\lambda\!\to\!0^{+}$, a phenomenon coined “implicit regularization”, which is well illustrated here.

Another interesting phenomenon is the non-monotonicity of the curve. This is actually an avatar of the so-called “double descent” phenomenon where the generalization “peaks” to a bad value and then decays again. This was observed and discussed recently in several papers , but similar observations appeared as early as 1996 in Opper and Kinzel . Indeed, we observed that the generalization error of the non-regularized square loss (in red) has a peak at $\alpha=1$ at which point the data matrix in the non-regularized square loss problem becomes invertible. It is interesting that for $\alpha>\alpha^{*}$ the generalization performance of the non-regularized square loss is better than the one of the maximum likelihood. This has been proven recently in , who showed that among all the convex non-regularized losses, the square loss is optimal.

Fig. 1 further depicts (in purple) the Bayes-optimal error eq. (19). We have also evaluated the performance of both the logistic and square loss at optimal value of the regularization parameter $\lambda$. This is where the symmetric case (left panel) differs crucially from the non-symmetric one (right panel). While in the high-dimensional limit of the symmetric case the optimal regularization $\lambda_{\rm opt}\!\to\!\infty$ and the corresponding error matches exactly the Bayes-optimal error, for the non-symmetric case $0<\lambda_{\rm opt}<\infty$ and the error for both losses is bounded away from the Bayes-optimal one for any $\alpha\!>\!0$.

We give a fully analytic argument in the Appendix E for the perhaps unexpected property of achieving the Bayes-optimal generalization at $\lambda_{\rm opt}\to\infty$ and $\rho=0.5$ for any loss that has a finite 2nd derivative at the origin. In simulations for finite value of $d$ we use a large but finite value of $\lambda$, details on the simulation are provided in the Appendix F.

Regularization and the interpolation peak —

In Fig. 2 we depict the dependence of the generalization error on the regularization $\lambda$ for the symmetric $\rho=0.5$ case for the square, hinge and logistic loss. The curves at small regularization show the interpolation peak/cusp at $\alpha=1$ for the square loss and $\alpha^{*}$ for all the losses that are zero whenever the data are linearly separable. We observe a smooth disappearance of the peak/cusp as regularization is added, similarly to what has been observed in other models that present the interpolation peak in the case of the square loss. Here we thus show that a similar phenomena arises with the logistic and hinge losses as well; this is of interest as this effect has been observed in deep neural networks using a logistic/cross-entropy loss . In fact, as the regularization increases, the error gets better in this model with equal-size cluster, and one reaches the Bayes-optimal values for large regularization.

Max-margin and weak regularization —

Fig. 3 illustrates the generic property that all non-regularized monotone non-increasing loss functions converge to the max-margin solution for linearly separable data . Fig. 3 depicts a very slow convergence towards this result as a function of regularization parameter $\lambda$ for the logistic loss. While for $\alpha>\alpha^{*}$ both the hinge and logistic losses performance is basically indistinguishable from the asymptotic one already at $\log\lambda\approx-3$, for $\alpha<\alpha^{*}$ the convergence of the logistic loss still did not happen even at $\log\lambda\approx-10$.

Cluster sizes and regularization —

In Fig. 4 we study in greater detail the dependence of the generalization error both on the regularization $\lambda$ and $\rho$ as $\rho\to 0.5$. We see that the optimality of $\lambda\to\infty$ holds only strictly at $\rho=0.5$ and at any $\rho$ only close to $0.5$ the error at $\lambda\to\infty$ is very large and there is a well delimited region of $\lambda$ for which the error is close to (but strictly above) the Bayes-optimal error. As $\rho\to 0.5$ this interval is getting longer and longer until it diverges at $\rho=0.5$. It needs to be stressed that this result is asymptotic, holding only when $n,d\to\infty$ while $n/d=\alpha$ is fixed. The finite size fluctuations cause that finite size system behaves rather as if $\rho$ was close but not equal to $0.5$, and at finite size if we set $\lambda$ arbitrarily large then we reach a high generalization error. We instead need to optimize the value of $\lambda$ for finite sizes either by cross-validation or otherwise.

Separability phase transition —

The position of the “interpolation” threshold when data become linearly separable has a well defined limit in the high-dimensional regime as a function of the ratio between the number of samples $n$ and the dimension $d$. The kink in generalization indeed occurs at a value $\alpha^{*}$ when the training loss of logistic and hinge losses goes to zero (while for the square loss the peak appears at $d=n$ when the system of $n$ linear equations with $d$ parameters becomes solvable). The position of $\alpha^{*}$, given by Proposition 2 is shown in Fig. 5 as a function of the cluster variance for different values of $\rho$. For very large cluster variance, the data become random and hence $\alpha=2$ for equal-sized cluster, as famously derived in classical work by . When $\rho<1/2$, however, it is easier to separate linearly the data points and the limiting value of $\alpha^{*}$ gets larger and differ from Cover’s. For finite $\Delta$, the two Gaussian distributions become distinguishable, and the data acquires structure. Consequently, the $\alpha^{*}$ is growing as the correlations make data easier to linearly separate again, similarly as described . This phenomenology of the separability phase transition, or equivalently of the existence of the maximum likelihood estimator, thus seems very generic.

Acknowledgements

We thank Pierfrancesco Urbani, Federica Gerace, and Bruno Loureiro for many clarifying discussions related to this project. This work is supported by the ERC under the European Union’s Horizon 2020 Research and Innovation Program 714608-SMiLe, by the French Agence Nationale de la Recherche under grant ANR-17-CE23-0023-01 PAIL and ANR-19-P3IA-0001 PRAIRIE, and by the US National Science Foundation under grants CCF-1718698 and CCF-1910410. We also acknowledge support from the chaire CFM-ENS “Science des données”. Part of this work was done when Yue Lu was visiting Ecole Normale as a CFM-ENS “Laplace” invited researcher. We thank Google Cloud for providing us access to their platform through the Research Credits Application program.

References

Appendix A Derivation of the generalization error formula

The generalization error is defined as the average fraction of mislabeled instances

where $y_{\text{new}}$ is the label of a new observation ${\bf x}_{\text{new}}$, and the estimator $\hat{y}_{\text{new}}$ is computed as

Eq. (A.2) holds for every vector ${\bf w}={\bf w}\left({\bf X},{\bf y}\right)$ and bias $b=b\left({\bf X},{\bf y}\right)$ computed on the training set $\left\{{\bf X},{\bf y}\right\}$. Using the fact that $y_{\text{new}},\hat{y}_{\text{new}}=\pm 1$, it is easy to show that (A.1) can be rewritten as

Let us consider the last term in (A.3). Using again $y_{\text{new}}=\pm 1$, we can move $y_{\text{new}}$ inside the argument of the sign function and rewrite

The term $y_{\text{new}}{\bf x}_{\text{new}}$ can be rewritten as

where ${\bf z}^{\prime}_{\text{new}}=y_{\text{new}}{\bf z}_{\text{new}}\sim\mathcal{N}({\bf 0},{\bf I}_{d})$ has the same distribution as ${\bf z}_{\text{new}}$, since $y_{\text{new}}$ and ${\bf z}_{\text{new}}$ are independent. Hence

The estimator ${\bf w}$ only depends on the training set, hence ${\bf w}$ and ${\bf z}^{\prime}_{\text{new}}$ are independent. We call their rescaled scalar product $\varsigma$, a random variable distributed as a standard normal

where we have used that $\sqrt{\frac{\Delta}{d}}\lVert{\bf w}\rVert>0$ to rescale the argument of the sign function. Finally, we obtain

where $Q(x)=\frac{1}{\sqrt{2\pi}}\int_{x}^{\infty}e^{-t^{2}/2}dt$ is the Gaussian tail function, and we have defined

In the large $d$ limit, the overlaps concentrate to deterministic quantities:

where $\rho\in(0,1)$ is the probability that $y_{\text{new}}=+1$.

Appendix B Derivation of the Bayes-optimal error

In order to compute the Bayes-optimal error, we consider the posterior distribution of a new label $y_{\text{new}}$, given the corresponding new data point $\bf x_{\text{new}}$ and the estimate $\bf v$ of the true centroid $\bf v^{*}$

where “$\propto$” takes into account the normalization over $y_{\text{new}}$. Similarly, the posterior on ${\bf v}$ given the training data is

where we remind that ${\bf v}$ has i.i.d. components taken in $\mathcal{N}(0,1)$, and “$\propto$” takes into account the normalization over ${\bf v}$. We would like to find an explicit expression for

where in the product over $\mu$ on the right-hand side we have used the notation $y_{0}=y_{\text{new}}$, ${\bf x}_{0}={\bf x}_{\text{new}}$. Let us call $I_{\text{v}}$ the integral over $\bf v$ in (B.5).

where in the last equality we have dropped the index $i$ from the components of $\bf v$ for simplicity, since they are all independent. Computing the integral over v, we obtain

Using the fact that $y_{\mu}{\bf x}_{\mu}=\frac{{\bf v}^{*}}{\sqrt{d}}+\sqrt{\Delta}{\bf z}_{\mu}$, ${\bf z}_{\mu}\sim\mathcal{N}(0,{\bf I}_{d})$ and ${\bf v}^{*}$ is the true realization of $\bf v$, the first term in (B.8) in the limit where $n,d\rightarrow\infty$ can be rewritten as

where $z^{\prime}_{\text{new}}\sim\mathcal{N}(0,1)$. Therefore, in the large $d$ limit we find that

It is useful to rewrite the generalization error as

where $Q(x)=\frac{1}{\sqrt{2\pi}}\int_{x}^{\infty}e^{-t^{2}/2}dt$ is the Gaussian tail function. If $y_{\text{new}}=-1$, (B.12) gives

Using the fact that $\rho={\rm p}_{y}(1)$ and $1-\rho={\rm p}_{y}(-1)$, we get that

It is worth noting that the optimal error in (B.15) can be achieved by the plug-in estimator

This result was already shown in for the case of symmetric clusters. The optimal bias is obtained from the minimization of the generalization error (A.13) with respect to $b$, at fixed $m,q$. This yields:

Substituting (B.16) in the definition of the overlaps (3) in the main text, we obtain that the values of $m$ and $q$ associated to the plugin estimator are

Hence, the generalization error of the plug-in estimator is

where we have used (B.9) in the last equality. The probability in (B.19) is the same as in (B.12). Hence, the plug-in estimator achieves the Bayes-optimal error.

Appendix C Details of proofs

In what follows, we provide more technical details for several key results stated in Section III. They serve as the basis of the proof of Proposition 1.

where in reaching the second equality we have used the fact that any ${\bf w}\in{\cal S}_{q,m}$ satisfies the equality $m=\frac{1}{d}{\bf w}^{\top}{\bf v}^{\ast}$. Introduce an auxiliary problem

where ${\bf g}=(g_{1},g_{2},\ldots,g_{d})^{\top}$ and ${\bf s}=(s_{1},s_{2},\ldots,s_{n})$ are two independent random vectors whose entries are drawn from the i.i.d. standard normal distribution, and $h_{i}=\sqrt{\Delta q}(y_{i}s_{i})+m+by_{i}$. As $y_{i}\in\left\{\pm 1\right\}$, independent of $s_{i}$, we note that $h_{i}$ has the same probability distribution as the quantity defined in (33) in the main text.

Gordon’s minimax inequalities allow us to make the following comparison: For any constants $c$ and $\delta>0$, we have

To connect this to the statements in Proposition 3, we note that

Combining this inequality with (C.1) gives us the first inequality in Proposition 3. To obtain the second inequality in the proposition, we use the fact that the unconstrained optimization problem in (22) for the global training loss ${\cal L}^{\ast}$ is convex. Following exactly the same strategy as used in , we can interchange the order of $\min$ and $\max$ in the dual formulation of (22), which then allows us to reach the two-sided inequality in (35).

C.2 Proof of Lemma 1

We first rewrite the optimization problem in (32) as

For the inner maximization, the constraint on the squared norm $\lVert{\bf u}\rVert^{2}$ weakly couples different coordinates of $\bm{u}$ together. To fully decouple these coordinates, we introduce a Lagrangian function

Since there is a one-to-one correspondence between the Lagrange multiplier $\gamma$ and the normalized squared norm $\mu=\lVert\bm{u}_{\gamma}\rVert^{2}/d$, it is thus equivalent to solve (C.2) in terms of

C.3 Proof of Proposition 1

We first establish (26) for the special case where the subset $\Omega$ is a singleton. In this case, we just need to show

for any fixed $q,m$ and $b$. Recall the characterization of ${\cal E}^{(d)}_{\lambda}(q,m,b)$ given in Lemma 1. The problem in (36) reaches its maximum at a point $\gamma^{\ast}_{d}$ where the derivative of the function to be maximized is equal to 0. In calculating this derivative, we need the quantity $\frac{du_{\gamma,i}}{d\gamma}$, which can be obtained as

Substituting these identities, we can characterize $v_{\gamma,i}$ via the implicit equation

and more importantly, (C.5) can be simplified as

Let $v_{\gamma}$ be a random variable defined via the implicit equation

where $h=\sqrt{\Delta q}s+m+by$ with $S\sim\mathcal{N}(0,1)$ and $y$ being a random variable independent of $s$ such that

uniformly over any compact subset of $\gamma$. It follows that $\gamma^{\ast}_{d}$ as defined in (C.7) converges to $\gamma^{\ast}$, which is the unique solution of (25). Moreover, we have

For any $\delta>0$, we can apply Proposition 3 to get

As the right-hand side tends to due to (C.9), we have (C.3).

Let $\Omega$ be an arbitrary compact subset of $\left\{(q,m,b):m^{2}\leq q\right\}$. We denote by $\Omega_{K}$ a finite subset of $\Omega$ consisting of $K$ points, i.e., $\Omega_{K}=\left\{(q_{k},m_{k},b_{k})\in\Omega:1\leq k\leq K\right\}$.

C.4 Proof of Proposition 2

It follows that the cumulant distribution function of $u_{\gamma}$ is given by

where $Q(\cdot)$ is the distribution function of a standard normal random variable. Writing (25) in terms of $u_{\gamma}$, we have

We further denote by $\widehat{\gamma}^{\ast}(q,\theta)$ the solution to (C.11). We can show that, for any fixed $\widetilde{\gamma}$ and $\theta$, the function $S(\widetilde{\gamma},q,\theta)$ is monotonically decreasing as we increase $q$. Moreover,

Clearly, $S^{\ast}(\widetilde{\gamma},\theta)$ is monotonic with respect to $\widetilde{\gamma}$, but it has a finite limit as $\widetilde{\gamma}\to\infty$, i.e.,

where $f(\cdot)$ is the probability density function of $\mathcal{N}(0,1)$. An implication of this limit being finite is that, although the Lagrange multiplier $\widehat{\gamma}^{\ast}(q,\theta)$ remains finite for any fixed $q$, it tends to $\infty$ as $q\to\infty$ if

This characterization can be interpreted as follows: If there exists a $\theta$ that satisfies (C.12), then as we move along the “ray” of constant slope $\theta=m/\sqrt{q}$, the training loss ${\cal E}_{\lambda=0}(q,m,b)$ will tend to . The critical threshold $\alpha^{\ast}$ can then be obtained by maximizing the right-hand side of (C.12), which gives us the final expression as stated in Proposition 2.

C.5 Derivation of Theorem 11 from Gordon’s characterization

In this section, we show that the fixed point equations in Theorem 11 can be mapped to Gordon’s characterization, namely (25) and (27) in the main text. First of all, we observe that (25) is trivially satisfied by the solution of system (4)-(9). Then, we consider the minimization of ${\cal E}_{\lambda}(q,m,b)$, derived in (C.9), with respect to $q,m,b$. This simply amounts to setting the derivatives to zero. Note that the partial derivatives of $v$ and $\gamma^{*}$ can be computed by taking the derivatives of both sides of (C.6) and (25) respectively. The minimization leads to the following system of equations:

where $s\sim\mathcal{N}(0,1)$, $y=+1$ with probability $\rho\in(0,1)$ and $y=-1$ otherwise. We observe that (C.15) is the same as (11) and (C.14) is equivalent to (4) and (6). Using again (6), we can rewrite (C.13) as

which leads to an identity if we substitute the definition of $\hat{\gamma}$ provided in (9).

Appendix D Evaluation of the fixed point equations

In this section we will compute the fixed-point equations for the square and hinge loss. The equations for the logistic loss cannot be computed analytically and require numerical integration.

where $h\sim\mathcal{N}(m+yb,\Delta q)$. Hence, we obtain

To compute the bias $b$, we have to solve

We can plug (D.3)-(D.5) in the equations for $m,q,\gamma$ to obtain

D.2 Hinge loss

Appendix E Bayes-optimality at λ=∞𝜆\lambda=\infty, for ρ=12𝜌12\rho=\tfrac{1}{2}

In this section we will show how the result on Bayes-optimality for balanced clusters at large regularization arises. First we start by considering the square loss. At $\rho=1/2$, it is straightforward to check from (11) that $b=0$ and the generalization error, given by (12) in the main text, is

where $m$ and $q$ are given by (D.9)-(D.10), evaluated at $\rho=\tfrac{1}{2}$. The Bayes-optimal error for this problem is given by (19) in the main text and reads

Therefore, in order to reach Bayes-optimality, we need a weight vector $\bf w$ with an overlap $m$ and a length $q$ such that

By using (D.3)-(D.4) evaluated at $\rho=\tfrac{1}{2}$, (E.3) can be rewritten as

Eq. (E.4) is verified by the fixed point equations only at $\lambda\to\infty$. Indeed in this limit we find that

so that any loss will behave like the square one. This is the origin of the peculiar behavior of Bayes optimally observed at $\lambda\to\infty$ for the symmetric case $\rho=1/2$. We observed numerically that this result is not valid anymore as soon as $\rho\neq 1/2$. This peculiar behaviour is shown in Fig. 6, which depicts the generalization error, computed from the solution of (4)-(11) in the main text, as a function of $\rho$ at zero, infinite and optimal regularization for the square and hinge losses.

Appendix F Details on the numerics

The solution $(q,m,b,\gamma)$ of the fixed point equations (4)-(9) can be obtained analytically only in the case of square loss. For the hinge and logistic loss, the equations have to be iterated until convergence. In our codes we used initialization $(q^{t=0},\gamma^{t=0},m^{t=0},b^{t=0})=(0.5,0.5,0.01,0)$. The stopping criterion for convergence consists in checking if the values of the generalization error at two consecutive iterations differ less than a threshold $eps$. In all figures, we used $eps\leq 10^{-8}$.

F.2 Simulations

In order to check the validity of the fixed point equations (4)-(9) we computed numerically the solution of the optimization problem defined in (2), and we averaged over multiple realizations of the noise. In the case of square loss, the solution is simply

In the case of logistic and hinge loss, the solution can be computed by a standard gradient descent algorithm. In particular, in Fig. 1 we used the Logistic Regression classifier provided by the scikitlearn package $linear\_model$ . In particular, we used the “lbfgs” solver, with L2-penalty, tolerance $tol=10^{-5}$ for the stopping criterion and maximum number of iterations $max\_iter=10^{-5}$. It is important to remind that all our analytic results are computed in the infinite-dimensional limit $d,n\rightarrow\infty$, while the ratio $\alpha=n/d$ remains finite. Therefore, all the simulations involve errors due to finite size effects. However, we found a very good agreement bewteen theory and simulations already at relatively small dimensionality ($d\leq 5000$). The only case in which finite size effects prevent simulations to match our theoretical predictions is the behaviour of the generalization error at large regularization $\lambda$, at $\rho=1/2$. Since at all finite dimensions $d$ the effective clusters size is $\rho\neq 1/2$, the result of reaching Bayes-optimality at $\lambda\rightarrow\infty$ cannot be obtained in simulations, since it holds strictly at $\rho=1/2$. However, we obtain greater and greater precision, i.e. the minimum of the generalization error moving towards higher values of $\lambda$ (see Fig. 4), as $d$ increases.