Learning curves of generic features maps for realistic datasets with a teacher-student model

Introduction

Teacher-student models are a popular framework to study the high-dimensional asymptotic performance of learning problems with synthetic data, and have been the subject of intense investigations spanning three decades . In the wake of understanding the limitations of classical statistical learning approaches , this direction is witnessing a renewal of interest . However, this framework is often assuming the input data to be Gaussian i.i.d., which is arguably too simplistic to be able to capture properties of realistic data. In this paper, we redeem this line of work by defining a Gaussian covariate model where the teacher and student act on different Gaussian correlated spaces with arbitrary covariance. We derive a rigorous asymptotic solution of this model generalizing the formulas found in the above mentioned classical works.

We then put forward a theory, supported by universality arguments and numerical experiments, that this model captures learning curves, i.e. the dependence of the training and test errors on the number of samples, for a generic class of feature maps applied to realistic datasets. These maps can be deterministic, random, or even learnt from the data. This analysis thus gives a unified framework to describe the learning curves of, for example, kernel regression and classification, the analysis of feature maps – random projections , neural tangent kernels , scattering transforms – as well as the analysis of transfer learning performance on data generated by generative adversarial networks . We also discuss limits of applicability of our results, by showing concrete situations where the learning curves of the Gaussian covariate model differ from the actual ones.

The labels $y^{\mu}$ are generated by a teacher function that is only using the vectors $\bm{u}^{\mu}$:

Our two main technical contributions are:

In Theorems 1 & 2, we give a rigorous closed-form characterisation of the properties of the estimator $\hat{\bm{w}}$ for the Gaussian covariate model (1.1), and the corresponding training and generalisation errors in the high-dimensional limit. We prove our result using Gaussian comparison inequalities ;

We show how the same expression can be obtained using the replica method from statistical physics . This is of additional interest given the wide range of applications of the replica approach in machine learning and computer science . In particular, this allows to put on a rigorous basis many results previously derived with the replica method.

The Gaussian covariate model (1.1) is exact in the case where $\bm{x}$ are Gaussian variables and the feature maps $(\bm{\varphi}_{s},\bm{\varphi}_{s})$ preserve the Gaussianity, for example linear features. In particular, this is the case for $\bm{u}\!=\!\bm{v}\!=\!\bm{x}$, which is the widely-studied vanilla teacher-student model . The interest of the model (1.1) is that it also captures a range of cases in which the feature maps $\bm{\varphi}_{t}$ and $\bm{\varphi}_{s}$ are deterministic, or even learnt from the data. The covariance matrices $\Psi$, $\Phi$, and $\Omega$ then represent different aspects of the data-generative process and learning model. The student (1.3) then corresponds to the last layer of the learning model. These observation can be distilled into the following conjecture:

(Gaussian equivalent model) For a wide class of data distributions $\{\bm{x}^{\mu}\}_{\mu=1}^{n}$, and features maps $\bm{u}=\bm{\varphi}_{t}(\bm{x}),\bm{v}=\bm{\varphi}_{s}(\bm{x})$, the generalisation and training errors of estimator (1.3) are asymptotically captured by the equivalent Gaussian model (1.1), where $[\bm{u},\bm{v}]$ are jointly Gaussian variables, and thus by the closed-form expressions of Theorem 1.

The second part of our main contributions are:

In Sec. 3.3 we show that the theoretical predictions from (C1) captures the learning curves in non-trivial cases, e.g. when input data are generated using a trained generative adversarial network, while extracting both the feature maps from a neural network trained on real data.

In Sec. 3.4, we show empirically that for ridge regression the asymptotic formula of Theorem 1 can be applied directly to real data sets, even though the Gaussian hypothesis is not satisfied. This universality-like property is a consequence of Theorem 3 and is illustrated in Fig. 1 (right) where the real learning curve of several features maps learning the odd-versus-even digit task on MNIST is compared to the theoretical prediction.

Rigorous results for teacher-student models: The Gaussian covariate model (1.1) contains the vanilla teacher-student model as a special case where one takes $\bm{u}$ and $\bm{v}$ identical, with unique covariance matrix $\Omega$. This special case has been extensively studied in the statistical physics community using the heuristic replica method . Many recent rigorous results for such models can be rederived as a special case of our formula, e.g. refs. . Numerous of these results are based on the same proof technique as we employed here: the Gordon’s Gaussian min-max inequalities . The asymptotic analysis of kernel ridge regression , of margin-based classification also follow from our theorem. See also Appendix A.6 for the details on these connections. Other examples include models of the double descent phenomenon . Closer to our work is the recent work of on the random feature model. For ridge regression, there are also precise predictions thanks to random matrix theory . A related set of results was obtained in for orthogonal random matrix models. The main technical novelty of our proof is the handling of a generic loss and regularisation, not only ridge, representing convex empirical risk minimization, for both classification and regression, with the generic correlation structure of the model (1.1).

Gaussian equivalence: A similar Gaussian conjecture has been discussed in a series of recent works, and some authors proved partial results in this direction . Ref. analyses a special case of the Gaussian model (corresponding to $\bm{\varphi}_{t}=\rm{id}$ here), and proves a Gaussian equivalence theorem (GET) for feature maps $\bm{\varphi}_{s}$ given by single-layer neural networks with fixed weights. They also show that for Gaussian data $\bm{x}\sim\mathcal{N}(\bm{0},\text{I}_{D})$, feature maps of the form $\bm{v}=\sigma(\text{W}\bm{x})$ (with some technical restriction on the weights) led to the jointly-Gaussian property for the two scalars $(\bm{v}\cdot\bm{w},\bm{u}\cdot\bm{\theta}_{0})$ for almost any vector $\bm{w}$. However, their stringent assumptions on random teacher weights limited the scope of applications to unrealistic label models. A related line of work discussed similar universality through the lens of random matrix theory . In particular, Seddik et al. showed that, in our notations, vectors $[\bm{u},\bm{v}]$ obtained from Gaussian inputs $\bm{x}\sim\mathcal{N}(\bm{0},\text{I}_{D})$ with Lipschitz feature maps satisfy a concentration property. In this case, again, one can expect the two scalars $\left(\bm{v}\cdot\bm{w},\bm{u}\cdot\bm{\theta}_{0}\right)$ to be jointly Gaussian with high-probability on $\bm{w}$. Remarkably, in the case of random feature maps, could go beyond this central-limit-like behavior and established the universality of the Gaussian covariate model (1.1) for the actual learned weights $\hat{\bm{w}}$.

Main technical results

We start by defining key quantities that we will use to characterize the estimator $\hat{\bm{w}}$. Let $\Omega=\text{S}^{\top}\text{diag}(\omega_{i})\text{S}$ be the spectral decomposition of $\Omega$. Let:

and define the joint empirical density $\hat{\mu}_{d}$ between $(\omega_{i},\bar{\theta}_{i})$:

Note that $\Phi^{\top}\bm{\theta}_{0}$ is the projection of the teacher weights on the student space, and therefore $\bar{\bm{\theta}}$ is the rotated projection on the basis of the student covariance, rescaled by the teacher variance. Together with the student eigenvalues $\omega_{i}$, these are relevant statistics of the model, encoded here in the joint distribution $\hat{\mu}_{d}$.

Consider the high-dimensional limit in which the number of samples $n$ and the dimensions $p,d$ go to infinity with fixed ratios:

The relevance of these assumptions in a supervised machine learning context is discussed in Appendix B.1. We are now in a position to state our result.

where prox stands for the proximal operator defined as

and where $(\nu,\lambda)$ are jointly Gaussian scalar variables:

and the overlap parameters $(V^{\star},q^{\star},m^{\star})$ are prescribed by the unique fixed point of the following set of self-consistent equations:

where we defined the scalar random functions $f_{g}(V,m,q)=\mbox{prox}_{Vg(.,f_{0}(\sqrt{\rho}s))}(\rho^{-1/2}ms+\sqrt{q-\rho^{-1}m^{2}}h)$ and $f^{\prime}_{g}(V,m,h)=\mbox{prox}^{\prime}_{Vg(.,f_{0}(\sqrt{\rho}s))}(\rho^{-1/2}ms+\sqrt{q-\rho^{-1}m^{2}}h)$ as the first derivative of the proximal operator.

Proof: This result is a consequence of Theorem 2, whose proof can be found in appendix B.

