Double Trouble in Double Descent : Bias and Variance(s) in the Lazy Regime

Introduction

Deep neural networks have achieved breakthroughs in a plethora of contexts, such as image classification , speech recognition , and automatic translation . Yet, theory lags far behind practice, and the key reasons underpinning the success of deep learning remain to be clarified.

One of the main puzzles is to understand the excellent generalization performance of heavily overparametrized deep neural networks able to fit random labels . Such interpolating estimators—that can reach zero training error— have attracted a growing amount of theoretical attention in the last few years, see e.g. . Indeed, classical learning theory suggests that generalization should first improve then worsen when increasing model complexity, following a U-shape curve characteristic of the bias-variance trade-off. Instead, deep neural networks as well as other machine learning models , follow a different curve, coined double descent.

This curve displays two regimes : the classical U-curve is superseded at high complexity by a modern interpolating regime where the test error decreases monotonically with overparametrization . Between these two regimes, i.e. at the interpolation threshold where training error vanishes, a peak occurs in absence of regularization, sometimes called the jamming peak due to similarities with a well-studied phenomenon in the Statistical Physics litterature . The reasons behind the performance of deep neural networks in the overparametrized regime are still poorly understood, even though some mechanisms are known to play an important role, such as the implicit regularization of stochastic gradient descent which allows to converge to the minimum norm solution, and the convergence to mean-field limits .

Here we present a detailed investigation of the double descent phenomenon, and its theoretical explanation in terms of bias and variance in the so-called lazy regime . This theoretically appealing scenario, where the weights stay close to their initial value during training, is called lazy learning as opposed to feature learning where the weights change enough to learn relevant features . Although replacing learnt features by random features may appear as a crude simplification, empirical results show that the loss in performance can be rather small in some cases . A burst of recent papers showed that in this regime, neural networks behave like kernel methods or equivalently random projection methods . This mapping makes the training analytically tractable, allowing, for example, to prove convergence to zero error solutions in overparametrized settings.

Optimization plays an important role in neural networks, and in particular for the double descent phenomenon, by inducing implicit regularization and fluctuations of the learnt estimator . Disentangling the variance stemming from the randomness of the optimization process from that the variance due to the randomness of the dataset is a crucial step towards a unified picture, as suggested in . In this paper, we address this issue and attempt to reconcile the behavior of bias and variance with the double descent phenomenon by providing a precise and quantitative theory in the lazy regime.

We focus on an analytically solvable model of random features (RF), introduced by , that can be viewed either as a randomized approximation to kernel ridge regression, or as a two-layer neural network whose first layer contains fixed random weights. The latter provides a simple model for lazy learning. Indeed, suppose that a neural network learns a function $f_{\theta}({\bf x})$ that relates labels (or responses) $y$ to inputs ${\bf x}$ via a set of weights ${\theta}$. The lazy regime is defined as the setting where the model can be linearized around the initial conditions $\theta_{0}$. Assuming that the initialization is such that $f_{\theta_{0}}\approx 0$One can alteratively define the estimator as $f_{\theta}-f_{\theta_{0}}$ ., one obtains:

In other words, the lazy regime corresponds to a linear fitting problem with a random feature vector $\left.\nabla_{\theta}f_{\theta}({\bf x})\right|_{\theta=\theta_{0}}$. In this setting our contributions are:

We demonstrate how to disentangle quantitatively the contributions to the test error of the bias and the various sources of variance of the estimator, stemming from the sampling of the dataset, from the additive noise corrupting the labels, and from the initialization of the random feature vectors.

We give a sharp asymptotic formula for the effect of ensembling (averaging the predictions of indepently initialized estimators) on these various terms. We show in particular how the over-fitting peak at the interpolation threshold can be attenuated by ensembling, as observed in real neural networks . We also compare the effect of ensembling, overparametrizing and optimally regularizing.

Several conclusions stem from the above analysis. First, the over-fitting near the interpolation threshold is entirely due to the variances due to the additive noise in the ground truth and the initialization of the random features. Second, the data sampling variance and the bias both display a phase transition at the interpolation threshold, and remain constant in the overparametrized regime. Hence, the benefit of ensembling and overparametrization beyond the interpolation threshold is solely due to a reduction of the noise and initialization variances.

