Generalisation error in learning with random features and the hidden manifold model^*

We study high-dimensional regression and classification for a synthetic dataset $\mathcal{D}={\left\{({\boldsymbol{x}}^{\mu },{y}^{\mu })\right\}}_{\mu =1}^{n}$ where each sample μ is created in the following three steps: (i) first, for each sample μ we create a vector ${\boldsymbol{c}}^{\mu }\in {\mathbb{R}}^{d}$ as

$$\begin{equation}{\boldsymbol{c}}^{\mu }\sim \mathcal{N}(0,{\mathrm{I}}_{d}).\end{equation} \tag{ 1.1 }$$

(ii) We then draw ${\boldsymbol{\theta }}^{0}\in {\mathbb{R}}^{d}$ from a separable distribution P_θ and draw independent labels ${\left\{{y}^{\mu }\right\}}_{\mu =1}^{n}$ from a (possibly probabilistic) rule f⁰:

$$\begin{equation}{y}^{\mu }={f}^{0}\left(\frac{1}{\sqrt{d}}{\boldsymbol{c}}^{\mu }\cdot {\boldsymbol{\theta }}^{0}\right)\in \mathbb{R}.\end{equation} \tag{ 1.2 }$$

(iii) The input data points ${\boldsymbol{x}}^{\mu }\in {\mathbb{R}}^{p}$ are created by a one-layer generative network with fixed and normalised weights $\mathrm{F}\in {\mathbb{R}}^{d\times p}$ and an activation function $\sigma :\mathbb{R}\to \mathbb{R}$, acting component-wise:

$$\begin{equation}{\boldsymbol{x}}^{\mu }=\sigma \left(\frac{1}{\sqrt{d}}{\mathrm{F}}^{\top }{\boldsymbol{c}}^{\mu }\right).\end{equation} \tag{ 1.3 }$$

We study the problem of supervised learning for the dataset $\mathcal{D}$ aiming at achieving a low generalisation error _g on a new sample x ^new, y^new drawn by the same rule as above, where:

$$\begin{equation}{{\epsilon}}_{g}=\frac{1}{{4}^{k}}{\mathbb{E}}_{{\boldsymbol{x}}^{\text{new}},{y}^{\text{new}}}\left[{\left({\hat{y}}_{\boldsymbol{w}}({\boldsymbol{x}}^{\text{new}})-{y}^{\mathrm{n}\mathrm{e}\mathrm{w}}\right)}^{2}\right]\end{equation} \tag{ 1.4 }$$

with k = 0 for regression task and k = 1 for classification task. Here, ${\hat{y}}_{\boldsymbol{w}}$ is the prediction on the new label y^new of the form:

$$\begin{equation}{\hat{y}}_{\boldsymbol{w}}(\boldsymbol{x})=\hat{f}\left(\boldsymbol{x}\cdot \hat{\boldsymbol{w}}\right).\end{equation} \tag{ 1.5 }$$

The weights $\hat{\boldsymbol{w}}\in {\mathbb{R}}^{p}$ are learned by minimising a loss function with a ridge regularisation term (for λ ⩾ 0) and defined as

$$\begin{equation}\hat{\boldsymbol{w}}=\underset{\boldsymbol{w}}{\mathrm{a}\mathrm{r}\mathrm{g}\;\;\mathrm{m}\mathrm{i}\mathrm{n}}\left[\sum\limits _{\mu =1}^{n}\ell ({y}^{\mu },{\boldsymbol{x}}^{\mu }\cdot \boldsymbol{w})+\frac{\lambda }{2}{\Vert}\boldsymbol{w}{{\Vert}}_{2}^{2}\right],\end{equation} \tag{ 1.6 }$$

where ℓ(⋅, ⋅) can be, for instance, a logistic, hinge, or square loss. Note that although our formula is valid for any f⁰ and $\hat{f}$, we take ${f}^{0}=\hat{f}=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}$, for the classification tasks and ${f}^{0}=\hat{f}=\text{id}$ for the regression tasks studied here. We now describe in more detail the above-discussed reasons why this model is of interest for machine learning.

Hidden manifold model: the dataset $\mathcal{D}$ can be seen as generated from the HMM introduced in [14]. From this perspective, although x ^μ lives in a p dimensional space, it is parametrised by a latent d-dimensional subspace spanned by the rows of the matrix F which are 'hidden' by the application of a scalar non-linear function σ. The labels y^μ are drawn from a generalised linear rule defined on the latent d-dimensional subspace via equation (1.2). In modern machine learning parlance, this can be seen as data generated by a one-layer generative neural network, such as those trained by generative adversarial networks or variational auto-encoders with random Gaussian inputs c^μ and a rather generic weight matrix F [19, 26–28].

Random features: the model considered in this paper is also an instance of the random features learning discussed in [20] as a way to speed up kernel-ridge-regression. From this perspective, the c ^μ s $\in {\mathbb{R}}^{d}$ are regarded as a set of d-dimensional i.i.d. Gaussian data points, which are projected by a feature matrix $\mathrm{F}={({\boldsymbol{f}}_{\rho })}_{\rho =1}^{p}\in {\mathbb{R}}^{d\times p}$ into a higher dimensional space, followed by a non-linearity σ. In the p → ∞ limit of infinite number of features, performing regression on $\mathcal{D}$ is equivalent to kernel regression on the c^μ s with a deterministic kernel $K({\mathbf{c}}^{{\mu }_{1}},{\mathbf{c}}^{{\mu }_{2}})={\mathbb{E}}_{\mathbf{f}}\left[\sigma (\mathbf{f}\cdot {\mathbf{c}}^{{\mu }_{1}}/\sqrt{d})\cdot \sigma (\mathbf{f}\cdot {\mathbf{c}}^{{\mu }_{2}}/\sqrt{d})\right]$ where $\mathbf{f}\in {\mathbb{R}}^{d}$ is sampled in the same way as the rows of F. Random features are also intimately linked with the lazy training regime, where the weights of a neural network stay close to their initial value during training. The training is lazy as opposed to a 'rich' one where the weights change enough to learn useful features. In this regime, neural networks become equivalent to a random feature model with correlated features [16–18, 29–31].

The main contribution of this work is a closed-form expression for the generalisation error _g , equation (2.1), that is valid in the high-dimensional limit where the number of samples n, and the two dimensions p and d are large, but their respective ratios are of order one, and for generic sequence of matrices F satisfying the following balance conditions:

$$\begin{equation}\frac{1}{\sqrt{p}}\sum\limits _{i=1}^{p}{w}_{i}^{{a}_{1}}{w}_{i}^{{a}_{2}}\dots {w}_{i}^{{a}_{s}}{\mathrm{F}}_{i{\rho }_{1}}{\mathrm{F}}_{i{\rho }_{2}}\dots {\mathrm{F}}_{i{\rho }_{q}}=O(1),\end{equation} \tag{ 1.7 }$$

where ${\left\{{\boldsymbol{w}}^{a}\right\}}_{a=1}^{r}$ are r independent samples from the Gibbs measure (2.7), and ρ₁, ρ₂, ..., ρ_q ∈ {1, ..., d}, a₁, a₂, ..., a_s ∈ {1, ..., r} are an arbitrary choice of subset of indices, with $s,q\in {\mathbb{Z}}_{+}$. The non-linearities ${f}^{0},\hat{f},\sigma$ and the loss function ℓ can be arbitrary. Our result for the generalisation error stems from the replica method and we conjecture it to be exact for convex loss functions ℓ. It can also be useful for non-convex loss functions but in those cases it is possible that the so-called replica symmetry breaking [21] needs to be taken into account to obtain an exact expression. In the present paper we hence focus on convex loss functions ℓ and leave the more general case for future work. The final formulas are simpler for nonlinearities σ that give zero when integrated over a centred Gaussian variable, and we hence focus on those cases.

An interesting application of our setting is ridge regression, i.e. taking $\hat{f}(x) = x$ with square loss, and random i.i.d. Gaussian feature matrices. For this particular case [13] proved an equivalent expression. Indeed, in this case there is an explicit solution of equation (1.6) that can be rigorously studied with random matrix theory. In a subsequent work [32] derived heuristically a formula for the special case of random i.i.d. Gaussian feature matrices for the maximum margin classification, corresponding to the hinge loss function in our setting, with the difference, however, that the labels y^μ are generated from the x^μ instead of the variable c^μ as in our case.

Our main technical contribution is thus to provide a generic formula for the model described in section 1.1 for any loss function and for fairly generic features F, including for instance deterministic ones.

The authors of [14] analysed the learning dynamics of a neural network containing several hidden units using a one-pass SGD for exactly the same model of data as here. In this online setting, the algorithm is never exposed to a sample twice, greatly simplifying the analysis as what has been learned at a given epoch can be considered independent of the randomness of a new sample. Another motivation of the present work is thus to study the sample complexity for this model (in our case only a bounded number of samples is available, and the one-pass SGD would be highly suboptimal).

An additional technical contribution of our work is to derive an extension of the equivalence between the considered data model and a model with Gaussian covariate, that has been observed and conjectured to hold rather generically in both [14, 32]. While we do not provide a rigorous proof for this equivalence, we show that it arises naturally using the replica method, giving further evidence for its validity.

Finally, the analysis of our formula for particular machine learning tasks of interest allows for an analytical investigation of a rich phenomenology that is also observed empirically in real-life scenarios. In particular

• The double descent behaviour, as termed in [33] and exemplified in [34], is exhibited for the non-regularized logistic regression loss. The peak of worst generalisation does not corresponds to p = n as for the square loss [13], but rather corresponds to the threshold of linear separability of the dataset. We also characterise the location of this threshold, generalising the results of [11] to our model.
• When using projections to approximate kernels, it has been observed that orthogonal features F perform better than random i.i.d. [35]. We show that this behaviour arises from our analytical formula, illustrating the 'unreasonable effectiveness of structured random orthogonal embeddings' [35].
• We compute the phase diagram for the generalisation error for the HMM and discuss the dependence on the various parameters, in particular the ratio between the ambient and latent dimensions.

Let F be a feature matrix satisfying the balance condition stated in equation (1.7). Then, in the high-dimensional limit where p, d, n → ∞ with α = n/p, γ = d/p fixed, the generalisation error, equation (1.4), of the model introduced in section (1.4) for σ such that its integral over a centered Gaussian variable is zero (so that κ₀ = 0 in equation (2.10)) is given by the following easy-to-evaluate integral:

$$\begin{equation}\underset{n\to \infty }{\mathrm{lim}}\enspace {{\epsilon}}_{g}={\mathbb{E}}_{\lambda ,\nu }\left[{({f}^{0}(\nu )-\hat{f}(\lambda ))}^{2}\right],\end{equation} \tag{ 2.1 }$$

where f⁰(.) is defined in (1.2), $\hat{f}(.)$ in (1.5) and (ν, λ) are jointly Gaussian random variables with zero mean and covariance matrix:

$$\begin{equation}{\Sigma}=\left(\begin{matrix} \rho & {M}^{\star } \\ {M}^{\star } & {Q}^{\star } \end{matrix}\right)\in {\mathbb{R}}^{2}\end{equation} \tag{ 2.2 }$$

with ${M}^{\star }={\kappa }_{1}{m}_{s}^{\star }$, ${Q}^{\star }={\kappa }_{1}^{2}{q}_{s}^{\star }+{\kappa }_{\star }^{2}{q}_{w}^{\star }$. The constants κ_⋆, κ₁ depend on the nonlinearity σ via equation (2.10), and ${q}_{s}^{\star },{q}_{w}^{\star },{m}_{s}^{\star }$, defined as:

$$\begin{equation}\rho =\frac{1}{d}{\Vert}{\boldsymbol{\theta }}^{0}{{\Vert}}^{2}\quad {q}_{s}^{\star }=\frac{1}{d}\mathbb{E}{\Vert}\mathrm{F}\hat{\boldsymbol{w}}{{\Vert}}^{2}\quad {q}_{w}^{\star }=\frac{1}{p}\mathbb{E}{\Vert}\hat{\boldsymbol{w}}{{\Vert}}^{2}\quad {m}_{s}^{\star }=\frac{1}{d}\mathbb{E}\left[(\mathrm{F}\hat{\boldsymbol{w}})\cdot {\boldsymbol{\theta }}^{0}\right].\end{equation} \tag{ 2.3 }$$

The values of these parameters correspond to the solution of the optimisation problem in equation (1.6), and can be obtained as the fixed point solutions of the following set of self-consistent saddle-point equations:

$$\begin{align} \begin{cases}{\hat{V}}_{s}=\frac{\alpha {\kappa }_{1}^{2}}{\gamma V}{\mathbb{E}}_{\xi }\left[{\int }_{\mathbb{R}}\mathrm{d}y\enspace \mathcal{Z}\left(y,{\omega }_{0}\right)\left(1-{\partial }_{\omega }\eta \left(y,{\omega }_{1}\right)\right)\right],\quad \\ {\hat{q}}_{s}=\frac{\alpha {\kappa }_{1}^{2}}{\gamma {V}^{2}}{\mathbb{E}}_{\xi }\left[{\int }_{\mathbb{R}}\mathrm{d}y\enspace \mathcal{Z}\left(y,{\omega }_{0}\right){\left(\eta \left(y,{\omega }_{1}\right)-{\omega }_{1}\right)}^{2}\right],\quad \\ {\hat{m}}_{s}=\frac{\alpha {\kappa }_{1}}{\gamma V}{\mathbb{E}}_{\xi }\left[{\int }_{\mathbb{R}}\mathrm{d}y\enspace {\partial }_{\omega }\mathcal{Z}\left(y,{\omega }_{0}\right)(\eta \left(y,{\omega }_{1}\right)-{\omega }_{1})\right],\quad \\ {\hat{V}}_{w}=\frac{\alpha {\kappa }_{\star }^{2}}{V}{\mathbb{E}}_{\xi }\left[{\int }_{\mathbb{R}}\mathrm{d}y\enspace \mathcal{Z}\left(y,{\omega }_{0}\right)\left(1-{\partial }_{\omega }\eta \left(y,{\omega }_{1}\right)\right)\right],\quad \\ {\hat{q}}_{w}=\frac{\alpha {\kappa }_{\star }^{2}}{{V}^{2}}{\mathbb{E}}_{\xi }\left[{\int }_{\mathbb{R}}\mathrm{d}y\enspace \mathcal{Z}\left(y,{\omega }_{0}\right){\left(\eta \left(y,{\omega }_{1}\right)-{\omega }_{1}\right)}^{2}\right],\quad \end{cases}& \begin{cases}{V}_{s}=\frac{1}{{\hat{V}}_{s}}\left(1-z{g}_{\mu }(-z)\right),\quad \\ {q}_{s}=\frac{{\hat{m}}_{s}^{2}+{\hat{q}}_{s}}{{\hat{V}}_{s}^{2}}\left[1-2z{g}_{\mu }(-z)+{z}^{2}{g}_{\mu }^{\prime }(-z)\right]\quad \\ \qquad -\frac{{\hat{q}}_{w}}{(\lambda +{\hat{V}}_{w}){\hat{V}}_{s}}\left[-z{g}_{\mu }(-z)+{z}^{2}{g}_{\mu }^{\prime }(-z)\right],\quad \\ {m}_{s}=\frac{{\hat{m}}_{s}}{{\hat{V}}_{s}}\left(1-z\enspace {g}_{\mu }(-z)\right),\quad \\ {V}_{w}=\frac{\gamma }{\lambda +{\hat{V}}_{w}}\left[\frac{1}{\gamma }-1+z{g}_{\mu }(-z)\right],\quad \\ {q}_{w}=\gamma \frac{{\hat{q}}_{w}}{{(\lambda +{\hat{V}}_{w})}^{2}}\left[\frac{1}{\gamma }-1+{z}^{2}{g}_{\mu }^{\prime }(-z)\right],\quad \\ \qquad -\gamma \frac{{\hat{m}}_{s}^{2}+{\hat{q}}_{s}}{(\lambda +{\hat{V}}_{w}){\hat{V}}_{s}}\left[-z{g}_{\mu }(-z)+{z}^{2}{g}_{\mu }^{\prime }(-z)\right],\quad \end{cases} \end{align} \tag{ 2.4 }$$

written in terms of the following auxiliary variables $\xi \sim \mathcal{N}(0,1)$, $z=\frac{\lambda +{\hat{V}}_{w}}{{\hat{V}}_{s}}$ and functions:

$$\begin{align} \begin{aligned} \eta (y,\omega )& =\underset{x\in \mathbb{R}}{\mathrm{a}\mathrm{r}\mathrm{g}\;\;\mathrm{m}\mathrm{i}\mathrm{n}}\left[\frac{{(x-\omega )}^{2}}{2V}+\ell (y,x)\right], \\ \mathcal{Z}(y,\omega )& =\int \frac{\mathrm{d}x}{\sqrt{2\pi {V}^{0}}}\enspace {\mathrm{e}}^{-\frac{1}{2{V}^{0}}{\left(x-\omega \right)}^{2}}\delta \left(y-{f}^{0}(x)\right) \end{aligned}\end{align} \tag{ 2.5 }$$

where $V={\kappa }_{1}^{2}{V}_{s}+{\kappa }_{\star }^{2}{V}_{w}$, ${V}^{0}=\rho -\frac{{M}^{2}}{Q}$, $Q={\kappa }_{1}^{2}{q}_{s}+{\kappa }_{\star }^{2}{q}_{w}$, M = κ₁ m_s , ${\omega }_{0}=\left(M/\sqrt{Q}\right)\xi$ and ${\omega }_{1}=\sqrt{Q}\xi$. In the above, we assume that the matrix $\mathrm{F}{\mathrm{F}}^{\top }\in {\mathbb{R}}^{d\times d}$ associated to the feature map F has a well behaved spectral density, and denote g_μ its Stieltjes transform.

The training loss on the dataset $\mathcal{D}={\left\{{\boldsymbol{x}}^{\mu },{y}^{\mu }\right\}}_{\mu =1}^{n}$ can also be obtained from the solution of the above equations as

$$\begin{equation}\underset{n\to \infty }{\mathrm{lim}}\enspace {{\epsilon}}_{t}=\frac{\lambda }{2\alpha }{q}_{w}^{\star }+{\mathbb{E}}_{\xi ,y}\left[\mathcal{Z}\left(y,{\omega }_{0}^{\star }\right)\ell \left(y,\eta (y,{\omega }_{1}^{\star })\right)\right]\end{equation} \tag{ 2.6 }$$

where as before $\xi \sim \mathcal{N}(0,1)$, $y\sim \text{Uni}(\mathbb{R})$ and $\mathcal{Z},\eta$ are the same as in equation (2.5), evaluated at the solution of the above saddle-point equations ${\omega }_{0}^{\star }=\left({M}^{\star }/\sqrt{{Q}^{\star }}\right)\xi$, ${\omega }_{1}^{\star }=\sqrt{{Q}^{\star }}\xi$.

Sketch of derivation—we now sketch the derivation of the above result. A complete and detailed account can be found in section C of the supplementary material. The derivation is based on the key observation that in the high-dimensional limit the asymptotic generalisation error only depends on the solution $\hat{\boldsymbol{w}}\in {\mathbb{R}}^{p}$ of equation (1.5) through the scalar parameters $({q}_{s}^{\star },{q}_{w}^{\star },{m}_{s}^{\star })$ defined in equation (2.3). The idea is therefore to rewrite the high-dimensional optimisation problem in terms of only these scalar parameters.

The first step is to note that the solution of equation (1.6) can be written as the average of the following Gibbs measure

$$\begin{equation}{\pi }_{\beta }(\boldsymbol{w}\vert \left\{{\boldsymbol{x}}^{\mu },{y}^{\mu }\right\})=\frac{1}{{\mathcal{Z}}_{\beta }}\enspace {\text{e}}^{-\beta \left[\sum\limits _{\mu =1}^{n}\ell \left({y}^{\mu },{\boldsymbol{x}}^{\mu }\cdot \boldsymbol{w}\right)+\frac{\lambda }{2}{\Vert}\boldsymbol{w}{{\Vert}}_{2}^{2}\right]},\end{equation} \tag{ 2.7 }$$

in the limit β → ∞. Of course, we have not gained much, since an exact calculation of π_β is intractable for large values of n, p and d. This is where the replica method comes in. It states that the distribution of the free energy density $f=-\mathrm{l}\mathrm{o}\mathrm{g}\enspace {\mathcal{Z}}_{\beta }$ (when seen as a random variable over different realisations of dataset $\mathcal{D}$) associated with the measure μ_β concentrates, in the high-dimensional limit, around a value f_β that depends only on the averaged replicated partition function ${\mathcal{Z}}_{\beta }^{r}$ obtained by taking r > 0 copies of ${\mathcal{Z}}_{\beta }$:

$$\begin{equation}{f}_{\beta }=\underset{r\to {0}^{+}}{\mathrm{lim}}\enspace \frac{\mathrm{d}}{\mathrm{d}r}\enspace \underset{p\to \infty }{\mathrm{lim}}\left[-\frac{1}{p}\left({\mathbb{E}}_{\left\{{\boldsymbol{x}}^{\mu },{y}^{\mu }\right\}}{\mathcal{Z}}_{\beta }^{r}\right)\right].\end{equation} \tag{ 2.8 }$$

Interestingly, ${\mathbb{E}}_{\left\{{\boldsymbol{x}}^{\mu },{y}^{\mu }\right\}}{\mathcal{Z}}_{\beta }^{r}$ can be computed explicitly for $r\in \mathbb{N}$, and the limit r → 0⁺ is taken by analytically continuing to r > 0 (see section C of the supplementary material). The upshot is that ${\mathcal{Z}}^{r}$ can be written as

$$\begin{equation}{\mathbb{E}}_{\left\{{\boldsymbol{x}}^{\mu },{y}^{\mu }\right\}}{\mathcal{Z}}_{\beta }^{r}\propto \int \mathrm{d}{q}_{s}\enspace \mathrm{d}{q}_{w}\enspace \mathrm{d}{m}_{s}\enspace {\mathrm{e}}^{p{{\Phi}}_{\beta }^{(r)}({m}_{s},{q}_{s},{q}_{w})}\end{equation} \tag{ 2.9 }$$

where Φ_β —known as the replica symmetric potential—is a concave function depending only on the following scalar parameters:

$$\begin{equation*}{q}_{s}=\frac{1}{d}{\Vert}\mathrm{F}\boldsymbol{w}{{\Vert}}^{2},\qquad {q}_{w}=\frac{1}{p}{\Vert}\boldsymbol{w}{{\Vert}}^{2},\qquad {m}_{s}=\frac{1}{d}(\mathrm{F}\boldsymbol{w})\cdot {\boldsymbol{\theta }}^{0}\end{equation*}$$

for w ∼ π_β . In the limit of p → ∞, this integral concentrates around the extremum of the potential ${{\Phi}}_{\beta }^{(0)}$ for any β. Since the optimisation problem in equation (1.5) is convex, by construction as β → ∞ the overlap parameters $({q}_{s}^{\star },{q}_{w}^{\star },{m}_{s}^{\star })$ satisfying this optimisation problem are precisely the ones of equation (2.3) corresponding to the solution $\hat{\boldsymbol{w}}\in {\mathbb{R}}^{p}$ of equation (1.5).

In summary, the replica method allows to circumvent the hard-to-solve high-dimensional optimisation problem equation (1.5) by directly computing the generalisation error in equation (1.4) in terms of a simpler scalar optimisation. Doing gradient descent (GD) in ${{\Phi}}_{\beta }^{(0)}$ and taking β → ∞ lead to the saddle-point equation (2.4).

The backbone of the replica derivation sketched above and detailed in section C of the supplementary material is a central limit theorem type result coined as the Gaussian equivalence theorem (GET) from [14] used in the context of the 'replicated' Gibbs measure obtained by taking r copies of (2.7). In this approach, we need to assume that the 'balance condition' (1.7) applies with probability one when the weights w are sampled from the replicated measure. We shall use this assumption in the following, checking its self-consistency via agreement with simulations.

It is interesting to observe that, when applying the GET in the context of our replica calculation, the resulting asymptotic generalisation error stated in section 2.1 is equivalent to the asymptotic generalisation error of the following linear model:

$$\begin{equation}{\boldsymbol{x}}^{\mu }={\kappa }_{0}\mathbf{1}+{\kappa }_{1}\frac{1}{\sqrt{d}}{\mathrm{F}}^{\top }{\boldsymbol{c}}^{\mu }+{\kappa }_{\star }\enspace {\boldsymbol{z}}^{\mu },\end{equation} \tag{ 2.10 }$$

with ${\kappa }_{0}=\mathbb{E}\left[\sigma (z)\right]$, ${\kappa }_{1}\equiv \mathbb{E}\left[z\sigma (z)\right]$, ${\kappa }_{\star }^{2}\equiv \mathbb{E}\left[\sigma {(z)}^{2}\right]-{\kappa }_{0}^{2}-{\kappa }_{1}^{2}$, and ${\boldsymbol{z}}^{\mu }\sim \mathcal{N}(\mathbf{0},{\mathrm{I}}_{p})$. We have for instance, $({\kappa }_{0},{\kappa }_{1},{\kappa }_{\star })\approx \left(0,\frac{2}{\sqrt{3\pi }},0.2003\right)$ for σ = erf and $({\kappa }_{0},{\kappa }_{1},{\kappa }_{\star })=\left(0,\sqrt{\frac{2}{\pi }},\sqrt{1-\frac{2}{\pi }}\right)$ for σ = sign, two cases explored in the next section. This equivalence constitutes a result with an interest in its own, with applicability beyond the scope of the generalised linear task equation (1.6) studied here.

Equation (2.10) is precisely the mapping obtained by [13], who proved its validity rigorously in the particular case of the square loss and Gaussian random matrix F using random matrix theory. The same equivalence arises in the analysis of kernel random matrices [39, 40] and in the study of online learning [14]. The replica method thus suggests that the equivalence actually holds in a much larger class of learning problem, as conjectured as well in [32], and numerically confirmed in all our numerical tests. It also potentially allows generalisation of the analysis in this paper for data coming from a learned generative adversarial network, along the lines of [28, 41].

Figure 1 illustrates the remarkable agreement between the result of the generalisation formula, equation (2.1) and simulations both on the data equation (1.3) with σ(x) = sign(x) non-linearity, and on the Gaussian equivalent data equation (2.10) where the non-linearity is replaced by rescaling by a constant plus noise. The agreement is flawless as implied by the theory in the high-dimensional limit, testifying that the used system size d = 200 is sufficiently large for the asymptotic theory to be relevant. We observed similar good agreement between the theory and simulation in all the cases we tested, in particular in all those presented in the following.

Zoom In Zoom Out Reset image size

Among the surprising observations in modern machine learning is the fact that one can use learning methods that achieve zero training error, yet their generalisation error does not deteriorate as more and more parameters are added into the neural network. The study of such 'interpolators' have attracted a growing attention over the last few years [9, 12, 13, 33, 34, 42–44], as it violates basic intuition on the bias-variance trade-off [45]. Indeed classical learning theory suggests that generalisation should first improve then worsen when increasing model complexity, following a U-shape curve. Many methods, including neural networks, instead follow a so-called 'double descent curve' [33] that displays two regimes: the 'classical' U-curve found at low number of parameters is followed at high number of parameters by an interpolation regime where the generalisation error decreases monotonically. Consequently neural networks do not drastically overfit even when using much more parameters than data samples [46], as actually observed already in the classical work [45]. Between the two regimes, a 'peak' occurs at the interpolation threshold [9, 22, 34, 47]. It should, however, be noted that existence of this 'interpolation' peak is an independent phenomenon from the lack of overfitting in highly over-parametrized networks, and indeed in a number of the related works these two phenomena were observed separately [9, 22, 45, 47]. Scaling properties of the peak and its relation to the jamming phenomena in physics are in particular studied in [43].

Among the simple models that allow to observe this behaviour, random projections—that are related to lazy training and kernel methods—are arguably the most natural one. The double descent has been analysed in detail in the present model in the specific case of a square loss on a regression task with random Gaussian features [13]. Our analysis allows to show the generality and the robustness of the phenomenon to other tasks, matrices and losses. In figure 2 we compare the double descent as present in the square loss (blue line) with the one of logistic loss (red line) for random Gaussian features. We plot the value of the generalisation error at small values of the regularisation λ (full line), and for optimal value of λ (dashed line) for a fixed ratio between the number of samples and the dimension n/d as a function of the number of random features per sample p/n. We also plot the value of the training error (lower panel) for a small regularisation value, showing that the peaks indeed occur when the training loss goes to zero. For the square loss the peak appears at 1/α = p/n = 1 when the system of n linear equations with p parameters becomes solvable. For the logistic loss the peak instead appears at a value 1/α* where the data $\mathcal{D}$ become linearly separable and hence the logistic loss can be optimised down to zero. These values 1/α* depends on the value n/d, and this dependence is plotted in figure 5. For very large dimension d, i.e. n/d → 0 the data matrix X is close to i.i.d. random matrix and hence the α*(n/d = 0) = 2 as famously derived in classical work by Cover [48]. For n/d > 0 the α* is growing (1/α* decreasing) as correlations make data easier to linearly separate, similarly as in [11].

Zoom In Zoom Out Reset image size

Figure 2 also shows that better error can be achieved with the logistic loss compared to the square loss, both for small and optimal regularisations, except in a small region around the logistic interpolation peak. In the Kernel limit, i.e. p/n → ∞, the generalization error at optimal regularisation saturates at _g (p/n → ∞) ≃ 0.17 for square loss and at _g (p/n → ∞) ≃ 0.16 for logistic loss. Figure 3 then depicts a 3D plot of the generalisation error also illustrating the position of the interpolation peak.

Zoom In Zoom Out Reset image size

Kernel methods are a very popular class of machine learning techniques, achieving state-of-the-art performance on a variety of tasks with theoretical guarantees [49–51]. In the context of neural network, they are the subject of a renewal of interest in the context of the neural tangent Kernel [17]. Applying kernel methods to large-scale 'big data' problems, however, poses many computational challenges, and this has motivated a variety of contributions to develop them at scale, see, e.g. [50, 52–54]. Random features [20] are among the most popular techniques to do so.

Here, we want to compare the performance of random projection with respect to structured ones, and in particular orthogonal random projections [35] or deterministic matrices such as real Fourier DCT and Hadamard matrices used in fast projection methods [55–57]. Up to normalisation, these matrices have the same spectral density. Since the asymptotic generalisation error only depends on the spectrum of FF^⊤, all these matrices share the same theoretical prediction when properly normalised, see figure 4. In our computation, left- and right-orthogonal invariance is parametrised by letting F = U^⊤DV for $\mathrm{U}\in {\mathbb{R}}^{d\times d}$, $\mathrm{V}\in {\mathbb{R}}^{p\times p}$ two orthogonal matrices drawn from the Haar measure, and $\mathrm{D}\in {\mathbb{R}}^{d\times p}$ a diagonal matrix of rank min(d, p). In order to compare the results with the Gaussian case, we fix the diagonal entries ${d}_{k}=\mathrm{max}(\sqrt{\gamma },1)$ of D such that an arbitrary projected vector has the same norm, on average, to the Gaussian case.

Zoom In Zoom Out Reset image size

Figure 4 shows that random orthogonal embeddings always outperform Gaussian random projections, in line with empirical observations, and that they allow to reach the kernel limit with fewer number of projections. Their behaviour is, however, qualitatively similar to the one of random i.i.d. projections. We also show in figure 5 that orthogonal projections allow to separate the data more easily than the Gaussian ones, as the phase transition curve delimiting the linear separability of the logistic loss get shifted to the left.

Zoom In Zoom Out Reset image size

In this subsection we consider the HMM where p-dimensional x data lie on a d-dimensional manifold, we have mainly in mind d < p. The labels y are generated using the coordinates on the manifold, equation (1.2).

In figure 6 we plot the generalisation error of classification with the square loss for various values of the regularisation λ. We fix the ratio between the dimension of the sub-manifold and the dimensionality of the input data to d/p = 0.1 and plot the learning curve, i.e. the error as a function of the number of samples per dimension. Depending on the value of the regularisation, we observe that the interpolation peak, which is at α = 1 at very small regularisation (here the over-parametrised regime is on the left-hand side), decreases for larger regularisation λ. A similar behaviour has been observed for other models in the past, see e.g. [47]. Finally, figure 6 depicts the error for optimised regularisation parameter in the black dashed line. For large number of samples we observe the generalisation error at optimal regularisation to saturate in this case at _g (α → ∞) → 0.0325. A challenge for future work is to see whether better performance can be achieved on this model by including hidden variables into the neural network.