Or in words: $m^{\star}$ is the correlation between the estimator projected in the teacher space, while $q^{\star}$ is the reweighted norm of the estimator by the covariance $\Omega$. The parameter $V^{*}$ also has a concrete interpretation : it parametrizes the deformation that must be applied to a Gaussian field specified by the solution of the fixed point equations to obtain the asymptotic behaviour of $\hat{\mathbf{z}}$. It prescribes the degree of non-linearity given to the linear output by the chosen loss function. This is coherent with the robust regression viewpoint, where one introduces non-square losses to deal with the potential non-linearity of the generative model. $\hat{V}^{*}$ plays a similar role for the estimator $\hat{\mathbf{w}}$ through the proximal operator of the regularisation, see Theorem 4 and 5 in the Appendix. Two cases are of particular relevance for the experiments that follow. The first is the case of ridge regression, in which $f_{0}(x)=\hat{f}(x)$ and both the loss $g$ and the performance measure $\hat{g}$ are taken to be the mean-squared error $\text{mse}(y,\hat{y})=\frac{1}{2}(y-\hat{y})^{2}$, and the asymptotic errors are given by the simple closed-form expression:

while the training error $\mathcal{E}^{\star}_{\rm{train}}$ depends on the choice of $g$ - which we will take to be the logistic loss $g(y,x)=\log\left(1+e^{-xy}\right)$ in all of the binary classification experiments.

where $\bm{s},\bm{h}\sim\mathcal{N}(0,\text{I}_{n})$ and $\bm{g}\sim\mathcal{N}(0,\text{I}_{d})$ are random vectors independent of the other quantities, $\bm{t}=\Omega^{-1/2}\Phi^{\top}\bm{\theta}_{0}$, $\bm{y}=\bm{f}_{0}\left(\sqrt{\rho}\bm{s}\right)$, and $(V^{*},\hat{V}^{*},q^{*},\hat{q}^{*},m^{*},\hat{m}^{*})$ is the unique solution to the fixed point equations presented in Lemma 12 of appendix B. Those fixed point equations are the generalization of (2.8) to generic, non-separable loss function and regularization. The formal concentration of measure result can then be stated in the following way:

Applications of the Gaussian model

We now discuss how the theorems above are applied to characterise the learning curves for a range of concrete cases. We present a number of cases – some rather surprising – for which Conjecture 1 seems valid, and point out some where it is not. An out-of-the-box iterator for all the cases studied hereafter is provided in the GitHub repository for this manuscript at https://github.com/IdePHICS/GCMProject.

If we choose random feature maps $\bm{\varphi}_{s}(\bm{x})=\sigma\left(\text{F}\bm{x}\right)$ for a random matrix F and a chosen scalar function $\sigma$ acting component-wise, we obtain the random kitchen sink model . This model has seen a surge of interest recently, and a sharp asymptotic analysis was provided in the particular case of uncorrelated Gaussian data $\bm{x}\sim\mathcal{N}(\bm{0},\text{I}_{D})$ and $\bm{\varphi}_{t}(\bm{x})=\bm{x}$ in for ridge regression and generalised by for generic convex losses. Both results can be framed as a Gaussian covariate model with:

In this case, the averages over $\mu$ in eq. (2.8) can be directly expressed in terms of the Stieltjes transform associated with the spectral density of $\text{F}\text{F}^{\top}$. Note, however, that our present framework can accommodate more involved random sinks models, such as when the teacher features are also a random feature model or multi-layer random architectures.

2 Kernel methods with Gaussian data

Another direct application of our formalism is to kernel methods. Kernel methods admit a dual representation in terms of optimization over feature space . The connection is given by Mercer’s theorem, which provides an eigen-decomposition of the kernel and of the target function in the feature basis, effectively mapping kernel regression to a teacher-student problem on feature space. The classical way of studying the performance of kernel methods is then to directly analyse the performance of convex learning in this space. In our notation, the teacher and student feature maps are equal, and we thus set $p=d,\Psi=\Phi=\Omega=\rm{diag}(\omega_{i})$ where $\omega_{i}$ are the eigenvalues of the kernel and we take the teacher weights $\bm{\theta}_{0}$ to be the decomposition of the target function in the kernel feature basis.

There are many results in classical learning theory on this problem for the case of ridge regression (where the teacher is usually called "the source" and the eigenvalues of the kernel matrix the "capacity", see e.g. ). However, these are worst case approaches, where no assumption is made on the true distribution of the data. In contrast, here we follow a typical case analysis, assuming Gaussianity in feature space. Through Theorem 1, this allows us to go beyond the restriction of the ridge loss. An example for logistic loss is in Fig. 2.

For the particular case of kernel ridge regression, Th. 1 provides a rigorous proof of the formula conjectured in . App. A.6 presents an explicit mapping to their results. Hard-margin Support Vector Machines (SVMs) have also been studied using the heuristic replica method from statistical physics in . In our framework, this corresponds to the hinge loss $g(x,y)=\text{max}(0,1-yx)$ when $\lambda\to 0^{+}$. Our theorem thus puts also these works on rigorous grounds, and extends them to more general losses and regularization.

3 GAN-generated data and learned teachers

To approach more realistic data sets, we now consider the case in which the input data $\bm{x}\in\mathcal{X}$ is given by a generative neural network $\bm{x}=\mathcal{G}(\bm{z})$, where $\bm{z}$ is a Gaussian i.i.d. latent vector. Therefore, the covariates $[\bm{u},\bm{v}]$ are the result of the following Markov chain:

Fig. 3 depicts the resulting learning curves obtained by training the last layer of the student. Interestingly, the performance of the feature map at epoch (random initialisation) beats the performance of the learned features during early phases of training in this experiment. Another interesting behaviour is given by the separability threshold of the learned features, i.e. the number of samples for which the training loss becomes larger than in logistic regression. At epoch $50$ the learned features are separable at lower sample complexity $\alpha=n/d$ than at epoch $200$ - even though in the later the training and generalisation performances are better.

4 Learning from real data sets

Given that the learning curves of realistic-looking inputs can be captured by the Gaussian covariate model, it is fair to ask whether the same might be true for real data sets. To test this idea, we first need to cast the real data set into the teacher-student formalism, and then compute the covariance matrices $\Omega,\Psi,\Phi$ and teacher vector $\bm{\theta}_{0}$ required by model (1.1).

At this point, we have all we need to run the self-consistent equations (2.8). The issue with this approach is that there is not a unique teacher map $\bm{\varphi}_{t}$ and teacher vector $\bm{\theta}_{0}$ that fit the true labels. However, we can show that all interpolating linear teachers are equivalent:

(Universality of linear teachers) For any teacher feature map $\bm{\varphi}_{t}$, and for any $\bm{\theta}_{0}$ that interpolates the data so that $y^{\mu}=\bm{\theta}_{0}^{\top}\bm{u}^{\mu}\,\forall\mu$, the asymptotic predictions of model (1.1) are equivalent.

It follows from the fact that the teacher weights and covariances only appear in eq. (2.8) through $\rho=\frac{1}{p}\bm{\theta}_{0}^{\top}\Psi\bm{\theta}_{0}$ and the projection $\Phi^{\top}\bm{\theta}_{0}$. Using the estimation (3.4) and the assumption that it exists $y^{\mu}=\bm{\theta}_{0}^{\top}\bm{u}^{\mu}$, one can write these quantities directly from the labels $y^{\mu}$:

For linear interpolating teachers, results are thus independent of the choice of the teacher. ∎

Although this result might seen surprising at first sight, it is quite intuitive. Indeed, the information about the teacher model only enters the Gaussian covariate model (1.1) through the statistics of $\bm{u}^{\top}\bm{\theta}_{0}$. For a linear teacher $f_{0}(x)=x$, this is precisely given by the labels.

Fig. 1 (right) shows a similar experiment on the MNIST data set, but for different out-of-the-box feature maps, such as random features and the scattering transform , and we chose the number of random features $d=1953$ to match the number of features from the scattering transform. Note the characteristic double-descent behaviour , and the accurate prediction of the peak where the interpolation transition occurs. We note in Appendix D.1 that for both Figs. 4 and 1, for a number of samples $n$ closer to $n_{\rm tot}$ we start to see deviations between the real learning curve and the theory. This is to be expected since in the teacher-student framework the student can, in principle, express the same function as the teacher if it recovers its weights exactly. Recovering the teacher weights becomes possible with a large training set. In that case, its test error will be zero. However, in our setup the test error on real data remains finite even if more training data is added, leading to the discrepancy between teacher-student learning curve and real data, see Appendix D.1 for further discussion.

Why is the Gaussian model so effective for describing learning with data that are not Gaussian? The point is that ridge regression is sensitive only to second order statistics, and not to the full distribution of the data. It is a classical property (see Appendix E) that the training and generalisation errors are only a function of the spectrum of the empirical and population covariances, and of their products. Random matrix theory teaches us that such quantities are very robust, and their asymptotic behaviour is universal for a broad class of distributions of $[\bm{u},\bm{v}]$ . The asymptotic behavior of kernel matrices has indeed been the subject of intense scrutiny . Indeed, a universality result akin to Theorem 3 was noted in in the specific case of kernel methods. We thus expect the validity of model (1.1) for ridge regression, with a linear teacher, to go way beyond the Gaussian assumption.