Finally, we present numerical results on a classic deep learning scenario in the lazy learning regime to show that our findings, obtained for simple random features and i.i.d. data, are relevant to realistic setups involving correlated random features and realistic data.

The analytical results we present are obtained using a heuristic method from Statistical Physics called the Replica Method , which despite being non-rigorous has shown its remarkable efficacy in many machine learning problems and random matrix topics, see e.g. . While it is an open problem to provide a rigorous proof of our computations, we check through numerical simulations that our asymptotic predictions are extremely accurate at moderately small sizes.

Our work crucially builds on two recent contributions by Geiger et al. , and Mei and Montanari . The authors of carried out a series of experiments in order to shed light on the generalization properties of neural networks. The current work is inspired by their observations and scaling theory about the role of the variance due to the random initialization of the weights in the double-descent curve. They argued that the decrease of the test error in the limit of very wide networks is due to this source of variance, which vanishes inversely proportional to the width of the network. They then used ensembling to empirically support these findings in more realistic situations. Another related work is , which empirically disentangles the various sources of variance in the process of training deep neural networks.

On the analytical side, our paper builds on the results of , which provide a precise expression of the test error of the RF model in the high-dimensional limit where the number of random features, the dimension of the input data and the number of data points are sent to infinity with their relative ratios fixed. The double descent was also studied analytically for various types of linear models, both for regression and classification . An example of a practical method that uses ensembling in kernel methods is detailed in . Note that this work performs an average over the sampling of the random feature vectors in contrast to where the average is taken over the sampling of the data set.

The codes necessary to reproduce the results presented in this paper and obtain new ones are given at https://github.com/mariaref/Random_Features.git.

Model

This work is centered around the RF model first introduced in . Although simpler settings such as linear regression display the double descent phenomenology , this model is more appealing in several ways. First, the presence of two layers allows to freely disentangle the dimensionality of the input data from the number of parameters of the model. Second, it closely relates to the lazy regime of neural networks, as described above. Third, and most importantly for our specific study, the randomness of the fixed first layer weights mimics the randomness due to weight initialization in neural networks.

The generalization to non-linear functions can also be performed as in .

The second layer weights, i.e the elements of ${\bm{a}}$, are calculated by the means of ridge regression:

Note that as $P\to\infty$, this is equivalent to kernel ridge regression with respect to the following kernel:

Analytical results

In this section, we present our main result, which is an analytical expression for the terms appearing in the decomposition of the test error.

The test error can be decomposed into its bias and variance components:

Noise variance: The first term is the variance associated with the additive noise corrupting the labels of the dataset which is learnt, $\bm{\varepsilon}$:

Initialization variance: The second term encodes the fluctuations stemming from the random initialization of the random feature vectors, $\bm{\Theta}$:

Sampling variance: The third term measures the fluctuations due to the sampling of the training data, $\bm{X}$:

Bias: Finally, the fourth term is the bias, i.e. the error that remains once all the sources of variance have been averaged out. It can be understood as the approximation error of our model and takes the form:

Note that since we are performing deterministic ridge regression, the noise induced by SGD, which can play an important role outside the lazy regime for deep neural networks, cannot be captured.

Consider the high-dimensional limit where the input dimension $D$, the hidden layer dimension $P$ (which is equal to the number of parameter in our model) and the number of training points $N$ go to infinity with their ratios fixed:

where the terms $\{\Psi_{1},\Psi_{2}^{v},\Psi_{3}^{v},\Psi_{2}^{e},\Psi_{3}^{e},\Psi_{2}^{d}\}$, whose full analytical expressions are deferred to Appendix B, are computed following the steps below:

Mapping to a random matrix theory problem. The first step is to express the right-hand sides of Equations 12-17 as traces over random matrices. This is achieved by replacing our model with its asymptotically equivalent Gaussian covariate model , in which the non-linearity of the activation function is encoded as an extra noise term. This enables to take the expectation value with respect to the test sample $\bm{x}$.

Mapping to a statistical physics model. The random matrix theory problem resulting from the solution of ridge regression (4) involves inverse random matrices. In order to evaluate their expection value, we use the formula:

which is based on the Replica Trick . The Gaussian integrals over $\bm{\varepsilon},\bm{\Theta},\bm{X}$ can then be straightforwardly performed and lead to a Statistical Physics model for the auxiliary variables $\eta_{i}^{\alpha}$.