Zoom In Zoom Out Reset image size

Figure 7 then shows the generalisation error for the optimised regularisation λ with square loss as a function of the ratio between the latent and the data dimensions d/p. In the limit d/p ≫ 1 the data matrix becomes close to a random i.i.d. matrix and the labels are effectively random, thus only bad generalisation can be reached. Interestingly, as d/p decreases to small values even the simple classification with regularised square loss is able to 'disentangle' the hidden manifold structure in the data and to reach a rather low generalisation error. The figure quantifies how the error deteriorates when the ratio between the two dimensions d/p increases. Rather remarkably, for a low d/p a good generalisation error is achieved even in the over-parametrised regime, where the dimension is larger than the number of samples, p > n. In a sense, the square loss linear classification is able to locate the low-dimensional subspace and find good generalisation even in the over-parametrised regime as long as roughly d ≲ n. The observed results are in qualitative agreement with the results of learning with SGD in [14] where for very low d/p good generalisation error was observed in the HMM, but a rather bad one for d/p = 0.5.

Zoom In Zoom Out Reset image size

In this work we study a series of regression and classification tasks for a dataset ${\left\{{\boldsymbol{x}}^{\mu },{y}^{\mu }\right\}}_{\mu =1}^{n}$ with labels ${y}^{\mu }\in \mathbb{R}$ sampled identically from a generalised linear model:

$$\begin{equation}{y}^{\mu }\sim {P}_{y}^{0}\left({y}^{\mu }\left\vert \frac{{\boldsymbol{c}}^{\mu }\cdot {\boldsymbol{\theta }}^{0}}{\sqrt{d}}\right.\right),\end{equation} \tag{ A.1 }$$

where the output-channel ${P}_{y}^{0}\left(\cdot \right)$ is defined as:

$$\begin{equation}{P}_{y}^{0}\left({y}^{\mu }\left\vert \frac{{\boldsymbol{c}}^{\mu }\cdot {\boldsymbol{\theta }}^{0}}{\sqrt{d}}\right.\right)=\int \mathrm{d}{\xi }^{\mu }P\left({\xi }^{\mu }\right)\delta \left({y}^{\mu }-{f}^{0}\left(\frac{{\boldsymbol{c}}^{\mu }\cdot {\boldsymbol{\theta }}^{0}}{\sqrt{d}};{\xi }^{\mu }\right)\right)\end{equation} \tag{ A.2 }$$

for some noise ξ^μ and for data points ${\boldsymbol{x}}^{\mu }\in {\mathbb{R}}^{p}$ given by:

$$\begin{equation}{\boldsymbol{x}}^{\mu }=\sigma \left(\frac{1}{\sqrt{d}}\sum\limits _{\rho =1}^{d}{c}_{\rho }^{\mu }{\boldsymbol{f}}_{\rho }\right).\end{equation} \tag{ A.3 }$$

The vectors ${\boldsymbol{c}}^{\mu }\in {\mathbb{R}}^{d}$ is assumed to be identically drawn from $\mathcal{N}(0,{\mathrm{I}}_{d})$, and ${\boldsymbol{\theta }}^{0}\in {\mathbb{R}}^{d}$ from a separable distribution P_θ . The family of vectors ${\boldsymbol{f}}_{\rho }\in {\mathbb{R}}^{p}$ and the scalar function $\sigma :\mathbb{R}\to \mathbb{R}$ can be arbitrary.

Although our results are valid for the general model introduced above, the two examples we will be exploring in this work are the noisy linear channel (for regression tasks) and the deterministic sign channel (for classification tasks):

$$\begin{equation}{y}^{\mu }=\frac{{\boldsymbol{c}}^{\mu }\cdot {\boldsymbol{\theta }}^{0}}{\sqrt{d}}+\sqrt{{\Delta}}{\xi }^{\mu }\quad {\Leftrightarrow}\quad {P}_{y}^{0}\left(\boldsymbol{y}\left\vert \frac{{\boldsymbol{c}}^{\mu }\cdot {\boldsymbol{\theta }}^{0}}{\sqrt{d}}\right.\right)=\prod\limits _{\mu =1}^{n}\mathcal{N}\left({y}^{\mu }\hspace{-1.5pt};\frac{{\boldsymbol{c}}^{\mu }\cdot {\boldsymbol{\theta }}^{0}}{\sqrt{d}},{\Delta}\right)\end{equation} \tag{ A.4 }$$

$$\begin{equation}{y}^{\mu }=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(\frac{{\boldsymbol{c}}^{\mu }\cdot {\boldsymbol{\theta }}^{0}}{\sqrt{d}}\right)\quad {\Leftrightarrow}\quad {P}_{y}^{0}\left(\boldsymbol{y}\left\vert \frac{{\boldsymbol{c}}^{\mu }\cdot {\boldsymbol{\theta }}^{0}}{\sqrt{d}}\right.\right)=\prod\limits _{\mu =1}^{n}\delta \left({y}^{\mu }-\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(\frac{{\boldsymbol{c}}^{\mu }\cdot {\boldsymbol{\theta }}^{0}}{\sqrt{d}}\right)\right)\end{equation} \tag{ A.5 }$$

where ${\xi }^{\mu }\sim \mathcal{N}(0,1)$ and Δ > 0.

This dataset can be regarded from two different perspectives.

Hidden manifold model: the dataset ${\left\{{\boldsymbol{x}}^{\mu },{y}^{\mu }\right\}}_{\mu =1,\dots ,n}$ is precisely the HMM introduced in [14] to study the dynamics of online learning in a synthetic but structured dataset. From this perspective, although x ^μ lives in a p dimensional space, it is parametrised by a latent d < p-dimensional subspace spanned by the basis ${\left\{{\boldsymbol{f}}_{\rho }\right\}}_{\rho =1,\dots ,d}$ which is 'hidden' by the application of a scalar nonlinear function σ acting component-wise. The labels y ^μ are then drawn from a generalised linear rule defined on the latent d-dimensional space.

Random features model: the dataset ${\left\{{\boldsymbol{x}}^{\mu },{y}^{\mu }\right\}}_{\mu =1,\dots ,n}$ is tightly related to the random features model studied in [20] as a random approximation for kernel ridge regression. In this perspective, ${\boldsymbol{c}}^{\mu }\in {\mathbb{R}}^{d}$ is regarded as a collection of d-dimensional data points which are projected by a random feature matrix $\mathrm{F}={({\boldsymbol{f}}_{\rho })}_{\rho =1}^{p}\in {\mathbb{R}}^{d\times p}$ into a higher dimensional space, followed by a non-linearity σ. In the limit of infinite number of features d, p → ∞ with fixed ratio d/p, performing ridge regression of x ^μ is equivalent to kernel ridge regression with a limiting kernel depending on the distribution of the feature matrix F and on the non-linearity σ.

In this work, we study the problem of learning the rule from equation (A.1) from the dataset ${\left\{({\boldsymbol{x}}^{\mu },{y}^{\mu })\right\}}_{\mu =1,\dots ,n}$ introduced above with a generalised linear model:

$$\begin{equation}{\hat{y}}^{\mu }=\hat{f}\left({\boldsymbol{x}}^{\mu }\cdot \hat{\boldsymbol{w}}\right)\end{equation} \tag{ A.6 }$$

where the weights $\boldsymbol{w}\in {\mathbb{R}}^{p}$ are learned by minimising a loss function with a ridge regularisation term:

$$\begin{equation}\hat{\boldsymbol{w}}=\underset{\boldsymbol{w}}{\mathrm{min}}\left[\sum\limits _{\mu =1}^{n}\ell ({y}^{\mu },{\boldsymbol{x}}^{\mu }\cdot \boldsymbol{w})+\frac{\lambda }{2}{\Vert}\boldsymbol{w}{{\Vert}}_{2}^{2}\right].\end{equation} \tag{ A.7 }$$

for λ > 0.

It is worth stressing that our results hold for general ℓ, $\hat{f}$ and f⁰—including non-convex loss functions. However, for the purpose of the applications explored in this manuscript, we will be mostly interested in the cases $\hat{f}(x)={f}^{0}(x)=x$ for regression and $\hat{f}(x)={f}^{0}(x)=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(x)$ for classification, and we will focus on the following two loss functions:

$$\begin{equation}\ell ({y}^{\mu },{\boldsymbol{x}}^{\mu }\cdot \boldsymbol{w})=\begin{cases}\frac{1}{2}{({y}^{\mu }-{\boldsymbol{x}}^{\mu }\cdot \boldsymbol{w})}^{2},\quad & \;\text{square}\;\text{loss} \\ \mathrm{log}\left(1+{\mathrm{e}}^{-{y}^{\mu }\left({\boldsymbol{x}}^{\mu }\cdot \boldsymbol{w}\right)}\right),\quad & \;\text{logistic}\;\text{loss}. \end{cases}\end{equation} \tag{ A.8 }$$

Note that these loss functions are strictly convex. Therefore, for these losses, the regularised optimisation problem in (A.7) has a unique solution.

Given a new pair ( x ^new, y^new) drawn independently from the same distribution as ${\left\{({\boldsymbol{x}}^{\mu },{y}^{\mu })\right\}}_{\mu =1}^{n}$, we define the success of our fit through the generalisation error, defined as:

$$\begin{equation}{{\epsilon}}_{g}=\frac{1}{{4}^{k}}{\mathbb{E}}_{{\boldsymbol{x}}^{\mathrm{n}\mathrm{e}\mathrm{w}},{y}^{\mathrm{n}\mathrm{e}\mathrm{w}}}{\left({y}^{\mathrm{n}\mathrm{e}\mathrm{w}}-{\hat{y}}^{\mathrm{n}\mathrm{e}\mathrm{w}}\right)}^{2}\end{equation} \tag{ A.9 }$$

where ${\hat{y}}^{\mathrm{n}\mathrm{e}\mathrm{w}}=\hat{f}({\boldsymbol{x}}^{\mathrm{n}\mathrm{e}\mathrm{w}}\cdot \hat{\boldsymbol{w}})$, and for convenience we choose k = 0 for the regression tasks and k = 1 for the classification task, such that the generalisation error in this case counts misclassification. Note that for a classification problem, the generalisation error is just one minus the classification error.

Similarly, we define the training loss on the dataset ${\left\{{\boldsymbol{x}}^{\mu },{y}^{\mu }\right\}}_{\mu =1}^{n}$ as:

$$\begin{equation}{{\epsilon}}_{t}=\frac{1}{n}{\mathbb{E}}_{\left\{{\boldsymbol{x}}^{\mu },{y}^{\mu }\right\}}\left[\sum\limits _{\mu =1}^{n}\ell \left({y}^{\mu },{\boldsymbol{x}}^{\mu }\cdot \hat{\boldsymbol{w}}\right)+\frac{\lambda }{2}{\Vert}\hat{\boldsymbol{w}}{{\Vert}}_{2}^{2}\right].\end{equation} \tag{ A.10 }$$

Finally, all the results of this manuscript are derived in the high-dimensional limit, also known as thermodynamic limit in the physics literature, in which we take p, d, n → ∞ while keeping the ratios α = n/p, γ = d/p fixed.

Let $\mathrm{F}\in {\mathbb{R}}^{d\times p}$ be a fixed matrix, ${\boldsymbol{w}}^{a}\in {\mathbb{R}}^{p}$, 1 ⩽ a ⩽ r be a family of vectors, ${\boldsymbol{\theta }}^{0}\in {\mathbb{R}}^{d}$ be a fixed vector and $\sigma :\mathbb{R}\to \mathbb{R}$ be a scalar function acting component-wise on vectors.

Let $\boldsymbol{c}\in {\mathbb{R}}^{d}$ be a Gaussian vector $\mathcal{N}(0,{\mathrm{I}}_{d})$. The GET is a statement about the (joint) statistics of the following r + 1 random variables

$$\begin{equation}{\lambda }^{a}=\frac{1}{\sqrt{p}}{\boldsymbol{w}}^{a}\cdot \sigma (\boldsymbol{u})\in \mathbb{R},\qquad \nu =\frac{1}{\sqrt{d}}\boldsymbol{c}\cdot {\boldsymbol{\theta }}^{0}\in \mathbb{R}\end{equation} \tag{ B.1 }$$

in the asymptotic limit where d, p → ∞ with fixed p/d and fixed r. For simplicity, assume that σ(x) = −σ(−x) is an odd function. Further, suppose that in the previously introduced limit the following two balance conditions hold:

Condition 1:

$$\begin{equation}\frac{1}{\sqrt{d}}\sum\limits _{\rho =1}^{d}{\mathrm{F}}_{i\rho }{\mathrm{F}}_{j\rho }=O(1),\end{equation} \tag{ B.2 }$$

for any ρ.

Condition 2:

$$\begin{equation}{S}_{{\rho }_{1},\dots ,{\rho }_{q}}^{{a}_{1},\dots ,{a}_{k}}=\frac{1}{\sqrt{p}}\sum\limits _{i=1}^{p}{w}_{i}^{{a}_{1}}{w}_{i}^{{a}_{2}}\dots {w}_{i}^{{a}_{k}}{\mathrm{F}}_{i{\rho }_{1}}{\mathrm{F}}_{i{\rho }_{2}}\dots {\mathrm{F}}_{i{\rho }_{q}}=O(1),\end{equation} \tag{ B.3 }$$

for any integers k ⩾ 0, q > 0, for any choice of indices ρ₁, ρ₂, ..., ρ_q ∈ {1, ..., d} all distinct from each other, and for any choice of indices a₁, a₂, ..., a_k ∈ {1, ..., r}. Under the aforementioned conditions, the following theorem holds:

Theorem 1. In the limit d, p → ∞ with fixed p/d, the random variables {λ^a , ν} are jointly normal, with zero mean and covariances:

$$\begin{equation*}\mathbb{E}\left[{\lambda }^{a}{\lambda }^{b}\right]=\frac{{\kappa }_{\star }^{2}}{p}{\boldsymbol{w}}^{a}\cdot {\boldsymbol{w}}^{b}+\frac{{\kappa }_{1}^{2}}{d}{\boldsymbol{s}}^{a}\cdot {\boldsymbol{s}}^{b},\quad \quad \mathbb{E}\left[{\nu }^{2}\right]=\frac{1}{d}{\Vert}{\boldsymbol{\theta }}^{0}{{\Vert}}^{2}\end{equation*}$$

$$\begin{equation}\mathbb{E}\left[{\lambda }^{a}\nu \right]=\frac{{\kappa }_{1}}{d}{\boldsymbol{s}}^{a}\cdot {\boldsymbol{\theta }}^{0}\end{equation} \tag{ B.4 }$$

where:

$$\begin{equation}{\boldsymbol{s}}^{a}=\frac{1}{\sqrt{p}}\mathrm{F}{\boldsymbol{w}}^{a}\in {\mathbb{R}}^{d},\quad a=1,\dots ,r\end{equation} \tag{ B.5 }$$

and

$$\begin{equation}{\kappa }_{0}={\mathbb{E}}_{z}\left[\sigma (z)\right],\qquad {\kappa }_{1}={\mathbb{E}}_{z}\left[z\sigma (z)\right],\qquad {\kappa }_{\star }^{2}={\mathbb{E}}_{z}\left[\sigma {(z)}^{2}\right]-{\kappa }_{0}^{2}-{\kappa }_{1}^{2}\end{equation} \tag{ B.6 }$$

where $z\sim \mathcal{N}(0,1)$.

Note that the GET holds for a fixed family ${\left\{{\boldsymbol{w}}^{a}\right\}}_{a=1}^{r}$ and matrix $\mathrm{F}\in {\mathbb{R}}^{d\times p}$ satisfying the balance condition from equation (B.3). In the replica setting, we will need to apply the GET under an average over r samples (referred to here as replicas) of the Gibbs distribution μ_β , introduced in equation (2.7) on the main. We therefore shall require the assumption that the balance condition equation (B.3) holds for any sample of μ_β . We refer to this stronger version of the GET as the rGE. Although proving this result is out of the scope of the present work, we check its self-consistency extensively with numerical simulations.

Given the dataset ${\left\{{\boldsymbol{x}}^{\mu },{y}^{\mu }\right\}}_{\mu =1}^{n}$ defined in appendix A.1, we define the following Gibbs measure over ${\mathbb{R}}^{p}$:

$$\begin{equation}{\mu }_{\beta }(\boldsymbol{w}\vert \left\{{\boldsymbol{x}}^{\mu },{y}^{\mu }\right\})=\frac{1}{{\mathcal{Z}}_{\beta }}\enspace {\text{e}}^{-\beta \left[\sum\limits _{\mu =1}^{n}\ell \left({y}^{\mu },{\boldsymbol{x}}^{\mu }\cdot \boldsymbol{w}\right)+\frac{\lambda }{2}{\Vert}\boldsymbol{w}{{\Vert}}_{2}^{2}\right]}=\frac{1}{{\mathcal{Z}}_{\beta }}\underset{\equiv {P}_{y}(\boldsymbol{y}\vert \boldsymbol{w}\cdot {\boldsymbol{x}}^{\mu })}{\underbrace{\prod\limits _{\mu =1}^{n}{\text{e}}^{-\beta \ell \left({y}^{\mu },{\boldsymbol{x}}^{\mu }\cdot \boldsymbol{w}\right)}}}\underset{\equiv {P}_{w}(\boldsymbol{w})}{\underbrace{\prod\limits _{i=1}^{p}{\mathrm{e}}^{-\frac{\beta \lambda }{2}{w}_{i}^{2}}}}\end{equation} \tag{ C.1 }$$

for β > 0. When β → ∞, the Gibbs measure peaks at the solution of the optimisation problem in equation (A.7)—which, in the particular case of a strictly convex loss, is unique. Note that in the second equality we defined the factorised distributions P_y and P_w , showing that μ_β can be interpreted as a posterior distribution of w given the dataset { x ^μ , y^μ }, with P_y and P_w being the likelihood and prior distributions respectively.

An exact calculation of μ_β is intractable for large values of n, p and d. However, the interest in μ_β is that in the limit n, p, d → ∞ with d/p and n/p fixed, the free energy density associated to the Gibbs measure:

$$\begin{equation}{f}_{\beta }=-\underset{p\to \infty }{\mathrm{lim}}\enspace \frac{1}{p}{\mathbb{E}}_{\left\{{\boldsymbol{x}}^{\mu },{y}^{\mu }\right\}}\enspace \mathrm{log}\enspace {\mathcal{Z}}_{\beta }\end{equation} \tag{ C.2 }$$

can be computed exactly using the replica method, and at β → ∞ give us the optimal overlaps:

$$\begin{equation}{q}_{w}=\frac{1}{p}\mathbb{E}{\Vert}\hat{\boldsymbol{w}}{{\Vert}}^{2}\quad {q}_{x}=\frac{1}{d}\mathbb{E}{\Vert}\mathrm{F}\hat{\boldsymbol{w}}{{\Vert}}^{2}\quad {m}_{x}=\frac{1}{d}\mathbb{E}\left[{\boldsymbol{\theta }}^{0}\cdot \mathrm{F}\hat{\boldsymbol{w}}\right]\end{equation} \tag{ C.3 }$$

that—as we will see—fully characterise the generalisation error defined in equation (A.9).

The replica calculation of f_β is based on a large deviation principle for the free energy density. Let

$$\begin{equation}{f}_{\beta }(\left\{{\boldsymbol{x}\enspace }^{\mu },{y}^{\mu }\right\})=-\frac{1}{p}\enspace \mathrm{log}\enspace {\mathcal{Z}}_{\beta }\end{equation} \tag{ C.4 }$$

be the free energy density for one given sample of the problem, i.e. a fixed dataset ${\left\{{\boldsymbol{x}\enspace }^{\mu },{y}^{\mu }\right\}}_{\mu =1}^{n}$. We assume that the distribution P(f) of the free energy density, seen as a random variable over different samples of the problem, satisfies a large deviation principle, in the sense that, in the thermodynamic limit:

$$\begin{equation}P(f)\simeq {\mathrm{e}}^{p{\Phi}(f)},\end{equation} \tag{ C.5 }$$

with Φ a concave function reaching its maximum at the free energy density f = f_β , with Φ(f_β ) = 0. This hypothesis includes the notion of self-averageness which states that the free-energy density is the same for almost all samples in the thermodynamic limit.

The value of f_β can be computed by computing the replicated partition function

$$\begin{equation}{\mathbb{E}}_{\left\{{\boldsymbol{x}}^{\mu },{y}^{\mu }\right\}}{\mathcal{Z}}_{\beta }^{r}=\int \mathrm{d}f\enspace {\mathrm{e}}^{p[{\Phi}(f)-rf]},\end{equation} \tag{ C.6 }$$

and taking the limit

$$\begin{equation}{f}_{\beta }=\underset{r\to {0}^{+}}{\mathrm{lim}}\enspace \frac{\mathrm{d}}{\mathrm{d}r}\enspace \underset{p\to \infty }{\mathrm{lim}}\left[-\frac{1}{p}\left({\mathbb{E}}_{\left\{{\boldsymbol{x}}^{\mu },{y}^{\mu }\right\}}{\mathcal{Z}}_{\beta }^{r}\right)\right].\end{equation} \tag{ C.7 }$$

Although this procedure is not fully rigorous, experience from the statistical physics of disordered systems shows that it gives exact results, and in fact the resulting expression can be verified to match the numerical simulations.

Using the replica method we need to evaluate:

$$\begin{align} {\mathbb{E}}_{\left\{{\boldsymbol{x}}^{\mu },{y}^{\mu }\right\}}{\mathcal{Z}}_{\beta }^{r}& =\int \mathrm{d}{\boldsymbol{\theta }}^{0}{P}_{\theta }({\boldsymbol{\theta }}^{0})\int \prod\limits _{a=1}^{r}\mathrm{d}\boldsymbol{w}{P}_{w}\left({\boldsymbol{w}}^{a}\right) \\ & \quad \times \prod\limits _{\mu =1}^{n}\int \mathrm{d}{y}^{\mu }\underset{(\text{I})}{\underbrace{{\mathbb{E}}_{{\boldsymbol{c}}^{\mu }}\left[{P}_{y}^{0}\left({y}^{\mu }\left\vert \frac{{\boldsymbol{c}}^{\mu }\cdot {\boldsymbol{\theta }}^{0}}{\sqrt{d}}\right.\right)\prod\limits _{a=1}^{r}{P}_{y}\left({y}^{\mu }\vert {\boldsymbol{w}}^{a}\cdot \sigma \left(\frac{1}{\sqrt{d}}{\mathrm{F}}^{\top }{\boldsymbol{c}}^{\mu }\right)\right)\right]}} \end{align} \tag{ C.8 }$$

where P_w and P_y have been defined in (C.1). In order to compute this quantity, we introduce, for each point μ in the database, the r + 1 variables

$$\begin{equation}{\nu }_{\mu }=\frac{1}{\sqrt{d}}{\boldsymbol{c}}^{\mu }\cdot {\boldsymbol{\theta }}^{0},\end{equation} \tag{ C.9 }$$

$$\begin{equation}{\lambda }_{\mu }^{a}={\boldsymbol{w}}^{a}\cdot \sigma \left(\frac{1}{\sqrt{d}}{\mathrm{F}}^{\top }{\boldsymbol{c}}^{\mu }\right).\end{equation} \tag{ C.10 }$$

Choosing c ^μ at random induces a joint distribution $P({\nu }_{\mu },{\lambda }_{\mu }^{a})$. In the thermodynamic limit p, d → ∞ with fixed p/n, and for matrices F satisfying the balance condition in equation (B.3), the rGE introduced in appendix B.2 tells us that, for a given μ, the r + 1 variables ${\left\{{\nu }_{\mu },{\lambda }_{\mu }^{a}\right\}}_{a=1}^{r}$ are Gaussian random values with zero mean and covariance given by:

$$\begin{equation}{{\Sigma}}^{ab}=\left(\begin{matrix} \rho & {M}^{a} \\ {M}^{a} & {Q}^{ab} \end{matrix}\right)\in {\mathbb{R}}^{(r+1)\times (r+1)}.\end{equation} \tag{ C.11 }$$

The elements of the covariance matrix M^a and Q^ab are the rescaled version of the so-called overlap parameters:

$$\begin{equation}\rho =\frac{1}{d}{\Vert}{\boldsymbol{\theta }}^{0}{{\Vert}}^{2},\qquad {m}_{s}^{a}=\frac{1}{d}{\boldsymbol{s}}^{a}\cdot {\boldsymbol{\theta }}^{0},\qquad {q}_{s}^{ab}=\frac{1}{d}{\boldsymbol{s}}^{a}\cdot {\boldsymbol{s}}^{b},\qquad {q}_{w}^{ab}=\frac{1}{p}{\boldsymbol{w}}^{a}\cdot {\boldsymbol{w}}^{b},\end{equation} \tag{ C.12 }$$

where ${\boldsymbol{s}}^{a}=\frac{1}{\sqrt{p}}\mathrm{F}{\boldsymbol{w}}^{a}$. They are thus given by:

$$\begin{equation}{M}^{a}={\kappa }_{1}{m}_{s}^{a},\qquad {Q}^{ab}={\kappa }_{\star }^{2}{q}_{w}^{ab}+{\kappa }_{1}^{2}{q}_{s}^{ab}.\end{equation} \tag{ C.13 }$$

where ${\kappa }_{1}={\mathbb{E}}_{z}\left[z\sigma (z)\right]$ and ${\kappa }_{\star }^{2}={\mathbb{E}}_{z}\left[\sigma {(z)}^{2}\right]-{\kappa }_{1}^{2}$ as in equation (B.6). With this notation, the asymptotic joint probability is simply written as:

$$\begin{equation}P({\nu }_{\mu },{\left\{{\lambda }_{\mu }^{a}\right\}}_{a=1}^{r})=\frac{1}{\sqrt{\mathrm{det}\left(2\pi {\Sigma}\right)}}\enspace {\mathrm{e}}^{-\frac{1}{2}\sum\limits _{a,b=0}^{r}{z}_{\mu }^{a}{\left({{\Sigma}}^{-1}\right)}^{ab}{z}_{\mu }^{b}}\end{equation} \tag{ C.14 }$$

with ${z}_{\mu }^{0}={\nu }_{\mu }$ and ${z}_{\mu }^{a}={\lambda }_{\mu }^{a}$ for a = 1, ..., r. The average over the replicated partition function (C.8) therefore reads:

$$\begin{align} {\mathbb{E}}_{\left\{{\boldsymbol{x}}^{\mu },{y}^{\mu }\right\}}{\mathcal{Z}}_{\beta }^{r}& =\int \mathrm{d}{\boldsymbol{\theta }}^{0}{P}_{\theta }({\boldsymbol{\theta }}^{0})\int \prod\limits _{a=1}^{r}\mathrm{d}\boldsymbol{w}{P}_{w}\left({\boldsymbol{w}}^{a}\right)\prod\limits _{\mu =1}^{n}\int \mathrm{d}{y}^{\mu } \\ & \quad \times \int \mathrm{d}{\nu }_{\mu }{P}_{y}^{0}\left({y}^{\mu }\vert {\nu }_{\mu }\right)\int \prod\limits _{a=1}^{r}\mathrm{d}{\lambda }_{\mu }^{a}P({\nu }_{\mu },\left\{{\lambda }_{\mu }^{a}\right\})\prod\limits _{a=1}^{r}{P}_{y}\left({y}^{\mu }\vert \left\{{\lambda }_{\mu }^{a}\right\}\right). \end{align} \tag{ C.15 }$$

C.2.1. Rewriting as a saddle-point problem

Note that after taking the average over x , the integrals involved in the replicated partition function only couple through the overlap parameters. It is therefore useful to introduce the following Dirac δ-functions to unconstrain them, introducing the decomposition:

$$\begin{align} 1& ={d}^{-{(r+1)}^{2}}\int \mathrm{d}\rho \delta \left(\mathrm{d}\rho -{\Vert}{\boldsymbol{\theta }}^{0}{{\Vert}}^{2}\right)\int \prod\limits _{a=1}^{r}\mathrm{d}{m}_{s}^{a}\delta \left(\mathrm{d}{m}_{s}^{a}-{\boldsymbol{s}}^{a}\cdot {\boldsymbol{\theta }}^{\mathbf{0}}\right) \\ & \quad \times \int \prod\limits _{1\leqslant a\leqslant b\leqslant r}\mathrm{d}{q}_{s}^{ab}\delta \left(\mathrm{d}{q}_{s}^{ab}-{\boldsymbol{s}}^{a}\cdot {\boldsymbol{s}}^{b}\right)\int \prod\limits _{1\leqslant a\leqslant b\leqslant r}\mathrm{d}{q}_{w}^{ab}\delta \left(p{q}_{w}^{ab}-{\boldsymbol{w}}^{a}\cdot {\boldsymbol{w}}^{b}\right) \\ & ={d}^{-{(r+1)}^{2}}\int \frac{\mathrm{d}\rho \enspace \mathrm{d}\hat{\rho }}{2\pi }\enspace {\text{e}}^{-\text{i}\hat{\rho }\left(d\rho -{\Vert}{\boldsymbol{\theta }}^{0}{{\Vert}}^{2}\right)}\int \prod\limits _{a=1}^{r}\frac{\mathrm{d}{m}_{s}^{a}\enspace \mathrm{d}{\hat{m}}_{s}^{a}}{2\pi }\enspace {\text{e}}^{-\text{i}\sum\limits _{a=1}^{r}{\hat{m}}_{s}^{a}\left(\mathrm{d}{m}_{s}^{a}-{\boldsymbol{s}}^{a}\cdot {\boldsymbol{\theta }}^{\mathbf{0}}\right)} \\ & \quad \times \int \prod\limits _{1\leqslant a\leqslant b\leqslant r}\frac{\mathrm{d}{q}_{s}^{ab}\enspace \mathrm{d}{\hat{q}}_{s}^{ab}}{2\pi }\enspace {\text{e}}^{-\text{i}\sum\limits _{1\leqslant a\leqslant b\leqslant r}{\hat{q}}_{s}^{ab}\left(d{q}_{s}^{ab}-{\boldsymbol{s}}^{a}\cdot {\boldsymbol{s}}^{b}\right)}\int \prod\limits _{1\leqslant a\leqslant b\leqslant r}\frac{\mathrm{d}{q}_{w}^{ab}{\hat{q}}_{w}^{ab}}{2\pi }\enspace {\text{e}}^{-\text{i}\sum\limits _{1\leqslant a\leqslant b\leqslant r}{\hat{q}}_{w}^{ab}\left(p{q}_{w}^{ab}-{\boldsymbol{w}}^{a}\cdot {\boldsymbol{w}}^{b}\right)}. \end{align} \tag{ C.16 }$$

Introducing the above in equation (C.15) and exchanging the integration order allows to factorise the integrals over the d, p, n dimensions and rewrite:

$$\begin{equation}{\mathbb{E}}_{\left\{{\boldsymbol{x}}^{\mu },{y}^{\mu }\right\}}{\mathcal{Z}}_{\beta }^{r}=\int \frac{\mathrm{d}\rho \enspace \mathrm{d}\hat{\rho }}{2\pi }\int \prod\limits _{a=1}^{r}\frac{\mathrm{d}{m}_{s}^{a}\enspace \mathrm{d}{\hat{m}}_{s}^{a}}{2\pi }\int \prod\limits _{1\leqslant a\leqslant b\leqslant r}\frac{\mathrm{d}{q}_{s}^{ab}\enspace \mathrm{d}{\hat{q}}_{s}^{ab}}{2\pi }\frac{\mathrm{d}{q}_{w}^{ab}\enspace \mathrm{d}{\hat{q}}_{w}^{ab}}{2\pi }\enspace {\mathrm{e}}^{p{{\Phi}}^{(r)}}\end{equation} \tag{ C.17 }$$

where the integrals over the variables ${m}_{s}^{a}$, ${q}_{s}^{ab}$ and ${q}_{w}^{ab}$ run over $\mathbb{R}$, while those over ${\hat{m}}_{s}^{a}$, ${\hat{q}}_{s}^{ab}$ and ${\hat{q}}_{w}^{ab}$ run over $i\mathbb{R}$. The function Φ^(r), a function of all the overlap parameters, is given by:

$$\begin{align} {{\Phi}}^{(r)}& =-\gamma \rho \hat{\rho }-\gamma \sum\limits _{a=1}^{r}{m}_{s}^{a}{\hat{m}}_{s}^{a}-\sum\limits _{1\leqslant a\leqslant b\leqslant r}\left(\gamma {q}_{s}^{ab}{\hat{q}}_{s}^{ab}+{q}_{w}{\hat{q}}_{w}\right)+\alpha {{\Psi}}_{y}^{(r)}\left(\rho ,{m}_{s}^{a},{q}_{s}^{ab},{q}_{w}^{ab}\right) \\ & \quad +{{\Psi}}_{w}^{(r)}\left(\hat{\rho },{\hat{m}}_{s}^{a},{\hat{q}}_{s}^{ab},{\hat{q}}_{w}^{ab}\right) \end{align} \tag{ C.18 }$$

where we recall that α = n/p, γ = d/p, and we have introduced:

$$\begin{align} \begin{aligned} {{\Psi}}_{y}^{(r)}& =\mathrm{log}\int \mathrm{d}y\int \mathrm{d}\nu {P}_{y}^{0}\left(y\vert \nu \right)\int \prod\limits _{a=1}^{r}\left[\mathrm{d}{\lambda }^{a}{P}_{y}\left(y\vert {\lambda }^{a}\right)\right]P(\nu ,\left\{{\lambda }^{a}\right\}) \\ {{\Psi}}_{w}^{(r)}& =\frac{1}{p}\enspace \mathrm{log}\int \mathrm{d}{\boldsymbol{\theta }}^{0}{P}_{\theta }({\boldsymbol{\theta }}^{0}){\mathrm{e}}^{-\hat{\rho }{\Vert}{\boldsymbol{\theta }}^{0}{{\Vert}}^{2}}\int \prod\limits _{a=1}^{r}\mathrm{d}{\boldsymbol{w}}^{a}\enspace {P}_{w}({\boldsymbol{w}}^{a}){\mathrm{e}}^{\sum\limits _{1\leqslant a\leqslant b\leqslant r}\left[{\hat{q}}_{w}^{ab}{\boldsymbol{w}}^{a}\cdot {\boldsymbol{w}}^{b}+{\hat{q}}_{s}^{ab}{\boldsymbol{s}}^{a}\cdot {\boldsymbol{s}}^{b}\right]-\sum\limits _{a=1}^{r}{\hat{m}}_{s}^{a}{\boldsymbol{s}}^{a}\cdot {\boldsymbol{\theta }}^{0}}. \end{aligned}\end{align} \tag{ C.19 }$$

Note that ${\boldsymbol{s}}^{a}=\frac{1}{\sqrt{p}}\mathrm{F}{\boldsymbol{w}}^{a}$ is a function of w ^a , and must be kept under the w ^a integral. In the thermodynamic limit where p → ∞ with n/p and d/p fixed, the integral in equation (C.17) concentrates around the values of the overlap parameters that extremize Φ^(r), and therefore

$$\begin{equation}f=-\underset{r\to {0}^{+}}{\mathrm{lim}}\enspace \frac{1}{r}\underset{\begin{subarray}{c}\;\;\;\left\{\rho ,\hat{\rho },{m}_{s}^{a},{\hat{m}}_{s}^{a}\right\}\\ \left\{{q}_{s}^{ab},{\hat{q}}_{s}^{ab},{q}_{w}^{ab},{\hat{q}}_{w}^{ab}\right\}\end{subarray}}{\mathbf{e}\mathbf{x}\mathbf{t}\mathbf{r}}{{\Phi}}^{(r)}.\end{equation} \tag{ C.20 }$$

C.2.2. Replica symmetric ansatz

In order to proceed with the r → 0⁺ limit, we restrict the extremization above to the following replica symmetric ansatz:

$$\begin{align} \begin{aligned} {m}_{s}^{a}={m}_{s}& & {\hat{m}}^{a}={\hat{m}}_{s}& & \quad \quad \quad \text{for}\enspace a=1,\dots ,r\\ {q}_{s/w}^{aa}={r}_{s/w}& & \qquad {\hat{q}}_{s/w}^{aa}=-\frac{1}{2}{\hat{r}}_{s/w}& & \quad \quad \text{for}\enspace a=1,\dots ,r\\ {q}_{s/w}^{ab}={q}_{s/w}& & {\hat{q}}_{s/w}^{ab}={\hat{q}}_{s/w}& & \quad \text{for}\enspace 1\leqslant a< b\leqslant r.\end{aligned}\end{align} \tag{ C.21 }$$

Note that, in the particular case of a convex loss function with λ > 0, the replica symmetric ansatz is justified: the problem only admitting one solution, it a fortiori coincides with the replica symmetric one. For non-convex losses, solutions that are not replica symmetric (also known as replica symmetry breaking) are possible, and the energy landscape of the free energy needs to be carefully analysed. In the practical applications explored in this manuscript, we focus on convex losses with ridge regularisation, and therefore the replica symmetric assumption is fully justified.

Before proceeding with the limit in equation (C.20), we need to verify that the above ansatz is well defined—in other words, that we have not introduced a spurious order one term in Φ that would diverge. This means we need to check that $\underset{r\to {0}^{+}}{\mathrm{lim}}\enspace {\Phi}=0$.

First, with a bit of algebra one can check that, within our replica symmetric ansatz:

$$\begin{equation}\underset{r\to {0}^{+}}{\mathrm{lim}}\enspace {{\Psi}}_{y}^{(r)}=0.\end{equation} \tag{ C.22 }$$

Therefore,

$$\begin{equation}\underset{r\to {0}^{+}}{\mathrm{lim}}\enspace {{\Phi}}^{(r)}=-\gamma \rho \hat{\rho }+\gamma \enspace \mathrm{log}{\int }_{\mathbb{R}}\mathrm{d}{\theta }^{0}{P}_{\theta }\left({\theta }^{0}\right){\mathrm{e}}^{\hat{\rho }{{\theta }^{0}}^{2}}\end{equation} \tag{ C.23 }$$

where we have used the fact that P_θ is a factorised distribution to take the p → ∞ limit. In order for this limit to be 0, we need that $\hat{\rho }=0$, which also fixes ρ to be a constant given by the second moment of θ⁰:

$$\begin{equation}\rho ={\mathbb{E}}_{{\theta }^{0}}\left[{{\theta }^{0}}^{2}\right].\end{equation} \tag{ C.24 }$$

We now proceed with the limit in equation (C.20). Let us look first at Ψ_y . The non-trivial limit comes from the fact that det Σ and Σ⁻¹ are non-trivial functions of r. It is not hard to see, however, that Σ⁻¹ itself has replica symmetric structure, with components given by:

$$\begin{align} \begin{aligned} {\left({{\Sigma}}^{-1}\right)}^{00}& =\tilde{\rho }=\frac{R+(r-1)Q}{\rho (R+(r-1)Q)-r{M}^{2}}, \\ {\left({{\Sigma}}^{-1}\right)}^{aa}& =\tilde{R}=\frac{\rho R+(r-2)\rho Q-(r-1){M}^{2}}{(R-Q)(\rho R+(r-1)\rho Q-r{M}^{2})} \\ {\left({{\Sigma}}^{-1}\right)}^{a0}& =\tilde{M}=\frac{M}{r{M}^{2}-\rho R-(r-1)\rho Q}, \\ {\left({{\Sigma}}^{-1}\right)}^{ab}& =\tilde{Q}=\frac{{M}^{2}-\rho Q}{(R-Q)(\rho R+(r-1)\rho Q-r{M}^{2})} \end{aligned}\end{align} \tag{ C.25 }$$

where M, Q and R are the rescaled overlap parameters in the replica symmetric ansatz, that is:

$$\begin{equation}M={\kappa }_{1}{m}_{s},\qquad Q={\kappa }_{\star }^{2}{q}_{w}+{\kappa }_{1}^{2}{q}_{s},\qquad R={\kappa }_{\star }^{2}{r}_{w}+{\kappa }_{1}^{2}{r}_{s}.\end{equation} \tag{ C.26 }$$

This allows us to write:

$$\begin{align} {{\Psi}}_{y}^{(r)}& =\mathrm{log}\int \mathrm{d}y\int \mathrm{d}\nu {P}_{y}^{0}\left(y\vert \nu \right){\mathrm{e}}^{-\frac{\tilde{\rho }}{2}{\nu }^{2}}\int \prod\limits _{a=1}^{r}\mathrm{d}{\lambda }^{a}{P}_{y}\left(y\vert {\lambda }^{a}\right){\mathrm{e}}^{-\frac{\tilde{Q}}{2}\sum\limits _{a,b=1}^{n}{\lambda }^{a}{\lambda }^{b}-\frac{\tilde{R}-\tilde{Q}}{2}\sum\limits _{a=1}^{r}{\left({\lambda }^{a}\right)}^{2}-\tilde{M}\nu \sum\limits _{a=1}^{n}{\lambda }^{a}} \\ & \quad -\frac{1}{2}\enspace \mathrm{log}\enspace \mathrm{det}\left(2\pi {\Sigma}\right). \end{align} \tag{ C.27 }$$

In order to completely factor the integral above in the replica space, we use the Hubbard–Stratonovich transformation:

$$\begin{equation}{\mathrm{e}}^{-\frac{\tilde{Q}}{2}\sum\limits _{a,b=1}^{r}{\lambda }^{a}{\lambda }^{b}}={\mathbb{E}}_{\xi }\enspace {\mathrm{e}}^{\sqrt{-\tilde{Q}}\xi \sum\limits _{a=1}^{r}{\lambda }^{a}}\end{equation} \tag{ C.28 }$$

for $\xi \sim \mathcal{N}(0,1)$, such that

$$\begin{align} {{\Psi}}_{y}^{(r)}& =\mathrm{log}\int \mathrm{d}y\int \mathrm{d}\nu {P}_{y}^{0}\left(y\vert \nu \right){\mathrm{e}}^{-\frac{\tilde{\rho }}{2}{\nu }^{2}}{\mathbb{E}}_{\xi }{\left[\int \mathrm{d}\lambda {P}_{y}\left(y\vert \lambda \right){\mathrm{e}}^{-\frac{\tilde{R}-\tilde{Q}}{2}{\lambda }^{2}+\left(\sqrt{-\tilde{Q}}\xi -\tilde{M}\nu \right)\lambda }\right]}^{r} \\ & \quad -\frac{1}{2}\enspace \mathrm{log}\enspace \mathrm{det}\left(2\pi {\Sigma}\right). \end{align} \tag{ C.29 }$$

Taking into account the r dependence of the inverse elements and of the determinant, we can take the limit to get:

$$\begin{align} \underset{r\to {0}^{+}}{\mathrm{lim}}\enspace \frac{1}{r}{{\Psi}}_{y}^{(r)}& ={\mathbb{E}}_{\xi }{\int }_{\mathbb{R}}\mathrm{d}y\int \frac{\mathrm{d}\nu }{\sqrt{2\pi \rho }}{P}_{y}^{0}\left(y\vert \nu \right){\mathrm{e}}^{-\frac{1}{2\rho }{\nu }^{2}} \\ & \quad \times \mathrm{log}\int \frac{\mathrm{d}\lambda }{\sqrt{2\pi }}{P}_{y}\left(y\vert \lambda \right){\mathrm{e}}^{-\frac{1}{2}\frac{{\lambda }^{2}}{R-Q}+\left(\frac{\sqrt{Q-{M}^{2}/\rho }}{R-Q}\xi +\frac{M/\rho }{R-Q}\nu \right)\lambda } \\ & \quad -\frac{1}{2}\enspace \mathrm{log}\left(R-Q\right)-\frac{1}{2}\frac{Q}{R-Q}. \end{align} \tag{ C.30 }$$

Finally, making a change of variables and defining:

$$\begin{equation}{\mathcal{Z}}_{y}^{\cdot /0}(y\hspace{-2pt};\omega ,V)=\int \frac{\mathrm{d}x}{\sqrt{2\pi V}}\enspace {\mathrm{e}}^{-\frac{1}{2V}{(x-\omega )}^{2}}{P}_{y}^{\cdot /0}\left(y\vert x\right)\end{equation} \tag{ C.31 }$$

allows us to rewrite the limit of Ψ_y —which abusing notation we still denote Ψ_y —as:

$$\begin{equation}{{\Psi}}_{y}={\mathbb{E}}_{\xi }\left[{\int }_{\mathbb{R}}\mathrm{d}y{\mathcal{Z}}_{y}^{0}\left(y;\frac{M}{\sqrt{Q}}\xi ,\rho -\frac{{M}^{2}}{Q}\right)\mathrm{log}\enspace {\mathcal{Z}}_{y}\left(y;\sqrt{Q}\xi ,R-Q\right)\right].\end{equation} \tag{ C.32 }$$

One can follow a very similar approach for the limit of Ψ_w , although in this case the limit is much simpler, since there is no r dependence on the hat variables. The limit can be written as:

$$\begin{equation}{{\Psi}}_{w}=\underset{p\to \infty }{\mathrm{lim}}\enspace \frac{1}{p}{\mathbb{E}}_{\xi ,\eta ,{\theta }^{0}}\enspace \mathrm{log}{\int }_{{\mathbb{R}}^{d}}\mathrm{d}\boldsymbol{s}{P}_{s}(\boldsymbol{s}\hspace{-2pt};\eta ){\mathrm{e}}^{-\frac{{\hat{V}}_{s}}{2}{\Vert}\boldsymbol{s}{{\Vert}}^{2}+{\left(\sqrt{{\hat{q}}_{s}}\xi {\mathbf{1}}_{d}+{\hat{m}}_{s}{\boldsymbol{\theta }}^{0}\right)}^{\top }\boldsymbol{s}}\end{equation} \tag{ C.33 }$$

for $\xi ,\eta \sim \mathcal{N}(0,1)$, and we have defined:

$$\begin{equation}{P}_{s}(\boldsymbol{s};\eta )={\int }_{{\mathbb{R}}^{p}}\mathrm{d}\boldsymbol{w}{P}_{w}(\boldsymbol{w}){\mathrm{e}}^{-\frac{{\hat{V}}_{w}}{2}{\Vert}\boldsymbol{w}{{\Vert}}^{2}+\sqrt{{\hat{q}}_{w}}\eta {\mathbf{1}}_{p}^{\top }\boldsymbol{w}}\delta \left(\boldsymbol{s}-\frac{1}{\sqrt{p}}\mathrm{F}\boldsymbol{w}\right)\end{equation} \tag{ C.34 }$$

and we have defined the shorthands ${\hat{V}}_{w}={\hat{r}}_{w}+{\hat{q}}_{w}$ and ${\hat{V}}_{s}={\hat{r}}_{s}+{\hat{q}}_{s}$.

C.2.3. Summary of the replica symmetric free energy density

Summarising the calculation above, the replica symmetric free energy density reads:

$$\begin{align} f& =\mathbf{e}\mathbf{x}\mathbf{t}\mathbf{r}\left\{-\frac{\gamma }{2}{r}_{s}{\hat{r}}_{s}-\frac{\gamma }{2}{q}_{s}{\hat{q}}_{s}+\gamma {m}_{s}{\hat{m}}_{s}-\frac{1}{2}{r}_{w}{\hat{r}}_{w}-\frac{1}{2}{q}_{w}{\hat{q}}_{w}\right. \\ & \quad \left.-\alpha {{\Psi}}_{y}(R,Q,M)-{{\Psi}}_{w}\left({\hat{r}}_{s},{\hat{q}}_{s},{\hat{m}}_{s},{\hat{r}}_{w},{\hat{q}}_{w}\right)\right\} \end{align} \tag{ C.35 }$$

with $\alpha =\frac{n}{p}$, $\gamma =\frac{d}{p}$, and:

$$\begin{equation}Q={\kappa }_{1}^{2}{q}_{s}+{\kappa }_{\star }^{2}{q}_{w},\qquad R={\kappa }_{1}^{2}{r}_{s}+{\kappa }_{\star }^{2}{r}_{w}\quad \quad M={\kappa }_{1}{m}_{s}.\end{equation} \tag{ C.36 }$$

The so-called potentials (Ψ_y , Ψ_w ) are given by:

$$\begin{equation}{{\Psi}}_{w}=\underset{p\to \infty }{\mathrm{lim}}\enspace \frac{1}{p}{\mathbb{E}}_{\xi ,\eta ,{\theta }^{0}}\enspace \mathrm{log}{\int }_{{\mathbb{R}}^{d}}\mathrm{d}\boldsymbol{s}{P}_{s}(\boldsymbol{s};\eta ){\mathrm{e}}^{-\frac{{\hat{V}}_{s}}{2}{\Vert}\boldsymbol{s}{{\Vert}}^{2}+{\left(\sqrt{{\hat{q}}_{s}}\xi {\mathbf{1}}_{d}+{\hat{m}}_{s}{\boldsymbol{\theta }}^{0}\right)}^{\top }\boldsymbol{s}}\end{equation} \tag{ C.37 }$$

$$\begin{equation}{{\Psi}}_{y}={\mathbb{E}}_{\xi }\left[{\int }_{\mathbb{R}}\mathrm{d}y{\mathcal{Z}}_{y}^{0}\left(y;\frac{M}{\sqrt{Q}}\xi ,\rho -\frac{{M}^{2}}{Q}\right)\mathrm{log}{\mathcal{Z}}_{y}\left(y;\sqrt{Q}\xi ,R-Q\right)\right]\end{equation} \tag{ C.38 }$$

where:

$$\begin{align} \begin{aligned} {P}_{s}(\boldsymbol{s};\eta )& ={\int }_{{\mathbb{R}}^{p}}\mathrm{d}\boldsymbol{w}{P}_{w}(\boldsymbol{w}){\mathrm{e}}^{-\frac{{\hat{V}}_{w}}{2}{\Vert}\boldsymbol{w}{{\Vert}}^{2}+\sqrt{{\hat{q}}_{w}}\eta {\mathbf{1}}_{p}^{\top }\boldsymbol{w}}\delta \left(\boldsymbol{s}-\frac{1}{\sqrt{p}}\mathrm{F}\boldsymbol{w}\right) \\ {\mathcal{Z}}_{y}^{\cdot /0}(y;\omega ,V)& =\int \frac{\mathrm{d}x}{\sqrt{2\pi V}}\enspace {\mathrm{e}}^{-\frac{1}{2V}{(x-\omega )}^{2}}{P}_{y}^{\cdot /0}\left(y\vert x\right). \end{aligned}\end{align} \tag{ C.39 }$$

Note that as long as the limit in Ψ_w is well defined, the equation (C.35) holds for any P_θ and P_w . However, as discussed in appendix A.1, we are interested in ${\boldsymbol{\theta }}^{0}\sim \mathcal{N}(0,{\mathrm{I}}_{d})$ and ridge regularisation so that ${P}_{w}=\mathrm{exp}\left(-\frac{\beta \lambda }{2}{\Vert}\boldsymbol{w}{{\Vert}}^{2}\right)$. In this case, we simply have:

$$\begin{equation}P(\boldsymbol{s};\eta )=\frac{{\mathrm{e}}^{\frac{p}{2}\frac{{\eta }^{2}{\hat{q}}_{w}}{\beta \lambda +{\hat{V}}_{w}}}}{{(\beta \lambda +{\hat{V}}_{w})}^{p/2}}\mathcal{N}(\boldsymbol{s};\boldsymbol{\mu },{\Sigma})\end{equation} \tag{ C.40 }$$

with:

$$\begin{equation}\boldsymbol{\mu }=\frac{\sqrt{{\hat{q}}_{w}}\eta }{\beta \lambda +{\hat{V}}_{w}}\frac{\mathrm{F}{\mathbf{1}}_{p}}{\sqrt{p}}\in {\mathbb{R}}^{d},\qquad {\Sigma}=\frac{1}{\beta \lambda +{\hat{V}}_{w}}\frac{\mathrm{F}{\mathrm{F}}^{\top }}{p}\in {\mathbb{R}}^{d\times d}.\end{equation} \tag{ C.41 }$$

Therefore the argument of the logarithm in Ψ_w is just another Gaussian integral we can do explicitly:

$$\begin{equation}{\mathbb{E}}_{s}\enspace {\mathrm{e}}^{-\frac{{\hat{V}}_{s}}{2}{\Vert}\boldsymbol{s}{{\Vert}}^{2}+{\boldsymbol{b}}^{\top }\boldsymbol{s}}=\frac{{\mathrm{e}}^{\frac{p}{2}\frac{{\eta }^{2}{\hat{q}}_{w}}{\beta \lambda +{\hat{V}}_{w}}}}{{\left(\beta \lambda +{\hat{V}}_{w}\right)}^{p/2}}\frac{{\mathrm{e}}^{-\frac{1}{2}{\boldsymbol{\mu }}^{\top }{{\Sigma}}^{-1}\boldsymbol{\mu }+\frac{1}{2{\hat{V}}_{s}}{\Vert}\boldsymbol{b}+{{\Sigma}}^{-1}\boldsymbol{\mu }{{\Vert}}^{2}}}{\sqrt{\mathrm{det}\left({\mathrm{I}}_{d}+{\hat{V}}_{s}{\Sigma}\right)}}{\mathrm{e}}^{-\frac{1}{2{\hat{V}}_{s}}{\left(b+{{\Sigma}}^{-1}\mu \right)}^{\top }{\left({\mathrm{I}}_{d}+{\hat{V}}_{s}{\Sigma}\right)}^{-1}\left(b+{{\Sigma}}^{-1}\mu \right)}\end{equation} \tag{ C.42 }$$

where we have defined the shorthand $\boldsymbol{b}=\left(\sqrt{{\hat{q}}_{s}}\xi {\mathbf{1}}_{d}+{\hat{m}}_{s}{\boldsymbol{\theta }}^{0}\right)\in {\mathbb{R}}^{d}$. Inserting back in equation (C.37) and taking the log,

$$\begin{align} {{\Psi}}_{w}& =\underset{p\to \infty }{\mathrm{lim}}\enspace {\mathbb{E}}_{{\theta }^{0},\xi ,\eta }\left[\frac{1}{2}\frac{{\eta }^{2}{\hat{q}}_{w}}{\beta \lambda +{\hat{V}}_{w}}-\frac{1}{2}\enspace \mathrm{log}\left(\beta \lambda +{\hat{V}}_{w}\right)-\frac{1}{2p}\mathrm{t}\mathrm{r}\enspace \mathrm{log}\left({\mathrm{I}}_{d}+{\hat{V}}_{s}{\Sigma}\right)-\frac{1}{2p}{\boldsymbol{\mu }}^{\top }{{\Sigma}}^{-1}\boldsymbol{\mu }\right. \\ & \quad \left.+\frac{1}{2p{\hat{V}}_{s}}{\Vert}\boldsymbol{b}+{{\Sigma}}^{-1}\boldsymbol{\mu }{{\Vert}}^{2}-\frac{1}{2p{\hat{V}}_{s}}{\left(b+{{\Sigma}}^{-1}\mu \right)}^{\top }{\left({\mathrm{I}}_{d}+{\hat{V}}_{s}{\Sigma}\right)}^{-1}\left(b+{{\Sigma}}^{-1}\mu \right)\right]. \end{align} \tag{ C.43 }$$

The averages over η, ξ, θ ⁰ simplify this expression considerably:

$$\begin{align} & {\mathbb{E}}_{\eta }\left[{\boldsymbol{\mu }}^{\top }{{\Sigma}}^{-1}\boldsymbol{\mu }\right]=\frac{1}{p}\frac{{\hat{q}}_{w}}{{(\beta \lambda +{\hat{V}}_{w})}^{2}}{\left(\mathrm{F}{\mathbf{1}}_{p}\right)}^{\top }{{\Sigma}}^{-1}\left(\mathrm{F}{\mathbf{1}}_{p}\right)=d\frac{{\hat{q}}_{w}}{\beta \lambda +{\hat{V}}_{w}} \\ & {\mathbb{E}}_{\eta ,\xi ,{\theta }^{0}}{\Vert}b+{{\Sigma}}^{-1}\boldsymbol{\mu }{{\Vert}}^{2}=d({\hat{m}}_{s}^{2}+{\hat{q}}_{s})+\frac{1}{p}{\hat{q}}_{w}\enspace \mathrm{t}\mathrm{r}{\left(\mathrm{F}{\mathrm{F}}^{\top }\right)}^{-1} \\ & {\mathbb{E}}_{\eta ,\xi ,{\theta }^{0}}{\left(b+{{\Sigma}}^{-1}\boldsymbol{\mu }\right)}^{\top }{\left({\mathrm{I}}_{d}+{\hat{V}}_{s}{\Sigma}\right)}^{-1}\left(b+{{\Sigma}}^{-1}\boldsymbol{\mu }\right)=\frac{1}{p}{\hat{q}}_{w}\enspace \mathrm{t}\mathrm{r}\left[\mathrm{F}{\mathrm{F}}^{\top }{\left({\mathrm{I}}_{d}+{\hat{V}}_{s}{\Sigma}\right)}^{-1}\right] \\ & +({\hat{m}}_{s}^{2}+{\hat{q}}_{s})\mathrm{t}\mathrm{r}{\left({\mathrm{I}}_{d}+{\hat{V}}_{s}{\Sigma}\right)}^{-1}. \end{align} \tag{ C.44 }$$

Finally, we can combine the two terms:

$$\begin{equation}\mathrm{t}\mathrm{r}\frac{\mathrm{F}{\mathrm{F}}^{\top }}{p}+\mathrm{t}\mathrm{r}\left[\frac{\mathrm{F}{\mathrm{F}}^{\top }}{p}{\left({\mathrm{I}}_{d}+{\hat{V}}_{s}{\Sigma}\right)}^{-1}\right]=\frac{{\hat{V}}_{s}}{\beta \lambda +{\hat{V}}_{w}}\enspace \mathrm{t}\mathrm{r}{\left({\mathrm{I}}_{d}+{\hat{V}}_{s}{\Sigma}\right)}^{-1},\end{equation} \tag{ C.45 }$$

and write:

$$\begin{align} {{\Psi}}_{w}& =-\frac{1}{2}\enspace \mathrm{log}\left(\beta \lambda +{\hat{V}}_{w}\right)-\frac{1}{2}\enspace \underset{p\to \infty }{\mathrm{lim}}\frac{1}{p}\enspace \mathrm{t}\mathrm{r}\enspace \mathrm{log}\left({\mathrm{I}}_{d}+\frac{{\hat{V}}_{s}}{\beta \lambda +{\hat{V}}_{w}}\frac{\mathrm{F}{\mathrm{F}}^{\top }}{p}\right) \\ & \quad +\frac{{\hat{m}}_{s}^{2}+{\hat{q}}_{s}}{2{\hat{V}}_{s}}\left[\gamma -\underset{p\to \infty }{\mathrm{lim}}\enspace \frac{1}{p}\enspace \mathrm{t}\mathrm{r}{\left({\mathrm{I}}_{d}+\frac{{\hat{V}}_{s}}{\beta \lambda +{\hat{V}}_{w}}\frac{\mathrm{F}{\mathrm{F}}^{\top }}{p}\right)}^{-1}\right] \\ & \quad +\frac{1}{2}\frac{{\hat{q}}_{w}}{\beta \lambda +{\hat{V}}_{w}}\left[1-\gamma +\underset{p\to \infty }{\mathrm{lim}}\enspace \frac{1}{p}\enspace \mathrm{t}\mathrm{r}{\left({\mathrm{I}}_{d}+\frac{{\hat{V}}_{s}}{\beta \lambda +{\hat{V}}_{w}}\frac{\mathrm{F}{\mathrm{F}}^{\top }}{p}\right)}^{-1}\right]. \end{align} \tag{ C.46 }$$

Note that Ψ only depends on the spectral properties of the matrix $\frac{1}{p}\mathrm{F}{\mathrm{F}}^{\top }\in {\mathbb{R}}^{d\times d}$, and more specifically on its resolvent in the asymptotic limit. A case of particular interest is when FF^⊤ has a well defined spectral measure μ on the p, d → ∞ limit with γ = d/p fixed. In that case, we can write:

$$\begin{equation}\underset{p\to \infty }{\mathrm{lim}}\enspace \frac{1}{p}\enspace \mathrm{t}\mathrm{r}{\left({\mathrm{I}}_{d}+\frac{{\hat{V}}_{s}}{\beta \lambda +{\hat{V}}_{w}}\frac{\mathrm{F}{\mathrm{F}}^{\top }}{p}\right)}^{-1}=\gamma \frac{\beta \lambda +{\hat{V}}_{w}}{{\hat{V}}_{s}}{g}_{\mu }\left(-\frac{\beta \lambda +{\hat{V}}_{w}}{{\hat{V}}_{s}}\right)\end{equation} \tag{ C.47 }$$

where g_μ is the Stieltjes transform of μ, defined by:

$$\begin{equation}{g}_{\mu }(z)=\int \frac{\mathrm{d}\mu (t)}{t-z}.\end{equation} \tag{ C.48 }$$

Similarly, the logarithm term can be expressed as the logarithm potential of μ—although for the purpose of evaluating the generalisation error we will only need the derivative of these terms, and therefore only the Stieltjes transforms and its derivative.

In what follows, we will mostly focus on two kinds of projection matrices F:

Gaussian projections: for $\mathrm{F}\in {\mathbb{R}}^{d\times p}$ a random matrix with i.i.d. Gaussian entries with zero mean and variance 1, μ is given by the well-known Marchenko–Pastur law, and the corresponding Stieltjes transform is given by:

$$\begin{equation}{g}_{\mu }(z)=\frac{1-z-\gamma -\sqrt{{(z-1-\gamma )}^{2}-4\gamma }}{2z\gamma },\quad z< 0.\end{equation} \tag{ C.49 }$$

Orthogonally invariant projection: for F = U^⊤DV with $\mathrm{U}\in {\mathbb{R}}^{d\times d}$ and $\mathrm{V}\in {\mathbb{R}}^{p\times p}$ two orthogonal matrices and $\mathrm{D}\in {\mathbb{R}}^{d\times p}$ a rectangular diagonal matrix of rank min(d, p) and diagonal entries d_k , the empirical spectral density μ_p is given by:

$$\begin{equation}{\mu }_{d}(\lambda )=\frac{1}{d}\sum\limits _{k=1}^{\mathrm{min}(r,p)}\delta (\lambda -{\lambda }_{k})=\left(1-\mathrm{min}\left(1,\frac{1}{\gamma }\right)\right)\delta (\lambda )+\frac{1}{p}\sum\limits _{k=1}^{\mathrm{min}(d,p)}\delta (\lambda -{d}_{k}^{2}).\end{equation} \tag{ C.50 }$$

Therefore the choice of diagonal elements d_k fully characterise the spectrum of FF^⊤. In order for the orthogonally invariant case to be comparable to the Gaussian case, we fix d_k in such a way that the projected vector F w is of the same order in both cases, i.e.

$$\begin{equation}{d}_{k}^{2}=\begin{cases}\gamma \quad & \quad \text{for}\enspace \gamma > 1 \\ 1\quad & \quad \text{for}\enspace \gamma \leqslant 1 \end{cases}.\end{equation} \tag{ C.51 }$$

With this choice, the Stieltjes transform of μ reads:

$$\begin{equation}{g}_{\mu }(z)=\begin{cases}-\left(1-\frac{1}{\gamma }\right)\frac{1}{z}+\frac{1}{\gamma }\frac{1}{\gamma -z}\quad & \quad \text{for}\enspace \gamma > 1 \\ \frac{1}{1-z}\quad & \quad \text{for}\enspace \gamma \leqslant 1. \end{cases}\end{equation} \tag{ C.52 }$$

It is interesting to note that the average over the dataset ${\left\{{\boldsymbol{x}}^{\mu },{y}^{\mu }\right\}}_{\mu =1}^{n}$ of the replicated partition function ${\mathcal{Z}}_{\beta }^{r}$ in equation (C.15), obtained after the application of the GET, is identical to the replicated partition function of the same task over the following dual dataset ${\left\{{\tilde{\boldsymbol{x}}}^{\mu },{y}^{\mu }\right\}}_{\mu =1}^{n}$, where:

$$\begin{equation}{\tilde{\boldsymbol{x}}}^{\mu }={\kappa }_{0}{\mathbf{1}}_{p}+{\kappa }_{1}\frac{1}{\sqrt{d}}{\mathrm{F}}^{\top }{\boldsymbol{c}}^{\mu }+{\kappa }_{\star }{\boldsymbol{z}}^{\mu }\end{equation} \tag{ C.53 }$$