We present an additional experiment in Fig. 3. We compare the learning curves of logistic regression on a classification task on the real CIFAR10 images with the real labels versus the one on dcGAN-generated CIFAR10-like images and teacher generated labels from Sec. 3.3. While the Gaussian theory captures well the behaviour of the later, it fails on the former. A histogram of the distribution of the product $\bm{u}^{\top}\hat{\bm{w}}$ for a fixed number of samples illustrates well the deviation from the prediction of the theory with the real case, in particular on the tails of the distribution. The difference between GAN generated data (that fits the Gaussian theory) and real data is clear. Given that for classification problems there exists a number of choices of "sign" teachers and feature maps that give the exact same labels as in the data set, an interesting open question is: is there a teacher that allows to reproduce the learning curves more accurately? This question is left for future works.

Acknowledgements

We thank Romain Couillet, Cosme Louart, Loucas Pillaud-Vivien, Matthieu Wyart, Federica Gerace, Luca Saglietti and Yue Lu for discussions. We are grateful to Kabir Aladin Chandrasekher, Ashwin Pananjady and Christos Thrampoulidis for pointing out discrepancies in the finite size rates and insightful related discussions. We acknowledge funding from the ERC under the European Union’s Horizon 2020 Research and Innovation Programme Grant Agreement 714608-SMiLe, and from the French National Research Agency grants ANR-17-CE23-0023-01 PAIL.

References

Appendix A Main result from the replica method

In this appendix we derive the formula for the performance of the Gaussian covariate model from a heuristic replica analysis. The computation closely follows the recent developments in . We refer to for an introduction to this remarkable heuristic (but seemingly never failing) approach.

In our analysis, we are interested in the training and generalisation performance of a linear classifier $\hat{y}=f_{\bm{w}}(\bm{v})=\hat{f}\left(\bm{w}^{\top}\bm{v}\right)$ trained on $n$ independent samples $\mathcal{D}=\{(\bm{v}^{\mu},y^{\mu})\}_{\mu=1}^{n}$ from $p_{\bm{\theta}_{0}}$ by minimising the regularised empirical risk:

where $\lambda>0$ is the regularisation strength. We define the sample complexity $\alpha=n/d$ and the aspect ratio $\gamma=p/d$.

As it was proven in Theorem 4 of the main manuscript, the asymptotic performance of the estimator in eq. (A.3) is fully characterised by the following scalar parameters:

The replica method is precisely a heuristic tool allowing us to circumvent the high-dimensional estimation problem defined in eq. (A.3) and giving us direct access to $(m^{\star},q^{\star})$.

where $\mathcal{Z}_{\beta}$, known as the partition function, is a constant normalising the Gibbs measure $\mu_{\beta}$:

Note that $P_{g}$ and $P_{w}$ can be interpreted as a (unormalised) likelihood and prior distribution respectively. In the limit $\beta\to\infty$, the measure $\mu_{\beta}$ concentrates around solutions of the minimisation in eq. (A.3). The aim in the replica method is to compute the free energy density, defined as:

A.1 Replica computation of the free energy

The average in eq. (A.7) is not straightforward due to the logarithm term. The replica method consists of computing it using the following trick to get rid of the logarithm:

Averaging

Applying the trick above, the computation of the free energy density boils down to the evaluation of the averaged replicated partition function:

Note that the term in brackets defines the joint density over $(\nu_{\mu},\lambda^{a}_{\mu})$. It is easy to check that these are Gaussian random variables with zero mean and covariance matrix given by:

where the so-called overlap parameters $(\rho,m^{a},Q^{ab})$ are related to the weights $\bm{\theta}_{0},\bm{w}$:

We can therefore write the averaged replicated partition function as:

Rewriting as a saddle-point problem

The next step is to free the overlap parameters by introducing delta functions:

Inserting this in eq. (A.11) allow us to rewrite:

where we have absorbed a $-i$ factor in the integrals (this won’t matter since we will look to the saddle-point) and defined the potential:

where we recall that $\alpha=n/d$, $\gamma=p/d$ and:

In the high-dimensional limit where $d\to\infty$ while $\alpha=n/d$ and $\gamma=p/d$ stay finite, the integral in eq. (A.13) concentrate around the values of the overlaps that extremise $\Phi^{(r)}$, and therefore we can write:

Replica symmetric ansatz

In order to proceed with the $r\to 0^{+}$ limit, we restrict the extremisation above to the following replica symmetric ansatz:

Inserting this ansatz in eq. (A.14) allows us to explicitly take the $r\to 0^{+}$ limit for each term. The first three terms are straightforward to obtain. The limit of $\Psi_{y}^{(r)}$ is cumbersome, but it common to many replica computations for the generalised linear likelihood $P_{g}$. We refer the curious reader to Appendix C of or to Appendix IV of for details, and write the final result here:

Summary

The replica symmetric free energy density is simply given by:

A.2 Ridge regression and fixed weights

where we have included a convenient constant, and therefore:

taking the log and using $\log\det=\tr\log$, up to the limit:

Defining the shorthand $\text{A}=\left(\beta\lambda\text{I}_{d}+\hat{V}\Omega\right)^{-1}$, we can now take the averages over $\xi$ explicitly:

A.3 Taking the β→∞→𝛽\beta\to\infty limit

Finally, in order to take the $\beta\to\infty$ limit explicitly, we note that under the rescaling

The potential $\Psi_{w}$ has a trivial limit:

while $\Psi_{y}$ requires more attention. Since $\mathcal{Z}_{0}$ only depends on $(q,m)$, it is invariant under the rescaling. On the other hand, we have that:

where $\mathcal{M}$ is the Moreau envelope associated to the loss $g$:

The zero temperature therefore is simply given by:

A.4 Saddle-point equations

To solve the extremisation problem defined by eq. (A.34), we search for vanishing gradient points of the potential. This lead to a set of self-consistent saddle-point equations:

where $f_{g}(y,\omega,V)=-\partial_{\omega}\mathcal{M}_{Vg(y,\cdot)}(\omega)$, which can also be obtained from the proximal operator

using the envelope theorem $\mathcal{M}^{\prime}_{\tau f}(x)=\tau^{-1}\left(x-\mbox{prox}_{\tau f}(x)\right)$. A python implementation of the saddle-point equations for the losses discussed below is available in https://github.com/IdePHICS/GCMProject

A.5 Examples

We now discuss a couple of examples in which the equations above simplify.

For the linear task, the asymptotic training and generalisation errors read:

where $\rho=\frac{1}{p}\bm{\theta}_{0}^{\top}\Psi\bm{\theta}_{0}$ and $(V^{\star},q^{\star},m^{\star})$ are the fixed point of the following set of self-consistent equations:

Note that quite interestingly we have the following relationship between the training and generalisation error:

This give us an interesting interpretation of $V^{\star}$ as parametrising the variance gap between the generalisation and training errorWe thank Stéphane d’Ascoli for bringing this relation to our attention.. In particular, note that $V^{\star}$ only depends on the spectrum of the population covariance, since it is the solution of:

where $\nu_{\Omega}$ is the spectral density of $\Omega$.

Binary classification

where again $(q^{\star},m^{\star})$ are solutions of the self-consistent saddle-point equations. The teacher measure is given by:

The explicit form of the equation depends on the choice of the loss function, three of which are of particular interest:

As in the ridge case, for $g(y,x)=\frac{1}{2}(y-x)^{2}$ the saddle-point equations simplify considerably:

Similarly, the asymptotic training error also admits a simple expression:

Different from the previous cases, for logistic loss $g(y,x)=\log\left(1+e^{-yx}\right)$ the equations for $(\hat{V},\hat{q},\hat{m})$ cannot be integrated explicitly, since the proximal operator doesn’t admit a closed form solution. Instead, $f_{g}$ can be found by solving the following self-consistent equation:

Another useful case in which the proximal operator has a closed form solution is for the hinge loss $g(y,x)=\rm{max}(0,1-yx)$. In this case:

Again, the equations cannot be integrated explicitly. Note that in the limit $\lambda\to 0$, both the logistic and soft-margin solutions converge to the max-margin estimator

A.6 Relation to previous models

The feature map for random features learning can be written as:

These relations hold asymptotically, and rely on the Gaussian equivalence theorem (GET), see for a proof.

In , a similar Gaussian covariate model was used to study the performance of random feature regression on data generated from pre-trained generative models:

Let $\mathcal{H}$ be a Kernel Reproducing Hilbert space (RKHS) associated to a given kernel $K$ and $\mathcal{D}=\{\bm{x}^{\mu},y^{\mu}\}_{\mu=1}^{n}$ be a labelled data set with $\bm{x}\sim p_{x}$ independently, and set $\mathcal{X}=\rm{supp}(p_{x})$. In Kernel regression, the aim is to solve:

where $||\cdot||_{\mathcal{H}}$ is the norm induced by the scalar product in $\mathcal{H}$. An alternative representation of this problem is given by the feature decomposition of the kernel given by Mercer’s theorem:

where $\omega_{i}$ and $e_{i}(\bm{x})$ are the eigenvalues and eigenvectors associated with the kernel:

Note that $\{e_{i}(\bm{x})\}_{i}^{\infty}$ form an orthonormal basis of the space of square-integrable functions $L^{2}\left(\mathcal{X}\right)$ (with respect to the standard scalar product of $L^{2}$). It is also convenient to define the feature map $\varphi_{i}(\bm{x})=\sqrt{\omega_{i}}e_{i}(\bm{x})$, which is an orthonormal basis of $\mathcal{H}\subset L^{2}(\mathcal{X})$ (with respect to the scalar product induced by $K$). Therefore, if we assume that the labels $y^{\mu}=f_{0}(\bm{x}^{\mu})$ are generated from a ground truth target function (not necessarely part of $\mathcal{H}$), we can expand both $f$ and $f_{0}$ the feature basis:

Note that $f\in\mathcal{H}$ implies that for this sum to make sense $w_{i}$ needs to decay fast enough with respect to $\sqrt{w_{i}}$, but in general we can have $f_{0}\notin\mathcal{H}$ meaning that $\theta_{0}^{i}$ decays slower than $\sqrt{\omega_{i}}$ but still fast enough such that $f_{0}\in L^{2}(\mathcal{X})$. If the number of features is finite ($\omega_{i}=0$ for $i\geq d$) or if we introduce a cut-off $d\gg n,|\mathcal{X}|$, the representation in the feature basis in eq. (A.54) allow us to rewrite Kernel regression problem in eq. (A.51) simply as ridge regression in feature space:

Indeed, inserting this expression equation (A.38):

and making a change of variables $\hat{q}\leftarrow\hat{q}\frac{d^{2}}{n}$, $\hat{m}\leftarrow\hat{m}\frac{d}{n}$, $\hat{V}\leftarrow\hat{q}\frac{d}{n}$, $\rho\leftarrow d\rho$, $m\leftarrow dm$, $q\leftarrow dq$, $\lambda\leftarrow d\lambda$ we recover exactly the self-consistent equations of for the performance of kernel ridge regression directly from our equations. Moreover, our model allow to generalise this discussion to more involved kernel tasks such as kernel logistic regression and support vector machines.

Appendix B Rigorous proof of the main result

Intuitively, the variables $\mathbf{u}_{\rm{new}}^{\top}\bm{\theta}_{0}$ and $\mathbf{v}_{\rm{new}}^{\top}\mathbf{w}$ will play a key role in the analysis. Given an instance of $\bm{\theta}_{0}$ and $\mathbf{w}$, the tuple $\left(\frac{1}{\sqrt{p}}\mathbf{u}_{\rm{new}}^{\top}\bm{\theta}_{0},\frac{1}{\sqrt{d}}\mathbf{v}_{\rm{new}}^{\top}\mathbf{w}\right)$ is a bivariate Gaussian with covariance:

We thus define the following overlaps, that will play a fundamental role in the analysis:

Note that here, we will not introduce the spectral decomposition 2.2 as it will not simplify the expressions as in the $l_{2}$ case. The representations are mathematically equivalent nonetheless. Our main result is that the distribution of the estimator $\hat{\mathbf{w}}$ can be exactly computed in the weak sense from the solution to six scalar fixed point equations with a unique solution.

We start with a list of the necessary assumptions for the most generic version of the result to hold. We also briefly discuss how they are relevant in a supervised machine learning context.

The spectral distributions of the matrices $\Phi,\Psi$ and $\Omega$ converge to distributions such that the overlaps defined by equation (B.5) are well-defined. Additionally, the maximum singular values of the covariance matrices are bounded with high probability when $n,p,d\to\infty$.

The functions $r$ and $g$ are proper, lower semi-continuous, convex functions. Additionally, we assume that the cost function $r+g$ is coercive, i.e.:

The random elements of the function $f_{0}$ are independent of the matrices $U$ and $\mathcal{V}$. Additionally the following limit exists and is finite

When we send the dimensions $n,p,d$ to infinity, they grow with finite ratios $\alpha=n/d$, $\gamma=p/d$.

Additional assumptions for linear finite sample size rates : the teacher vector $\bm{\theta}_{0}$ has sub-Gaussian one dimensional marginals. The functions $r,g,\phi_{1},\phi_{2}$ are pseudo-Lipschitz of finite order. The eigenvalues of the covariance matrices are bounded with probability one.

Additional assumptions for exponential finite sample size rates: all of the above, and the loss function $g$ is separable and pseudo-Lipschitz of order 2, the regularisation is either a ridge or a Lipschitz function, the functions $\phi_{1},\phi_{2}$ are respectively separable, pseudo-Lipschitz of order 2, and a square or Lipschitz function.

The first assumption (A1) ensures that the teacher distribution is non-vanishing. The positive definiteness in (A2) means the covariance matrices of the blocks U and V are well-specified. Note that the cross-correlation matrix $\Phi$ can have singular values equal to zero. The assumption about the limiting spectral distribution is essentially a summability condition which is immediately verified if the limiting spectral distributions have compact support, a common case. The scaling assumptions from (A3) are natural as they imply that non-diverging inputs result in non-diverging outputs in the functions $f$ and $g$, as well as the sub-differentials. Similar scaling assumptions are encountered in proofs such as . They also allow to show Gaussian concentration of Moreau envelopes, as we will see in Lemma 5. The coercivity assumption is verified in most common machine learning setups : any convex loss with ridge regularisation, or any convex loss that is bounded below with a coercive regularisation (LASSO, elastic-net,…), see Corollary 11.15 from . Assumption (A4) is a classical assumption of teacher-student setups, where any correlation between the teacher and the student is modeled by the covariance matrices and not by the label generating function $f_{0}$. The summability condition ensures generalization error is well-defined for squared performance measures. Finally, (A5) is the typical high-dimensional limit used in statistical physics of learning, random matrix theory and a large recent body of work in high-dimensional statistical learning.

B.2 Main theorem

First, let’s define quantities and a scalar optimization problem that will be used to state the asymptotic behaviour of (1.2-1.3):

(Scalar potentials/replica free energy) Define the following functions of the scalar variables $\tau_{1}>0,\tau_{2}>0,\kappa\geqslant 0,\eta\geqslant 0,\nu,m$:

where $\mathbf{s},\mathbf{h}\sim\mathcal{N}(0,\text{I}_{n})$ and $\mathbf{g}\sim\mathcal{N}(0,\text{I}_{d})$ are random vectors independent of the other quantities, $\mathbf{t}=\Omega^{-1/2}\Phi^{\top}\bm{\theta}_{0}$, $\mathbf{y}=\bm{f}_{0}\left(\sqrt{\rho}\mathbf{s}\right)$, and $\mathcal{M}$ denotes the Moreau envelope of a target function.

From these quantities define the following potential:

Under Assumption (B.1), the previously defined quantities all admit finite limits when $n,p,d\to\infty$.

Proof: This follows directly from Lemma 5.

The next lemma characterizes important properties of the "potential" function $\mathcal{E}(\tau_{1},\tau_{2},\kappa,\eta,\nu,m)$:

(Geometry and minimizers of $\mathcal{E}$) The function $\mathcal{E}(\tau_{1},\tau_{2},\kappa,\eta,\nu,m)$ is jointly convex in $(m,\eta,\tau_{1})$ and jointly concave in $(\nu,\kappa,\tau_{2})$, and the optimization problem

has a unique solution $(\tau_{1}^{*},\tau_{2}^{*},\kappa^{*},\eta^{*},\nu^{*},m^{*})$ on $\mbox{dom}(\mathcal{E})$.

Proof: see Appendix 13. The optimality condition of problem (B.12) yields the set of self-consistent fixed point equations given in Lemma 12 of Appendix B. Finally, define the following variables:

where prox denotes the proximal operator. With these definitions, we can now state our main result:

(Training loss and generalisation error) Under Assumption (B.1), there exist constants $C,c,c^{\prime}>0$ such that, for any optimal solution $\hat{\mathbf{w}}$ to (1.3), the training loss and generalisation error defined by equation verify, for any $0<\epsilon<c^{\prime}$:

where $\mathcal{E}_{train}^{*}$ is defined as follows:

and the random variables $(\omega,\xi)$ are jointly Gaussian with covariance