Mean-Field Theory. The model for the $\eta_{i}^{\alpha}$ variables can then be solved by introducing as order parameters the $n\times n$ overlap matrices $Q^{\alpha\beta}=\frac{1}{P}\sum_{i=1}^{P}\eta^{\alpha}_{i}\eta^{\beta}_{i}$ and using replica theory , see Appendix C for the detailed computationIn order to obtain the asymptotic formulas for the $\Psi$’s we need to compute (what are called in the Statistical Physics jargon) fluctuations around mean-field theory..

The $\Psi$’s may also be estimated numerically at finite size by evaluating the traces of the random matrices appearing in the Gaussian covariate model at the end of step 1. Figure 2 shows that results thus obtained are in excellent agreement with the asymptotic expressions even at moderate sizes, e.g. $D=200$, proving the robustness of steps 2 and 3, which differ from the approach presented in .

The indices $v,e,d$ in $\{\Psi_{1},\Psi_{2}^{v},\Psi_{3}^{v},\Psi_{2}^{e},\Psi_{3}^{e},\Psi_{2}^{d}\}$ stand for vanilla, ensemble and divide and conquer. The vanilla terms are sufficient to obtain the test error of a single RF model and were computed in . The ensemble and divide and conquer terms allow to obtain the test error obtained when averaging the predictions of several different learners trained respectively on the same dataset and on different splits of the original dataset (see section 5). Figure 2 shows that the vanilla terms exhibit a radically different behavior from the others: at vanishing regularization, they diverge at $P=N$ then decrease monotonically, whereas the others display a kink followed by a plateau. This behavior will be key to the following analysis.

Analysis of Bias and Variances

The results of the previous section, allow to rewrite the decomposition of the test error as follows:

These contributions, together with the test error, are shown in figure 3 in the case of small (top) and large (bottom) regularization.

The peak at the interpolation threshold is completely due to noise and initialization variance, which both diverge at vanishing regularization. In contrast, the sampling variance and the bias remain finite and exhibit a phase transition at $P=N$, which is revealed by a kink at vanishing regularization. Adding regularization smooths out these singular behaviours: it removes the divergence and irons out the kink.

In the overparametrized regime, the sampling variance and the bias do not vary substantially (they remain constant for vanishing regularization). The decrease of the test error is entirely due to the decrease of the noise and initalization variances for $P\!>\!N$. In the limit $P/N\!\to\!\infty$, the initialization variance vanishes, whereas there remains an irreducible noise variance.

In conclusion, we find that the origin of the double descent curve lies in the behavior of noise and initialization variances. The benefit of overparametrizing stems only from reducing these two contributions.

These results are qualitatively similar to the empirical decomposition of for real neural networks. The divergence of the test error as $(P/N-1)^{-1}$ at the interpolation threshold is in agreement with the results of Note that for classification problems the singularity is different .. As for the decrease of the test error in the over-parametrized regime, we find consistently with the scaling arguments of that the initialization error asymptotically decays to zero inversely proportional to the width (see Appendix A.1 for more details). The interpretation of our results differ from those of where the authors relate the over-fitting peak occurring at $P=N$ to a divergence in both the variance and the bias terms. This is due to the fact the bias term, as defined in that paper, also includes the initialization varianceFor a given set of random features this is legitimate, but from the perspective of lazy learning the randomness in the features corresponds to the one due to initialization, which is an additional source of variance.. When the two are disentangled, it becomes clear that it is only the latter which is responsible for the divergence: the bias is, in fact, well-behaved at $P=N$.

At $P\!=\!N$, the problem becomes fully determined: the data is perfectly interpolated for vanishing $\lambda$. Two types of noise are overfit: (i) the stochastic noise corrupting the labels, yielding the divergence in noise variance, and (ii) the deterministic noise stemming from the non-linearity of the activation function which cannot be captured, yielding the divergence in initialization variance. However, by further increasing $P$, the noise is spread over more and more random features and is effectively averaged out. Consequently, the test error decreases again as $P$ increases.

When we make the problem deterministic by averaging out all sources of randomness, i.e. by considering the bias, we see that increasing $P$ beyond $N$ has no effect whatsoever. Indeed, the extra degrees of freedom, which lie in the null space of $\bm{Z}$, do not provide any extra expressivity: at vanishing regularization, they are killed by the pseudo-inverse to reach the minimum norm solution. For non-vanishing $\lambda$, a similar phenomenology is observed but the interpolation threshold is reached slightly after $P=N$ since the expressivity of the learner is lowered by regularization.