where ${\boldsymbol{z}}^{\mu }\sim \mathcal{N}(\mathbf{0},{\mathrm{I}}_{p})$, and the labels y^μ ∼ P_y are the same. Indeed, calling ${\tilde{\mathcal{Z}}}_{\beta }^{r}$ the replicated partition function for this equivalent dataset, and considering κ₀ we have:

$$\begin{align} {\mathbb{E}}_{\left\{{\tilde{\boldsymbol{x}}}^{\mu },{y}^{\mu }\right\}}{\tilde{\mathcal{Z}}}_{\beta }^{r}& =\int \mathrm{d}{\boldsymbol{\theta }}^{0}{P}_{\theta }({\boldsymbol{\theta }}^{0})\int \prod\limits _{a=1}^{r}\mathrm{d}\boldsymbol{w}{P}_{w}\left({\boldsymbol{w}}^{a}\right)\times \\ & \quad \times \prod\limits _{\mu =1}^{n}\int \mathrm{d}{y}^{\mu }\underset{(\text{I})}{\underbrace{{\mathbb{E}}_{{\boldsymbol{c}}^{\mu },{\boldsymbol{z}}^{\mu }}\left[{P}_{y}^{0}\left({y}^{\mu }\left\vert \frac{{\boldsymbol{c}}^{\mu }\cdot {\boldsymbol{\theta }}^{0}}{\sqrt{d}}\right.\right)\prod\limits _{a=1}^{r}{P}_{y}\left({y}^{\mu }\vert {\boldsymbol{w}}^{a}\cdot \left(\frac{{\kappa }_{1}}{\sqrt{d}}{\mathrm{F}}^{\top }{\boldsymbol{c}}^{\mu }+{\kappa }_{\star }{\boldsymbol{z}}^{\mu }\right)\right)\right]}}. \end{align} \tag{ C.54 }$$

Rewriting (I):

$$\begin{align} (\text{I})& =\int \mathrm{d}{\nu }_{\mu }{P}_{y}^{0}\left({y}^{\mu }\vert {\nu }_{\mu }\right)\int \prod\limits _{a=1}^{r}\mathrm{d}{\lambda }_{\mu }^{a}{P}_{y}\left({y}^{\mu }\vert {\lambda }_{\mu }^{a}\right)\times \\ & \quad \times \underset{\equiv P(\nu ,\lambda )}{\underbrace{{\mathbb{E}}_{{\boldsymbol{c}}^{\mu },{\boldsymbol{z}}^{\mu }}\left[\delta \left({\nu }_{\mu }-\frac{1}{\sqrt{d}}{\boldsymbol{c}}^{\mu }\cdot {\boldsymbol{\theta }}^{0}\right)\prod\limits _{a=1}^{r}\delta \left({\lambda }_{\mu }^{a}-\frac{{\kappa }_{1}}{\sqrt{d}}{\boldsymbol{w}}^{a}\cdot {\mathrm{F}}^{\top }{\boldsymbol{c}}^{\mu }+{\kappa }_{\star }{\boldsymbol{w}}^{a}\cdot {\boldsymbol{z}}^{\mu }\right)\right]}}. \end{align} \tag{ C.55 }$$

It is easy to show that taking (κ₀, κ₁) to match those from equation (B.6), the variables $\left({\nu }_{\mu },\left\{{\lambda }_{\mu }^{a}\right\}\right)$ are jointly Gaussian variables with correlation matrix given by Σ exactly as in equation (C.11). This establishes the equivalence

$$\begin{equation}{\tilde{\mathcal{Z}}}_{\beta }^{r}={\mathcal{Z}}_{\beta }^{r}\end{equation} \tag{ C.56 }$$

from which follows the equivalence between the asymptotic generalisation and test error of these two models.

Let { x ^new, y^new} be a new sample independently drawn from the same distribution of our data ${\left\{{\boldsymbol{x}}^{\mu },{y}^{\mu }\right\}}_{\mu =1}^{n}$. The generalisation error can then be written as:

$$\begin{align} {{\epsilon}}_{g}& =\frac{1}{{4}^{k}}{\mathbb{E}}_{{\boldsymbol{x}}^{\mathrm{n}\mathrm{e}\mathrm{w}},{y}^{\mathrm{n}\mathrm{e}\mathrm{w}}}{\left({y}^{\mathrm{n}\mathrm{e}\mathrm{w}}-\hat{f}\left(\sigma \left({\mathrm{F}}^{\top }{\boldsymbol{c}}^{\mathrm{n}\mathrm{e}\mathrm{w}}\right)\cdot \hat{\boldsymbol{w}}\right)\right)}^{2} \\ & =\frac{1}{{4}^{k}}\int \mathrm{d}y\int \mathrm{d}\nu {P}_{y}^{0}(y\vert \nu )\int \mathrm{d}\lambda {(y-\hat{f}(\lambda ))}^{2}{\mathbb{E}}_{{\boldsymbol{c}}^{\mathrm{n}\mathrm{e}\mathrm{w}}}\left[\delta \left(\nu -{\boldsymbol{c}}^{\mathrm{n}\mathrm{e}\mathrm{w}}\cdot {\boldsymbol{\theta }}^{0}\right)\delta \left(\lambda -\sigma \left({\mathrm{F}}^{\top }{\boldsymbol{c}}^{\mathrm{n}\mathrm{e}\mathrm{w}}\right)\cdot \hat{\boldsymbol{w}}\right)\right] \end{align} \tag{ D.2 }$$

where for convenience, we normalise k = 0 for the regression task and k = 1 for the classification task. Again, we apply the GET from appendix B to write the joint distribution over {ν, λ}:

$$\begin{equation}P(\nu ,\lambda )=\frac{1}{\sqrt{\mathrm{det}\left(2\pi {\Sigma}\right)}}\enspace {\mathrm{e}}^{-\frac{1}{2}{\boldsymbol{z}}^{\top }{{\Sigma}}^{-1}\boldsymbol{z}},\end{equation} \tag{ D.3 }$$

where $\boldsymbol{z}={(\nu ,\lambda )}^{\top }\in {\mathbb{R}}^{2}$ and Σ is given by

$$\begin{align} {\Sigma}=\left(\begin{matrix} \rho & {M}^{\star } \\ {M}^{\star } & {Q}^{\star } \end{matrix}\right),\quad & \rho =\frac{1}{d}{\Vert}{\boldsymbol{\theta }}^{\mathbf{0}}{{\Vert}}^{2} & {M}^{\star }=\frac{{\kappa }_{1}}{d}\left(\mathrm{F}\hat{\boldsymbol{w}}\right)\cdot {\boldsymbol{\theta }}^{0},\quad & {Q}^{\star }=\frac{{\kappa }_{1}^{2}}{d}{\Vert}\mathrm{F}\hat{\boldsymbol{w}}{{\Vert}}^{2}+\frac{{\kappa }_{\star }^{2}}{p}{\Vert}\hat{\boldsymbol{w}}{{\Vert}}^{2}. \end{align} \tag{ D.4 }$$

Inserting in equation (D.2) gives the desired representation of the generalisation error in terms of the optimal overlap parameters:

$$\begin{equation}{{\epsilon}}_{g}=\frac{1}{{4}^{k}}\int \mathrm{d}y\int \mathrm{d}\nu \enspace {P}_{y}^{0}(y\vert \nu )\int \mathrm{d}\lambda P(\nu ,\lambda ){(y-\hat{f}(\lambda ))}^{2}.\end{equation} \tag{ D.5 }$$

For linear labels y = c ⋅ θ ⁰ in the regression problem, we simply have:

$$\begin{equation}{{\epsilon}}_{g}=\rho +{Q}^{\star }-2{M}^{\star }\end{equation} \tag{ D.6 }$$

while for the corresponding classification problem with $y=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(\boldsymbol{c}\cdot {\boldsymbol{\theta }}^{0}\right)$:

$$\begin{equation}{{\epsilon}}_{g}=\frac{1}{\pi }\enspace {\mathrm{cos}}^{-1}\left(\frac{{M}^{\star }}{\sqrt{{Q}^{\star }}}\right)\end{equation} \tag{ D.7 }$$

which, as expected, only depend on the angle between $\mathrm{F}\hat{\boldsymbol{w}}$ and θ ⁰.

Similarly to the generalisation error, the asymptotic of the training loss, defined for the training data ${\left\{{\boldsymbol{x}}^{\mu },{y}^{\mu }\right\}}_{\mu =1}^{n}$ as:

$$\begin{equation}{{\epsilon}}_{t}=\frac{1}{n}{\mathbb{E}}_{\left\{{\boldsymbol{x}}^{\mu },{y}^{\mu }\right\}}\left[\sum\limits _{\mu =1}^{n}\ell \left({y}^{\mu },{\boldsymbol{x}}^{\mu }\cdot \hat{\boldsymbol{w}}\right)+\frac{\lambda }{2}{\Vert}\hat{\boldsymbol{w}}{{\Vert}}_{2}^{2}\right],\end{equation} \tag{ D.8 }$$

can also be written only in terms of the overlap parameters. Indeed, it is closely related to the free energy density defined in equation (C.2). A close inspection on this definition tells us that:

$$\begin{equation}\underset{n\to \infty }{\mathrm{lim}}\enspace {{\epsilon}}_{t}=\underset{\beta \to \infty }{\mathrm{lim}}\enspace {\partial }_{\beta }{f}_{\beta }.\end{equation} \tag{ D.9 }$$

Taking the derivative of the free energy with respect to the parameter β and recalling that p = αn, we can then get:

$$\begin{equation}\underset{n\to \infty }{\mathrm{lim}}\enspace {{\epsilon}}_{t}=\frac{\lambda }{2\alpha }\underset{p\to \infty }{\mathrm{lim}}\enspace {\mathbb{E}}_{\left\{{\boldsymbol{x}}^{\mu },{y}^{\mu }\right\}}\left[\frac{{\Vert}\hat{\boldsymbol{w}}{{\Vert}}_{2}^{2}}{p}\right]-\underset{\beta \to \infty }{\mathrm{lim}}\enspace {\partial }_{\beta }{{\Psi}}_{y}.\end{equation} \tag{ D.10 }$$

For what concerns the contribution of the regulariser, we simply note that in the limit of p → ∞, the average concentrates around the overlap parameter ${q}_{w}^{\star }$. Instead, for what concerns the contribution of the loss function, we can start by explicitly taking the derivative with respect to β of Ψ_y in equation (C.32), i.e.:

$$\begin{equation}{\partial }_{\beta }{{\Psi}}_{y}=-{\mathbb{E}}_{\xi }\left[{\int }_{\mathbb{R}}\mathrm{d}y\frac{{\mathcal{Z}}_{y}^{0}\left(y,{\omega }_{0}^{\star }\right)}{{\mathcal{Z}}_{y}\left(y,{\omega }_{1}^{\star }\right)}\int \frac{\mathrm{d}x}{\sqrt{2\pi {V}_{1}^{\star }}}\enspace {\mathrm{e}}^{-\frac{1}{2{V}_{1}^{\star }}{(x-{\omega }_{1}^{\star })}^{2}-\beta \ell \left(y,x\right)}\ell \left(y,x\right)\right],\end{equation} \tag{ D.11 }$$

with ${\mathcal{Z}}_{y}^{\cdot /0}$ defined in equation (C.31). At this point, as explained more in details in appendix D.4, we can notice that in the limit of β → ∞, it holds:

$$\begin{equation}\underset{\beta \to \infty }{\mathrm{lim}}\enspace {\partial }_{\beta }\enspace {{\Psi}}_{y}=-{\mathbb{E}}_{\xi }\left[{\int }_{\mathbb{R}}\mathrm{d}y{\mathcal{Z}}_{y}^{0}\left(y,{\omega }_{0}^{\star }\right)\ell \left(y,\eta \left(y,{\omega }_{1}^{\star }\right)\right)\right],\end{equation} \tag{ D.12 }$$

with $\eta \left(y,{\omega }_{1}^{\star }\right)$ given in equation (D.21). Combining the two results together we then finally get:

$$\begin{equation}\underset{n\to \infty }{\mathrm{lim}}\enspace {{\epsilon}}_{t}\to \frac{\lambda }{2\alpha }{q}_{w}^{\star }+{\mathbb{E}}_{\xi }\left[{\int }_{\mathbb{R}}\mathrm{d}y{\mathcal{Z}}_{y}^{0}\left(y,{\omega }_{0}^{\star }\right)\ell \left(y,\eta \left(y,{\omega }_{1}^{\star }\right)\right)\right].\end{equation} \tag{ D.13 }$$

As we showed above, both the generalisation error and the training loss are completely determined by the β → ∞ solution of the extremization problem in equation (C.35). For strictly convex losses ℓ, there is a unique solution to this problem, that can be found by considering the derivatives of the replica potential. This leads to a set of self-consistent saddle-point equations that can be solved iteratively:

$$\begin{align} \begin{cases}{\hat{r}}_{s}=-2\frac{\alpha {\kappa }_{1}^{2}}{\gamma }{\partial }_{{r}_{s}}{{\Psi}}_{y}\left(R,Q,M\right)\quad \\ {\hat{q}}_{s}=-2\frac{\alpha {\kappa }_{1}^{2}}{\gamma }{\partial }_{{q}_{s}}{{\Psi}}_{y}\left(R,Q,M\right)\quad \\ {\hat{m}}_{s}=\frac{\alpha {\kappa }_{1}}{\gamma }{\partial }_{{m}_{s}}{{\Psi}}_{y}\left(R,Q,M\right)\quad \\ {\hat{r}}_{w}=-2\alpha {\kappa }_{\star }^{2}{\partial }_{{r}_{w}}{{\Psi}}_{y}\left(R,Q,M\right)\quad \\ {\hat{q}}_{w}=-2\alpha {\kappa }_{\star }^{2}{\partial }_{{q}_{w}}{{\Psi}}_{y}\left(R,Q,M\right)\quad \end{cases}& \begin{cases}{r}_{s}=-\frac{2}{\gamma }{\partial }_{{\hat{r}}_{s}}{{\Psi}}_{w}\left({\hat{r}}_{s},{\hat{q}}_{s},{\hat{m}}_{s},{\hat{r}}_{w},{\hat{q}}_{w}\right)\quad \\ {q}_{s}=-\frac{2}{\gamma }{\partial }_{{\hat{q}}_{s}}{{\Psi}}_{w}\left({\hat{r}}_{s},{\hat{q}}_{s},{\hat{m}}_{s},{\hat{r}}_{w},{\hat{q}}_{w}\right)\quad \\ {m}_{s}=\frac{1}{\gamma }{\partial }_{{\hat{m}}_{s}}{{\Psi}}_{w}\left({\hat{r}}_{s},{\hat{q}}_{s},{\hat{m}}_{s},{\hat{r}}_{w},{\hat{q}}_{w}\right)\quad \\ {r}_{w}=-2{\partial }_{{\hat{r}}_{w}}{{\Psi}}_{w}\left({\hat{r}}_{s},{\hat{q}}_{s},{\hat{m}}_{s},{\hat{r}}_{w},{\hat{q}}_{w}\right)\quad \\ {q}_{w}=-2{\partial }_{{\hat{q}}_{w}}{{\Psi}}_{w}\left({\hat{r}}_{s},{\hat{q}}_{s},{\hat{m}}_{s},{\hat{r}}_{w},{\hat{q}}_{w}\right).\quad \end{cases} \end{align} \tag{ D.14 }$$

In the case of an F with well-defined spectral density μ, we can be more explicit and write:

$$\begin{equation}\begin{cases}{V}_{s}=\frac{1}{{\hat{V}}_{s}}\left(1-z{g}_{\mu }(-z)\right)\quad \\ {q}_{s}=\frac{{\hat{m}}_{s}^{2}+{\hat{q}}_{s}}{{\hat{V}}_{s}^{2}}\left[1-2z{g}_{\mu }(-z)+{z}^{2}{g}_{\mu }^{\prime }(-z)\right]-\frac{{\hat{q}}_{w}}{(\beta \lambda +{\hat{V}}_{w}){\hat{V}}_{s}}\left[-z{g}_{\mu }(-z)+{z}^{2}{g}_{\mu }^{\prime }(-z)\right]\quad \\ {m}_{s}=\frac{{\hat{m}}_{s}}{{\hat{V}}_{s}}\left(1-z{g}_{\mu }(-z)\right)\quad \\ {V}_{w}=\frac{\gamma }{\beta \lambda +{\hat{V}}_{w}}\left[\frac{1}{\gamma }-1+z{g}_{\mu }(-z)\right]\quad \\ {q}_{w}=\gamma \frac{{\hat{q}}_{w}}{{(\beta \lambda +{\hat{V}}_{w})}^{2}}\left[\frac{1}{\gamma }-1+{z}^{2}{g}_{\mu }^{\prime }(-z)\right]-\gamma \frac{{\hat{m}}_{s}^{2}+{\hat{q}}_{s}}{(\beta \lambda +{\hat{V}}_{w}){\hat{V}}_{s}}\left[-z{g}_{\mu }(-z)+{z}^{2}{g}_{\mu }^{\prime }(-z)\right]\quad \end{cases}\end{equation} \tag{ D.15 }$$

where:

$$\begin{equation}{V}_{s/w}={r}_{s/w}-{q}_{r/w}\quad {\hat{V}}_{s/w}={\hat{r}}_{s/w}+{\hat{q}}_{r/w}\quad z=\frac{\beta \lambda +{\hat{V}}_{w}}{{\hat{V}}_{s}}.\end{equation} \tag{ D.16 }$$

We can also simplify slightly the derivatives of Ψ_y without loosing generality by applying Stein's lemma, yielding:

$$\begin{equation}\begin{cases}{\hat{V}}_{s}=-\frac{\alpha {\kappa }_{1}^{2}}{\gamma }{\mathbb{E}}_{\xi }\left[{\int }_{\mathbb{R}}\mathrm{d}y{\mathcal{Z}}_{y}^{0}\left(y\hspace{-2pt};\frac{M}{\sqrt{Q}}\xi ,\rho -\frac{{M}^{2}}{Q}\right){\partial }_{\omega }{f}_{y}\left(y\hspace{-2pt};\sqrt{Q}\xi ,R-Q\right)\right]\quad \\ {\hat{q}}_{s}=\frac{\alpha {\kappa }_{1}^{2}}{\gamma }{\mathbb{E}}_{\xi }\left[{\int }_{\mathbb{R}}\mathrm{d}y{\mathcal{Z}}_{y}^{0}\left(y\hspace{-2pt};\frac{M}{\sqrt{Q}}\xi ,\rho -\frac{{M}^{2}}{Q}\right){f}_{y}{\left(y\hspace{-2pt};\sqrt{Q}\xi ,R-Q\right)}^{2}\right]\quad \\ {\hat{m}}_{s}=\frac{\alpha {\kappa }_{1}}{\gamma }{\mathbb{E}}_{\xi }\left[{\int }_{\mathbb{R}}\mathrm{d}y{\mathcal{Z}}_{y}^{0}\left(y\hspace{-2pt};\frac{M}{\sqrt{Q}}\xi ,\rho -\frac{{M}^{2}}{Q}\right){f}_{y}^{0}\left(y\hspace{-2pt};\frac{M}{\sqrt{Q}}\xi ,\rho -\frac{{M}^{2}}{Q}\right){f}_{y}\left(y\hspace{-2pt};\sqrt{Q}\xi ,R-Q\right)\right]\quad \\ {\hat{V}}_{w}=-\alpha {\kappa }_{\star }^{2}{\mathbb{E}}_{\xi }\left[{\int }_{\mathbb{R}}\mathrm{d}y{\mathcal{Z}}_{y}^{0}\left(y\hspace{-2pt};\frac{M}{\sqrt{Q}}\xi ,\rho -\frac{{M}^{2}}{Q}\right){\partial }_{\omega }{f}_{y}\left(y\hspace{-2pt};\sqrt{Q}\xi ,R-Q\right)\right]\quad \\ {\hat{q}}_{w}=\alpha {\kappa }_{\star }^{2}{\mathbb{E}}_{\xi }\left[{\int }_{\mathbb{R}}\mathrm{d}y{\mathcal{Z}}_{y}^{0}\left(y\hspace{-2pt};\frac{M}{\sqrt{Q}}\xi ,\rho -\frac{{M}^{2}}{Q}\right){f}_{y}{\left(y\hspace{-2pt};\sqrt{Q}\xi ,R-Q\right)}^{2}\right]\quad \end{cases}\end{equation} \tag{ D.17 }$$

with ${f}_{y}^{\cdot /0}(y;\omega ,V)={\partial }_{\omega }\enspace \mathrm{log}\enspace {\mathcal{Z}}_{y}^{\cdot /0}$. For a given choice of spectral density μ (corresponding to a choice of projection F), label rule ${P}_{y}^{0}$ and loss function ℓ, the auxiliary functions $({\mathcal{Z}}^{0},\mathcal{Z})$ can be computed, and from them the right-hand side of the update equations above. The equations can then be iterated until the convergence to the fixed point minimising the free energy at fixed (α, γ, β). For convex losses and β → ∞, the fixed point of these equations gives the overlap corresponding to the estimator solving equation (A.7).

Although the saddle-point equations above can be iterated explicitly for any β > 0, it is envisageable to take the limit β → ∞ explicitly, since β is an auxiliary parameter we introduced, and that was not present in the original problem defined in equation (A.7).

Since the overlap parameters depend on β only implicitly through ${\mathcal{Z}}_{y}$ and its derivatives, we proceed with the following ansatz for their β → ∞ scaling:

$$\begin{align} \begin{aligned} {V}_{s/w}^{\infty }=\beta {V}_{s/w}& & {q}_{s/w}^{\infty }={q}_{s/w}& & \qquad {m}_{s}^{\infty }={m}_{s}\\ {\hat{V}}_{s/w}^{\infty }=\frac{1}{\beta }{\hat{V}}_{s/w}& & \qquad {\hat{q}}_{s/w}^{\infty }=\frac{1}{{\beta }^{2}}{\hat{q}}_{s/w}& & \qquad {\hat{m}}_{s}^{\infty }={\hat{m}}_{s}.\end{aligned}\end{align} \tag{ D.18 }$$

This ansatz can be motivated as follows. Recall that:

$$\begin{equation}{\mathcal{Z}}_{y}(y;\omega ,V)=\int \frac{\mathrm{d}x}{\sqrt{2\pi V}}\enspace {\text{e}}^{-\beta \left[\frac{{(x-\omega )}^{2}}{2\beta V}+\ell (x,y)\right]}=\int \frac{\mathrm{d}x}{\sqrt{2\pi V}}\enspace {\text{e}}^{-\beta \mathcal{L}(x)}.\end{equation} \tag{ D.19 }$$

Therefore, letting $V={\mu }_{1}^{2}{V}_{s}+{\mu }_{\star }^{2}{V}_{w}$ scale as V^∞ = βV, at β → ∞:

$$\begin{equation}{\mathcal{Z}}_{y}(y;\omega ,V)\underset{\beta \to \infty }{=}{\text{e}}^{-\beta \mathcal{L}(\eta )}\end{equation} \tag{ D.20 }$$

where:

$$\begin{equation}\eta (y;\omega ,V)=\underset{x\in \mathbb{R}}{\mathrm{a}\mathrm{r}\mathrm{g}\;\;\mathrm{m}\mathrm{i}\mathrm{n}}\left[\frac{{(x-\omega )}^{2}}{2{V}^{\infty }}+\ell (x,y)\right].\end{equation} \tag{ D.21 }$$

For convex losses ℓ with λ > 0, this one-dimensional minimisation problem has a unique solution that can be easily evaluated. Intuitively, this ansatz translates the fact the variance of our estimator goes to zero as a power law at β → ∞, meaning the Gibbs measure concentrates around the solution of the optimisation problem equation (A.7). The other scalings in equation (D.19) follow from analysing the dependence of the saddle-point equations in V.

The ansatz in equation (D.18) allow us to take the β → ∞ and rewrite the saddle-point equations as:

$$\begin{equation}\begin{cases}{\hat{V}}_{s}^{\infty }=\frac{\alpha {\mu }_{1}^{2}}{\gamma }{\mathbb{E}}_{\xi }\left[{\int }_{\mathbb{R}}\mathrm{d}y{\mathcal{Z}}_{y}^{0}\left(\frac{1-{\partial }_{\omega }\eta }{{V}^{\infty }}\right)\right]\quad \\ {\hat{q}}_{s}^{\infty }=\frac{\alpha {\mu }_{1}^{2}}{\gamma }{\mathbb{E}}_{\xi }\left[{\int }_{\mathbb{R}}\mathrm{d}y{\mathcal{Z}}_{y}^{0}{\left(\frac{\eta -\omega }{{V}^{\infty }}\right)}^{2}\right]\quad \\ {\hat{m}}_{s}^{\infty }=\frac{\alpha {\mu }_{1}}{\gamma }{\mathbb{E}}_{\xi }\left[{\int }_{\mathbb{R}}\mathrm{d}y{\partial }_{\omega }{\mathcal{Z}}_{y}^{0}\left(\frac{\eta -\omega }{{V}^{\infty }}\right)\right]\quad \\ {\hat{V}}_{w}^{\infty }=\alpha {\mu }_{\star }^{2}{\mathbb{E}}_{\xi }\left[{\int }_{\mathbb{R}}\mathrm{d}y{\mathcal{Z}}_{y}^{0}\left(\frac{1-{\partial }_{\omega }\eta }{{V}^{\infty }}\right)\right]\quad \\ {\hat{q}}_{w}^{\infty }=\alpha {\mu }_{\star }^{2}{\mathbb{E}}_{\xi }\left[{\int }_{\mathbb{R}}\mathrm{d}y{\mathcal{Z}}_{y}^{0}{\left(\frac{\eta -\omega }{{V}^{\infty }}\right)}^{2}\right]\quad \end{cases}\end{equation} \tag{ D.22 }$$

$$\begin{align} \\ & \begin{cases}{V}_{s}^{\infty }=\frac{1}{{\hat{V}}_{s}^{\infty }}\left(1-z{g}_{\mu }(-z)\right)\quad \\ {q}_{s}^{\infty }=\frac{{\left({\hat{m}}_{s}^{\infty }\right)}^{2}+{\hat{q}}_{s}^{\infty }}{{\left({\hat{V}}_{s}^{\infty }\right)}^{2}}\left[1 - 2z{g}_{\mu }(-z) + {z}^{2}{g}_{\mu }^{\prime }(-z)\right] - \frac{{\hat{q}}_{w}^{\infty }}{(\lambda +{\hat{V}}_{w}){\hat{V}}_{s}}\left[-z{g}_{\mu }(-z) + {z}^{2}{g}_{\mu }^{\prime }(-z)\right]\quad \\ {m}_{s}^{\infty }=\frac{{\hat{m}}_{s}^{\infty }}{{\hat{V}}_{s}^{\infty }}\left(1-z{g}_{\mu }(-z)\right)\quad \\ {V}_{w}^{\infty }=\frac{\gamma }{\lambda +{\hat{V}}_{w}^{\infty }}\left[\frac{1}{\gamma }-1+z{g}_{\mu }(-z)\right]\quad \\ {q}_{w}^{\infty }=\gamma \frac{{\hat{q}}_{w}^{\infty }}{{(\lambda +{\hat{V}}_{w}^{\infty })}^{2}}\left[\frac{1}{\gamma }-1+{z}^{2}{g}_{\mu }^{\prime }(-z)\right]-\gamma \frac{{\left({\hat{m}}_{s}^{\infty }\right)}^{2}+{\hat{q}}_{s}^{\infty }}{(\lambda +{\hat{V}}_{w}^{\infty }){\hat{V}}_{s}^{\infty }}\left[-z{g}_{\mu }(-z)+{z}^{2}{g}_{\mu }^{\prime }(-z)\right]\quad \end{cases} \end{align} \tag{ D.23 }$$

where ${\mathcal{Z}}_{y}^{0}(y;\omega ,V)$ is always evaluated at $(\omega ,V)=\left(\frac{{M}^{\infty }}{\sqrt{{Q}^{\infty }}}\xi ,\rho -\frac{{{M}^{\infty }}^{2}}{{Q}^{\infty }}\right)$, η(y; ω, V) at $\left(\omega ,V\right)=\left(\sqrt{{Q}^{\infty }}\xi ,{V}^{\infty }\right)$ and $z=\frac{\lambda +{\hat{V}}_{w}^{\infty }}{{\hat{V}}_{s}^{\infty }}$.

In this section we exemplify our general result in two particular cases for which the integrals in the right-hand side of equation (D.22) can be analytically performed: the ridge regression task with linear labels and a classification problem with square loss and ridge regularisation term. The former example appears in figure 4 (left) and the later in figure 2 (blue curve), figures 6 and 7 of the main.

Ridge regression with linear labels: consider the task of doing ridge regression $\ell (y,x)=\frac{1}{2}{\left(y-x\right)}^{2}$, λ > 0 on the linear patterns $\boldsymbol{y}=\frac{1}{\sqrt{d}}\mathrm{C}{\boldsymbol{\theta }}^{0}+\sqrt{{\Delta}}\boldsymbol{z}$, with $\boldsymbol{z}\sim \mathcal{N}(\mathbf{0},{\mathrm{I}}_{n})$ and ${\boldsymbol{\theta }}^{\star }\sim \mathcal{N}(\mathbf{0},{\mathrm{I}}_{d})$. In this case, we have:

$$\begin{equation}\eta (y;\omega ,V)=\frac{\omega +yV}{1+V}\end{equation} \tag{ D.24 }$$

and the saddle-point equations for the hat overlap read:

$$\begin{equation*}{\hat{V}}_{s}^{\infty }=\frac{\alpha }{\gamma }\frac{{\kappa }_{1}^{2}}{1+{V}^{\infty }}\quad {\hat{q}}_{s}^{0}=\frac{\alpha {\kappa }_{1}^{2}}{\gamma }\frac{1+{\Delta}+{Q}^{\infty }-2{M}^{\infty }}{{\left(1+{V}^{\infty }\right)}^{2}}\quad {\hat{m}}_{s}=\frac{\alpha }{\gamma }\frac{{\kappa }_{1}}{1+{V}^{\infty }}\end{equation*}$$

$$\begin{equation}{\hat{V}}_{w}^{\infty }=\frac{\alpha {\kappa }_{\star }^{2}}{1+{V}^{\infty }}\qquad {\hat{q}}_{w}^{\infty }=\alpha {\kappa }_{\star }^{2}\frac{1+{\Delta}+{Q}^{\infty }-2{M}^{\infty }}{{\left(1+{V}^{\infty }\right)}^{2}}.\end{equation} \tag{ D.25 }$$

This particular example corresponds precisely to the setting studied in [32].

Classification with square loss and ridge regularisation: consider a classification task with square loss $\ell (y,x)=\frac{1}{2}{\left(y-x\right)}^{2}$ and labels generated as $\boldsymbol{y}=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(\frac{1}{\sqrt{d}}\mathrm{C}{\boldsymbol{\theta }}^{0}\right)$, with ${\boldsymbol{\theta }}^{0}\sim \mathcal{N}(\mathbf{0},{\mathrm{I}}_{d})$. Then the saddle-point equations are simply:

$$\begin{equation*}{\hat{V}}_{s}^{\infty }=\frac{\alpha }{\gamma }\frac{{\kappa }_{1}^{2}}{1+{V}^{\infty }}\qquad {\hat{q}}_{s}^{\infty }=\frac{\alpha }{\gamma }{\kappa }_{1}^{2}\frac{1+{Q}^{\infty }-\frac{2\sqrt{2}{M}^{\infty }}{\sqrt{\pi }}}{{\left(1+{V}^{\infty }\right)}^{2}}\qquad {\hat{m}}_{s}=\frac{\alpha }{\gamma }\sqrt{\frac{2}{\pi }}\frac{{\kappa }_{1}}{1+{V}^{\infty }}\end{equation*}$$

$$\begin{equation}{\hat{V}}_{w}^{\infty }=\frac{\alpha {\kappa }_{\star }^{2}}{1+{V}^{\infty }}\qquad {\hat{q}}_{w}^{\infty }=\alpha {\kappa }_{\star }^{2}\frac{1+{Q}^{\infty }-\frac{2{M}^{\infty }}{\sqrt{\pi }}}{{\left(1+{V}^{\infty }\right)}^{2}}\end{equation} \tag{ D.26 }$$