Proof: see Appendix B.6. Note that the regularisation may be removed to evaluate the training loss. A more generic result, aiming directly at the estimator $\hat{\mathbf{w}}$, can also be stated:

Proof: see Appendix B.6. Concentration still holds for a larger class of functions $\phi_{1,2}$, but exponential rates are lost. This is discussed in Appendix B.1.

B.3 Theoretical toolbox

Here we remind a few known results that are used throughout the proof. We also provide proofs of useful, straightforward consequences of theses results that do not appear explicitly in the literature for completeness.

We start with the Convex Gaussian Min-max Theorem, as presented in , which is a tight version of an inequality initially derived in .

Following , we will say that any reformulation of a target problem matching the form of (B.19) is an acceptable primary optimization problem (PO), and the corresponding form (B.20) is an acceptable auxiliary problem (AO). The main idea of this approach is to study the asymptotic properties of the (PO) by studying the simpler (AO).

B.3.2 Proximal operators and Moreau envelopes : differentials and useful functions

As reminded in , the Moreau envelope is jointly convex in $(\tau,\mathbf{x})$ and differentiable almost everywhere, with gradients:

We remind that $\mbox{prox}_{\tau f}(\mathbf{x})$ is the unique point which solves the strongly convex optimization problem defining the Moreau envelope, i.e.:

We also remind the definition of order k pseudo-Lipschitz function.

We now give some further properties that will be helpful throughout the proof.

Proof of Lemma 2: For any $\mathbf{x},\mathbf{y}$ in $\mbox{dom}(f)$, we have, using the pseudo-Lipschitz property:

where the second line follows immediately with the same constant $L(k)$ owing to the firm-nonexpansiveness of the proximal operator. Furthermore

due to the pseudo-Lipschitz property, one has

This, along with the firm-nonexpansiveness of $\rm{Id}-\mbox{prox}$, concludes the proof. ∎

$h_{1}$ is nondecreasing, and $h_{2},h_{3},h_{4}$ are nonincreasing.

which implies that $h_{2}$ is non-increasing. Using the Moreau decomposition, see e.g. , we have:

Proof of Lemma 4 : the subdifferential of a proper convex function is a monotone operator, thus:

additionally, $\mbox{prox}_{\gamma f}(\mathbf{x})=\left(\rm{Id}+\gamma\partial f\right)^{-1}(\mathbf{x})$, hence:

B.3.3 Useful concentration of measure elements

We begin by reminding the Gaussian-Poincaré inequality, see e.g. .

We now use this previous result to show Gaussian concentration of Moreau envelopes of appropriately scaled convex functions.

where the second line is integrable under a multivariate Gaussian measure. Then, using Proposition 1, we get:

Using Proposition 12.27 and Corollary 4.3 from , $\mathbf{g}\to\mathbf{z}-\mbox{prox}_{\tau f}(\mathbf{g})$ is firmly non-expansive and:

Chebyshev’s inequality then gives, for any $\epsilon>0$:

Gaussian concentration of pseudo-Lipschitz functions of finite order can also be proven using the Gaussian Poincaré inequality to yield a bound similar to the one obtained for Moreau envelopes. We thus give the result without proof:

We now cite an exponential concentration lemma for separable, pseudo-Lipschitz functions of order 2, taken from .

where it is understood that $f(\mathbf{g})=\sum_{i=1}^{n}f(g_{i})$.

B.4 Determining a candidate primary problem, auxiliary problem and its solution.

We start with a reformulation of the problem (1.2-1.3) in order to obtain an acceptable primary problem in the framework of Theorem 6. Partitioning the Gaussian distribution, we can rewrite the matrices U and $\mathcal{V}$ in the following way, introducing the standard normal vector:

We can then rewrite the vectors $\mathbf{u},\mathbf{v}$ and matrices $U,\mathcal{V}$ as:

where the matrices $A$ and $B$ have independent standard normal entries and are independent of $\bm{\theta}_{0}$. The learning problem then becomes equivalent to :

We are then interested in the optimal cost of the following problem

Introducing the auxiliary variable $\mathbf{z}$:

In the remainder of the proof, the preceding cost function will be denoted

such that the problem reads $\min_{\mathbf{w},\mathbf{z}}\mathbf{C}(\mathbf{w},\mathbf{z})$. Theorem 6 requires working with compact feasibility sets. Adopting similar approaches to the ones from , the next lemma shows that the optimization problem (B.61) can be equivalently recast as one over compact sets.

(Compactness of feasibility set) Let $\mathbf{w}^{*},\mathbf{z}^{*},\bm{\lambda}^{*}$ be optimal in (B.61). Then there exists positive constants $C_{\mathbf{w}},C_{\mathbf{z}}$ and $C_{\bm{\lambda}}$ such that

Proof of Lemma 8: consider the initial minimisation problem:

Now consider the equivalent formulation of problem (B.64):

The optimality condition in $\bm{\lambda}$ gives:

According to assumption (A2), the operator norms of the matrices involving the covariance matrices are bounded with high probability and using known results on random matrices, see e.g. , the operator norms of $\frac{1}{\sqrt{d}}A$ and $\frac{1}{\sqrt{d}}B$ are bounded by finite constants with high probability when the dimensions go to infinity. Thus there exists a constant $C_{\mathbf{z}}$ also independent of d such that:

Finally, the scaling condition from assumption (A3) directly shows that there exists a constant $C_{\bm{\lambda}}$ such that

The rest of this section can then be summarized by the following lemma, the proof of which shows how to find an acceptable (PO) for problem (B.71), the corresponding (AO) and how to reduce the (AO) to a scalar optimization problem. At this point we will assume the teacher vector $\bm{\theta}_{0}$ is deterministic, and relax this assumption in paragraph B.7. For this reason we do not add it to the initial list of assumptions in section B.1.

(Scalar equivalent problem) In the framework of Theorem 6, acceptable (AO)s of problem (B.71) can be reduced to the following scalar optimization problems

Proof of Lemma 9: We need to find an i.i.d. Gaussian matrix independent from the rest of the problem in order to use Theorem 6. We thus decompose the mixing matrix A by taking conditional expectations w.r.t. $\mathbf{y}$, which amounts to conditioning on a linear subset of the Gaussian space generated by A. Dropping the feasibility sets for confort of notation in the following lines:

Two cases must now be considered, $\bm{\theta}_{0}\notin\mbox{Ker}(\phi^{\top})$ and $\bm{\theta}_{0}\in\mbox{Ker}(\phi^{\top})$. Another possible case is $\Phi=0_{p\times d}$, however it leads to the same steps as the case $\bm{\theta}_{0}\in\mbox{Ker}(\Phi^{\top})$.

A follow-up of the previous equations shows that the feasibility set now reads :

We now turn to the simplification of this problem. The variable $\bm{\lambda}$ only appears in linear terms, we can thus directly optimize over its direction, introducing the positive scalar variable $\kappa=\norm{\bm{\lambda}}_{2}/\sqrt{d}$:

The previous expression may not be convex-concave because of the term $\norm{\mathbf{p}}_{2}\mathbf{h}$. However, it was shown in that the order of the min and max can still be inverted in this case, because of the convexity of the original problem. As the proof would be very similar, we do not reproduce it. Inverting the max-min order and performing the linear optimization on $\mathbf{p}$ with $\eta=\norm{\mathbf{p}}_{2}/\sqrt{d}$:

using the following representation of the norm, as in , for any vector $t$, $\norm{t}_{2}=\min_{\tau>0}\frac{\tau}{2}+\frac{\norm{t}_{2}^{2}}{2\tau}$:

performing the minimisation over $\mathbf{z}$ and recognizing the Moreau envelope of $g(.,\mathbf{y})$:

At this point we have a convex-concave problem. Inverting the min-max order, $\bm{\mu}$ appears in a well defined strictly convex least-square problem.

Isolating the terms depending on $\mathbf{w}_{\perp}$, we get a strictly convex least-square problem, remembering that $\Omega\in\mathcal{S}_{d}^{++}$:

Another strictly convex least-square problem appears on $\bm{\eta}$, the solution and optimal value of which read

At this point we have expressed feasible solutions of $\bm{\eta},\mathbf{w}_{\perp}$ as functions of the remaining variables. For any feasible solution in those variables, $\mathbf{w}$ and $\mathbf{w}^{\prime}$ are the same. Replacing in the (AO) and a completion of squares leads to

where the Moreau envelopes of $f$ and $g$ are respectively defined w.r.t. the variables $\mathbf{w^{\prime\prime}}$ and $\mathbf{z}$. At this point we have reduced the initial high-dimensional minimisation problem (B.4) to a scalar problem over six parameters. Another follow-up of the feasibility set shows that there exist positive constants $C_{m},C_{\kappa},C_{\eta}$ independent of $n,p,d$ such that $0\leqslant\kappa\leqslant C_{\kappa}$, $0\leqslant\eta\leqslant C_{\eta}$ and $0\leqslant m\leqslant C_{m}$.