On the effect of ensembling

In order to further study the effect of the variances on the test error, we follow and study the impact of ensembling. In the lazy regime of deep neural networks, the initial values of the weights only affect the gradient at initialization, which corresponds to the vector of random features. Hence, we can study the effect of ensembling in the lazy regime by averaging the predictions of RF models with independently drawn random feature vectors.

Consider a set of $K\!>\!1$ RF networks whose first layer weights are drawn independently. These networks are trained independently on the same training set. In the analogy outlined above, they correspond to $K$ independent inizializations of the neural network. At the end of training, one obtains $K$ estimators $\{\hat{f}_{\bm{\Theta}^{k}}\}$ ($k=1,...,K$). When a new sample ${\bm{x}}$ is presented to the system, the output is taken to be the average over the outputs of the $K$ networks, as illustrated in figure 4. By expanding the square and taking the expectation with respect to the random initalizations, the test error can then be written as:

We see that ensembling amounts to a linear interpolation between the vanilla terms $\Psi_{2}^{v},\Psi_{3}^{v}$, for $K=1$, and the ensemble terms $\Psi_{2}^{e},\Psi_{3}^{e}$ for $K\to\infty$.

The effect of ensembling on the double descent curve is shown in figure 5. As $K$ increases, the overfitting peak at the interpolation threshold is diminished. This observation is very similar to the empirical findings of in the context of real neural networks. Our analytic expression agrees with the numerical results obtained by training RF models, even at moderate size $D=200$.

Note that a related procedure is the divide and conquer approach, where the dataset is partitioned into $K$ splits of equal size and each one of the $K$ differently initialized learners is trained on a distinct split. This approach was studied for kernel learning in , and is analyzed within our framework in Appendix A.2.

The bias-variance decomposition of the test error makes the suppression of the divergence explicit. The bias and variances contribution read for the averaged estimator:

As we have shown, ensembling and overparametrizing have similar effects in the lazy regime. But which is more powerful: ensembling $K$ models, or using a single model with $K$ times more features? The answer is given in figure 7 for $K=2$ where we plot our analytical results while varying the number of data points, $N$. Two observations are particularly interesting. First, overparametrization shifts the interpolation threshold, opening up a region where ensembling outperforms overparametrizing. Second, overparametrization yields a higher asymptotic improvement in the large dataset limit $N/D\to\infty$, but the gap between overparametrizing and ensembling is reduced as $P/D$ increases. At $P\!\gg\!D$, where we are already close to the kernel limit, both methods yield a similar improvement. Note that from the point of view of efficiency, ridge regression involves the inversion of a $P\times P$ matrix, therefore ensembling is significantly more efficient.

In all the results presented above, we keep the regularization constant $\lambda$ fixed. However, by appropriately choosing the value of $\lambda$ at each value of $P/N$, the performance is improved. As figure 7 (left) reveals, the optimal value of $\lambda$ decreases with $K$ since the minimum of the test error shifts to the left when increasing $K$. In other words, ensembling is best when the predictors one ensembles upon are individually under-regularized, as was observed previously for kernel learning in . Figure 8 (right) shows that an infinitely ensembled model ($K\to\infty$) always performs better than an optimally regularized single model ($K=1$).

Numerical experiments on neural networks

Finally, we investigate whether the phenomenology described here holds for realistic neural networks learning real data in the lazy regime. We follow here the protocol used in and train a 5-layer fully-connected network on the CIFAR-10 dataset. We keep only the first ten PCA components of the images, and divide the images in two classes according to the parity of the labels. We perform $10^{5}$ steps of full-batch gradient descent with the Adam optimizer and a learning rate of $0.1$, and scale the weights as prescribed in .

We gradually go from the usual feature learning regime to the lazy learning regime using the trick introduced in , which consists in scaling the output of the network by a factor $\alpha$ and replacing the learning function $f_{\theta}(\bm{x})$ by $\alpha(f_{\theta}(\bm{x})\!-\!f_{\theta_{0}}(\bm{x}))$. For $\alpha\!\gg\!1$, one must have that $\theta\!-\!\theta_{0}\!\sim\!{1}/{\alpha}$ in order for the learning function to remain of order one. In other words, the weights are forced to stay close to their initialization, hence the name lazy learning.