In this case, the min-max problem (B.85) becomes:

where the compactness of the feasibility set is preserved almost surely from the almost sure boundedness of the eigenvalues of $\Omega$. We can thus write the corresponding auxiliary optimization problem, reintroducing the normalization by d:

introducing the convex conjugate of $f$ with dual parameter $\bm{\eta}$:

Using the square root trick with parameters $\tau_{1},\tau_{2}$:

performing the optimizations on $\mathbf{z},\bm{\eta}$ and recognizing the Moreau envelopes, the problem becomes:

B.5 Study of the scalar equivalent problem : geometry and asymptotics.

Here we study the geometry, solutions and asymptotics of the scalar optimization problem (B.4). We will focus on the case $\bm{\theta}_{0}\notin\mbox{Ker}(\Phi^{\top})$ as the other case simply shows that no learning is performed (see the remark at the end of this section). The following lemma characterizes the continuity and geometry of the cost function $\mathcal{E}_{n}$.

(Geometry of $\mathcal{E}_{n}$) Recall the function:

Then $\mathcal{E}_{n}(\tau_{1},\tau_{2},\kappa,\eta,\nu,m)$ is continuous on its domain, jointly convex in $(m,\eta,\tau_{1})$ and jointly concave in $(\kappa,\nu,\tau_{2})$.

Proof of Lemma 10 : $\mathcal{E}_{n}(\tau_{1},\tau_{2},\kappa,\eta,\nu,m)$ is a linear combination of linear and quadratic terms with Moreau envelopes, which are all continuous on their domain. Remembering the formulation

The squared term in $m,\eta,\mathbf{z}$ can be written as

(Asymptotics of $\mathcal{E}_{n}$) Recall the following quantities:

and $\mathcal{E}(\tau_{1},\tau_{2},\kappa,\eta,\nu,m)$ is continuously differentiable on its domain, jointly convex in $(m,\eta,\tau_{1})$ and jointly concave in $(\kappa,\nu,\tau_{2})$.

using the Gaussian tail. The Borel-Cantelli lemma and summability of this tail gives

Concentration of the Moreau envelopes of both $f$ and $g$ follows directly from lemma 5. We thus have the pointwise convergence:

Since pointwise convergence preserves convexity, $\mathcal{E}(\tau_{1},\tau_{2},\kappa,\eta,\nu,m)$ is jointly convex in $(m,\eta,\tau_{1})$ and jointly concave in $(\kappa,\nu,\tau_{2})$. Now recall the expression of $\mathcal{E}$

The feasibility sets of $\kappa,\eta,m$ are compact from Lemma 8 and the subsequent follow-up of the feasibility sets. Then, using Proposition 12.32 from , for fixed $(\tau_{2},\kappa,\eta,\nu,m)$, we have:

which is a finite quantity since $g(.,\mathbf{y})$ is a proper, convex function verifying the scaling assumptions B.1. Then, since $\kappa>0$, we have:

Similarly, for fixed $(\tau_{1},\kappa,\eta,\nu,m)$ and noting that composing $f$ with the positive definite matrix $\Omega^{-1/2}$ does not change its convexity, or it being proper and lower semi-continuous, we get:

which is also a bounded quantity from the scaling assumptions made on $f$. Since $\beta>0$, we then have:

Finally, the limit $\lim_{\nu\to+\infty}\mathcal{E}_{n}(\tau_{1},\tau_{2},\kappa,\eta,\nu,m)$ needs to be checked for both $+\infty$ and $-\infty$ since there is no restriction on the sign of $\nu$. From the definition of the Moreau envelope, we can write:

Thus, for any fixed $(\tau_{1},\tau_{2},m,\kappa,\eta)$:

which immediately gives $\lim_{\nu\to-\infty}\mathcal{E}_{n}=-\infty$. Turning to the other limit, remembering that $\mathcal{E}_{n}$ is continuously differentiable on its domain, we have:

Thus $\lim_{\nu\to+\infty}\frac{\partial\mathcal{E}_{n}(\tau_{1},\tau_{2},\kappa,\eta,\nu,m)}{\partial\nu}\to-\infty$. Since $\mathcal{E}_{n}$ is continuously differentiable in $\nu$ on $[0,+\infty[$, and from the short argument led above, we have shown

Using similar arguments as in the proof of Lemma 8, we can now reduce the feasibility set of $\tau_{1},\tau_{2},\nu$ to a compact one. Then, using the fact that convergence of convex functions on compact sets implies uniform convergence , we obtain

which is the desired result. ∎ At this point, it is necessary to characterize the set of solutions of the asymptotic minimisation problem (B.12). We start with the explicit form of the optimality condition associated to any solution.

(Fixed point equations) The zero-gradient condition of the optimization problem (B.12) prescribes the following set of fixed point equations for any feasible solution:

This set of equations can be converted to the replica notations using the table (C.2).

Proof of Lemma 12: Using arguments similar to the ones in the proof of Lemma 5, Moreau envelopes and their derivatives verify the necessary conditions of the dominated convergence theorem. Additionally, uniform convergence of the sequence of derivatives can be verified in a straightforward manner as all involved functions are firmly non-expansive and integrated w.r.t. Gaussian measures. We can therefore invert the limits and derivatives, and invert expectations and derivatives. We can now write explicitly the optimality condition for the scalar problem (B.126), using the expressions for derivatives of Moreau envelopes from Appendix B.3. Some algebra and replacing with prescriptions obtained from each partial derivative leads to the set of equations above. ∎ Remark : Here we see that the potential function (B.126) can be further studied using the fixed point equations (12) and the relation (B.24). For any optimal $(\tau_{1},\tau_{2},\kappa,\eta,\nu,m)$, it holds that

Finally, we give a strict-convexity and strict-concavity property of the asymptotic potential $\mathcal{E}$ which will be helpful to prove Lemma 1.

(Strict convexity and strict concavity near minimisers) Consider the asymptotic potential function $\mathcal{E}(\tau_{1},\tau_{2},\kappa,\eta,\nu,m)$. Then for any fixed $(\eta,m,\tau_{1})$ in their feasibility sets, the function

is jointly strictly concave in $(\tau_{2},\kappa,\nu)$. Additionally, consider the set $\mathcal{S}_{\partial_{\nu,\tau_{2}}}$ defined by:

then for any fixed $\tau_{2},\kappa,\nu$ in $\mathcal{S}_{\partial_{\nu,\tau_{2}}}$, the function $(\eta,m,\tau_{1})\to\mathcal{E}(\tau_{1},\tau_{2},\kappa,\eta,\nu,m)$ is jointly strictly convex in $(\eta,m,\tau_{1})$ on $\mathcal{S}_{\partial_{\nu,\tau_{2}}}$

Proof of Lemma 13: We will use the following first order characterization of strictly convex functions: $f\thickspace\mbox{is strictly convex}\iff\langle\mathbf{x}-\mathbf{y}|\nabla f(\mathbf{x})-\nabla f(\mathbf{y})\rangle>0\thickspace\forall\mathbf{x}\neq\mathbf{y}\in\mbox{dom}(f)$. To simplify notations, we will write, for any fixed $(m,\eta,\tau_{1})$

As we will see in the next section, this will lead to estimators solely based on noise.

B.6 Back to the original problem : proof of Theorem 4 and 5

We begin this part by considering that the "necessary assumptions for exponential rates" from the set of assumptions B.1 are verified. In the end we will discuss how relaxing these assumptions modifies the convergence speed. We closely follow the analysis introduced in and further developed in . The main difference resides in checking the concentration properties of generic Moreau envelopes depending on the regularity of the target function instead of specific instances such as the LASSO. Since the dimensions $n,p,d$ are linked by multiplicative constants, we can express the rates with any of the three. Recall the original reformulation of the problem defining the student.

Consider the finite size scalar optimization problem

where the feasibility set of $(\tau_{1},\tau_{2},\kappa,\eta,\nu,m)$ is compact and $\tau_{1}>0,\tau_{2}>0$. Then any optimal values $\kappa^{*},\tau_{2}^{*}$ verify:

For any $\epsilon>0$, there exist constants $C,c,\gamma$ such that:

Proof of Lemma 15: for any fixed $(\tau_{1},\tau_{2},\kappa,\nu,\eta,m)$, we can determine the rates of convergence of all the random quantities in $\mathcal{E}_{n}$. The linear terms involving $\frac{1}{d}\mathbf{g}^{T}\mathbf{v}$ are sub-Gaussian with sub-Gaussian norm bounded by $C/d$ for some constant $C>0$. Thus we can find constants, $C,c>0$ such that, for any $\epsilon>0$ :

The term involving $\mathbf{v}^{T}\Omega\mathbf{v}$ is deterministic in this setting. We will see in section B.7 how a random $\bm{\theta}_{0}$ affects the convergence rates. The term involving $\frac{1}{d}\mathbf{g}^{T}\mathbf{g}$ is a weighted sum of sub-exponential random variables, the tail of which can be determined using Bernstein’s inequality, see e.g. Corollary 2.8.3, which gives a sub-Gaussian tail for small deviations and a sub-exponential tail for large deviations. Parametrizing the deviation $\epsilon$ with a scalar variable $c^{\prime}$, we thus get the following bound : for any $\epsilon>0$, there exists constants $C,c,c^{\prime}>0$ such that:

For any $\epsilon>0$, there exist constants $C,c,\gamma$ such that:

For any $\epsilon>0$, there exists constants $\gamma,c,C>0$ such that the event

has probability at most $\frac{C}{\epsilon}e^{-cn\epsilon^{2}}$.

For any $\epsilon>0$, there exists constants $\gamma,c,C>0$ such that the event

has probability at most $\frac{C}{\epsilon^{2}}e^{-cn\epsilon^{4}}$.

which proves Theorem 5 using the fact that $\hat{\mathbf{w}},\hat{\mathbf{z}}$ are minimizers of the initial cost function. Theorem 4 is a consequence of Theorem 5. If the restriction on $f,g,\phi_{1},\phi_{2}$ are relaxed to any pseudo-Lipschitz functions of finite orders, the exponential rates involving them are lost and become linear following Lemma 5.

B.7 Relaxing the deterministic teacher assumption

The entirety of the previous proof has been done with a deterministic vector $\bm{\theta}_{0}$. Now, if $\bm{\theta}_{0}$ is assumed to be a random vector independent of all other quantities, as prescribed in the set of assumptions B.1, we can "freeze" the variable $\bm{\theta}_{0}$ by conditioning on it. The whole proof can then be understood as studying the value of the cost conditioned on the value of $\bm{\theta}_{0}$. Note that, in the Gaussian case, correlations between the teacher and student are expressed through the covariance matrices, thus leaving the possibility to parametrise the teacher with a vector $\bm{\theta}_{0}$ indeed independent of all the rest. To lift the conditioning in the end, one only needs to average out on the distribution of $\bm{\theta}_{0}$, the summability conditions of which are prescribed in the set of assumptions B.1. Thus, random teacher vectors can be treated simply by taking an additional expectation in the expressions of Theorem 5, provided $\bm{\theta}_{0}$ is independent of the matrices $A,B$ and the randomness in $f_{0}$. As mentioned at the end of the previous section, the finite size rates will be determined by the assumptions made on the teacher vector and decay of the eigenvalues of the covariance matrices. We do not investigate in detail the limiting assumptions under which exponential rates still hold regarding the randomness of the teacher or tails of the eigenvalue distributions of covariance matrices.

B.8 The ’vanilla’ teacher-student scenario

In this section, we give the explicit forms of the fixed points equations and optimal asymptotic estimators in the case where the teacher and the student are sampled from the same distribution, i.e. $\Omega=\Phi=\Psi=\Sigma$ where $\Sigma$ is a positive definite matrix with sub-Gaussian eigenvalue decay. This setup was rigorously studied in for the LASSO and heuristically in for the ridge regularized logistic regression. In this case, the fixed point equations become

and the asymptotic optimal estimators read:

Appendix C Equivalence replica-Gordon

In this Appendix, we show that the rigorous result of Theorem 5 can be used to prove the replica prediction in the case of a separable loss, a ridge penalty. For simplicity, we restrict ourselves to the case of random teacher weights with $\bm{\theta}_{0}\sim\mathcal{N}(0,\text{I}_{p})$. We provide an exact analytical matching between the replica prediction and the one obtained with Gordon’s theorem. We start by an explicit derivation of the form presented in Corollary 1 from the main result (1).

using Lemma 5 with a separable function, the expectation over the Moreau envelope converges to:

where $s$ and $h$ are standard normal random variables and $y=f_{0}(\sqrt{\rho}s)$. The corresponding optimality conditions then reads:

simplifying these equations using Stein’s lemma, we get:

C.2 Matching with Replica equations

In this section, we show that the fixed point equations obtained from the asymptotic optimality condition of the scalar minimization problem 1 match the ones obtained using the replica method. In what follows we will use the same notations as in , and an explicit, clear match with the notations from the proof of the main theorem will be shown. The replica computation, similar to the one from , leads to the following fixed point equations, in the replica notations:

where $f_{g}(y,\omega,V)=-\partial_{\omega}\mathcal{M}_{Vg(y,\cdot)}(\omega)$ and $\mathcal{Z}_{0}$ is given by:

To be explicit with the notation, let’s open the equations up. Take for instance the one for $\hat{m}$. Opening all the integrals:

where in $(a)$ we integrated over $y$ explicitly. A direct comparison between the two sets of equations suggests the following mapping to navigate between the replica derivation and the proof using Gaussian comparison theorems. We denote replica quantities with Rep indices:

The first three equations match the replica prediction, the last three can be exactly matched using the following change of variable and Gaussian integration:

Appendix D Details on the simulations

In this Appendix we give full details on the numerics used to generate the plots in the main manuscript. An implementation of all the pipelines described below is available at https://github.com/IdePHICS/GCMProject.

First, we discuss in detail how we conducted the numerical simulations in Figs. 4 and 4 in the main manuscript.

In Fig. 4, the feature $\bm{\varphi}_{s}^{t}$ is taken from a learned neural network at different epochs $t\in\{0,5,50,200\}$ of training. For this experiment, we chose the following architecture implemented in \pythPytorch: {python} Sequential( (0): Linear(in_features=784, out_features=2352, bias=False) (1): ReLU() (2): Linear(in_features=2352, out_features=2352, bias=False) (3): ReLU() (4): Linear(in_features=2352, out_features=1, bias=False) )

In both experiments, we ran ridge regression at fixed regularisation $\lambda>0$ by sub-sampling $n$ samples from the data set $\mathcal{D}=\{\bm{v}^{\mu},y^{\mu}\}_{\mu=1}^{n_{\text{tot}}}$, $\bm{v}^{\mu}=\bm{\varphi}_{s}\left(\bm{x}^{\mu}\right)$, with the estimator given by the closed-form expression:

In principle, the teacher weights need to be estimated by inverting $\bm{y}=\text{U}\bm{\theta}_{0}$. However, as explained in $n_{\text{tot}}$ in Sec. 3.4, one can avoid doing so by noting the teacher weights only appear in the self-consistent equations 2.8 through $\rho=\frac{1}{k}\bm{\theta}_{0}^{\top}\Psi\bm{\theta}_{0}$ and $\Phi^{\top}\bm{\theta}_{0}$. Therefore, all teacher vector $\bm{\theta}_{0}$ and feature map $\varphi_{s}$ that linearly interpolate the data set $\{\bm{x}^{\mu},y^{\mu}\}_{\mu=1}^{n_{\text{tot}}}$ are equivalent, since we can write:

which is independent from $\left(\bm{\varphi}_{t},\bm{\theta}_{0}\right)$. In particular, note that for our binary labels $y^{\mu}\in\{+1,-1\}$, we have $\rho=1$. In both Fig. 4 and 4 of the main, we estimated the covariance $\Omega$ as in eq. (D.3) by applying the feature maps $\bm{\varphi}_{s}$ described above to the whole data set, took $\rho=1$ (since in both we have binary labels) and used eq. (D.4) to estimate $\Phi^{\top}\bm{\theta}_{0}$. This was then fed to our iterator package (https://github.com/IdePHICS/GCMProject) to compute the curves. For the kernel curve, we used the random features approximation of eq. (D.1) with a $d=20\times 1953$ dimensional feature space to estimate the covariance $\Omega$. We have checked that this indeed provide a good approximation of $K$ for the sample range considered, see Fig.5.

As we have discussed above, a key ingredient of our theoretical analysis is the estimation of the population covariances. For real data, this relies on the empirical covariance of the whole data set with $n_{\text{tot}}$ samples. We expect this approximation to be good only for $n\ll n_{\text{tot}}$ samples, as it is the case for the ranges plotted in Figs. 4 and 4. Indeed, as $n\approx n_{\text{tot}}$ we start observing deviations between the theoretical prediction and the simulations. In Fig. 6 (right) we show an example of a NTK kernel regression task on 8 vs 9 MNIST digit classification, for which $n_{\text{tot}}=7000$. Note that while the theoretical prediction reach perfect generalisation at $n\approx n_{\text{tot}}$, the simulated error approaches a plateau. Alternatively, instead of varying the sample range, in Fig. 6 (left) we show how the matching betweem theory and simulation degrades by varying $n_{\text{tot}}$ on a fixed sample range for a MNIST odd vs. even task.

As it was discussed in Sec. 3.4 of the main manuscript, the universality argument sketched above is only valid in the case of a linear student. For instance, applying the same construction to a binary classification task with $f_{0}(x)=\hat{f}(x)=\text{sign}(x)$ lead to a mismatch between theory and experiments, as exemplified in Fig. 3 of the main for a logistic regression task on CIFAR10 gray-scale images. Interestingly, this is even the case for binary classification with the square loss $g(x,y)=\frac{1}{2}(x-y)^{2}$, in which the estimator $\hat{\bm{w}}$ is the same as for ridge regression. In other words, by simply changing the predictor $\hat{f}(x)=\text{sign}(x)$, we have a breakdown of universality, as shown in Fig. 7.

D.2 Binary classification on GAN generated data

and which has been trained on the full CIFAR10 data set. It therefore takes a $100$-dimensional latent vector and returns a $D=32\times 32\times 3=3072$ CIFAR10-looking image. The GAN was trained on the original CIFAR10 data set without data augmentation for 50 epochs. Both the discriminator and the generator were trained using Adam, with Adam parameters $\beta_{1}=0.5$ and $\beta_{2}=0.999$. In practice, the advantage of working with a GAN is that we have a generative process to sample as many independent data points as we need, both for the simulations and for the estimation of the population covariances.

The experiment shown in Fig. 3 follow a similar pipeline as the one described in Sec. D.1. The student feature maps $\bm{\varphi}_{s}^{t}$ are obtained by removing the last layer of a trained a 3-layer student network with architecture: {python} Sequential( (0): Linear(in_features=1024, out_features=2304, bias=False) (1): ReLU() (2): Linear(in_features=2304, out_features=2304, bias=False) (3): ReLU() (4): Linear(in_features=2304, out_features=1, bias=False) ) Training was performed on a data set composed of $n=30000$ independent samples drawn from the dcGAN described above, with labels $y^{\mu}\in\{+1,-1\}$ assigned by the learned teacher $y^{\mu}=\text{sign}\left(\bm{u}^{\top}\bm{\theta}_{0}\right)$, $\bm{u}^{\mu}=\bm{\varphi}_{t}\left(\bm{x}^{\mu}\right)$. The network was trained for $300$ epochs using Adam optimiser on the MSE loss and \pythpyTorch’s default Kaiming initialisation, and snapshops of the weights were extracted at epochs $t\in\{0,5,50,200\}$. Finally, logistic regression was performed on the learned features $\bm{v}=\bm{\varphi}_{s}(\bm{x})$ on fresh pair of dcGAN generated samples and labels using the out-of-the-box \pythLogisticRegression solver from \pythScikit-learn. The points and error bars in Fig. 3 were computed by averaging over $10$ independent runs. The same pipeline was used for Fig. 3, but for $\bm{\varphi}_{s}=\rm{id}$ and on dcGAN generated CIFAR10 gray-scale images.

As before, the self-consistent eqs. 2.8 require the population covariances $(\Omega,\Phi,\Psi)$ and the teacher weights $\bm{\theta}_{0}$. For synthetic GAN data, the population covariances of the feature maps $\left(\bm{\varphi}_{t},\bm{\varphi}^{t}_{s}\right)$ used in the simulations can be estimated as well as needed with a Monte Carlo sampling algorithm. For the curves shown in Figs. 3 and 3, the covariances were estimated with $n=10^{6}$ samples with a precision of the order of $10^{-5}$. Together with the teacher weights $\bm{\theta}_{0}$ used to generate the labels, this provides everything needed to compute the theoretical learning curves from the self-consistent equations.

Appendix E Ridge regression with linear teachers

In this Appendix we discuss briefly random matrix theory, and consider heuristic reasons behind the validity of our asymptotic result beyond Gaussian covariates $(\bf{u},\bf{v})$ in the context of ridge regression, with linear teacher. As is well known, the computation of the training and test MSE for ridge regression can be written as a random matrix theory problem. We do not attempt a rigorous approach, but rather to motivate with simple arguments, many of them actually well known, the observed universality and its limits.

We assume the existence of a linear teacher generating the labels ${\bf y}=\text{U}\bm{\theta}_{0}$, and recall the student performs ridge regression on the data matrix $\mathcal{V}$.

The Gaussian model we consider can therefore be rewritten with a Rademacher vector $\bm{\theta}$ provided we change the correlation matrix ${\Psi}$ (as well as the cross-correlation $\Phi$) accordingly.

where we have defined the empirical covariance matrices

Given this vector, one can now readily write the expected value of the training and test losses as follows:

where we have denoted the population correlation matrices $\Psi\equiv\Sigma_{u,u},\Phi\equiv\Sigma_{u,v},\Omega\equiv\Sigma_{v,v}$ for readability and a direct comparison with their empirical counterpart. The traces appears by the left and right multiplication by the random vector $\bm{\theta}$.

At this point, the entire problem has been mapped to a random matrix theory exercise: assuming data are indeed Gaussian, one can use RMT to compute the six traces that appears in (E.5,E.6). Indeed, this is the canonical approach used in most rigorous works for the ridge regression task in the teacher-student framework, instance in . Remarkably, the replica (and the rigorous Gordon counterpart) allow to find the same result without the explicit use of RMT.

We now discuss, heuristically, why these results are valid even though the distribution of $[\bf u,v]$ is not actually Gaussian, and in some instances even for real data. Indeed, that both $\left(\mathcal{E}_{\rm{train.}},\mathcal{E}_{\rm{gen.}}\right)$ do not depend explicitly on the distribution of the data, but —assuming some concentration (or self-averaging)— only on:

The spectrum of the population covariances $\Sigma_{u,u},\Sigma_{u,v},\Sigma_{v,v}$.

The spectrum of the empirical covariances $\text{S}_{u,u},\text{S}_{u,v},\text{S}_{v,v}$.

The expectation of the trace of products between empirical and population covariances.

We expect that asymptotically the prediction from the theory will thus be valid for much more generic distributions $[{\bf u,\bf v}]\sim P_{{\bf u,v}}$, provided they share the same population covariances $\Sigma_{u,u},\Sigma_{u,v},\Sigma_{v,v}$ (which we call $\Psi,\Phi,\Omega$). To see this, we need to check how this change in distribution would affect points (1),(2) and (3). Fixing the population covariances, the first bullet point (1) is automatically taken into account. Point (2) and (3) are, however, less trivial: in order to have universality we need that a) the spectrum of the empirical covariances of the non-Gaussian distribution to converge the one obtained with the Gaussian one; and b) the trace of products between the empirical and the population covariances also to converge to the universal values computed from Gaussians data.

These two last points have been investigated in RMT [anderson2010introduction], and it is a classical result that such quantities are universal and converge to the Gaussian-predicted values for many distribution, way beyond the Gaussian assumption (in which case the spectral densities are known as the Wigner and Wishart model, or Marcenko-Pastur distribution [marchenko1967distribution]): this powerful universality of RMT is at the origin of the applicability of the model beyond Gaussian data. For instance, showed that these assumptions are verified for any data generated as $\bf u=\Sigma_{u,u}^{1/2}\bm{\omega}$, assuming the components of the vector $\bm{\omega}$ are drawn i.i.d. from any distribution (with some assumption on the larger moments). While this is still restrictive, stronger results can be shown, and [63, 65, chafai2018convergence, 49] extended them (also loosening the independence assumption) for a very generic class of distributions of correlated random vectors $\bf u$.

Let us give a concrete example. For simplicity, consider the restricted case where ${\bf u}={\bf v}$, i.e. the teacher acts on the same space as the student. In this case, eqs. (E.5,E.6) simplify (this is essentially the analysis in ) to:

In the expression of the training loss eq. (E.7), we see terms such as

In the expression of the generalisation loss eq. (E.8), however, terms such as

appears. These can be computed using classical RMT results on the concentration of the inverse of the covariance [hachem2007deterministic, 64]. The strongest result we are aware of for such problems is from the remarkable work of . This universality of random matrix theory is thus at the origin of the surprisingly successful application of our Gaussian theory to real data with arbitrary feature maps. Of course the discussion here is limited to the case where ${\bf u}={\bf v}$ and a concrete mathematical statement would require the generalisation of these arguments to the more generic case of eqs.(E.5,E.6), which are closer to the work of . We leave this discussion to future works.

A similar universality has been discussed for kernel methods in very recent works, but for the slightly different setting in which data is drawn from a Mixture of Gaussians (in which case there is no teacher, the label depends on which Gaussian has been chosen). The universality observed here for ridge regression with linear student, albeit different, is of a similar nature, and it would be interesting to discuss the link between these two approaches.