Results are shown in figure 9. Close to the lazy regime ($\alpha\!=\!100$, right panel), a very similar behavior as the RF model is observed. The test error curveNote that we are considering a binary classification task here: the error is defined as the fraction of misclassified images. obtained when ensembling $K=20$ independently initialized networks becomes roughly flat after the interpolation threshold (which here is signalled by the peak in the test accuracy). As we move away from the lazy regime ($\alpha\!=\!10$, left panel), the same curve develops a dip around the interpolation threshold and increases beyond $P\!>\!N$ as observed previously in . This may arguably be associated to the beneficial effect of feature learning, as discussed in where the transition from lazy to feature learning was investigated.

We thank Matthieu Wyart and Lenka Zdeborová for discussions related to this project. This work is supported by the French Agence Nationale de la Recherche under grant ANR-17-CE23-0023-01 PAIL and ANR-19-P3IA-0001 PRAIRIE, and by the Simons Foundation ($\#$454935, Giulio Biroli). We also acknowledge support from the chaire CFM-ENS “Science des données”.

Upon completion of this paper, we became aware of two related parallel works presenting bias-variance tradeoffs for RF models. The variance due to the sampling of the dataset was considered in , whereas focused on the variance due to the randomness of the random feature vectors.

References

Appendix A Further analytical results

Figure 10 (left) shows that the various terms entering the decomposition of the generalization error approach their asymptotic values at a rate $(P/N)^{-1}$. This scaling law is consistent with that found in for real neural networks, where $P$ is replaced by the width of the layers of the network. As for the divergence of the noise and initialization variances observed at the interpolation threshold, figure 10 (right) shows that they also follow an inverse power law $(P/N-1)^{-1}$ at vanishing regularization.

A.2 Divide and Conquer approach

As mentioned in the main text, another way to average the predictions of differently initialized learners is the divide and conquer approach . In this framework, the data set is divided into $K$ splits of size $N/K$. Each of the $K$ differently initalized learner is trained on a distinct split. This approach is extremely useful for kernel learning , where the computational burden is in the inversion of the Gram matrix which is of size $N\times N$. In the random projection approach considered here, it does not offer any computational gain, however it is interesting how it affects the generalization error.

Within our framework, the generalization error can easily be calculated as:

In Figure 11, we see that the kernel limit error of the divide and conquer approach, i.e. the asymptotic value of the error at $P/N\to\infty$, is different from the usual kernel limit error, since the effective dataset is two times smaller at $K=2$. The denoising effect of the divide and conquer approach is illustrated by the fact that its kernel limit error is higher at high SNR, but lower at low SNR. This is of practical relevance, and is much related to the beneficial effect of bagging in noisy dataset scenarios. The divide and conquer approach, which only differs from bagging by the fact that the different partitions of the dataset are disjoint, was shown to reach bagging-like performance in various setups such as decision trees and neural networks .

A.3 Is it always better to be overparametrized ?

A common thought is that the double descent curve always reaches its minimum in the over-parametrized regime, leading to the idea that the corresponding model ”cannot overfit”. In this section, we show that this is not always the case. Three factors tend to shift the optimal generalization to the underparametrized regime: (i) increasing the numbers of learners from which we average the predictions, $K$, (ii) decreasing the signal-to-noise ratio (SNR), $F/\tau$, and (iii) decreasing the size of the dataset, $N/D$. In other words, when ensembling on a small, noisy dataset, one is better off using an underparametrized model.

These three effects are shown in figure 12. In the left panel, we see that as we increase $K$, the minimum of generalization error jumps to the underparametrized regime $P<N$ for a high enough value of $K$. In the central/right panels, a similar effect occurs when decreasing the SNR or decreasing $N/D$.

Appendix B Statement of the Main Result

First, we state precisely the assumptions under which our main result is valid. Note, that these are the same as in .

Assumption 2: We work in the high-dimensional limit, i.e. in the limit where the input dimension $D$, the hidden layer dimension $P$ and the number of training points $N$ go to infinity with their ratios fixed. That is:

This condition implies that, in the computation of the risk $\mathcal{R}$, we can neglect all the terms of order $\mathcal{O}(1)$ in favour of the terms of order $\mathcal{O}(D)$.

Assumption 3: The labels are given by a linear ground truth, or teacher function:

Note that as explained in , it is easy to add a non linear component to the teacher, but the latter would not be captured by the model (the student) in the regime $N/D=\mathcal{O}(1)$, and would simply amount to an extra noise term.

B.2 Results

Here we give the explicit form of the quantities appearing in our main result. In these expressions, the index $a\in\{v,e,d\}$ distinguishes the vanilla, ensembling and divide and conquer terms.

where we defined the scalars $P_{XX},P_{WX},P_{WW}$ as follows:

and the $2\times 2$ matrices $M_{X},M_{W},N_{X}$ as follows:

Appendix C Replica Computation

In order to obtain the main result for the generalisation error, we perform the averages over all the sources of randomness in the system in the following order: over the dataset $X$, then over the noise $W$, and finally over the random feature layers $\Theta$. Here are some useful formulaes used throughout the computations:

C.1.2 Replica representation of an inverse matrix

To obtain gaussian integrals we will use the ”replica” representation the element $(ij)$ of a matrix $M$ of size $D$:

Indeed, using the gaussian integral representation of the inverse of $M$,

Using the replica identity, we rewrite this as

Renaming the integration variable of the integral on the left as $\eta^{1}$ and the $n-1$ others as $\eta^{\alpha},\alpha\in\{2,n\}$, we obtain expression (33).

C.2 The Random Feature model

In what follows, we will explicitly leave the indices of all the quantities used. We use the notation, called Einstein summation convention in physics, in which all repeated indices are summed but the sum is not explicitly written. Indices $i\in\{1...D\}$ are used to refer to the input dimension, $h\in\{1...P\}$ to refer to the hidden layer dimension and $\mu\in\{1...N\}$ to refer to the number of data points.

In the random features model, the predictor can be computed explicitly:

Hence the generalization error can be computed as:

C.2.2 Ensembling over K𝐾K learners

When ensembling over $K$ learners with independently sampled random feature vectors, the predictor becomes:

The generalisation error is then given by:

C.2.3 Equivalent Gaussian Covariate Model

It was shown in that the random features model is equivalent, in the high-dimensional limit of Assumption 2, to a Gaussian covariate model in which the activation function $\sigma$ is replaced as:

This powerful mapping allows to express the quantities $\bm{U},\bm{V}$. We will not repeat their calculations here: the only difference here is $\bm{U}^{kl}$, which carries extra indices $k,l$ due to the different initialization of the random features $\bm{\Theta}^{(k)}$. In our case,

Hence we can rewrite the generalization error as

where $\Psi_{1}$,$\Psi_{2}^{v}$,$\Psi_{2}^{e}$,$\Psi_{3}^{v}$,$\Psi_{3}^{e}$ are given by:

C.3 Computation of the vanilla terms

In the vanilla terms, the two inverse matrices that appear are the same. Hence we use twice the replica identity (33), introducing $2n$ replicas which all play the same role:

The first step is to perform the averages, i.e. the Gaussian integrals, over the dataset $X$, the deterministic noise $W$ induced by the non-linearity of the activation function and the random features $\Theta$.

Replacing the activation function by its Gaussian covariate equivalent model and using (52), the term $\Psi_{3}$ can be expanded as:

Now, we introduce $\lambda_{i}^{\alpha}:=\frac{1}{\sqrt{P}}\eta_{h}^{\alpha}\bm{\Theta}_{hi}$, and enforce this relation using the Fourier representation of the delta-function:

The average over the dataset $\bm{X}_{\mu i}$ has the form of (32) with:

Note that due to with a slight abuse of notation we got rid of indices $\mu$, which all sum up trivially to give a global factor $N$.

C.3.2 Averaging over the deterministic noise

The expectation over the deterministic noise $\bm{W}_{h}$ is a Gausssian integral of the form (32) with:

Note that the prefactor involves, constant, linear and quadratic terms in $\bm{W}$ since:

C.3.3 Averaging over the random feature vectors

The expectation over the random feature vectors $\bm{\Theta}_{hi}$ is a Gausssian integral of the form (32) with:

C.3.4 Expression of the action and the prefactor

To complete the computation we integrate with respect to $\hat{\lambda}^{\alpha}_{i}$, using again formulae (32):

This yields the final expression of the term:

with the prefactor $P_{\Psi_{3}^{v}}$ and the action $S^{v}$ defined as:

C.3.5 Expression of the action and the prefactor in terms of order parameters

Here we see that we have a factor $D\to\infty$ in the exponential part, which can be estimated using the saddle point method. Before doing so, we introduce the following order parameters using the Fourier representation of the delta-function:

This allows to rewrite the prefactor only in terms of $Q,R$: for example,

To do this, there are two key quantities we need to calculate: $\lambda G^{-1}_{X}\lambda$ and $\eta G^{-1}_{W}\eta$. To calculate both, we note that $G_{X}$ ang $G_{W}$ are both of the form $\mathbf{I}+\mathbf{X}$, therefore there inverse may be calculated using their series representation. The result is:

The integrals over $\eta,\lambda$ become simple Gaussian integrals with covariance matrices given by $\hat{Q},\hat{R}$, yielding:

The next step is to take the saddle point with respect to the auxiliary variables $\hat{Q}$ and $\hat{R}$ in order to eliminate them:

C.3.6 Saddle point equations

The aim is now to use the saddle point method in order to evaluate the integrals over the order parameters. Thus, one looks for $R$ and $Q$ solutions to the equations:

To solve the above, it is common to make a replica symmetric ansatz. In this case, we assume that the solutions to the saddle points equations take the form:

C.3.7 Fluctuations around the saddle point

Therefore we must go beyond the saddle point contribution to obtain a non zero result, i.e. we have to examine the quadratic fluctuations around the saddle point. To do so we preform a second-order expansion of the action (77) as a function of $Q$ and $R$:

Computing the second derivative of (77), it is easy to show that:

The last equality follows from the fact that for a matrix of size $n\times n$ of the form $M_{\alpha\beta}=a\delta_{\alpha\beta}+b(1-\delta_{\alpha\beta})$, we have

C.3.8 Expression of the vanilla terms

Using the above procedure, one can compute the terms $\Psi_{1},\Psi_{2}^{v},\Psi_{3}^{v}$ of (51): for each of these terms, the action is the same as in (77), and the prefactors can be obtained as:

C.4 Computation of the ensembling terms

In the ensembling terms, the two inverse matrices are different, hence one has to introduce two distinct replica variables. We distinguish them by the use of an extra index $a\in\{1,2\}$, denoted in brackets in order not to be confused with the replica indices $\alpha$.

Calculations of the Gaussian integrals follow through in a very similar way as for the vanilla terms. The matrices appearing in the process are:

Starting with the computation of $\Psi_{3}^{e}$ in order to illustrate the method used, the prefactor $P_{\Psi_{3}^{e}}$ and the action are $S^{e}$ are given by:

C.4.2 Expression of the action and the prefactor in terms of order parameters

This time, because of the two different systems, the order parameters carry an additional index $a$, which turns them into $2\times 2$ block matrices:

The systems being decoupled, we make the following ansatz for the order parameters:

In virtue of the simple structure of the above matrices, the replica indices $\alpha$ trivialize and we may replace the matrices $Q$ and $R$ by the $2\times 2$ matrices:

where products are now over $2\times 2$ matrices. Then, one has:

C.4.3 Expression of the ensembling terms

Evaluating the fluctuations around the saddle point follows through in the same way as for the vanilla terms, with the following expressions of the prefactors:

C.5 Computation of the divide and conquer term

Here, we are interested in computing the term $\Psi_{2}^{d}$. This term differs from the previous ones in that there are now two independent data matrices $\bm{X}^{(1)}$ and $\bm{X}^{(2)}$. The calculations for the action and the prefactor are very similar to calculations performed for the ensembling terms $\Psi_{2}^{e},\Psi_{3}^{e}$, with the addition that $\bm{X}$ now also carries an index $(a)$. Firstly let us write $\Psi_{2}^{d}$ as a trace over random matrices:

Calculations follow through in the same way as in the previous sections. Using the replica formula (90), and performing the integrals over the Gaussian variables the following quantities appear:

The saddle point ansatz for $Q$ and $R$ is the same as the one for the ensembling terms (see (135)). The procedure to evaluate $\Psi_{2}^{d}$ is also the same as the one for $\Psi_{2}^{e}$ except the Hessian is taken with respect to $S^{d}$. The final result is given below.