# High-dimensional Asymptotics of Denoising Autoencoders

Hugo Cui  Statistical Physics of Computation laboratory
École Polytechnique Fédérale de Lausanne, 1015 Lausanne Switzerland  Lenka
Zdeborová  Statistical Physics of Computation laboratory
École Polytechnique Fédérale de Lausanne, 1015 Lausanne Switzerland

###### Abstract

We address the problem of denoising data from a Gaussian mixture using a two-
layer non-linear autoencoder with tied weights and a skip connection. We
consider the high-dimensional limit where the number of training samples and
the input dimension jointly tend to infinity while the number of hidden units
remains bounded. We provide closed-form expressions for the denoising mean-
squared test error. Building on this result, we quantitatively characterize
the advantage of the considered architecture over the autoencoder without the
skip connection that relates closely to principal component analysis. We
further show that our results accurately capture the learning curves on a
range of real data sets.

## Introduction

Machine learning techniques have a long history of success in denoising tasks.
The recent breakthrough of diffusion-based generation [1, 2] has further
revived the interest in denoising networks, demonstrating how they can also be
leveraged, beyond denoising, for generative tasks. However, this rapidly
expanding range of applications stands in sharp contrast to the relatively
scarce theoretical understanding of denoising neural networks, even for the
simplest instance thereof – namely Denoising Auto Encoders (DAEs) [3].

Theoretical studies of autoencoders have hitherto almost exclusively focused
on data compression tasks using Reconstruction Auto Encoders (RAEs), where the
goal is to learn a concise latent representation of the data. A majority of
this body of work addresses linear autoencoders [4, 5, 6, 7]. The authors of
[8, 9] analyze the gradient-based training of non-linear autoencoders with
online stochastic gradient descent or in population, thus implicitly assuming
the availability of an infinite number of training samples. Furthermore, two-
layer RAEs were shown to learn to essentially perform Principal Component
Analysis (PCA) [10, 11, 12], i.e. to learn a linear model. Ref. [13] shows
that this is also true for infinite-width architectures. Learning in DAEs has
been the object of theoretical investigations only in the linear case [14],
while the case of non-linear DAEs remains theoretically largely unexplored.

### Main contributions

The present work considers the problem of denoising data sampled from a
Gaussian mixture by learning a two-layer DAE with a skip connection and tied
weights via empirical risk minimization. Throughout the manuscript, we
consider the high-dimensional limit where the number of training samples $n$
and the dimension $d$ are large ($n,d\to\infty$) while remaining comparable,
i.e. $\alpha\equiv\nicefrac{{n}}{{d}}=\Theta(1)$. Our main contributions are:

  * •

Leveraging the replica method, we provide sharp, closed-form formulae for the
mean squared denoising test error (MSE) for DAEs, as a function of the sample
complexity $\alpha$ and the problem parameters. We also provide a sharp
characterization for other learning metrics including the weights norms, skip
connection strength, and cosine similarity between the weights and the cluster
means. These formulae encompass as a corollary the case of RAEs. We show that
these formulae also describe quantitatively rather well the denoising MSE for
real data sets, including MNIST [15] and FashionMNIST [16].

  * •

We find that PCA denoising (namely denoising by projecting the noisy data
along the principal component of the training samples) is widely sub-optimal
compared to the DAE, leading to a MSE superior by a difference of $\Theta(d)$,
thereby establishing that DAEs do not simply learn to perform PCA.

  * •

Building on the formulae, we quantify the role of each component of the DAE
architecture (skip connection and the bottleneck network) in its overall
performance. We find that the two components have complementary effects in the
denoising process –namely preserving the data nuances and removing the noise–
and discuss how the training of the DAE results from a tradeoff between these
effects.

The code used in the present manuscript can be found in the following
repository.

### Related works

#### Theory of autoencoders

– Various aspects of RAEs have been studied, for example, memorization [17],
or latent space alignment [18]. However, the largest body of work has been
dedicated to the analysis of gradient-based algorithms when training RAEs.
Ref. [5] established that minimizing the training loss leads to learning the
principal components of the data. Authors of [11, 12] have analyzed how a
linear RAE learns these components during training. These studies were later
extended to non-linear networks by [19, 8, 9], at the sacrifice of further
assuming an infinite number of training samples to be available –either by
considering online stochastic gradient descent, or the population loss. Refs.
[20, 13] are able to address a finite sample complexity, but in exchange, have
to consider infinite-width architectures, which [13] further shows, also tend
to a large extent to learn to perform PCA.

#### Exact asymptotics from the replica method

– The replica method [21, 22, 23, 24] has proven a very valuable gateway to
access sharp asymptotic characterizations of learning metrics for high-
dimensional machine learning problems. Past works have addressed –among
others– single-[25, 26, 27, 28] and multi-index models [29], or kernel methods
[30, 31, 32, 33]. While the approach has traditionally addressed convex
problems, for which its prediction can be proven e.g. using the convex Gordon
minimax theorem [34], the replica method allows to average over all the global
minimizers of the loss, and therefore also accommodates non-convex settings.
Refs. [35, 36] are two recent examples of its application to non-convex
losses. In the present manuscript, we leverage this versatility to study the
minimization of the empirical risk of DAEs, whose non-convexity represents a
considerable hurdle to many other types of analyses.

##  I Setting

#### Data model

We consider the problem of denoising data $\bm{x}\in\mathbb{R}^{d}$ corrupted
by Gaussian white noise of variance $\Delta$,

$\tilde{\bm{x}}=\sqrt{1-\Delta}\bm{x}+\sqrt{\Delta}\bm{\xi},$ (1)

where we denoted $\tilde{\bm{x}}$ the noisy data point, and
$\bm{\xi}\sim\mathcal{N}(0,\mathbb{I}_{d})$ the additive noise. The rescaling
of the clean data point by a factor $\sqrt{1-\Delta}$ is a practical choice
that entails no loss of generality, and allows to easily interpolate between
the noiseless case ($\Delta=0$) and the case where the signal-to-noise ratio
vanishes ($\Delta=1$). Furthermore, it allows us to seamlessly connect with
works on diffusion-based generative models, where the rescaling naturally
follows from the way the data is corrupted by an Ornstein-Uhlenbeck process
[1, 2]. In the present work, we assume the clean data $\bm{x}$ to be drawn
from a Gaussian mixture distribution $\mathbb{P}$ with $K$ clusters

$\displaystyle\bm{x}\sim\sum\limits_{k=1}^{K}\rho_{k}\mathcal{N}(\bm{\mu}_{k},\bm{\Sigma}_{k}).$
(2)

The $k-$th cluster is thus centered around $\bm{\mu}_{k}\in\mathbb{R}^{d}$,
has covariance $\bm{\Sigma}_{k}\succeq 0$, and relative weight $\rho_{k}$.

#### DAE model

An algorithmic way to retrieve the clean data $\bm{x}$ from the noisy data
$\tilde{\bm{x}}$ is to build a neural network taking the latter as an input
and yielding the former as an output. A particularly natural choice for such a
network is an autoencoder architecture [3]. The intuition is that the narrow
hidden layer of an autoencoder forces the network to learn a succinct latent
representation of the data, which is robust against noise corruption of the
input. In this work, we analyze a two-layer DAE. We further assume that the
weights are tied. Additionally, mirroring modern denoising architectures like
U-nets [37] or [38, 39, 40, 41], we also allow for a (trainable) skip-
connection:

$\displaystyle
f_{b,\bm{w}}(\tilde{\bm{x}})=b\times\tilde{\bm{x}}+\frac{\bm{w}^{\top}}{\sqrt{d}}\sigma\left(\frac{\bm{w}\tilde{\bm{x}}}{\sqrt{d}}\right).$
(3)

The DAE (3) is therefore parametrized by the scalar skip connection strength
$b\in\mathbb{R}$ and the weights $\bm{w}\in\mathbb{R}^{p\times d}$, with $p$
the width of the DAE hidden layer. The normalization of the weight $\bm{w}$ by
$\sqrt{d}$ in (3) is the natural choice which ensures for high dimensional
settings $d\gg 1$ that the argument of the non-linearity $\sigma(\cdot)$ stays
$\Theta(1)$. Like [8], we focus on the case with $p\ll d$. The assumption of
weight-tying affords a more concise theoretical characterization and thus
clearer discussions. Note that it is also a strategy with substantial
practical history, dating back to [3], as it prevents the DAE from functioning
in the linear region of its non-linearity $\sigma(\cdot)$. This choice of
architecture is also motivated by a particular case of Tweedie’s formula [42]
(see eq. (82) in Appendix B), which will be the object of further discussion
in Section III.

We also consider two other simple architectures

$\displaystyle
u_{\bm{v}}(\tilde{\bm{x}})=\frac{\bm{v}^{\top}}{\sqrt{d}}\sigma\left(\frac{\bm{v}\tilde{\bm{x}}}{\sqrt{d}}\right),$
$\displaystyle r_{c}(\tilde{\bm{x}})=c\times\tilde{\bm{x}},$ (4)

which correspond to the building blocks of the complete DAE architecture
$f_{b,\bm{w}}$ (3) (hereafter referred to as the full DAE). Note that indeed
$f_{b,\bm{w}}=r_{b}+u_{\bm{w}}$. The part $u_{\bm{v}}(\cdot)$ is a DAE without
skip connection (hereafter called the bottleneck network component), while
$r_{c}(\cdot)$ correspond to a simple single-parameter trainable rescaling of
the input (hereafter called the rescaling component).

To train the DAE (3), we assume the availability of a training set
$\mathcal{D}=\\{\tilde{\bm{x}}^{\mu},\bm{x}^{\mu}\\}_{\mu=1}^{n}$, with $n$
clean samples $\bm{x}^{\mu}$ drawn i.i.d from $\mathbb{P}$ (2) and the
corresponding noisy samples $\tilde{\bm{x}}^{\mu}=\bm{x}^{\mu}+\bm{\xi}^{\mu}$
(with the noises $\bm{\xi}^{\mu}$ assumed mutually independent). The DAE is
trained to recover the clean samples $\bm{x}^{\mu}$ from the noisy samples
$\tilde{\bm{x}}^{\mu}$ by minimizing the empirical risk

$\displaystyle\hat{\mathcal{R}}(b,\bm{w})=\sum\limits_{\mu=1}^{n}\left\lVert\bm{x}^{\mu}-f_{b,\bm{w}}(\tilde{\bm{x}}^{\mu})\right\lVert^{2}+g(\bm{w}),$
(5)

where $g:\mathbb{R}^{p\times d}\to\mathbb{R}_{+}$ is an arbitrary convex
regularizing function. We denote by $\hat{b},\hat{\bm{w}}$ the minimizers of
the empirical risk (5) and by $\hat{f}\equiv f_{\hat{b},\hat{\bm{w}}}$ the
corresponding trained DAE (3). For future discussion, we also consider
training independently the components (4) via empirical risk minimization, by
which we mean replacing $f_{b,\bm{w}}$ by $u_{\bm{v}}$ or $r_{c}$ in (5). We
similarly denote $\hat{\bm{v}}$ (resp. $\hat{c}$) the learnt weight of the
bottleneck network (resp. rescaling) component and $\hat{u}\equiv
u_{\hat{\bm{v}}}$ (resp. $\hat{r}\equiv r_{\hat{c}}$). Note that generically,
$\hat{\bm{v}}\neq\hat{\bm{w}}$ and $\hat{c}\neq\hat{b}$, and therefore
$\hat{f}\neq\hat{u}+\hat{r}$, since $\hat{b},\hat{\bm{w}}$ result from their
joint optimization as parts of the full DAE $f_{b,\bm{w}}$, while $\hat{c}$
(or $\hat{\bm{v}}$) are optimized independently. As we discuss in Section III,
training the sole rescaling $r_{c}$ does not afford an expressive enough
denoiser, while an independently learnt bottleneck network component
$u_{\bm{v}}$ essentially only learns to implement PCA. However, when jointly
trained as components of the full DAE $f_{b,\bm{w}}$ (3), the resulting
denoiser $\hat{f}$ is a genuinely non-linear model which yields a much lower
test error than PCA, and learns to leverage flexibly its two components to
balance the preservation of the data nuances and the removal of the noise.

#### Learning metrics

The performance of the DAE (3) trained with the loss (5) is quantified by its
reconstruction (denoising) test MSE, defined as

$\displaystyle\mathrm{mse}_{\hat{f}}\equiv\mathbb{E}_{\mathcal{D}}\mathbb{E}_{\bm{x}\sim\mathbb{P}}\mathbb{E}_{\bm{\xi}\sim\mathcal{N}(0,\mathbb{I}_{d})}\left\lVert\bm{x}-f_{\hat{b},\hat{\bm{w}}}\left(\sqrt{1-\Delta}\bm{x}+\sqrt{\Delta}\bm{\xi}\right)\right\lVert^{2}.$
(6)

The expectations run over a fresh test sample $\bm{x}$ sampled from the
Gaussian mixture $\mathbb{P}$ (2), and a new additive noise $\bm{\xi}$
corrupting it. Note that an expectation over the train set $\mathcal{D}$ is
also included to make $\mathrm{mse}_{\hat{f}}$ a metric that does not depend
on the particular realization of the train set. The denoising test MSEs
$\mathrm{mse}_{\hat{u}},\mathrm{mse}_{\hat{r}}$ are defined similarly as the
denoising test errors of the independently learnt components (4). Aside from
the denoising MSE (6), another question of interest is how much the DAE
manages to learn the structure of the data distribution, as described by the
cluster means $\bm{\mu}_{k}$. This is measured by the cosine similarity matrix
$\bm{\theta}\in\mathbb{R}^{p\times K}$, where for $i\in\llbracket
1,p\rrbracket$ and $k\in\llbracket 1,K\rrbracket$,

$\displaystyle\theta_{ik}\equiv\mathbb{E}_{\mathcal{D}}\left[\frac{\hat{\bm{w}}_{i}^{\top}\bm{\mu}_{k}}{\lVert\hat{\bm{w}}_{i}\lVert\lVert\bm{\mu}_{k}\lVert}\right].$
(7)

In other words, $\theta_{ik}$ measures the alignment of the $i-$th row
$\hat{\bm{w}}_{i}$ of the trained weight matrix $\hat{\bm{w}}$ with the mean
of the $k-$th cluster $\bm{\mu}_{k}$.

#### High-dimensional limit

We analyze the optimization problem (5) in the high-dimensional limit where
the input dimension $d$ and number of training samples $n$ jointly tend to
infinity, while their ratio $\alpha=\nicefrac{{n}}{{d}}$ stays $\Theta(1)$.
The hidden layer width $p$, the noise level $\Delta$, the number of clusters
$K$ and the norm of the cluster means $\lVert\bm{\mu}_{k}\lVert$ are also
assumed to remain $\Theta(1)$. This corresponds to a rich limit, where the
number of parameters of the DAE is not large compared to the number of samples
like in [20, 13], and therefore cannot trivially fit the train set, or simply
memorize it [17]. Conversely, the number of samples $n$ is not infinite like
in [8, 9, 19], and therefore importantly allows to study the effect of a
finite train set on the representation learnt by the DAE.

##  II Asymptotic formulae for DAEs

We now state the main result of the present work, namely the closed-form
asymptotic formulae for the learning metrics $\mathrm{mse}_{\hat{f}}$ (6) and
$\theta$ (7) for a DAE (3) learnt with the empirical loss (5), derived using
the replica method in its replica-symmetric formulation (see Appendix A for
the derivation).

######  Assumption II.1.

The covariances $\\{\bm{\Sigma}\\}_{k=1}^{K}$ admit a common set of
eigenvectors $\\{\bm{e}_{i}\\}_{i=1}^{d}$. We further note
$\\{\lambda^{k}_{i}\\}_{i=1}^{d}$ the eigenvalues of $\bm{\Sigma}_{k}$. The
eigenvalues $\\{\lambda^{k}_{i}\\}_{i=1}^{d}$ and the projection of the
cluster means on the eigenvectors $\\{\bm{e}_{i}^{\top}\bm{\mu}_{k}\\}_{i,k}$
are assumed to admit a well-defined joint distribution $\nu$ as $d\to\infty$ –
namely, for $\gamma=(\gamma_{1},...,\gamma_{K})\in\mathbb{R}^{K}$ and
$\tau=(\tau_{1},...,\tau_{K})\in\mathbb{R}^{K}$:

$\displaystyle\frac{1}{d}\sum\limits_{i=1}^{d}\prod\limits_{k=1}^{K}\delta\left(\lambda_{i}^{k}-\gamma_{k}\right)\delta\left(\sqrt{d}\bm{e}_{i}^{\top}\bm{\mu}_{k}-\tau_{k}\right)\xrightarrow[]{d\to\infty}\nu\left(\gamma,\tau\right).$
(8)

Moreover, the marginals $\nu_{\gamma}$ (resp. $\nu_{\tau}$) are assumed to
have a well-defined first (resp. second) moment.

######  Assumption II.2.

$g(\cdot)$ is a $\ell_{2}$ regularizer with strength $\lambda$, i.e.
$g(\cdot)=\nicefrac{{\lambda}}{{2}}\lVert\cdot\lVert_{F}^{2}$.

######  Result II.3.

(Closed-form asymptotics for DAEs trained with empirical risk minimization)
Under Assumptions II.1 and II.2, in the high-dimensional limit $n,d\to\infty$
with fixed ratio $\alpha$, the denoising test MSE $\mathrm{mse}_{\hat{f}}$ (6)
admits the expression

$\displaystyle\mathrm{mse}_{\hat{f}}-\mathrm{mse}_{\circ}=$
$\displaystyle\sum\limits_{k=1}^{K}\rho_{k}\mathbb{E}_{z}\Tr[q\sigma\left({\scriptstyle{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}m_{k}+\sqrt{\Delta
q+{\tiny{\color[rgb]{0,0,0}(1-\Delta)}}q_{k}}z}\right)^{\otimes 2}]$ (9)
$\displaystyle-2\sum\limits_{k=1}^{K}\rho_{k}\mathbb{E}_{u,v}\left[{\scriptstyle\sigma\left({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}m_{k}+\sqrt{q_{k}({\tiny{\color[rgb]{0,0,0}1-\Delta}})}u+\sqrt{\Delta
q}v\right)^{\top}\left((1-\hat{b}{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}})(m_{k}+\sqrt{q_{k}}u)-\hat{b}\sqrt{\Delta
q}v\right)}\right]+o(1),$

where the averages bear over independent Gaussian variables
$z,u,v\sim\mathcal{N}(0,\mathbb{I}_{p})$. We denoted

$\displaystyle\mathrm{mse}_{\circ}=d\Delta\hat{b}^{2}+\left(1-{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}\hat{b}\right)^{2}\left[\sum\limits_{k=1}^{K}\rho_{k}\left(\int
d\nu_{\tau}(\tau)\tau_{k}^{2}+d\int
d\nu_{\gamma}(\gamma)\gamma_{k}\right)\right].$ (10)

The learnt skip connection strength $\hat{b}$ is

$\displaystyle\hat{b}=\frac{\left(\sum\limits_{k=1}^{K}\rho_{k}\int
d\nu_{\gamma}(\gamma)\gamma_{k}\right){\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}}{\left(\sum\limits_{k=1}^{K}\rho_{k}\int
d\nu_{\gamma}(\gamma)\gamma_{k}\right){\tiny{\color[rgb]{0,0,0}(1-\Delta)}}+\Delta}+o(1).$
(11)

The cosine similarity $\theta$ (7) admits the compact formula for
$i\in\llbracket 1,p\rrbracket$ and $k\in\llbracket 1,K\rrbracket$

$\displaystyle\theta_{ik}=\frac{(m_{k})_{i}}{\sqrt{q_{ii}\int
d\nu_{\tau}(\tau)\tau_{k}^{2}}},$ (12)

where we have introduced the summary statistics

$\displaystyle
q=\underset{d\to\infty}{\lim}\mathbb{E}_{\mathcal{D}}\left[\frac{\hat{\bm{w}}\hat{\bm{w}}^{\top}}{d}\right],$
$\displaystyle
q_{k}=\underset{d\to\infty}{\lim}\mathbb{E}_{\mathcal{D}}\left[\frac{\hat{\bm{w}}\bm{\Sigma}_{k}\hat{\bm{w}}^{\top}}{d}\right],$
$\displaystyle
m_{k}=\underset{d\to\infty}{\lim}\mathbb{E}_{\mathcal{D}}\left[\frac{\hat{\bm{w}}\bm{\mu}_{k}}{\sqrt{d}}\right].$
(13)

Thus $q,q_{k}\in\mathbb{R}^{p\times p},~{}m_{k}\in\mathbb{R}^{p}$. The summary
statistics $q,q_{k},m_{k}$ can be determined as solutions of the system of
equations

$\displaystyle\begin{cases}\hat{q}_{k}=\alpha\rho_{k}\mathbb{E}_{\xi,\eta}V_{k}^{-1}\left(\mathrm{prox}_{y}^{k}-q_{k}^{\frac{1}{2}}\eta-
m_{k}\right)^{\otimes 2}V_{k}^{-1}\\\
\hat{V}_{k}=-\alpha\rho_{k}q_{k}^{-\frac{1}{2}}\mathbb{E}_{\xi,\eta}V_{k}^{-1}\left(\mathrm{prox}_{y}^{k}-q_{k}^{\frac{1}{2}}\eta-
m_{k}\right)\eta^{\top}\\\
\hat{m}_{k}=\alpha\rho_{k}\mathbb{E}_{\xi,\eta}V_{k}^{-1}\left(\mathrm{prox}_{y}^{k}-q_{k}^{\frac{1}{2}}\eta-
m_{k}\right)\\\
\hat{q}=\frac{\alpha}{\Delta}\sum\limits_{k=1}^{K}\rho_{k}\mathbb{E}_{\xi,\eta}V^{-1}\left(\mathrm{prox}^{k}_{x}-\sqrt{\Delta}q^{\frac{1}{2}}\xi\right)^{\otimes
2}V^{-1}\\\
\hat{V}=-\alpha\sum\limits_{k=1}^{K}\rho_{k}\mathbb{E}_{\xi,\eta}\Bigg{[}\frac{1}{\sqrt{\Delta}}q^{-\frac{1}{2}}V^{-1}\left(\mathrm{prox}^{k}_{x}-\sqrt{\Delta}q^{\frac{1}{2}}\xi\right)\xi^{\top}-\sigma\left(\sqrt{1-\Delta}\mathrm{prox}_{y}^{k}+\mathrm{prox}_{x}^{k}\right)^{\otimes
2}\Bigg{]}\end{cases}$ $\displaystyle\begin{cases}q_{k}=\int
d\nu(\gamma,\tau)\gamma_{k}\left(\lambda\mathbb{I}_{p}+\hat{V}+\sum\limits_{j=1}^{K}\gamma_{j}\hat{V}_{j}\right)^{-2}\left(\hat{q}+\sum\limits_{j=1}^{K}\gamma_{j}\hat{q}_{j}+\sum\limits_{1\leq
j,l\leq K}\tau_{j}\tau_{l}\hat{m}_{j}\hat{m}_{l}^{\top}\right)\\\ V_{k}=\int
d\nu(\gamma,\tau)\gamma_{k}\left(\lambda\mathbb{I}_{p}+\hat{V}+\sum\limits_{j=1}^{K}\gamma_{j}\hat{V}_{j}\right)^{-1}\\\
m_{k}=\int
d\nu(\gamma,\tau)\tau_{k}\left(\lambda\mathbb{I}_{p}+\hat{V}+\sum\limits_{j=1}^{K}\gamma_{j}\hat{V}_{j}\right)^{-1}\sum\limits_{j=1}^{K}\tau_{j}\hat{m}_{j}\\\
q=\int
d\nu(\gamma,\tau)\left(\lambda\mathbb{I}_{p}+\hat{V}+\sum\limits_{j=1}^{K}\gamma_{j}\hat{V}_{j}\right)^{-2}\left(\hat{q}+\sum\limits_{j=1}^{K}\gamma_{j}\hat{q}_{j}+\sum\limits_{1\leq
j,l\leq K}\tau_{j}\tau_{l}\hat{m}_{j}\hat{m}_{l}^{\top}\right)\\\ V=\int
d\nu(\gamma,\tau)\left(\lambda\mathbb{I}_{p}+\hat{V}+\sum\limits_{j=1}^{K}\gamma_{j}\hat{V}_{j}\right)^{-1}\end{cases}$
(14)

In (II.3), $\hat{q}_{k},\hat{V}_{k},\hat{V},V\in\mathbb{R}^{p\times p}$ and
$\hat{m}_{k}\in\mathbb{R}^{p}$, and the averages bear over finite-dimensional
i.i.d Gaussians $\xi,\eta\sim\mathcal{N}(0,\mathbb{I}_{p})$. Finally,
$\mathrm{prox}^{k}_{x},\mathrm{prox}^{k}_{y}$ are given as the solutions of
the finite-dimensional optimization

$\displaystyle\mathrm{prox}^{k}_{x},~{}$
$\displaystyle\mathrm{prox}^{k}_{y}=\underset{x,y\in\mathbb{R}^{p}}{\mathrm{arginf}}\Bigg{\\{}\Tr[V_{k}^{-1}\left(y-q_{k}^{\frac{1}{2}}\eta-
m_{k}\right)^{\otimes
2}]+\frac{1}{\Delta}\Tr[V^{-1}\left(x-\sqrt{\Delta}q^{\frac{1}{2}}\xi\right)^{\otimes
2}]$
$\displaystyle+\Tr[q\sigma({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}y+x)^{\otimes
2}]-2\sigma({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}y+x)^{\top}((1-{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}\hat{b})y-\hat{b}x)\Bigg{\\}}.$
(15)

In fact, Assumptions II.1 and II.2 are not strictly necessary, and can be
simultaneously relaxed to address arbitrary convex regularizer $g(\cdot)$ and
generically non-commuting $\\{\bm{\Sigma}_{k}\\}_{k=1}^{K}$ – but at the price
of more intricate formulae. For this reason, we choose to discuss here Result
II.3, and defer a discussion and detailed derivation of the generic
asymptotically exact formulae to Appendix A, see eq. (61). Let us mention that
a sharp asymptotic characterization of the train MSE can also be derived; for
conciseness, we do not present it here and refer the interested reader to
equation (A.6) in Appendix A. Result II.3 encompasses as special cases the
asymptotic characterization of the components $\hat{r},\hat{u}$ (4):

######  Corollary II.4.

(MSE of components) The test MSE of $\hat{r}$ (4) is given by
$\mathrm{mse}_{\hat{r}}=\mathrm{mse}_{\circ}$ (10). Furthermore, the learnt
value of its single parameter $\hat{c}$ is given by (11). The test MSE, cosine
similarity and summary statistics of the bottleneck network $\hat{u}$ (4)
follow from Result II.3 by setting $\hat{b}=0$.

The implications of Corollary II.4 shall be further discussed in Section III,
and a full derivation is provided in Appendix E. Finally, remark that in the
noiseless limit $\Delta=0$, the denoising task reduces to a reconstruction
task, with the autoencoder being tasked with reproducing the clean data as an
output when taking the same clean sample as an input. Therefore Result II.3
also includes RAEs (by definition, without skip connection) as a special case.

######  Corollary II.5.

(RAEs) In the $n,d\to\infty$ limit, the MSE, cosine similarity and summary
statistics for an RAE follow from Result II.3 by setting $x=0$ in (II.3),
removing the first term in the brackets in the equation of $\hat{V}$ (II.3)
and taking the limit $\Delta,\hat{q},\hat{b}\to 0$.

Corollary II.5 will be the object of further discussion in Section III. A
detailed derivation is presented in Appendix F. Note that Corollary II.5
provides a characterization of RAEs as a function of the sample complexity
$\alpha$, where previous studies on non-linear RAEs rely on the assumption of
an infinite number of available training samples [13, 8, 9].

Equations (11) and (13) of Result II.3 thus characterize the statistics of the
learnt parameters $\hat{b},\hat{\bm{w}}$ of the trained DAE (3). These summary
statistics are, in turn, sufficient to fully characterize the learning metrics
(6) and (7) via equations (9) and (12). We thus have reduced the high-
dimensional optimization (5) and the high-dimensional average over the train
set $\mathcal{D}$ involved in the definition of the metrics (6) and (7) to a
simpler system of equations over $4+6K$ variables (II.3) which can be solved
numerically. It is important to note that all the summary statistics involved
in (II.3) are finite-dimensional as $d\to\infty$, and therefore Result II.3 is
a fully asymptotic characterization, in the sense that it does not involve any
high-dimensional object. In the next paragraphs, we give two examples of
applications of Result II.3, to a simple binary isotropic mixture, and to real
data sets.

Figure 1:
$\alpha=1,K=2,\rho_{1,2}=\nicefrac{{1}}{{2}},\bm{\Sigma}_{1,2}=0.09\times\mathbb{I}_{d},p=1,\lambda=0.1,\sigma(\cdot)=\tanh(\cdot)$;
the cluster mean $\bm{\mu}_{1}=-\bm{\mu}_{2}$ was taken as a random Gaussian
vector of norm $1$. (left) In blue, the difference in MSE between the full DAE
$\hat{f}$ (3) and the rescaling component $\hat{r}$ (4). Solid lines
correspond to the sharp asymptotic characterization of Result II.3. Dots
represent numerical simulations for $d=700$, training the DAE using the
Pytorch implementation of full-batch Adam, with learning rate $\eta=0.05$ over
$2000$ epochs, averaged over $N=10$ instances. Error bars represent one
standard deviation. For completeness, the MSE of the oracle denoiser is given
as a baseline in green, see Section III. The performance of a linear DAE
($\sigma(x)=x$) is represented in dashed red. (right) Cosine similarity
$\theta$ (7) (green), squared weight norm
$\nicefrac{{\lVert\hat{\bm{w}}\lVert^{2}_{F}}}{{d}}$ (red) and skip connection
strength $\hat{b}$ (blue). Solid lines correspond to the formulae (12)(13) and
(11) of Result II.3; dots are numerical simulations. For completeness, the
cosine similarity of the first principal component of the clean train data
$\\{\bm{x}^{\mu}\\}_{\mu=1}^{n}$ is plotted in dashed black.

#### Example 1: Isotropic homoscedastic mixture

We give as a first example the case of a synthetic binary Gaussian mixture
with
$K=2,\bm{\mu}_{1}=-\bm{\mu}_{2},\bm{\Sigma}_{1,2}=0.09\times\mathbb{I}_{d},\rho_{1,2}=\nicefrac{{1}}{{2}}$,
using a DAE with $\sigma=\tanh$ and $p=1$. Since this simple case exhibits the
key phenomenology discussed in the present work, we refer to it in future
discussions. The MSE $\mathrm{mse}_{\hat{f}}$ (9) evaluated from the solutions
of the self-consistent equations (II.3) is plotted as the solid blue line in
Fig. 1 (left) and compared to numerical simulations corresponding to training
the DAE (3) with the Pytorch implementation of the Adam optimizer [43] (blue
dots), for sample complexity $\alpha=1$ and $\ell_{2}$ regularization (weight
decay) $\lambda=0.1$. The agreement between the theory and simulation is
compelling. The green solid line and corresponding green dots in Fig. 1
(right) correspond to the replica prediction (12) and simulations for the
cosine similarity $\theta$ (7), and again display very good agreement.

A particularly striking observation is that due to the non-convexity of the
loss (5), there is a priori no guarantee that an Adam-optimized DAE should
find a global minimum, as described by the Result II.3, rather than a local
minimum. The compelling agreement between theory and simulations in Fig. 1
temptingly suggests that the loss landscape of DAEs (3) trained with the loss
(5) for the data model (2) should in some way be benign. Authors of [12] have
shown, for linear RAEs, that there exists a unique global and local minimum
for the square loss and no regularizer. Ref. [14] offers further insight for a
linear DAE in dimension $d=1$, and shows that, aside from the global minima,
the loss landscape only includes an unstable saddle point from which the
dynamics easily escapes. Extending these works and intuitions to non-linear
DAEs is an exciting research topic for future work.

---

Figure 2: Difference in MSE between the full DAE (3) and the rescaling
component (4) for the MNIST data set (middle), of which for simplicity only
$1$s and $7$s were kept, and FashionMNIST (right), of which only boots and
shoes were kept. In blue, the theoretical predictions resulting from using
Result II.3 with the empirically estimated covariances and means, see Appendix
D for further details. In red, numerical simulations of a DAE ($p=1$,
$\sigma=\tanh$) trained with $n=784$ training points, using the Pytorch
implementation of full-batch Adam, with learning rate $\eta=0.05$ and weight
decay $\lambda=0.1$ over $2000$ epochs, averaged over $N=10$ instances. Error
bars represent one standard deviation. (left) illustration of the denoised
images: (top left) original image, (top right) noisy image, (bottom left) DAE
$\hat{f}$ (3), (bottom right) rescaling $\hat{r}$ (4).

#### Example 2: MNIST, FashionMNIST

It is reasonable to ask whether Result II.3 is restricted to Gaussian mixtures
(2). The answer is negative – in fact, Result II.3 also describes well a
number of real data distributions. We provide such an example for FashionMNIST
[16] (from which, for simplicity, we only kept boots and shoes) and MNIST [15]
(1s and 7s), in Fig. 2. For each data set, samples sharing the same label were
considered to belong to the same cluster. The mean and covariance thereof were
estimated numerically, and combined with Result II.3. The resulting denoising
MSE predictions $\mathrm{mse}_{\hat{f}}$ are plotted as solid lines in Fig. 2,
and agree very well with numerical simulations of DAEs optimized over the real
data sets using the Pytorch implementation of Adam [43]. A full description of
this experiment is given in Appendix D.

The observation that the MSEs of real data sets are to such degree of accuracy
captured by the equivalent Gaussian mixture strongly hints at the presence of
Gaussian universality [44]. This opens a gateway to future research, as
Gaussian universality has hitherto been exclusively addressed in
classification and regression (rather than denoising) settings, see e.g. [44,
45, 46]. Denoising tasks further constitute a particularly intriguing setting
for universality results, as Gaussian universality would signify that only
second-order statistics of the data can be reconstructed using a shallow
autoencoder.

##  III The role and importance of the skip connection.

Result II.3 for the full DAE $\hat{f}$ (3) and Corollary II.4 for its
components $\hat{r},~{}\hat{u}$ (4) allow to disentangle the contribution of
each part, and thus to pinpoint their respective roles in the DAE
architecture. We sequentially present a comparison of $\hat{f}$ with
$\hat{r}$, and $\hat{f}$ with $\hat{u}$. We remind that $\hat{f}$, $\hat{r}$
and $\hat{u}$ result from independent optimizations over the same train set
$\mathcal{D}$, and that while $f_{b,\bm{w}}=u_{\bm{w}}+h_{b}$,
$\hat{f}\neq\hat{u}+\hat{r}$.

#### Full DAE and the rescaling component

We start this section by observing that for noise levels $\Delta$ below a
certain threshold, the full DAE $\hat{f}$ yields better MSE than the learnt
rescaling $\hat{r}$, as can be seen by the negative value of
$\mathrm{mse}_{\hat{f}}-\mathrm{mse}_{\hat{r}}$ in Fig. 1 and Fig. 2. The
improvement is more sizeable at intermediate noise levels $\Delta$, and is
observed for a growing region of $\Delta$ as the sample complexity $\alpha$
increases, see Fig. 3 (a). This lower MSE further translates into visible
qualitative changes in the result of denoising. As can be seen from Fig. 2
(left), the full DAE $\hat{f}$ (3) (bottom left) yields denoised images with
sensibly higher definition and overall contrast, while a simple rescaling
$\hat{r}$ (bottom right) leads to a still largely blurred image.

We provide one more comparison: for the isotropic binary mixture (see Fig. 1),
the DAE test error $\mathrm{mse}_{\hat{f}}$ in fact approaches the
information-theoretic lowest achievable MSE $\mathrm{mse}^{\star}$ as the
sample complexity $\alpha$ increases. To see this, note that
$\mathrm{mse}^{\star}$ is given by the application of Tweedie’s formula [42],
that requires perfect knowledge of the cluster means $\bm{\mu}_{k}$ and
covariances $\bm{\Sigma}_{k}$ – it is, therefore, an oracle denoiser.
A sharp asymptotic characterization of the oracle denoiser is provided in
Appendix B. As can be observed from Fig. 3 (a), the DAE MSE (3) approaches the
oracle test error $\mathrm{mse}^{\star}$ as the number of available training
samples $n$ grows, and is already sensibly close to the optimal value for
$\alpha=8$.

Figure 3: (left) Solid lines: difference in MSE between the full DAE $\hat{f}$
(3), with $\sigma=\tanh$, $p=1$, and the rescaling $\hat{r}$ (4). Dashed: the
same curve for the oracle denoiser. Different colours represent different
sample complexities $\alpha$ (solid lines). (right) Difference in MSE between
the bottleneck network $\hat{u}$ (4) and the complete DAE $\hat{f}$ (3). In
blue, the theoretical prediction (16); in red, numerical simulations for the
bottleneck network (4) ($\sigma=\tanh$, $p=1$) trained with the Pytorch
implementation of full-batch Adam, with learning rate $\eta=0.05$ and weight
decay $\lambda=0.1$ over $2000$ epochs, averaged over $N=5$ instances, for
$d=700$. In green, the MSE (minus the MSE of the complete DAE (3)) achieved by
PCA. Error bars represent one standard deviation. The model and parameters are
the same as in Fig. 1.

(a) (b) (c) (d) (e) (f)

Figure 4: Illustration of the denoised image for the various networks and
algorithms. (a) original image (b) noisy image, for $\sqrt{\Delta}=0.2$ (c)
trained rescaling $\hat{r}$ (4) (d) full DAE $\hat{f}$ (3) (e) bottleneck
network $\hat{u}$ (4) (f) PCA. The DAE and training parameters are the same as
Fig. 2, see also Appendix D.

#### DAEs with(out) skip connection

We now turn our attention to comparing the full DAE $\hat{f}$ (3) to the
bottleneck network component $\hat{u}$ (4). It follows from Result II.3 and
Corollary II.4 that $\hat{u}$ (4) leads to a higher MSE than the full DAE
$\hat{f}$ (3), with the gap being $\Theta(d)$. More precisely,

$\displaystyle\frac{1}{d}\left(\mathrm{mse}_{\hat{u}}-\mathrm{mse}_{\hat{f}}\right)=\frac{\left(\int
d\nu_{\gamma}(\gamma)\sum\limits_{k=1}^{K}\rho_{k}\gamma_{k}\right)^{2}(1-\Delta)}{\left(\int
d\nu_{\gamma}(\gamma)\sum\limits_{k=1}^{K}\rho_{k}\gamma_{k}\right){\tiny{\color[rgb]{0,0,0}(1-\Delta)}}+\Delta}.$
(16)

The theoretical prediction (16) compares excellently with numerical
simulations; see Fig. 3 (right). Strikingly, we find that PCA denoising yields
an MSE almost indistinguishable from $\hat{u}$, see Fig. 3, strongly
suggesting that $\hat{u}$ essentially learns, also in the denoising setting,
to project the noisy data $\tilde{\bm{x}}$ along the principal components of
the training set. The last two images of Fig. 4 respectively correspond to
$\hat{u}$ and PCA, which can indeed be observed to lead to visually near-
identical results. This echoes the findings of [10, 11, 8, 9, 13] in the case
of RAEs that bottleneck networks are limited by the PCA reconstruction
performance – a conclusion that we also recover from Corollary II.5, see
Appendix F. Crucially however, it also means that compared to the full DAE
$\hat{f}$ (3), PCA is sizeably suboptimal, since
$\mathrm{mse}_{\mathrm{PCA}}\approx\mathrm{mse}_{\hat{u}}=\mathrm{mse}_{\hat{f}}+\Theta(d)$.

This last observation has an important consequence: in contrast to previously
studied RAEs [10, 12, 11, 9, 13], the full DAE $\hat{f}$ does not simply learn
to perform PCA. In contrast to bottleneck RAE networks [8, 9, 13], the non-
linear DAE hence does not reduce to a linear model after training. The non-
linearity is important to improve the denoising MSE, see Fig. 1. We stress
this finding: trained alone, the bottleneck network $\hat{u}$ only learns to
perform PCA; trained jointly with the rescaling component as part of the full
DAE $f_{b,\bm{w}}$ (3), it learns a richer, non-linear representation. The
full DAE (3) thus offers a genuinely non-linear learning model and opens
exciting research avenues for the theory of autoencoders, beyond linear (or
effectively linear) cases. In the next paragraph, we explore further the
interaction between the rescaling component and the bottleneck network.

#### A tradeoff between the rescaling and the bottleneck network

Result (II.3), alongside Corollary (II.4) and the discussion in Section III
provide a firm theoretical basis for the well-known empirical intuition
(discussed e.g. in [38]) that skip-connections allow to better propagate
information from the input to the output of the DAE, thereby contributing to
preserving intrinsic characteristics of the input. This effect is clearly
illustrated in Fig. 4, where the resulting denoised image of an MNIST $7$ by
$\hat{r}$, $\hat{f}$, $\hat{u}$ and PCA are presented. While the bottleneck
network $\hat{u}$ perfectly eliminates the background noise and produces an
image with a very good resolution, it essentially collapses the image to the
cluster mean, and yields, like PCA, the average MNIST $7$. As a consequence,
the denoised image bears little resemblance with the original image – in
particular, the horizontal bar of the $7$ is lost in the process. Conversely,
the rescaling $\hat{r}$ preserves the nuances of the original image, but the
result is still largely blurred and displays overall poor contrast. Finally,
the complete DAE (3) manages to preserve the characteristic features of the
original data, while enhancing the image resolution by slightly overlaying the
average $7$ thereupon.

The optimization of the DAE (5) is therefore described by a tradeoff between
two competing effects – namely the preservation of the input nuances by the
skip connection, and the enhancement of the resolution/noise removal by the
bottleneck network. This allows us to discuss the curious non-monotonicity of
the cosine similarity $\theta$ as a function of the noise level $\Delta$, see
Fig. 1 (left). While it may at first seem curious that the DAE seemingly does
not manage to learn the data structure better for low $\Delta$ than for
intermediate $\Delta$ (where the cosine similarity $\theta$ is observed to be
higher), this is actually due to the afore-dicussed tradeoff. Indeed, for
small $\Delta$, the data is still substantially clean, and there is therefore
no incentive to enhance the contrast by using the cluster means –which are
consequently not learnt. This phase is thus characterized by a large skip
connection strength $\hat{b}$, and small cosine similarity $\theta$ and weight
norm $\lVert\hat{\bm{w}}\lVert_{F}$. Conversely, at high noise levels
$\Delta$, the nuances of the data are already lost because of the noise. Hence
the DAE does not rely on the skip connection component (whence the small
values of $\hat{b}$), and the only way to produce reasonably denoised data is
to collapse to the cluster mean using the network component (whence a large
$\lVert\hat{\bm{w}}\lVert_{F}$).

## Conclusion

We consider the problem of denoising a high-dimensional Gaussian mixture, by
training a DAE via empirical risk minimization, in the limit where the number
of training samples and the dimension are proportionally large. We provide a
sharp asymptotic characterization of a number of summary statistics of the
trained DAE weight, average MSE, and cosine similarity with the cluster means.
These results contain as a corollary the case of RAEs. Building on these
findings, we isolate the role of the skip connection and the bottleneck
network in the DAE architecture and characterize the tradeoff between those
two components in terms of preservation of the data nuances and noise removal
– thereby providing some theoretical insight into a longstanding practical
intuition in machine learning.

We believe the present work also opens exciting research avenues. First, our
real data experiments hint at the presence of Gaussian universality. While
this topic has gathered considerable attention in recent years, only
classification/regression supervised learning tasks have been hitherto
addressed. Which aspects of universality carry over to denoising tasks, and
how they differ from the current understanding of supervised
regression/classification, is an important question. Second, the DAE with skip
connection (3) provides an autoencoder model which does not just simply learn
the principal components of the training set. It, therefore, affords a
genuinely non-linear network model where richer learning settings can be
investigated.

### Acknowledgements

We thank Eric Vanden-Eijnden for his wonderful lecture at the les Houches
Summer School of July 2022, which inspired parts of this work. We thank Maria
Refinetti and Sebastian Goldt for insightful discussions at the very early
stages of this project.

## References

  * [1] Yang Song, Jascha Narain Sohl-Dickstein, Diederik P. Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole.  Score-based generative modeling through stochastic differential equations.  International Conference on Learning Representations, 2021.
  * [2] Jonathan Ho, Ajay Jain, and Pieter Abbeel.  Denoising diffusion probabilistic models.  Advances in Neural Information Processing Systems, 33:6840–6851, 2020.
  * [3] Pascal Vincent, H. Larochelle, Isabelle Lajoie, Yoshua Bengio, and Pierre-Antoine Manzagol.  Stacked denoising autoencoders: Learning useful representations in a deep network with a local denoising criterion.  J. Mach. Learn. Res., 11:3371–3408, 2010.
  * [4] Reza Oftadeh, Jiayi Shen, Zhangyang Wang, and Dylan A. Shell.  Eliminating the invariance on the loss landscape of linear autoencoders.  In International Conference on Machine Learning, 2020.
  * [5] Daniel Kunin, Jonathan M. Bloom, Aleksandrina Goeva, and Cotton Seed.  Loss landscapes of regularized linear autoencoders.  In International Conference on Machine Learning, 2019.
  * [6] Xuchan Bao, James Lucas, Sushant Sachdeva, and Roger B Grosse.  Regularized linear autoencoders recover the principal components, eventually.  Advances in Neural Information Processing Systems, 33:6971–6981, 2020.
  * [7] Gauthier Gidel, Francis R. Bach, and Simon Lacoste-Julien.  Implicit regularization of discrete gradient dynamics in deep linear neural networks.  In Neural Information Processing Systems, 2019.
  * [8] Maria Refinetti and Sebastian Goldt.  The dynamics of representation learning in shallow, non-linear autoencoders.  In International Conference on Machine Learning, pages 18499–18519. PMLR, 2022.
  * [9] A. E. Shevchenko, Kevin Kögler, Hamed Hassani, and Marco Mondelli.  Fundamental limits of two-layer autoencoders, and achieving them with gradient methods.  ArXiv, abs/2212.13468, 2022.
  * [10] Carl Eckart and G. Marion Young.  The approximation of one matrix by another of lower rank.  Psychometrika, 1:211–218, 1936.
  * [11] Hervé Bourlard and Yves Kamp.  Auto-association by multilayer perceptrons and singular value decomposition.  Biological Cybernetics, 59:291–294, 1988.
  * [12] Pierre Baldi and Kurt Hornik.  Neural networks and principal component analysis: Learning from examples without local minima.  Neural Networks, 2:53–58, 1989.
  * [13] Phan-Minh Nguyen.  Analysis of feature learning in weight-tied autoencoders via the mean field lens.  ArXiv, abs/2102.08373, 2021.
  * [14] Arnu Pretorius, Steve Kroon, and Herman Kamper.  Learning dynamics of linear denoising autoencoders.  In International Conference on Machine Learning, 2018.
  * [15] Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner.  Gradient-based learning applied to document recognition.  Proc. IEEE, 86:2278–2324, 1998.
  * [16] Han Xiao, Kashif Rasul, and Roland Vollgraf.  Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms.  ArXiv, abs/1708.07747, 2017.
  * [17] Adityanarayanan Radhakrishnan, Mikhail Belkin, and Caroline Uhler.  Overparameterized neural networks implement associative memory.  Proceedings of the National Academy of Sciences of the United States of America, 117:27162 – 27170, 2019.
  * [18] Saachi Jain, Adityanarayanan Radhakrishnan, and Caroline Uhler.  A mechanism for producing aligned latent spaces with autoencoders.  ArXiv, abs/2106.15456, 2021.
  * [19] Thanh Van Nguyen, Raymond K. W. Wong, and Chinmay Hegde.  On the dynamics of gradient descent for autoencoders.  In International Conference on Artificial Intelligence and Statistics, 2019.
  * [20] Thanh Van Nguyen, Raymond K. W. Wong, and Chinmay Hegde.  Benefits of jointly training autoencoders: An improved neural tangent kernel analysis.  IEEE Transactions on Information Theory, 67:4669–4692, 2019.
  * [21] Giorgio Parisi.  Towards a mean field theory for spin glasses.  Phys. Lett, 73(A):203–205, 1979.
  * [22] Giorgio Parisi.  Order parameter for spin glasses.  Phys. Rev. Lett, 50:1946–1948, 1983.
  * [23] Lenka Zdeborová and Florent Krzakala.  Statistical physics of inference: thresholds and algorithms.  Advances in Physics, 65:453 – 552, 2015.
  * [24] Marylou Gabrié.  Mean-field inference methods for neural networks.  Journal of Physics A: Mathematical and Theoretical, 53, 2019.
  * [25] Elizabeth Gardner and Bernard Derrida.  Optimal storage properties of neural network models.  Journal of Physics A: Mathematical and general, 21(1):271, 1988\.
  * [26] Manfred Opper and David Haussler.  Generalization performance of bayes optimal classification algorithm for learning a perceptron.  Physical Review Letters, 66(20):2677, 1991.
  * [27] Jean Barbier, Florent Krzakala, Nicolas Macris, Léo Miolane, and Lenka Zdeborová.  Optimal errors and phase transitions in high-dimensional generalized linear models.  Proceedings of the National Academy of Sciences of the United States of America, 116:5451 – 5460, 2017.
  * [28] Benjamin Aubin, Florent Krzakala, Yue M. Lu, and Lenka Zdeborová.  Generalization error in high-dimensional perceptrons: Approaching bayes error with convex optimization.  Advances in Neural Information Processing Systems, 33:12199–12210, 2020.
  * [29] Benjamin Aubin, Antoine Maillard, Jean Barbier, Florent Krzakala, Nicolas Macris, and Lenka Zdeborová.  The committee machine: computational to statistical gaps in learning a two-layers neural network.  Journal of Statistical Mechanics: Theory and Experiment, 2019, 2018\.
  * [30] Rainer Dietrich, Manfred Opper, and Haim Sompolinsky.  Statistical mechanics of support vector networks.  Physical review letters, 82(14):2975, 1999.
  * [31] Blake Bordelon, Abdulkadir Canatar, and Cengiz Pehlevan.  Spectrum dependent learning curves in kernel regression and wide neural networks.  In International Conference on Machine Learning, 2020.
  * [32] Federica Gerace, Bruno Loureiro, Florent Krzakala, Marc Mézard, and Lenka Zdeborová.  Generalisation error in learning with random features and the hidden manifold model.  Journal of Statistical Mechanics: Theory and Experiment, 2021, 2020\.
  * [33] Hugo Cui, Bruno Loureiro, Florent Krzakala, and Lenka Zdeborová.  Error rates for kernel classification under source and capacity conditions.  ArXiv, abs/2201.12655, 2022.
  * [34] Christos Thrampoulidis, Ehsan Abbasi, and Babak Hassibi.  Precise error analysis of regularized $m$-estimators in high dimensions.  IEEE Transactions on Information Theory, 64(8):5592–5628, 2018\.
  * [35] Jacob A. Zavatone-Veth, William Tong, and Cengiz Pehlevan.  Contrasting random and learned features in deep bayesian linear regression.  Physical review. E, 105 6-1:064118, 2022.
  * [36] Hugo Cui, Florent Krzakala, and Lenka Zdeborová.  Optimal learning of deep random networks of extensive-width.  ArXiv, abs/2302.00375, 2023.
  * [37] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18, pages 234–241. Springer, 2015.
  * [38] Xiao-Jiao Mao, Chunhua Shen, and Yubin Yang.  Image restoration using convolutional auto-encoders with symmetric skip connections.  ArXiv, abs/1606.08921, 2016.
  * [39] Tong Tong, Gen Li, Xiejie Liu, and Qinquan Gao.  Image super-resolution using dense skip connections.  In Proceedings of the IEEE international conference on computer vision, pages 4799–4807, 2017.
  * [40] Jiwon Kim, Jung Kwon Lee, and Kyoung Mu Lee.  Accurate image super-resolution using very deep convolutional networks.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 1646–1654, 2016.
  * [41] Jiwon Kim, Jung Kwon Lee, and Kyoung Mu Lee.  Deeply-recursive convolutional network for image super-resolution.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 1637–1645, 2016.
  * [42] Bradley Efron.  Tweedie’s formula and selection bias.  Journal of the American Statistical Association, 106:1602 – 1614, 2011.
  * [43] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  CoRR, abs/1412.6980, 2014.
  * [44] Sebastian Goldt, Bruno Loureiro, Galen Reeves, Florent Krzakala, Marc Mézard, and Lenka Zdeborová.  The gaussian equivalence of generative models for learning with shallow neural networks.  In Mathematical and Scientific Machine Learning, 2020.
  * [45] Hong Hu and Yue M. Lu.  Universality laws for high-dimensional learning with random features.  IEEE Transactions on Information Theory, 69:1932–1964, 2020.
  * [46] Andrea Montanari and Basil N Saeed.  Universality of empirical risk minimization.  In Conference on Learning Theory, pages 4310–4312. PMLR, 2022.
  * [47] Bruno Loureiro, Cédric Gerbelot, Hugo Cui, Sebastian Goldt, Florent Krzakala, Marc Mézard, and Lenka Zdeborová.  Learning curves of generic features maps for realistic datasets with a teacher-student model.  Journal of Statistical Mechanics: Theory and Experiment, 2022, 2021\.
  * [48] Bruno Loureiro, Gabriele Sicuro, Cédric Gerbelot, Alessandro Pacco, Florent Krzakala, and Lenka Zdeborová.  Learning gaussian mixtures with generalised linear models: Precise asymptotics in high-dimensions.  In Neural Information Processing Systems, 2021.

##  Appendix A Derivation of Result II.3

In this section, we detail the derivation of Result (II.3).

###  A.1 Derivation technique

In order to access sharp asymptotic characterizations for the metrics and
summary statistics $\mathrm{mse}$ (6),$\theta$ (7), $\hat{b}$ and
$\nicefrac{{\lVert\hat{\bm{w}}\lVert}}{{d}}$, the first step is to observe
that any test function $\varphi(\hat{\bm{w}},\hat{b})$ of the learnt network
parameters $\hat{\bm{w}},\hat{b}$ can be written as an average over a limit
probability distribution

$\displaystyle\mathbb{E}_{\mathcal{D}}\left[\varphi(\hat{\bm{w}},\hat{b})\right]=\underset{\beta\to\infty}{\lim}\mathbb{E}_{\mathcal{D}}\frac{1}{Z}\int
dwdb\varphi(\bm{w},b)e^{-\beta\mathcal{\hat{R}}(\bm{w},b)}$ (17)

where $\mathcal{\hat{R}}(\cdot,\cdot)$ is the empirical risk (5) and the
corresponding probability density

$\displaystyle\mathbb{P}_{\beta}(\bm{w},b)=\frac{e^{-\beta\mathcal{\hat{R}}(\bm{w},b)}}{Z}$
(18)

is known as the Boltzmann measure in statistical physics, with the parameter
$\beta$ typically referred to as the inverse temperature. The normalization
$Z$ is called the partition function of $\mathbb{P}_{\beta}$:

$\displaystyle Z=\int d\bm{w}dbe^{-\beta\mathcal{\hat{R}}(\bm{w},b)}.$ (19)

The idea of the replica method [21, 22, 23, 24], building on the expression
(17), is to compute the moment generating function (also known as the free
entropy) as

$\displaystyle\mathbb{E}_{\mathcal{D}}\ln Z=\underset{s\to
0}{\lim}\frac{\ln\left[\mathbb{E}_{\mathcal{D}}Z^{s}\right]}{s}.$ (20)

With the moment generating function, it then becomes possible to compute
averages like (17). The backbone of the derivation, which we detail in the
following subsections, therefore lies in the computation of the replicated
partition function $\mathbb{E}_{\mathcal{D}}Z^{s}$.

###  A.2 Replicated partition function

For subsequent clarity, we shall introduce the variable

$\displaystyle\bm{\eta}^{\mu}\equiv\bm{x}^{\mu}-\bm{\mu}_{k}$ (21)

if $\bm{x}^{\mu}$ belongs to the $k-$th cluster. By definition of the Gaussian
mixture model, this is just a Gaussian variable:
$\bm{\eta}^{\mu}\sim\mathcal{N}(0,\bm{\Sigma}_{k})$, and the average over the
data set $\\{\bm{x}^{\mu}\\}$ reduces to averaging over the cluster index $k$
and $\\{\bm{\eta}^{\mu}\\}$. For notational brevity, we shall name in this
appendix ${\bm{\xi}}^{\mu}$ what is called $\sqrt{\Delta}{\bm{\xi}}^{\mu}$ in
the main text (i.e., absorbing the variance in the definition of
${\bm{\xi}}^{\mu}$). The replicated partition function then reads

$\displaystyle\mathbb{E}_{\mathcal{D}}Z^{s}=$
$\displaystyle\int\prod\limits_{a=1}^{s}d\bm{w}_{a}db_{a}e^{-\beta\sum\limits_{a=1}^{s}g(\bm{w}^{a})}$
$\displaystyle\times\prod\limits_{\mu=1}^{n}\sum\limits_{k=1}^{K}\rho_{k}\mathbb{E}_{\bm{\eta},{\bm{\xi}}}\underbrace{e^{-\frac{\beta}{2}\sum\limits_{a=1}^{s}\left\lVert\bm{\mu}_{k}+\bm{\eta}-\left[b_{a}\times({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}(\bm{\mu}_{k}+\bm{\eta})+{\bm{\xi}})+\frac{\bm{w}_{a}^{\top}}{\sqrt{d}}\sigma\left(\frac{\bm{w}_{a}({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}(\bm{\mu}_{k}+\bm{\eta})+{\bm{\xi}})}{\sqrt{d}}\right)\right]\right\rVert^{2}}}_{(\star)}.$
(22)

One can expand the exponent $(\star)$as

$\displaystyle
e^{{\scriptscriptstyle-\frac{\beta}{2}\sum\limits_{a=1}^{s}\left[\Tr[\frac{\bm{w}_{a}\bm{w}_{a}^{\top}}{d}\sigma\left(\frac{\bm{w}_{a}({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}(\bm{\mu}_{k}+\bm{\eta})+{\bm{\xi}})}{\sqrt{d}}\right)^{\otimes
2}]-2\sigma\left(\frac{\bm{w}_{a}({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}(\bm{\mu}_{k}+\bm{\eta})+{\bm{\xi}})}{\sqrt{d}}\right)^{\top}\frac{\bm{w}_{a}((1-b_{a}{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}})(\bm{\mu}_{k}+\bm{\eta})-b_{a}{\bm{\xi}})}{\sqrt{d}}\right]}}$
$\displaystyle\times
e^{{\scriptscriptstyle-\frac{\beta}{2}\sum\limits_{a=1}^{s}\left[(1-{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b_{a})^{2}({\color[rgb]{0,0,0}\lVert\bm{\mu}_{k}\lVert^{2}}+\lVert\bm{\eta}\rVert^{2}+2\bm{\mu}_{k}^{\top}\bm{\eta})+b_{a}^{2}\lVert{\bm{\xi}}\rVert^{2}+2b_{a}({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b_{a}-1)(\bm{\mu}_{k}^{\top}{\bm{\xi}}+\bm{\eta}^{\top}{\bm{\xi}})\right]}}.$
(23)

Therefore

$\displaystyle\mathbb{E}_{s,\bm{\eta},{\bm{\xi}}}(\star)=\sum\limits_{k=1}^{K}\rho_{k}e^{-\frac{\beta}{2}{\color[rgb]{0,0,0}\lVert\bm{\mu}_{k}\lVert^{2}}\sum\limits_{a=1}^{s}({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b_{a}-1)^{2}}\int\frac{d\bm{\eta}d{\bm{\xi}}}{(2\pi\sqrt{\Delta})^{d}\sqrt{\det\bm{\Sigma}_{k}}}$
$\displaystyle\underbrace{e^{{\scriptscriptstyle-\frac{1}{2}\bm{\eta}^{\top}\left(\bm{\Sigma}_{k}^{-1}+\beta\sum\limits_{a=1}^{s}({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b_{a}-1)^{2}\mathbb{I}_{d}\right)\bm{\eta}-\frac{1}{2}\lVert{\bm{\xi}}\rVert^{2}\left(\frac{1}{\Delta}+\beta\sum\limits_{a=1}^{s}b_{a}^{2}\right)-\beta\sum\limits_{a=1}^{s}\left((1-{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b_{a})^{2}\bm{\mu}_{k}^{\top}\bm{\eta}+b_{a}({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b_{a}-1)\bm{\mu}_{k}^{\top}{\bm{\xi}}+b_{a}({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b_{a}-1)\bm{\eta}^{\top}{\bm{\xi}}\right)}}}_{P_{\mathrm{eff.}}(\bm{\eta},{\bm{\xi}})}$
$\displaystyle\times
e^{{\scriptscriptstyle-\frac{\beta}{2}\sum\limits_{a=1}^{s}\left[\Tr[\frac{\bm{w}_{a}\bm{w}_{a}^{\top}}{d}\sigma\left(\frac{\bm{w}_{a}({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}(\bm{\mu}_{k}+\bm{\eta})+{\bm{\xi}})}{\sqrt{d}}\right)^{\otimes
2}]-2\sigma\left(\frac{\bm{w}_{a}({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}(\bm{\mu}_{k}+\bm{\eta})+{\bm{\xi}})}{\sqrt{d}}\right)^{\top}\frac{\bm{w}_{a}((1-b_{a}{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}})(\bm{\mu}_{k}+\bm{\eta})-b_{a}{\bm{\xi}})}{\sqrt{d}}\right]}}.$
(24)

The effective prior over the noises $\bm{\eta},{\bm{\xi}}$ is Gaussian with
means $\bm{\mu}_{\xi},\bm{\mu}_{\eta}$ and covariance $C$:

$\displaystyle
P_{\mathrm{eff.}}(\bm{\eta},{\bm{\xi}})=\mathcal{N}\left(\bm{\eta},{\bm{\xi}};\bm{\mu}_{\eta},\bm{\mu}_{\xi},\left[\begin{array}[]{cc}\bm{C}^{-1}_{\eta\eta}&C^{-1}_{\xi\eta}\mathbb{I}_{d}\\\
C^{-1}_{\xi\eta}\mathbb{I}_{d}&C^{-1}_{\xi\xi}\mathbb{I}_{d}\end{array}\right]^{-1}\right).$
(27)

Identifying the covariance components leads to

$\displaystyle
C^{-1}_{\xi\xi}=\frac{1}{\Delta}+\beta\sum\limits_{a=1}^{s}b_{a}^{2},$
$\displaystyle\bm{C}^{-1}_{\eta\eta}=\bm{\Sigma}_{k}^{-1}+\beta\sum\limits_{a=1}^{s}({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b_{a}-1)^{2}\mathbb{I}_{d},$
$\displaystyle
C^{-1}_{\xi\eta}=\beta\sum\limits_{a=1}^{s}b_{a}({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b_{a}-1).$
(28)

Note that all the sums indexed by $a$ introduce terms of $\mathcal{O}(s)$,
which can safely be neglected. Identifying the means leads to the following
system of equations:

$\displaystyle\begin{cases}C^{-1}_{\xi\xi}\bm{\mu}_{\xi}+C^{-1}_{\xi\eta}\bm{\mu}_{\eta}=\beta\sum\limits_{a=1}^{s}b_{a}(1-{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b_{a})\bm{\mu}_{k}\\\
\bm{C}^{-1}_{\eta\eta}\bm{\mu}_{\eta}+C^{-1}_{\xi\eta}\bm{\mu}_{\xi}=-\beta\sum\limits_{a=1}^{s}(1-{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b_{a})^{2}\bm{\mu}_{k}\end{cases},$
(29)

which is solved as

$\displaystyle\bm{\mu}_{\xi}=\underbrace{\beta\left(\bm{C}^{-1}_{\eta\eta}C^{-1}_{\xi\xi}-(C^{-1}_{\xi\eta})^{2}\mathbb{I}_{d}\right)\left(\sum\limits_{a=1}^{s}\left(\bm{C}^{-1}_{\eta\eta}b_{a}(1-{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b_{a})+C^{-1}_{\xi\eta}(1-{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b_{a})^{2}\mathbb{I}_{d}\right)\right)}_{\equiv
g^{\mu}_{\xi}}\bm{\mu}_{k},$
$\displaystyle\bm{\mu}_{\eta}=\underbrace{-\beta\left(\bm{C}^{-1}_{\eta\eta}C^{-1}_{\xi\xi}-(C^{-1}_{\xi\eta})^{2}\mathbb{I}_{d}\right)\left(\sum\limits_{a=1}^{s}\left(C^{-1}_{\xi\xi}(1-{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b_{a})^{2}+C^{-1}_{\xi\eta}b_{a}(1-{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b_{a})\mathbb{I}_{d}\right)\right)}_{\equiv
g^{\mu}_{\eta}}\bm{\mu}_{k}.$ (30)

Again, observe that $g^{\mu}_{\eta},g^{\mu}_{\xi}=\mathcal{O}(s)$ and will be
safely neglected later in the derivation. Therefore

$\displaystyle\mathbb{E}_{s,\bm{\eta},\xi}(\star)=\sum\limits_{k=1}^{K}\rho_{k}\frac{e^{{\scriptscriptstyle-\frac{\beta}{2}{\color[rgb]{0,0,0}\lVert\bm{\mu}_{k}\lVert^{2}}\sum\limits_{a=1}^{s}({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b_{a}-1)^{2}+\frac{1}{2}C^{-1}_{\xi\xi}\lVert\bm{\mu}_{\xi}\rVert^{2}+\frac{1}{2}\bm{\mu}_{\eta}^{\top}\bm{C}^{-1}_{\eta\eta}\bm{\mu}_{\eta}+C^{-1}_{\xi\eta}\bm{\mu}_{\xi}^{\top}\bm{\mu}_{\eta}}}}{{\scriptstyle\sqrt{\det\bm{\Sigma}_{k}}\Delta^{\frac{d}{2}}\det\left(\bm{C}^{-1}_{\eta\eta}C^{-1}_{\xi\xi}-(C^{-1}_{\xi\eta})^{2}\mathbb{I}_{d}\right)^{\frac{1}{2}}}}$
$\displaystyle\times\underbrace{\left\langle
e^{{\scriptscriptstyle-\frac{\beta}{2}\sum\limits_{a=1}^{s}\left[\Tr[\frac{\bm{w}_{a}\bm{w}_{a}^{\top}}{d}\sigma\left(\frac{\bm{w}_{a}({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}(\bm{\mu}_{k}+\bm{\eta})+{\bm{\xi}})}{\sqrt{d}}\right)^{\otimes
2}]-2\sigma\left(\frac{\bm{w}_{a}({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}(\bm{\mu}_{k}+\bm{\eta})+{\bm{\xi}})}{\sqrt{d}}\right)^{\top}\frac{\bm{w}_{a}((1-b_{a}{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}})(\bm{\mu}_{k}+\bm{\eta})-b_{a}{\bm{\xi}})}{\sqrt{d}}\right]}}\right\rangle_{P_{\mathrm{eff.}}(\bm{\eta},{\bm{\xi}})}}_{(a)}.$
(31)

We are now in a position to introduce the local fields:

$\displaystyle\lambda_{a}^{\xi}=\frac{\bm{w}_{a}({\bm{\xi}}-\bm{\mu}_{\xi})}{\sqrt{d}},$
$\displaystyle\lambda_{a}^{\eta}=\frac{\bm{w}_{a}(\bm{\eta}-\bm{\mu}_{\eta})}{\sqrt{d}},$
$\displaystyle h^{k}_{a}=\frac{\bm{w}_{a}\bm{\mu}_{k}}{\sqrt{d}},$ (32)

with statistics

$\displaystyle\langle\lambda_{a}^{\eta}\lambda_{b}^{\eta}\rangle\approx\frac{\bm{w}_{a}\bm{\Sigma}_{k}\bm{w}_{b}^{\top}}{d},$
$\displaystyle\langle\lambda_{a}^{\xi}\lambda_{b}^{\eta}\rangle=C_{\xi\eta}\frac{\bm{w}_{a}\bm{w}_{b}^{\top}}{d}\approx
0,$
$\displaystyle\langle\lambda_{a}^{\xi}\lambda_{b}^{\xi}\rangle\approx\Delta\frac{\bm{w}_{a}\bm{w}_{b}^{\top}}{d}.$
(33)

We used the leading order of the covariances
$\bm{C}_{\xi\xi,\eta\eta,\eta\xi}$. One therefore has to introduce the summary
statistics:

$\displaystyle
Q_{ab}=\frac{\bm{w}_{a}\bm{w}_{b}^{\top}}{d}\in\mathbb{R}^{p\times p},$
$\displaystyle
S^{k}_{ab}=\frac{\bm{w}_{a}\bm{\Sigma}_{k}\bm{w}_{b}^{\top}}{d}\in\mathbb{R}^{p\times
p},$ $\displaystyle
m_{a}^{k}=\frac{\bm{w}_{a}\bm{\mu}_{k}}{\sqrt{d}}\in\mathbb{R}^{p}.$ (34)

Note that the local fields $\lambda_{a}^{\eta,\xi}$ thus follow a Gaussian
distribution:

$\displaystyle(\lambda^{\xi}_{a},\lambda^{\eta}_{a})_{a=1}^{s}\sim\mathcal{N}\left(0,\underbrace{\begin{bmatrix}S^{k}&0\\\
0&\Delta Q\end{bmatrix}}_{\equiv\Omega_{k}}\right).$ (35)

Going back to the computation, $(a)$ can be rewritten as

$\displaystyle\Big{\langle}$ $\displaystyle
e^{{\scriptscriptstyle-\frac{\beta}{2}\sum\limits_{a=1}^{s}\Tr\left[Q_{aa}\sigma\left(m_{a}^{k}({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}(1+g^{\mu}_{\eta})+g^{\mu}_{\xi})+{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}\lambda_{a}^{\eta}+\lambda_{a}^{\xi}\right)^{\otimes
2}\right]}}$ $\displaystyle\times
e^{{\scriptscriptstyle-\frac{\beta}{2}\sum\limits_{a=1}^{s}\left[-2\sigma\left(m^{k}_{a}({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}(1+g^{\mu}_{\eta})+g^{\mu}_{\xi})+{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}\lambda_{a}^{\eta}+\lambda_{a}^{\xi}\right)^{\top}\left(m^{k}_{a}((1-{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b_{a})(1+g^{\mu}_{\eta})-b_{a}g^{\mu}_{\xi})+(1-{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b_{a})\lambda^{\eta}_{a}-b_{a}\lambda^{\xi}_{a}\right)\right]}}\Big{\rangle}_{\\{\lambda^{\xi}_{a},\lambda^{\eta}_{a}\\}_{a=1}^{s}}.$
(36)

###  A.3 Reformulating as a saddle-point problem

Introducing Dirac functions enforcing the definitions of $Q_{ab},m_{a}$ brings
the replicated function in the following form:

$\displaystyle\mathbb{E}Z^{s}=$
$\displaystyle\int\prod\limits_{a=1}^{s}db_{a}\prod\limits_{a,b}dQ_{ab}d\hat{Q}_{ab}\prod\limits_{k=1}^{K}\prod\limits_{a=1}^{s}dm_{a}d\hat{m}_{a}\prod\limits_{k=1}^{K}\prod\limits_{a,b}dS^{k}_{ab}d\hat{S}^{k}_{ab}$
$\displaystyle\underbrace{e^{-d\sum\limits_{a\leq
b}Q_{ab}\hat{Q}_{ab}-d\sum\limits_{k=1}^{K}\sum\limits_{a\leq
b}S^{k}_{ab}\hat{S}^{k}_{ab}-\sum\limits_{k=1}^{K}\sum\limits_{a}dm^{k}_{a}\hat{m}^{k}_{a}}}_{e^{\beta
sd\Psi_{t}}}$
$\displaystyle\underbrace{\int\prod\limits_{a=1}^{s}dw_{a}e^{-\beta\sum\limits_{a}g(\bm{w}^{a})+\sum\limits_{a\leq
b}\hat{Q}_{ab}\bm{w}_{a}\bm{w}_{b}^{\top}+\sum\limits_{k=1}^{K}\sum\limits_{a\leq
b}\hat{S}^{k}_{ab}\bm{w}_{a}\bm{\Sigma}_{k}\bm{w}_{b}^{\top}+\sum\limits_{k=1}^{K}\sum\limits_{a}\hat{m}^{k}_{a}\sqrt{d}\bm{w}_{a}\bm{\mu}_{k}}}_{e^{\beta
sd\Psi_{w}}}$
$\displaystyle\underbrace{\left[\sum\limits_{k=1}^{K}\rho_{k}\frac{e^{{\scriptscriptstyle\mathcal{O}(s^{2})e^{-\frac{\beta}{2}{\color[rgb]{0,0,0}\lVert\bm{\mu}_{k}\lVert^{2}}\sum\limits_{a=1}^{s}({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b_{a}-1)^{2}}}}}{{\scriptstyle\sqrt{\det\bm{\Sigma}_{k}}\Delta^{\frac{d}{2}}\det\left(\bm{C}^{-1}_{\eta\eta}C^{-1}_{\xi\xi}-\mathcal{O}(s^{2})\right)^{\frac{1}{2}}}}(a)\right]^{\alpha
d}}_{e^{\alpha sd^{2}\Psi_{\mathrm{quad.}}+\beta sd\Psi_{y}}}.$ (37)

We introduced the trace, entropic and energetic potentials
$\Psi_{t},\Psi_{w},\Psi_{y}$. Since all the integrands scale exponentially (or
faster) with $d$, this integral can be computed using a saddle-point method.
To proceed further, note that the energetic term encompasses two types of
terms, scaling like $sd^{2}$ and $sd$. More precisely,

$\displaystyle\left[\sum\limits_{k=1}^{K}\rho_{k}\frac{e^{-\frac{\beta}{2}{\color[rgb]{0,0,0}\lVert\bm{\mu}_{k}\lVert^{2}\sum\limits_{a=1}^{s}({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b_{a}-1)^{2}}}}{{\scriptstyle\sqrt{\det\bm{\Sigma}_{k}}\Delta^{\frac{d}{2}}\det\left(\bm{C}^{-1}_{\eta\eta}C^{-1}_{\xi\xi}\right)^{\frac{1}{2}}}}(a)\right]^{\alpha
d}$
$\displaystyle=\\!\\!\Bigg{[}1\\!\\!-\\!\\!sd\sum\limits_{k=1}^{K}\rho_{k}\frac{1}{sd}\ln\left({\scriptstyle\sqrt{\det\bm{\Sigma}_{k}}\Delta^{\frac{d}{2}}\det\left(\bm{C}^{-1}_{\eta\eta}C^{-1}_{\xi\xi}\right)^{\frac{1}{2}}}\right)\\!\\!+\\!\\!s\sum\limits_{k=1}^{K}\rho_{k}\frac{1}{s}\ln(a)-\\!\\!\frac{\beta}{2}\sum\limits_{k=1}^{K}\rho_{k}\lVert\bm{\mu}_{k}\lVert^{2}\sum\limits_{a=1}^{s}({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b_{a}\\!\\!-\\!\\!1)^{2}\Bigg{]}^{\alpha
d}$ $\displaystyle\approx e^{-s\alpha
d^{2}\sum\limits_{k=1}^{K}\rho_{k}\frac{1}{sd}\ln\left({\scriptstyle\sqrt{\det\bm{\Sigma}_{k}}\Delta^{\frac{d}{2}}\det\left(\bm{C}^{-1}_{\eta\eta}C^{-1}_{\xi\xi}\right)^{\frac{1}{2}}}\right)+s\alpha
d\sum\limits_{k=1}^{K}\rho_{k}\frac{1}{s}\ln(a)}.$ (38)

Note that the saddle-point method involves an intricate extremization over
$s\times s$ matrices $Q,S_{k}$. To make further progress, we assume the
extremum is realized at the Replica Symmetric (RS) fixed point

$\displaystyle\forall a,\qquad b_{a}=b$ $\displaystyle\forall a,\qquad
m_{a}^{k}=m_{k},\qquad\hat{m}^{k}_{a}=\hat{m}_{k}$ $\displaystyle\forall
a,b,\qquad
Q_{ab}=(r-q)\delta_{ab}+q,\qquad\hat{Q}_{ab}=-\left(\frac{\hat{r}}{2}+\hat{q}\right)\delta_{ab}+\hat{q}$
$\displaystyle\forall a,b,\qquad
S^{k}_{ab}=(r_{k}-q_{k})\delta_{ab}+q_{k},\qquad\hat{S}^{k}_{ab}=-\left(\frac{\hat{r}_{k}}{2}+\hat{q}_{k}\right)\delta_{ab}+\hat{q}_{k}$
(39)

This is a standard assumption known as the RS ansatz [21, 22, 23, 24].

#### Quadratic potential

The quadratic potential $\Psi_{\mathrm{quad.}}$, which correspond to the
leading order term in $d$ in the exponent, needs to be extremized first in the
framework of our saddle-point analysis. Its expression can be simplified as

$\displaystyle\sum\limits_{k=1}^{K}\rho_{k}\ln\left({\scriptstyle\sqrt{\det\bm{\Sigma}_{k}}\Delta^{\frac{d}{2}}\det\left(\bm{C}^{-1}_{\eta\eta}C^{-1}_{\xi\xi}\right)^{\frac{1}{2}}}\right)$
$\displaystyle=\frac{1}{2}\sum\limits_{k=1}^{K}\rho_{k}\ln\det\left[(1+\Delta\beta
sb^{2})\left(\mathbb{I}_{d}+\bm{\Sigma}_{k}\beta
s(\sqrt{1-\Delta}b-1)^{2}\right)\right]$ $\displaystyle=\frac{\beta
s}{2}\sum\limits_{k=1}^{K}\rho_{k}\Tr[\Delta
b^{2}\mathbb{I}_{d}+\bm{\Sigma}_{k}(\sqrt{1-\Delta}b-1)^{2}]$
$\displaystyle=\frac{\beta ds}{2}\left[\Delta
b^{2}+(\sqrt{1-\Delta}b-1)^{2}\frac{1}{d}\sum\limits_{k=1}^{K}\rho_{k}\Tr\bm{\Sigma}_{k}\right],$
(40)

which is extremized for

$\displaystyle
b=\frac{\frac{1}{d}\left(\sum\limits_{k=1}^{K}\rho_{k}\Tr\bm{\Sigma}_{k}\right){\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}}{\frac{1}{d}\left(\sum\limits_{k=1}^{K}\rho_{k}\Tr\bm{\Sigma}_{k}\right){\tiny{\color[rgb]{0,0,0}(1-\Delta)}}+\Delta}.$
(41)

This fixes the skip connection strength $b$.

#### Entropic potential

We now turn to the entropic potential $\Psi_{w}$. It is convenient to
introduce the variance order parameters

$\displaystyle\hat{V}\equiv\hat{r}+\hat{q},$
$\displaystyle\hat{V}_{k}\equiv\hat{r}_{k}+\hat{q}_{k}.$ (42)

The entropic potential can then be expressed as

$\displaystyle e^{\beta sd\Psi_{w}}$
$\displaystyle=\int\prod\limits_{a=1}^{s}dw_{a}e^{-\beta\sum\limits_{a}g(\bm{w}^{a})-\frac{1}{2}\sum\limits_{a=1}^{s}\Tr[\hat{V}\bm{w}_{a}\bm{w}_{a}^{\top}]+\frac{1}{2}\sum\limits_{a,b}\Tr[\hat{q}\bm{w}_{a}\bm{w}_{b}^{\top}]+\hat{m}\sum\limits_{k=1}^{K}\sum\limits_{a=1}^{s}\sqrt{d}\hat{m}_{k}^{\top}\bm{w}_{a}^{\top}\bm{\mu}_{k}}$
$\displaystyle\qquad\qquad\qquad\times
e^{-\frac{1}{2}\sum\limits_{k=1}^{K}\sum\limits_{a=1}^{s}\Tr[\hat{V}_{k}\bm{w}_{a}\bm{\Sigma}_{k}\bm{w}_{a}^{\top}]+\frac{1}{2}\sum\limits_{k=1}^{K}\sum\limits_{a,b}\Tr[\hat{q}_{k}\bm{w}_{a}\bm{\Sigma}_{k}\bm{w}_{b}^{\top}]}$
$\displaystyle=\int D\bm{\Xi}_{0}\prod\limits_{k=1}^{K}D\bm{\Xi}_{k}$
$\displaystyle\qquad\left[\int dwe^{-\beta
g(\bm{w})-\frac{1}{2}\Tr[\hat{V}\bm{w}\bm{w}^{\top}+\sum\limits_{k=1}^{K}\hat{V}_{k}\bm{w}\bm{\Sigma}_{k}\bm{w}^{\top}]+\left(\sum\limits_{k=1}^{K}\sqrt{d}\hat{m}_{k}\bm{\mu}^{\top}+\sum\limits_{k=1}^{K}\bm{\Xi}_{k}\odot(\hat{q}_{k}\otimes\bm{\Sigma}_{k})^{\frac{1}{2}}+\bm{\Xi}_{0}\odot(\hat{q}\otimes\mathbb{I}_{d})^{\frac{1}{2}}\right)\odot\bm{w}}\right]^{s}$
$\displaystyle=\int\underbrace{D\bm{\Xi}_{0}\prod\limits_{k=1}^{K}D\bm{\Xi}_{k}}_{\equiv
D\bm{\Xi}}$ $\displaystyle\qquad\left[\int dwe^{-\beta
g(\bm{w})-\frac{1}{2}\bm{w}\odot\left[\hat{V}\otimes\mathbb{I}_{d}+\sum\limits_{k=1}^{K}\hat{V}_{k}\otimes\bm{\Sigma}_{k}\right]\odot\bm{w}+\left(\sum\limits_{k=1}^{K}\sqrt{d}\hat{m}_{k}\bm{\mu}^{\top}+\sum\limits_{k=1}^{K}\bm{\Xi}_{k}\odot(\hat{q}_{k}\otimes\bm{\Sigma}_{k})^{\frac{1}{2}}+\bm{\Xi}_{0}\odot(\hat{q}\otimes\mathbb{I}_{d})^{\frac{1}{2}}\right)\odot\bm{w}}\right]^{s}.$
(43)

Therefore

$\displaystyle\beta\Psi_{w}=\frac{1}{d}\int D\bm{\Xi}\ln\left[\int dwe^{-\beta
g(\bm{w})-\frac{1}{2}\bm{w}\odot\left[\hat{V}\otimes\mathbb{I}_{d}+\sum\limits_{k=1}^{K}\hat{V}_{k}\otimes\bm{\Sigma}_{k}\right]\odot\bm{w}+\left(\sum\limits_{k=1}^{K}\sqrt{d}\hat{m}_{k}\bm{\mu}^{\top}+\sum\limits_{k=1}^{K}\bm{\Xi}_{k}\odot(\hat{q}_{k}\otimes\bm{\Sigma}_{k})^{\frac{1}{2}}+\bm{\Xi}_{0}\odot(\hat{q}\otimes\mathbb{I}_{d})^{\frac{1}{2}}\right)\odot\bm{w}}\right].$
(44)

For a matrix $\bm{\Xi}\in\mathbb{R}^{p\times d}$ and tensors
$\bm{A},\bm{B}\in\mathbb{R}^{p\times d}\otimes\mathbb{R}^{p\times d}$, we
denoted $(\bm{\Xi}\odot\bm{A})_{kl}=\sum_{ij}\bm{\Xi}^{ij}\bm{A}_{ij,kl}$ and
$(\bm{A}\odot\bm{B})_{ij,kl}=\sum_{rs}\bm{A}_{ij,rs}\bm{B}_{rs,kl}$.

#### Energetic potential

In order to compute the energetic potential $\Psi_{y}$, one must first compute
the inverse of the covariance $\Omega_{k}$, given by

$\displaystyle\Omega_{k}^{-1}=\begin{bmatrix}S_{k}^{-1}&0\\\
0&\frac{1}{\Delta}Q^{-1}\end{bmatrix}.$ (45)

It is straightforward to see that $S_{k}^{-1},Q^{-1}$ share the same block
structure as $S_{k},Q$, and we shall name for clarity

$\displaystyle\left(S_{k}^{-1}\right)_{ab}=\left(\tilde{r}_{k}-\tilde{q}_{k}\right)\delta_{ab}+\tilde{q}_{k},$
$\displaystyle\left(Q^{-1}\right)_{ab}=\left(\tilde{r}-\tilde{q}\right)\delta_{ab}+\tilde{q}.$
(46)

From the Sherman-Morisson lemma, and noting

$\displaystyle V_{k}\equiv r_{k}-q_{k}\in\mathbb{R}^{p\times p},$
$\displaystyle V\equiv r-q\in\mathbb{R}^{p\times p},$ (47)

one reaches

$\displaystyle\begin{cases}\tilde{r}=V^{-1}-(V+sq)^{-1}qV^{-1}\\\
\tilde{q}=-(V+sq)^{-1}qV^{-1}\end{cases},$
$\displaystyle\begin{cases}\tilde{r}_{k}=V_{k}^{-1}-(V_{k}+sq_{k})^{-1}q_{k}V_{k}^{-1}\\\
\tilde{q}_{k}=-(V_{k}+sq_{k})^{-1}q_{k}V_{k}^{-1}\end{cases}.$ (48)

We will also need the expression of the determinants

$\displaystyle\ln\det S_{k}=s\det V_{k}+s\Tr[V^{-1}_{k}q_{k}],$ $\displaystyle
Q=s\det V+s\Tr[V^{-1}q].$ (49)

One can finally proceed in analyzing the term $(a)$:

$\displaystyle(a)$
$\displaystyle=\int_{\mathbb{R}^{p}}\frac{\prod\limits_{a=1}^{s}d\lambda^{\eta}_{a}d\lambda^{\xi}_{a}}{(2\pi)^{ps}\sqrt{\det
Q\det S_{k}}\sqrt{\Delta}^{ps}}$ $\displaystyle\qquad\times
e^{-\frac{1}{2}\sum\limits_{a=1}^{s}(\lambda^{\eta}_{a})^{\top}(\tilde{r}_{k}-\tilde{q}_{k})\lambda_{a}^{\eta}-\frac{1}{2}\sum\limits_{a,b}(\lambda^{\eta}_{a})^{\top}\tilde{q}_{k}\lambda_{b}-\frac{1}{2\Delta}\sum\limits_{a=1}^{s}(\lambda^{\xi}_{a})^{\top}(\tilde{r}-\tilde{q})\lambda_{a}^{\xi}-\frac{1}{2\Delta}\sum\limits_{a,b}(\lambda^{\xi}_{a})^{\top}\tilde{q}\lambda_{b}}e^{-\frac{\beta}{2}\sum\limits_{a=1}^{s}(*)}$
$\displaystyle=\int_{\mathbb{R}^{p}}\frac{D\xi D\eta}{(2\pi)^{ps}\sqrt{\det
Q\det S_{k}}\sqrt{\Delta}^{ps}}$
$\displaystyle\qquad\times\left[\int_{\mathbb{R}^{p}}d\lambda_{\eta}d\lambda_{\xi}e^{-\frac{1}{2}\lambda_{\eta}^{\top}V_{k}^{-1}\lambda_{\eta}+\lambda_{\eta}^{\top}V_{k}^{-1}q_{k}^{\frac{1}{2}}\eta-\frac{1}{2\Delta}\lambda_{\xi}^{\top}V^{-1}\lambda_{\xi}+\frac{1}{\sqrt{\Delta}}\lambda_{\xi}^{\top}V^{-1}q^{\frac{1}{2}}\xi-\frac{\beta}{2}(*)}\right]^{s}$
$\displaystyle=\int_{\mathbb{R}^{p}}D\xi
D\eta\left[\int\mathcal{N}\left(\lambda_{\xi};\sqrt{\Delta}q^{\frac{1}{2}}\xi,\Delta
V\right)\mathcal{N}\left(\lambda_{\eta};q^{\frac{1}{2}}_{k}\xi,V_{k}\right)e^{-\frac{\beta}{2}(*)}\right]^{s}$
$\displaystyle=1+s\int_{\mathbb{R}^{p}}D\xi
D\eta\ln\left[\int\mathcal{N}\left(\lambda_{\xi};\sqrt{\Delta}q^{\frac{1}{2}}\xi,\Delta
V\right)\mathcal{N}\left(\lambda_{\eta};q^{\frac{1}{2}}_{k}\xi,V_{k}\right)e^{-\frac{\beta}{2}(*)}\right]$
(50)

where we noted with capital $D$ an integral over a $p-$ dimensional Gaussian
distribution $\mathcal{N}(0,\mathbb{I}_{p})$. Therefore

$\displaystyle\beta\Psi_{y}=\sum\limits_{k=1}^{K}\rho_{k}\int_{\mathbb{R}^{p}}D\xi
D\eta\ln\left[\int\mathcal{N}\left(\lambda_{\xi};\sqrt{\Delta}q^{\frac{1}{2}}\xi,\Delta
V\right)\mathcal{N}\left(\lambda_{\eta};q^{\frac{1}{2}}_{k}\xi,V_{k}\right)e^{-\frac{\beta}{2}(*)}\right].$
(51)

###  A.4 zero-temperature limit

In this subsection, we take the zero temperature $\beta\rightarrow\infty$
limit. Rescaling

$\displaystyle\frac{1}{\beta}V\leftarrow V,$
$\displaystyle\beta\hat{V}\leftarrow\hat{V},$
$\displaystyle\beta^{2}\hat{q}\leftarrow\hat{q},$
$\displaystyle\frac{1}{\beta}V_{k}\leftarrow V_{k},$
$\displaystyle\beta\hat{V}_{k}\leftarrow\hat{V}_{k},$
$\displaystyle\beta^{2}\hat{q}_{k}\leftarrow\hat{q}_{k},$
$\displaystyle\beta\hat{m}_{k}\leftarrow\hat{m}_{k},$ (52)

one has that

$\displaystyle\Psi_{w}$
$\displaystyle=\frac{1}{2d}\Tr[\left(\hat{V}\otimes\mathbb{I}_{d}+\sum\limits_{k=1}^{K}\hat{V}_{k}\otimes\bm{\Sigma}_{k}\right)^{-1}\odot\left(\hat{q}\otimes\mathbb{I}_{d}+\sum\limits_{k=1}^{K}\hat{q}_{k}\otimes\bm{\Sigma}_{k}+d\left(\sum\limits_{k=1}^{K}\hat{m}_{k}\bm{\mu}_{k}^{\top}\right)^{\otimes
2}\right)]-\frac{1}{d}\mathbb{E}_{\bm{\Xi}}\mathcal{M}_{r}(\bm{\Xi}),$ (53)

where we introduced the Moreau enveloppe

$\displaystyle\mathcal{M}_{r}(\bm{\Xi})$
$\displaystyle=\underset{\bm{w}}{\inf}\left\\{\frac{1}{2}\left\lVert{\scriptstyle\left(\hat{V}\otimes\mathbb{I}_{d}+\sum\limits_{k=1}^{K}\hat{V}_{k}\otimes\bm{\Sigma}_{k}\right)^{\frac{1}{2}}\\!\\!\\!\\!\\!\\!\odot\bm{w}-\left(\hat{V}\otimes\mathbb{I}_{d}+\sum\limits_{k=1}^{K}\hat{V}_{k}\otimes\bm{\Sigma}_{k}\right)^{-\frac{1}{2}}\\!\\!\\!\\!\\!\odot\left((\hat{q}\otimes\mathbb{I}_{d})^{\frac{1}{2}}\odot\bm{\Xi}_{0}+\sum\limits_{k=1}^{K}(\hat{q}_{k}\otimes\bm{\Sigma}_{k})^{\frac{1}{2}}\odot\bm{\Xi}_{k}+\sqrt{d}\sum\limits_{k=1}^{K}\hat{m}_{k}\bm{\mu}_{k}^{\top}\right)}\right\lVert^{2}\\!\\!\\!\\!+g(\bm{w})\right\\}.$
(54)

Note that while this optimization problem is still cast in a space of
dimension $p\times d$, for some regularizers $g(\cdot)$ including the usual
$\ell_{1,2}$ penalties, the Moreau enveloppe admits a compact analytical
expression.

The energetic potential $\Psi_{y}$ also simplifies in this limit to

$\displaystyle\Psi_{y}=-\sum\limits_{k=1}^{K}\rho_{k}\mathbb{E}_{\eta,\xi}\mathcal{M}_{k}(\eta,\xi),$
(55)

Where

$\displaystyle\mathcal{M}_{k}(\xi,\eta)=\frac{1}{2}\underset{x,y}{\inf}\Bigg{\\{}$
$\displaystyle\Tr[V_{k}^{-1}\left(y-q_{k}^{\frac{1}{2}}\eta-
m_{k}\right)^{\otimes
2}]+\frac{1}{\Delta}\Tr[V^{-1}\left(x-\sqrt{\Delta}q^{\frac{1}{2}}\xi\right)^{\otimes
2}]$
$\displaystyle\Tr[q\sigma({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}y+x)^{\otimes
2}]-2\sigma({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}y+x)^{\top}((1-{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b)y-bx)\Bigg{\\}}.$
(56)

It is immediate to see that the trace potential $\Psi_{t}$ can be expressed as

$\displaystyle\Psi_{t}=\frac{\Tr[\hat{V}q]-\Tr[\hat{q}V]}{2}+\frac{1}{2}\sum\limits_{k=1}^{K}(\Tr[\hat{V}_{k}q_{k}]-\Tr[\hat{q}_{k}V_{k}])-\sum\limits_{k=1}^{K}m_{k}\hat{m}_{k}.$
(57)

Putting these results together, the total free entropy reads

$\displaystyle\Phi=$
$\displaystyle\underset{q,m,V,\hat{q},\hat{m},\hat{V},\\{q_{k},m_{k},V_{k},\hat{q}_{k},\hat{m}_{k},\hat{V}_{k}\\}_{k=1}^{K}}{\mathrm{extr}}\frac{\Tr[\hat{V}q]-\Tr[\hat{q}V]}{2}+\frac{1}{2}\sum\limits_{k=1}^{K}(\Tr[\hat{V}_{k}q_{k}]-\Tr[\hat{q}_{k}V_{k}])-\sum\limits_{k=1}^{K}m_{k}\hat{m}_{k}$
$\displaystyle-\alpha\sum\limits_{k=1}^{K}\rho_{k}\mathbb{E}_{\xi,\eta}\mathcal{M}_{k}(\xi,\eta)-\frac{1}{d}\mathbb{E}_{\bm{\Xi}}\mathcal{M}_{r}(\bm{\Xi})$
$\displaystyle+\frac{1}{2d}\Tr[\left(\hat{V}\otimes\mathbb{I}_{d}+\sum\limits_{k=1}^{K}\hat{V}_{k}\otimes\bm{\Sigma}_{k}\right)^{-1}\odot\left(\hat{q}\otimes\mathbb{I}_{d}+\sum\limits_{k=1}^{K}\hat{q}_{k}\otimes\bm{\Sigma}_{k}+d\left(\sum\limits_{k=1}^{K}\hat{m}_{k}\bm{\mu}_{k}^{\top}\right)^{\otimes
2}\right)]$ (58)

where we remind that

$\displaystyle\mathcal{M}_{k}(\xi,\eta)=\frac{1}{2}\underset{x,y}{\inf}\Bigg{\\{}$
$\displaystyle\Tr[V_{k}^{-1}\left(y-q_{k}^{\frac{1}{2}}\eta-
m_{k}\right)^{\otimes 2}]+\Tr[V^{-1}\left(x-q^{\frac{1}{2}}\xi\right)^{\otimes
2}]$
$\displaystyle+\Tr[q\sigma({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}y+x)^{\otimes
2}]-2\sigma({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}y+x)^{\top}((1-{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b)y-bx)\Bigg{\\}}.$
(59)

and

$\displaystyle\mathcal{M}_{r}(\bm{\Xi})$
$\displaystyle=\underset{\bm{w}}{\inf}\left\\{\frac{1}{2}\left\lVert{\scriptstyle\left(\hat{V}\otimes\mathbb{I}_{d}+\sum\limits_{k=1}^{K}\hat{V}_{k}\otimes\bm{\Sigma}_{k}\right)^{\frac{1}{2}}\\!\\!\\!\\!\\!\odot\bm{w}-\left(\hat{V}\otimes\mathbb{I}_{d}+\sum\limits_{k=1}^{K}\hat{V}_{k}\otimes\bm{\Sigma}_{k}\right)^{-\frac{1}{2}}\\!\\!\\!\\!\\!\odot\left((\hat{q}\otimes\mathbb{I}_{d})^{\frac{1}{2}}\odot\bm{\Xi}_{0}+\sum\limits_{k=1}^{K}(\hat{q}_{k}\otimes\bm{\Sigma}_{k})^{\frac{1}{2}}\odot\bm{\Xi}_{k}+\sqrt{d}\sum\limits_{k=1}^{K}\hat{m}_{k}\bm{\mu}_{k}^{\top}\right)}\right\lVert^{2}\\!\\!+r(\bm{w})\right\\}$
(60)

###  A.5 Self-consistent equations

The extremization problem (61) can be recast as a system of self-consistent
equations by requiring the gradient with respect to each summary statistic
involved in (61) to be zero. This translates into:

$\displaystyle\begin{cases}\hat{q}_{k}=\alpha\rho_{k}\mathbb{E}_{\xi,\eta}V_{k}^{-1}\left(\mathrm{prox}_{y}^{k}-q_{k}^{\frac{1}{2}}\eta-
m_{k}\right)^{\otimes 2}V_{k}^{-1}\\\
\hat{V}_{k}=-\alpha\rho_{k}q_{k}^{-\frac{1}{2}}\mathbb{E}_{\xi,\eta}V_{k}^{-1}\left(\mathrm{prox}_{y}^{k}-q_{k}^{\frac{1}{2}}\eta-
m_{k}\right)\eta^{\top}\\\
\hat{m}_{k}=\alpha\rho_{k}\mathbb{E}_{\xi,\eta}V_{k}^{-1}\left(\mathrm{prox}_{y}^{k}-q_{k}^{\frac{1}{2}}\eta-
m_{k}\right)\\\
\hat{q}=\alpha\sum\limits_{k=1}^{K}\rho_{k}\mathbb{E}_{\xi,\eta}V^{-1}\left(\mathrm{prox}^{k}_{x}-q^{\frac{1}{2}}\xi\right)^{\otimes
2}V^{-1}\\\
\hat{V}=-\alpha\sum\limits_{k=1}^{K}\rho_{k}\mathbb{E}_{\xi,\eta}\Bigg{[}q^{-\frac{1}{2}}V^{-1}\left(\mathrm{prox}^{k}_{x}-q^{\frac{1}{2}}\xi\right)\xi^{\top}-\sigma\left(\sqrt{1-\Delta}\mathrm{prox}_{y}+\mathrm{prox}_{x}\right)^{\otimes
2}\Bigg{]}\end{cases},$ (61)
$\displaystyle\begin{cases}V_{k}=\frac{1}{d}\mathbb{E}_{\bm{\Xi}}\left[\left(\mathrm{prox}_{r}\odot\left(\hat{q}_{k}\otimes\bm{\Sigma}_{k}\right)^{-\frac{1}{2}}\odot\left(\mathbb{I}_{p}\otimes\bm{\Sigma}_{k}\right)\right)\bm{\Xi}_{k}^{\top}\right]\\\
q_{k}=\frac{1}{d}\mathbb{E}_{\bm{\Xi}}\left[\mathrm{prox}_{r}\bm{\Sigma}_{k}\mathrm{prox}_{r}^{\top}\right]\\\
m_{k}=\frac{1}{\sqrt{d}}\mathbb{E}_{\xi,\eta}\left[\mathrm{prox}_{r}\bm{\mu}_{k}\right]\\\
V=\frac{1}{d}\mathbb{E}_{\bm{\Xi}}\left[\left(\mathrm{prox}_{r}\odot\left(\hat{q}\otimes\mathbb{I}_{d}\right)^{-\frac{1}{2}}\right)\bm{\Xi}_{0}^{\top}\right]\\\
q=\frac{1}{d}\mathbb{E}_{\bm{\Xi}}\left[\mathrm{prox}_{r}\mathrm{prox}_{r}^{\top}\right]\end{cases}.$
(62)

We remind that the skip connection $b$ is fixed to

$\displaystyle
b=\frac{\frac{1}{d}\left(\sum\limits_{k=1}^{K}\rho_{k}\Tr\bm{\Sigma}_{k}\right){\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}}{\frac{1}{d}\left(\sum\limits_{k=1}^{K}\rho_{k}\Tr\bm{\Sigma}_{k}\right){\tiny{\color[rgb]{0,0,0}(1-\Delta)}}+\Delta}.$
(63)

###  A.6 Sharp asymptotic formulae for the learning metrics

The previous subsections have allowed to obtain sharp asymptotic
characterization for a number of summary statistics of the probability (18).
These statistics are in turn sufficient to asymptotically characterize the
learning metrics discussed in the main text. We successively address the test
MSE (6), cosine similarity (7) and training MSE.

#### MSE

The denoising MSE reads

$\displaystyle\mathrm{mse}$
$\displaystyle=\mathbb{E}_{k,\bm{\eta},{\bm{\xi}}}\left\lVert\bm{\mu}_{k}+\bm{\eta}-\left(\hat{b}\sqrt{1-\Delta}(\bm{\mu}_{k}+\bm{\eta})+\hat{b}\sqrt{\Delta}{\bm{\xi}}+\frac{\hat{\bm{w}}^{\top}}{\sqrt{d}}\sigma\left(\frac{\hat{\bm{w}}(\sqrt{1-\Delta}(\bm{\mu}_{k}+\bm{\eta})+\sqrt{\Delta}{\bm{\xi}})}{\sqrt{d}}\right)\right)\right\lVert^{2}$
$\displaystyle=\mathrm{mse}_{\circ}+\sum\limits_{k=1}^{K}\rho_{k}\mathbb{E}_{z}^{\mathcal{N}(0,1)}\left[\Tr[q\sigma\left({\scriptstyle{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}m_{k}+\sqrt{\Delta
q+{\tiny{\color[rgb]{0,0,0}(1-\Delta)}}q_{k}}}z\right)^{\otimes 2}]\right]$
$\displaystyle-2\sum\limits_{k=1}^{K}\rho_{k}\mathbb{E}_{u,v}^{\mathcal{N}(0,1)}\left[{\scriptstyle\sigma\left({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}m_{k}+\sqrt{q_{k}({\tiny{\color[rgb]{0,0,0}1-\Delta}})}u+\sqrt{q}\sqrt{\Delta}v\right)}\right]^{\top}\left({\scriptstyle(1-b{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}})m_{k}+(1-b{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}})\sqrt{q_{k}}u-b\sqrt{q}\sqrt{\Delta}v}\right)$
(64)

where

$\displaystyle\mathrm{mse}_{\circ}=d\Delta
b^{2}+(1-{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b)^{2}\left(\sum\limits_{k=1}^{K}\rho_{k}\left({\color[rgb]{0,0,0}\lVert\bm{\mu}_{k}\lVert^{2}}+\Tr\bm{\Sigma}_{k}\right)\right)$
(65)

#### Cosine similarity

By the very definition of the summary statistics $m_{k},q$:

$\displaystyle\theta_{ik}\equiv\frac{\hat{\bm{w}}_{i}\bm{\mu}_{k}}{\lVert\hat{\bm{w}}\lVert\times\lVert\bm{\mu}_{k}\lVert}=\frac{(m_{k})_{i}}{\sqrt{q_{ii}}\lVert\bm{\mu}_{k}\lVert}$
(66)

Figure 5: (left) Training MSE for the full DAE (3) ($p=1,~{}\sigma=\tanh$).
Solid lines represent the sharp asymptotic formula (A.6); dots correspond to
simulation, training the DAE with the Pytorch implementation of full-batch
Adam, over $T=2000$ epochs using learning rate $\eta=0.05$ and weight decay
$\lambda=0.1$. The data was averaged over $N=5$ instances; error bars are
smaller than the point size. (right) Generalization gap
$\mathrm{mse}_{\hat{f}}-\epsilon_{t}$. Solid lines correspond to the
asymptotic prediction of Result II.3 (for the test MSE) and of (A.6) (for the
train MSE), while dots correspond to simulations. Error bars represent one
standard deviation. The Gaussian mixture is the isotropic binary mixture,
whose parameters are specified in the caption of Fig. 1 in the main text.

#### Training error

The replica formalism also allows to compute a sharp asymptotic
characterization for the training MSE

$\displaystyle\epsilon_{t}\equiv\mathbb{E}_{\mathcal{D}}\left[\frac{1}{nd}\sum\limits_{\mu=1}^{n}\left\lVert\bm{x}^{\mu}-\hat{f}(\tilde{\bm{x}}^{\mu})\right\lVert^{2}\right].$
(67)

We have skipped the discussion of this learning metric in the main text for
the sake of conciseness, but will now detail the derivation. Note that the
training error $\epsilon_{t}$ can be deduced from the risk (5) and the average
regularization as

$\displaystyle\frac{1}{2}\epsilon_{t}=\mathbb{E}_{\mathcal{D}}\left[\hat{\mathcal{R}}(\hat{\bm{w}},\hat{b})\right]-\mathbb{E}_{\mathcal{D}}\left[r(\hat{\bm{w}})\right].$
(68)

Note that in turn, from the definition of the free entropy, the average risk
(training loss) can be computed as

$\displaystyle\mathbb{E}_{\mathcal{D}}\left[\hat{\mathcal{R}}(\hat{\bm{w}},\hat{b})\right]=-\underset{\beta\to\infty}{\lim}\frac{1}{\beta}\frac{\partial\ln
Z}{\partial\beta}=-\underset{\beta\to\infty}{\lim}\frac{1}{\beta}\frac{\partial\Phi}{\partial\beta}.$
(69)

Doing so reveals

$\displaystyle\epsilon_{t}$
$\displaystyle=-2\underset{\beta\to\infty}{\lim}\frac{\partial(\Psi_{y}-d\Psi_{\mathrm{quad}}-\frac{\beta}{2}m(\sqrt{1-\Delta}b-1)^{2}))}{\partial\beta}$
$\displaystyle=\mathrm{mse}_{\circ}+\underset{\beta\to\infty}{\lim}\sum\limits_{k=1}^{K}\rho_{k}\int_{\mathbb{R}^{p}}D\xi
D\eta\frac{\int\mathcal{N}\left(\lambda_{\xi};q^{\frac{1}{2}}\xi,V_{k}\right)\mathcal{N}\left(\lambda_{\eta};q^{\frac{1}{2}}_{k}\xi,V\right)(*)e^{-\frac{\beta}{2}(*)}}{\int\mathcal{N}\left(\lambda_{\xi};q^{\frac{1}{2}}\xi,V_{k}\right)\mathcal{N}\left(\lambda_{\eta};q^{\frac{1}{2}}_{k}\xi,V\right)e^{-\frac{\beta}{2}(*)}}.$
(70)

Taking the $\beta\to\infty$ limit finally leads to

$\displaystyle\epsilon_{t}=\mathrm{mse}_{\circ}+\sum\limits_{k=1}^{K}\rho_{k}\mathbb{E}_{\xi,\eta}$
$\displaystyle\Bigg{[}\Tr[q\sigma({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}\mathrm{prox}_{y}+\mathrm{prox}_{x})^{\otimes
2}]$
$\displaystyle-2\sigma({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}\mathrm{prox}_{y}+\mathrm{prox}_{x})^{\top}((1-{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}\hat{b})\mathrm{prox}_{y}-\hat{b}\mathrm{prox}_{x})\Bigg{]}$
(71)

The sharp asymptotic formula (A.6) agrees well with numerical simulations,
Fig. 5. As is intuitive, the generalization gap
$\mathrm{mse}_{\hat{f}}-\epsilon_{t}$ grows with the noise level $\Delta$, as
the learning problem grows harder.

###  A.7 Generic asymptotically exact formulae

This last result ends the derivation of the generic version of Result II.3,
i.e. not assuming Assumptions II.1 and II.2.

Importantly, note that the characterization (61) is, like the formulae in e.g.
[47, 48], asymptotically exact, but not fully asymptotic, as the equations
(61) still involve high-dimensional quantities. In practice however, for
standard regularizers $r(\cdot)$ like $\ell_{1}$ or $\ell_{2}$, the proximal
$\mathrm{prox}_{r}$ admits a simple closed-form expression, and the average
over $\bm{\Xi}$ can be carried out analytically. The only high-dimensional
operation left is then taking the trace of linear combinations of
$\bm{\Sigma}_{k}$ matrices (or the inverse of such combinations), which is
numerically less demanding, and analytically much simpler, than the high
dimensional optimization (5) and averages (6) involved in the original
problem. We give in the next paragraph an example of how Assumption II.1 can
be relaxed, for a binary mixture with arbitrary (i.e. not necessarily jointly
diagonalizable) covariances $\bm{\Sigma}_{1,2}$.

#### Example: Anisotropic, heteroscedastic binary mixture

As an illustration, we provide the equivalent of (61) for a binary mixture,
with generically anisotropic and distinct covariances, for $p=1$, thereby
breaking free from Assumption II.1. We index the clusters by $+,-$ rather than
$1,2$ for notational convenience. We remind that

$\displaystyle\mathcal{M}_{\pm}(\xi,\eta)=\frac{1}{2}\underset{x,y}{\inf}\left\\{{\scriptstyle\frac{1}{\Delta
V}(x-\sqrt{q}\sqrt{\Delta}\xi)^{2}+\frac{1}{V_{\pm}}(y-\sqrt{q_{\pm}}\eta\mp
m)^{2}+q\sigma({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}y+x)^{2}-2\sigma({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}y+x)((1-{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}b)y-bx)}\right\\}.$
(72)

Then the self-consistent equations (61) simplify to

$\displaystyle\begin{cases}q=\frac{1}{d}\Tr[\left((\lambda+\hat{V})\mathbb{I}_{d}+\hat{V}_{+}\bm{\Sigma}_{+}+\hat{V}_{-}\bm{\Sigma}_{-}\right)^{-2}\left(\hat{q}\mathbb{I}_{d}+\hat{q}_{+}\bm{\Sigma}_{+}+\hat{q}_{-}\bm{\Sigma}_{-}+\hat{m}^{2}d\bm{\mu}\bm{\mu}^{\top}\right)]\\\
m=\hat{m}\Tr[\left((\lambda+\hat{V})\mathbb{I}_{d}+\hat{V}_{+}\bm{\Sigma}_{+}+\hat{V}_{-}\bm{\Sigma}_{-}\right)^{-1}\bm{\mu}\bm{\mu}^{\top}]\\\
V=\frac{1}{d}\Tr[\left(\lambda+\hat{V})\mathbb{I}_{d}+\hat{V}_{+}\bm{\Sigma}_{+}+\hat{V}_{-}\bm{\Sigma}_{-}\right)^{-1}]\\\
q_{\pm}=\frac{1}{d}\Tr[{\scriptstyle\left((\lambda+\hat{V})\mathbb{I}_{d}+\hat{V}_{+}\bm{\Sigma}_{+}+\hat{V}_{-}\bm{\Sigma}_{-}\right)^{-1}\bm{\Sigma}_{\pm}\left((\lambda+\hat{V})\mathbb{I}_{d}+\hat{V}_{+}\bm{\Sigma}_{+}+\hat{V}_{-}\bm{\Sigma}_{-}\right)^{-1}\left(\hat{q}\mathbb{I}_{d}+\hat{q}_{+}\bm{\Sigma}_{+}+\hat{q}_{-}\bm{\Sigma}_{-}+\hat{m}^{2}d\bm{\mu}\bm{\mu}^{\top}\right)}]\\\
V_{\pm}=\frac{1}{d}\Tr[\bm{\Sigma}_{\pm}\left((\lambda+\hat{V})\mathbb{I}_{d}+\hat{V}_{+}\bm{\Sigma}_{+}+\hat{V}_{-}\bm{\Sigma}_{-}\right)^{-1}]\end{cases}$
$\displaystyle\begin{cases}\hat{m}=\alpha\left[\rho\frac{1}{V_{+}}\mathbb{E}_{\xi,\eta}\left(\mathrm{prox}_{y}^{+}-\sqrt{q_{+}}\eta-m\right)-(1-\rho)\frac{1}{V_{-}}\mathbb{E}_{\xi,\eta}\left(\mathrm{prox}_{y}^{-}-\sqrt{q_{-}}\eta+m\right)\right]\\\
\hat{q}=\frac{\alpha}{V^{2}\Delta}\left[\rho\mathbb{E}_{\xi,\eta}(\mathrm{prox}_{x}^{+}-\sqrt{q}\sqrt{\Delta}\xi)^{2}+(1-\rho)\mathbb{E}_{\xi,\eta}(\mathrm{prox}_{x}^{-}-\sqrt{q}\sqrt{\Delta}\xi)^{2}\right]\\\
\hat{V}=-\alpha\left[\rho\mathbb{E}_{\xi,\eta}\Bigg{[}\frac{\xi(\mathrm{prox}_{x}^{+}-\sqrt{q}\sqrt{\Delta}\xi)}{\sqrt{q}\sqrt{\Delta}V}-\sigma({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}\mathrm{prox}_{y}^{+}+\mathrm{prox}_{x}^{+})^{2}\right]\\\
\qquad+(1-\rho)\mathbb{E}_{\xi,\eta}\left[\frac{\xi(\mathrm{prox}_{x}^{-}-\sqrt{q}\sqrt{\Delta}\xi)}{\sqrt{q}\sqrt{\Delta}V}-\sigma({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}\mathrm{prox}_{y}^{-}+\mathrm{prox}_{x}^{-})^{2}\right]\Bigg{]}\\\
\hat{q}_{+}=\alpha\rho\frac{1}{V_{+}^{2}}\mathbb{E}_{\xi,\eta}\left(\mathrm{prox}_{y}^{+}-\sqrt{q_{+}}\eta-m\right)^{2}\\\
\hat{q}_{-}=\alpha(1-\rho)\frac{1}{V_{-}^{2}}\mathbb{E}_{\xi,\eta}\left(\mathrm{prox}_{y}^{-}-\sqrt{q_{-}}\eta+m\right)^{2}\\\
\hat{V}_{+}=-\alpha\rho\frac{1}{\sqrt{q_{+}}V_{+}}\mathbb{E}_{\xi,\eta}\left(\mathrm{prox}_{y}^{+}-\sqrt{q_{+}}\eta-m\right)\eta\\\
\hat{V}_{-}=-\alpha(1-\rho)\frac{1}{\sqrt{q_{-}}V_{-}}\mathbb{E}_{\xi,\eta}\left(\mathrm{prox}_{y}^{-}-\sqrt{q_{-}}\eta+m\right)\eta\end{cases}$
(73)

These equations are, as previously discussed, asymptotically exact as
$d\to\infty$. While they still involve traces over high-dimensional matrices
$\bm{\Sigma}_{k}$, this is a very simple operation. We have therefore reduced
the original high-dimensional optimization problem (5) to the much simpler one
of computing traces like (A.7). Crucially, while these traces cannot be
generally simplified without Assumption II.1, they provide the benefit of a
simple and compact expression which bypasses the need of jointly
diagonalizable covariances, and which can thus be readily evaluated for any
covariance, like real-data covariances. This versatility is leveraged in
Appendix D.

###  A.8 End of the derivation: Result II.3

We finally provide the last step in the derivation of the fully asymptotic
characterization of Result II.3. Assuming
$r(\cdot)=\nicefrac{{\lambda}}{{2}}\lVert\cdot\lVert_{F}^{2}$ (Assumption
II.2), the Moreau envelope $\mathcal{M}_{r}$ assumes a very simple expression,
and no longer needs to be written as a high-dimensional optimization problem.
The resulting free energy can be compactly written as:

$\displaystyle\Phi=$
$\displaystyle\underset{q,m,V,\hat{q},\hat{m},\hat{V},\\{q_{k},m_{k},V_{k},\hat{q}_{k},\hat{m}_{k},\hat{V}_{k}\\}_{k=1}^{K}}{\mathrm{extr}}\frac{\Tr[\hat{V}q]-\Tr[\hat{q}V]}{2}+\frac{1}{2}\sum\limits_{k=1}^{K}(\Tr[\hat{V}_{k}q_{k}]-\Tr[\hat{q}_{k}V_{k}])-\sum\limits_{k=1}^{K}m_{k}\hat{m}_{k}$
$\displaystyle-\alpha\sum\limits_{k=1}^{K}\rho_{k}\mathbb{E}_{\xi,\eta}\mathcal{M}_{k}(\xi,\eta)$
$\displaystyle+\underbrace{\frac{1}{2d}\Tr[\left(\lambda\mathbb{I}_{p}\odot\mathbb{I}_{d}+\hat{V}\otimes\mathbb{I}_{d}+\sum\limits_{k=1}^{K}\hat{V}_{k}\otimes\bm{\Sigma}_{k}\right)^{-1}\\!\\!\\!\\!\odot\left(\hat{q}\otimes\mathbb{I}_{d}+\sum\limits_{k=1}^{K}\hat{q}_{k}\otimes\bm{\Sigma}_{k}+d\left(\sum\limits_{k=1}^{K}\hat{m}_{k}\bm{\mu}_{k}^{\top}\right)^{\otimes
2}\right)]}_{(*)}$ (74)

$(*)$ is the only segment which still involves high-dimensional matrices, and
is therefore not yet full asymptotic. Assuming Assumption II.1, $(*)$ can be
massaged into

$\displaystyle(*)$
$\displaystyle=\frac{1}{2d}\sum\limits_{i=1}^{d}\Tr[\left(\lambda+\hat{V}+\lambda^{k}_{i}\hat{V}_{k}\right)^{-1}\left(\hat{q}+\sum\limits_{k=1}^{K}\lambda^{k}_{i}\hat{q}_{k}+\sum\limits_{1\leq
k,j\leq
K}e_{i}^{\top}\mu_{j}\hat{m}_{j}\hat{m}_{k}^{\top}\bm{\mu}_{k}^{\top}e_{i}\right)]$
$\displaystyle\overset{d\to\infty}{=}\frac{1}{2}\int
d\nu(\gamma,\tau)\Tr[\left(\lambda+\hat{V}+\gamma_{k}\hat{V}_{k}\right)^{-1}\left(\hat{q}+\sum\limits_{k=1}^{K}\gamma_{k}\hat{q}_{k}+\sum\limits_{1\leq
k,j\leq K}\tau_{j}\tau_{k}\hat{m}_{j}\hat{m}_{k}^{\top}\right)]$ (75)

The fully asymptotic free energy thus becomes

$\displaystyle\Phi=$
$\displaystyle\underset{q,m,V,\hat{q},\hat{m},\hat{V},\\{q_{k},m_{k},V_{k},\hat{q}_{k},\hat{m}_{k},\hat{V}_{k}\\}_{k=1}^{K}}{\mathrm{extr}}\frac{\Tr[\hat{V}q]-\Tr[\hat{q}V]}{2}+\frac{1}{2}\sum\limits_{k=1}^{K}(\Tr[\hat{V}_{k}q_{k}]-\Tr[\hat{q}_{k}V_{k}])-\sum\limits_{k=1}^{K}m_{k}\hat{m}_{k}$
$\displaystyle-\alpha\sum\limits_{k=1}^{K}\rho_{k}\mathbb{E}_{\xi,\eta}\mathcal{M}_{k}(\xi,\eta)$
$\displaystyle+\frac{1}{2}\int
d\nu(\gamma,\tau)\Tr[\left(\lambda+\hat{V}+\gamma_{k}\hat{V}_{k}\right)^{-1}\left(\hat{q}+\sum\limits_{k=1}^{K}\gamma_{k}\hat{q}_{k}+\sum\limits_{1\leq
k,j\leq K}\tau_{j}\tau_{k}\hat{m}_{j}\hat{m}_{k}^{\top}\right)]$ (76)

Requiring the free energy to be extremized implies

$\displaystyle\begin{cases}\hat{q}_{k}=\alpha\rho_{k}\mathbb{E}_{\xi,\eta}V_{k}^{-1}\left(\mathrm{prox}_{y}^{k}-q_{k}^{\frac{1}{2}}\eta-
m_{k}\right)^{\otimes 2}V_{k}^{-1}\\\
\hat{V}_{k}=-\alpha\rho_{k}q_{k}^{-\frac{1}{2}}\mathbb{E}_{\xi,\eta}V_{k}^{-1}\left(\mathrm{prox}_{y}^{k}-q_{k}^{\frac{1}{2}}\eta-
m_{k}\right)\eta^{\top}\\\
\hat{m}_{k}=\alpha\rho_{k}\mathbb{E}_{\xi,\eta}V_{k}^{-1}\left(\mathrm{prox}_{y}^{k}-q_{k}^{\frac{1}{2}}\eta-
m_{k}\right)\\\
\hat{q}=\alpha\sum\limits_{k=1}^{K}\rho_{k}\mathbb{E}_{\xi,\eta}V^{-1}\left(\mathrm{prox}^{k}_{x}-q^{\frac{1}{2}}\xi\right)^{\otimes
2}V^{-1}\\\
\hat{V}=-\alpha\sum\limits_{k=1}^{K}\rho_{k}\mathbb{E}_{\xi,\eta}\Bigg{[}q^{-\frac{1}{2}}V^{-1}\left(\mathrm{prox}^{k}_{x}-q^{\frac{1}{2}}\xi\right)\xi^{\top}-\sigma\left(\sqrt{1-\Delta}\mathrm{prox}_{y}^{k}+\mathrm{prox}_{x}^{k}\right)^{\otimes
2}\Bigg{]}\end{cases}$ (77) $\displaystyle\begin{cases}q_{k}=\int
d\nu(\gamma,\tau)\gamma_{k}\left(\lambda\mathbb{I}_{p}+\hat{V}+\sum\limits_{j=1}^{K}\gamma_{j}\hat{V}_{j}\right)^{-2}\left(\hat{q}+\sum\limits_{j=1}^{K}\gamma_{j}\hat{q}_{j}+\sum\limits_{1\leq
j,l\leq K}\tau_{j}\tau_{l}\hat{m}_{j}\hat{m}_{l}^{\top}\right)\\\ V_{k}=\int
d\nu(\gamma,\tau)\gamma_{k}\left(\lambda\mathbb{I}_{p}+\hat{V}+\sum\limits_{j=1}^{K}\gamma_{j}\hat{V}_{j}\right)^{-1}\\\
m_{k}=\int
d\nu(\gamma,\tau)\tau_{k}\left(\lambda\mathbb{I}_{p}+\hat{V}+\sum\limits_{j=1}^{K}\gamma_{j}\hat{V}_{j}\right)^{-1}\sum\limits_{j=1}^{K}\tau_{j}\hat{m}_{j}\\\
q=\int
d\nu(\gamma,\tau)\left(\lambda\mathbb{I}_{p}+\hat{V}+\sum\limits_{j=1}^{K}\gamma_{j}\hat{V}_{j}\right)^{-2}\left(\hat{q}+\sum\limits_{j=1}^{K}\gamma_{j}\hat{q}_{j}+\sum\limits_{1\leq
j,l\leq K}\tau_{j}\tau_{l}\hat{m}_{j}\hat{m}_{l}^{\top}\right)\\\ V=\int
d\nu(\gamma,\tau)\left(\lambda\mathbb{I}_{p}+\hat{V}+\sum\limits_{j=1}^{K}\gamma_{j}\hat{V}_{j}\right)^{-1}\end{cases},$
(78)

which is exactly (II.3). This closes the derivation of Result II.3. ∎

###  A.9 Additional comparisons

Figure 6: Same parameters as Fig. 1 in the main text.
$\alpha=1,K=2,\rho_{1,2}=\nicefrac{{1}}{{2}},\bm{\Sigma}_{1,2}=0.3\times\mathbb{I}_{d},p=1$;
the cluster mean $\mu_{1}=-\mu_{2}$ was taken as a random Gaussian vector of
norm $1$. Difference in MSE between the full DAE $\hat{f}$ (3) and the
rescaling network $\hat{r}$ (4) for $\sigma(\cdot)=\tanh$ (blue) or
$\sigma(x)=x$ (red). Solid lines correspond to the sharp asymptotic
characterization (II.3), which is a particular case of Result II.3. Dots
represent numerical simulations for $d=700$, training the DAE using the
Pytorch implementation of full-batch Adam [43], with learning rate $\eta=0.05$
over $2000$ epochs, averaged over $N=2$ instances. Error bars represent one
standard deviation. In dashed green, numerical simulations for the same
architecture, but untied weights. The learning parameters were left unchanged.

We close this appendix by providing further discussion on some points,
including the role of the non-linearity $\sigma(\cdot)$, the influence of the
weight-tying assumption, and of the (an-)isotropy and
(homo/hetero)-scedasticity of the clusters.

#### The activation function

The figures in the main text were generated using $\sigma=\tanh$. Several
studies for RAEs, however, have highlighted that an auto-encoder would
optimally seek to place itself in the linear region of its non-linearity, so
as to learn the principal components of the training data, at least in the
non-linear untied weights case, see e.g. [3, 8]. In light of these findings,
it is legitimate to wonder whether setting a linear activation $\sigma(x)=x$
would not yield a better MSE. Fig. 6 shows that it is not the case. For the
binary isotropic mixture (II.3), the linear activations yield a worse
performance than the tanh activation.

#### Weight-tying

We have assumed that the weights of the DAE (3) were tied. While this
assumption originates mainly for technical reasons in the derivation, note
that it has been also used and discussed in practical setting [3] as a way to
prevent the DAE from going to the linear region of its non-linearity, by
making the norm of the first layer weights very small, and that of the second
layer very large to compensate [8]. A full extension of Result II.3 to the
case of the untied weight is a lengthy theoretical endeavour, which we leave
for future work. However, note that Fig. 6 shows, in the binary isotropic
case, that untying the weights actually leads to a worse MSE than the DAE with
tied weights. This is a very interesting observation, as a DAE with untied
weights obviously has the expressive power to implement a DAE with tied
weights. Therefore, this effect should be mainly due to the untied landscape
presenting local minima trapping the optimization. A full landscape analysis
would present a very important future research topic. Finally, remark that
this stands in sharp contrast to the observation of [8] for non-linear RAEs,
where weight-tying worsened the performance, as it prevents the RAE from
implementing a linear principal component analysis.

An important conclusion from these two observations is that, in contrast to
RAEs, the DAE (3) does not just implement a linear principal component
analysis — weight tying, and the presence of a non-linearity, which are
obstacles for RAEs in reaching the PCA MSE, lead for DAEs (3) to a better MSE.
The model (3) therefore constitutes a model of a genuinely non-linear
algorithm, where the non-linearity is helpful and is not undone during
training. Further discussion can be found in Section III of the main text.

Figure 7:
$\alpha=1,K=2,\rho_{1,2}=\nicefrac{{1}}{{2}},p=1,\lambda=0.1,\sigma=\tanh$;
the cluster mean $\mu_{1}=-\mu_{2}$ was taken as a random Gaussian vector of
norm $1$. The cluster covariances were independently drawn from the Wishart-
Laguerre ensemble, with aspect ratios $\nicefrac{{5}}{{6}}$ and
$\nicefrac{{5}}{{7}}$, before being normalized by $10$ to match the trace of
the isotropic case in Fig. 1. In particular, the eigenvectors are independent,
so the mixture is totally anisotropic with different main variance directions
and heteroscedastic. (left) In blue, the difference in MSE between the full
DAE $\hat{f}$ (3) and the rescaling network $\hat{r}$ (4). Solid lines
correspond to the sharp asymptotic characterization of Result II.3 in its
generic formulation (61). Dots represent numerical simulations for $d=500$,
training the DAE using the Pytorch implementation of full-batch Adam [43],
with learning rate $\eta=0.05$ over $2000$ epochs, averaged over $N=5$
instances Error bars represent one standard deviation. (right) Cosine
similarity $\theta$ (7) (green), squared weight norm
$\nicefrac{{\lVert\hat{\bm{w}}\lVert^{2}}}{{d}}$ (red) and skip connection
strength $\hat{b}$ (blue). Solid lines correspond to the formulas (12)(13) and
(11) of Result II.3; dots are numerical simulations. For completeness, the
cosine similarity of the first principal component of the train set is plotted
in dashed black.

#### Heteroscedasticity, anisotropy

Fig. 1 and 3 in the main text all represent the isotropic, homoscedastic case.
This simple model actually encapsulates all the interesting phenomenology. In
this paragraph, we discuss for completeness the generically heteroscedastic,
anisotropic case. We consider a binary mixture, with covariances
$\bm{\Sigma}_{1},\bm{\Sigma}_{2}$ independently drawn from the Wishart-
Laguerre ensemble, with aspect ratios $\nicefrac{{5}}{{7}}$ and
$\nicefrac{{5}}{{6}}$. Therefore, the clusters are anisotropic, and
eigendirections associated with the largest eigenvalues are more "stretched"
(i.e. induce higher cluster variance). Furthermore, since the set of
eigenvectors of $\bm{\Sigma}_{1},\bm{\Sigma}_{2}$ are independent, the two
clusters are stretched in different directions. To ease the comparison with
Fig. 1 (isotropic, homoscedastic), these Wishart matrices were further divided
by $10$, so that the trace is approximately $0.09$ like in Fig. 1 – i.e., the
clusters have the same average extension. Fig. 7 presents the resulting
metrics and summary statistics. Again, the agreement between the theory II.3
and numerical simulation using Pytorch is compelling. Qualitatively, the
observations made in Section III still hold true in the anisotropic
heteroscedastic case:

  * •

The curve of the cosine similarity still displays a non-monotonic behaviour,
signalling the preferential activation by the DAE of its skip connection at
low $\Delta$ and of its network component at high $\Delta$.

  * •

The skip connection strength $\hat{b}$ is higher at small $\Delta$ and
decreases, in connection to the previous remark.

  * •

The norm of the weight matrix $\hat{\bm{w}}$ overall increases with the noise
level $\Delta$, signalling an increasing usage by the DAE of its bottleneck
network component, again in accordance to the previous remark.

Therefore, the generic case does not introduce qualitative changes compared to
the isotropic case.

Figure 8: Same parameters as Fig. 1 in the main text.
$\alpha=1,K=2,\rho_{1,2}=\nicefrac{{1}}{{2}},\bm{\Sigma}_{1,2}=0.09\times\mathbb{I}_{d},p=1,\sigma=\tanh$;
the cluster mean $\mu_{1}=-\mu_{2}$ was taken as a random Gaussian vector of
norm $1$. Difference in MSE between the full DAE $\hat{f}$ (3) and the
rescaling network $\hat{r}$ (4) for $\lambda=0.1$ (blue) or $\lambda=0.001$
(red). Lines correspond to the sharp asymptotic characterization (II.3) of
Result II.3. Dots represent numerical simulations for $d=700$, training the
DAE using the Pytorch implementation of full-batch Adam [43], with learning
rate $\eta=0.05$ over $2000$ epochs, averaged over $N=2$ instances. (right)
Cosine similarities, for $\lambda=0.1$ (green) and $\lambda=0.001$ (orange);
lines correspond to the asymptotic formulae of Result II.3, dots represent the
numerical simulations.

#### Strength of the weight decay $\lambda$

A $\ell_{2}$ regularization (weight decay) $\lambda=0.1$ was adopted in the
main text. In fact, the value of the strength of the weight decay was not
found to sensibly influence the curves, and again, the qualitative phenomena
discussed in Section III are observed for any value. Fig. 8 shows, for an
isotropic binary mixture, the MSE difference
$\mathrm{mse}_{\hat{f}}-\mathrm{mse}_{\hat{r}}$ (6) and the cosine similarity
$\theta$ (7) for regularization strength $\lambda=0.1$ and $\lambda=0.001$. As
can be observed, even reducing the regularization a hundredfold does not
change at all the qualitative picture – in particular, the non-monotonicity of
the cosine similarity, discussed in Section III – , and very little the
quantitative values.

##  Appendix B Tweedie baseline

###  B.1 Oracle denoiser

Tweedie’s formula [42] provides the best estimator of a clean sample $x$ given
a noisy sample $\tilde{x}$ corrupted by Gaussian noise, as

$\displaystyle
t(\tilde{\bm{x}})=\frac{\tilde{\bm{x}}}{\sqrt{1-\Delta}}+\frac{\Delta}{\sqrt{1-\Delta}}\gradient\mathbb{P}|_{\tilde{\bm{x}}}.$
(79)

Note that this gives the Bayes-optimal estimator for the MSE, assuming perfect
knowledge of the clean data distribution $\mathbb{P}$ from which the so-called
score $\gradient\mathbb{P}$ has to be computed. Of course, this knowledge is
inaccessible in general, where the only information on $\mathbb{P}$ is
provided in the form of the train set $\mathcal{D}$. Therefore, Tweedie’s
formula does not give a learning algorithm, but allows to give an oracle
lower-bound on the achievable MSEs, as any learning algorithm will yield a
higher MSE than this information-theoretic baseline.

For a generic Gaussian mixture (2), Tweedie’s formula reduces to the
expression

$\displaystyle
t(\bm{x})=\frac{1}{\sqrt{1-\Delta}}\frac{\sum\limits_{k=1}^{K}\rho_{k}\left[\tilde{\bm{\Sigma}}_{k}(\tilde{\bm{\Sigma}}_{k}+\Delta\mathbb{I}_{d})^{-1}\bm{x}+(\tilde{\bm{\Sigma}}_{k}+\Delta\mathbb{I}_{d})^{-1}\tilde{\bm{\mu}}_{k})\right]\mathcal{N}\left(\bm{x};\tilde{\bm{\mu}}_{k},\tilde{\bm{\Sigma}}_{k}+\Delta\mathbb{I}_{d}\right)}{\sum\limits_{k=1}^{K}\rho_{k}\mathcal{N}\left(\bm{x};\tilde{\bm{\mu}}_{k},\tilde{\bm{\Sigma}}_{k}+\Delta\mathbb{I}_{d}\right)},$
(80)

where we noted

$\displaystyle\tilde{\bm{\mu}}_{k}=\sqrt{1-\Delta}\bm{\mu}_{k},$
$\displaystyle\tilde{\bm{\Sigma}}_{k}=(1-\Delta)\bm{\Sigma}_{k}.$ (81)

###  B.2 Oracle test MSE

In general, except in the complete homoscedastic case where all clusters have
the same covariance, there is no closed-form asymptotic expression for the MSE
achieved by the Tweedie denoiser. In the binary homoscedastic case
$\bm{\mu}_{1}=-\bm{\mu}_{2}\equiv\bm{\mu}$,
$\bm{\Sigma}_{1}=\bm{\Sigma}_{2}=\sigma^{2}\mathbb{I}_{d}$ (see Fig. 1),
Tweedie’s formula (80) reduces to the compact form

$\displaystyle
t(\bm{x})=\frac{\sqrt{1-\Delta}\sigma^{2}}{\sigma^{2}(1-\Delta)+\Delta}\tilde{\bm{x}}+\frac{\Delta}{\sigma^{2}(1-\Delta)+\Delta}\tanh\left(\frac{\tilde{\bm{x}}^{\top}\bm{\mu}\sqrt{1-\Delta}}{\sigma^{2}(1-\Delta)+\Delta}\right)\times\bm{\mu}$
(82)

Note that this is of the same form as the DAE architecture (3). The associated
MSE reads

$\displaystyle\mathrm{mse}^{\star}$
$\displaystyle=\mathbb{E}_{\bm{x},\bm{\xi}}\left\lVert\bm{x}-\left[b\times({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}\bm{x}+\bm{\xi})+\frac{\Delta\bm{\mu}}{\sigma^{2}{\tiny{\color[rgb]{0,0,0}(1-\Delta)}}+\Delta}\sigma\left(\frac{{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}\bm{\mu}^{\top}({\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}\bm{x}+\bm{\xi})}{\sigma^{2}{\tiny{\color[rgb]{0,0,0}(1-\Delta)}}+\Delta}\right)\right]\right\rVert^{2}.$
(83)

A sharp asymptotic formula can be found to be

$\displaystyle\mathrm{mse}^{\star}$
$\displaystyle=\mathrm{mse}_{\mathrm{\circ}}+\frac{\sigma_{e}^{4}}{(\sigma^{2}{\tiny{\color[rgb]{0,0,0}(1-\Delta)}}+\Delta)^{2}}\frac{1}{2}\sum\limits_{s=\pm
1}\mathbb{E}_{z}^{\mathcal{N}(0,1)}\left[\sigma\left(\frac{{\scriptstyle{\tiny{\color[rgb]{0,0,0}(1-\Delta)}}s+{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}z\times\sqrt{\Delta+{\tiny{\color[rgb]{0,0,0}(1-\Delta)}}\sigma^{2}}}}{\sigma^{2}{\tiny{\color[rgb]{0,0,0}(1-\Delta)}}+\Delta}\right)^{2}\right]$
$\displaystyle-\frac{2\Delta}{\sigma^{2}{\tiny{\color[rgb]{0,0,0}(1-\Delta)}}+\Delta}\frac{1}{2}\sum\limits_{s=\pm
1}\mathbb{E}_{z}^{\mathcal{N}(0,1)}\left[\sigma\left(\frac{{\scriptstyle{\tiny{\color[rgb]{0,0,0}(1-\Delta)}}s+{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}}z\times\sqrt{\Delta+{\tiny{\color[rgb]{0,0,0}(1-\Delta)}}\sigma^{2}}}}{\sigma^{2}{\tiny{\color[rgb]{0,0,0}(1-\Delta)}}+\Delta}\right)\right]\times\left((1-b{\tiny{\color[rgb]{0,0,0}\sqrt{1-\Delta}}})s\right).$
(84)

This is the theoretical characterization plotted in Figs. 1 and 3 in the main
text.

##  Appendix C Bayes-optimal MSE

###  C.1 Assuming imperfect prior knowledge

The Tweedie denoiser discussed in Appendix B is an oracle denoiser, in the
sense that perfect knowledge of all the parameters
$\\{\bm{\mu}_{k},\bm{\Sigma}_{k}\\}_{k=1}^{K}$ of the Gaussian mixture
distribution $\mathbb{P}$ (2) is assumed. Therefore, though the comparison of
the DAE (3) MSE with the oracle Tweedie MSE $\mathrm{mse}^{\star}$ does
provide useful intuition (see Fig. 3), it importantly does not allow to
disentangle which part of the MSE of the DAE is due to the limited
expressivity of the architecture, and which is due to the limited
availablility of training samples – which entails imperfect knowledge of the
parameters $\\{\bm{\mu}_{k},\bm{\Sigma}_{k}\\}_{k=1}^{K}$ of the Gaussian
mixture distribution $\mathbb{P}$ (2). In this appendix, we derive the MSE of
the Bayes optimal estimator (with respect to the MSE), assuming only knowledge
of the distribution from which the parameters
$\\{\bm{\mu}_{k},\bm{\Sigma}_{k}\\}_{k=1}^{K}$ are drawn. For simplicity and
conciseness, we limit ourselves to the binary isotropic and homoscedastic
mixture (see e.g. Fig. 1), for which
$\bm{\mu}_{1}=-\bm{\mu}_{2}\equiv\bm{\mu}^{\star}$ and
$\bm{\Sigma}_{1}=\bm{\Sigma}_{2}=\sigma\mathbb{I}_{d}$. For definiteness, we
assume $\bm{\mu}\sim\mathbb{P}_{\mu}$, with
$\mathbb{P}_{\mu}=\mathcal{N}(0,\nicefrac{{\mathbb{I}_{d}}}{{d}})$, so that
with high probability $\lVert\bm{\mu}\lVert=1$. We shall moreover assume for
ease of discussion that $\sigma$ is perfectly known. Thus, the centroid
$\bm{\mu}$ is the only unknown parameter.

###  C.2 Bayes-optimal denoiser without knowledge of cluster means

The corresponding Bayes-optimal denoiser then reads

$\displaystyle\mathfrak{b}(\tilde{\bm{x}})$
$\displaystyle\equiv\mathbb{E}\left[\bm{x}|\tilde{\bm{x}},\sigma,\mathcal{D}\right]=\int
d\bm{\mu}\underbrace{\mathbb{E}\left[\bm{x}|\tilde{\bm{x}},\bm{\mu},\sigma,\mathcal{D}\right]}_{t(\bm{\mu},\tilde{\bm{x}})}\mathbb{P}(\bm{\mu}|\sigma,\mathcal{D})$
(85)

Where we have identified the oracle denoiser (82), and emphasized the
dependence on $\bm{\mu}$. Note that this is a slight abuse of notation, as the
true oracle $t(\bm{\mu}^{\star},\tilde{\bm{x}})$ involves the ground truth
centroid. Further note that
$\mathbb{P}(\bm{\mu}|\sigma,\mathcal{D})=\mathbb{P}(\bm{\mu}|\mathcal{D})$,
since knowledge of $\sigma$ does not bring information about $\bm{\mu}$.
Furthermore, note that as the noisy part of the dataset
$\tilde{\mathcal{D}}\equiv\\{\tilde{\bm{x}}^{\mu}\\}_{\mu=1}^{n}$ brings only
redundant information, one further has that

$\displaystyle\mathbb{P}(\bm{\mu}|\mathcal{D})=\mathbb{P}(\bm{\mu}|\mathcal{\check{D}}),$
(86)

where we noted $\mathcal{\check{D}}=\\{\bm{x}^{\mu}\\}_{\mu=1}^{n}$ the clean
part of the dataset. Thus

$\displaystyle\mathfrak{b}(\tilde{\bm{x}})$
$\displaystyle=\frac{1}{\mathbb{P}(\mathcal{\check{D}})}\int
d\bm{\mu}~{}~{}t(\bm{\mu},\tilde{\bm{x}})\mathbb{P}(\mathcal{\check{D}}|\bm{\mu})\mathbb{P}_{\mu}(\bm{\mu})$
$\displaystyle=\frac{1}{\mathbb{P}(\mathcal{\check{D}})}\int
d\bm{\mu}~{}~{}t(\bm{\mu},\tilde{\bm{x}})\mathbb{P}_{\mu}(\bm{\mu})\prod\limits_{\mu=1}^{n}\frac{1}{2}\frac{1}{(2\pi\sigma^{2})^{\nicefrac{{d}}{{2}}}}\left[e^{-\frac{1}{2\sigma^{2}}\lVert\bm{x}^{\mu}-\bm{\mu}\lVert^{2}}+e^{-\frac{1}{2\sigma^{2}}\lVert\bm{x}^{\mu}+\bm{\mu}\lVert^{2}}\right]$
$\displaystyle=\frac{1}{\mathbb{P}(\mathcal{\check{D}})}\int
d\bm{\mu}~{}~{}t(\bm{\mu},\tilde{\bm{x}})\frac{1}{(\nicefrac{{2\pi}}{{d}})^{\nicefrac{{d}}{{2}}}}e^{-\frac{d}{2}\lVert\bm{\mu}\lVert^{2}}e^{-\frac{n}{2\sigma^{2}}\lVert\bm{\mu}\lVert^{2}}\prod\limits_{\mu=1}^{n}\frac{e^{-\frac{1}{2\sigma^{2}}\lVert\bm{x}^{\mu}\lVert^{2}}}{(2\pi\sigma^{2})^{\nicefrac{{d}}{{2}}}}\cosh\left(\frac{\bm{\mu}^{\top}\bm{x}^{\mu}}{\sigma^{2}}\right)$
$\displaystyle=\frac{1}{\sqrt{d}\mathbb{P}(\mathcal{\check{D}})}\int
d\bm{\mu}~{}~{}t(\nicefrac{{\bm{\mu}}}{{\sqrt{d}}},\tilde{\bm{x}})\frac{1}{(\nicefrac{{2\pi}}{{d}})^{\nicefrac{{d}}{{2}}}}e^{-\frac{1}{2}\lVert\bm{\mu}\lVert^{2}}e^{-\frac{\alpha}{2\sigma^{2}}\lVert\bm{\mu}\lVert^{2}}\prod\limits_{\mu=1}^{n}\frac{e^{-\frac{1}{2\sigma^{2}}\lVert\bm{x}^{\mu}\lVert^{2}}}{(2\pi\sigma^{2})^{\nicefrac{{d}}{{2}}}}\cosh\left(\frac{\bm{\mu}^{\top}\bm{x}^{\mu}}{\sigma^{2}\sqrt{d}}\right)$
$\displaystyle=\frac{1}{Z}\int
d\bm{\mu}~{}~{}t(\nicefrac{{\bm{\mu}}}{{\sqrt{d}}},\tilde{\bm{x}})\mathbb{P}_{b}(\bm{\mu}),$
(87)

where we noted

$\displaystyle\mathbb{P}_{b}(\bm{\mu})\equiv
e^{-\frac{1}{2\tilde{\sigma}^{2}}\lVert\bm{\mu}\lVert^{2}}\prod\limits_{\mu=1}^{n}\cosh\left(\frac{\bm{\mu}^{\top}\bm{x}^{\mu}}{\sigma^{2}\sqrt{d}}\right),$
(88)

with the effective variance

$\displaystyle\hat{\sigma}^{2}\equiv\frac{\sigma^{2}}{\sigma^{2}+\alpha}.$
(89)

Finally, the partition function $Z$ is defined as the normalisation of
$\mathbb{P}_{b}$, i.e.

$\displaystyle Z\equiv\int d\bm{\mu}~{}~{}\mathbb{P}_{b}(\bm{\mu}).$ (90)

One is now in a position to compute the MSE associated with the Bayes
estimator $\mathfrak{b}(\tilde{\bm{x}})$:

$\displaystyle\mathrm{mse}_{\mathfrak{b}}$
$\displaystyle=\mathbb{E}_{\mathcal{D}}\mathbb{E}_{\bm{x}\sim\mathcal{P}}\mathbb{E}_{\bm{\xi}\sim\mathcal{N}(0,\mathbb{I}_{d})}\left\lVert\bm{x}-\left\langle
t(\nicefrac{{\bm{\mu}}}{{\sqrt{d}}},\sqrt{1-\Delta}\bm{x}+\sqrt{\Delta}\bm{\xi})\right\rangle_{\bm{\mu}\sim\mathbb{P}_{b}}\right\lVert^{2}$
$\displaystyle=\mathrm{mse}_{\circ}-\mathbb{E}_{\mathcal{D}}\sum\limits_{s=\pm
1}\mathbb{E}_{\bm{\xi},\bm{\eta}\sim\mathcal{N}(0,\mathbb{I}_{d})}\frac{\Delta}{\sigma^{2}(1-\Delta)+\Delta}\Bigg{\langle}\left({\scriptstyle\frac{(1-b\sqrt{1-\Delta})s\bm{\mu}^{\top}\bm{\mu}^{\star}}{d}+\frac{\bm{\mu}^{\top}((1-b\sqrt{1-\Delta})\sigma\bm{\eta}-b\sqrt{\Delta}\bm{\xi}}{\sqrt{d}}}\right)$
$\displaystyle\qquad\tanh\left({\scriptstyle\frac{\sqrt{1-\Delta}\frac{s\bm{\mu}^{\top}\bm{\mu}^{\star}}{2}+\sqrt{1-\Delta}\sigma\frac{\bm{\mu}^{\top}\bm{\eta}}{\sqrt{d}}+\sqrt{\Delta}\frac{\bm{\mu}^{\top}\bm{\xi}}{\sqrt{d}}}{\sigma^{2}(1-\Delta)+\Delta}}\right)\Bigg{\rangle}_{\bm{\mu}\sim\mathbb{P}_{b}}$
$\displaystyle+\mathbb{E}_{\mathcal{D}}\frac{1}{2}\sum\limits_{s=\pm
1}\mathbb{E}_{\bm{\xi},\bm{\eta}\sim\mathcal{N}(0,\mathbb{I}_{d})}\Bigg{\langle}\left(\frac{\Delta}{\sigma^{2}(1-\Delta)+\Delta}\right)^{2}\frac{\bm{\mu}_{1}^{\top}\bm{\mu}_{2}}{d}\tanh\left({\scriptstyle\frac{\sqrt{1-\Delta}\frac{s\bm{\mu}_{1}^{\top}\bm{\mu}^{\star}}{d}+\sqrt{1-\Delta}\sigma\frac{\bm{\mu}_{1}^{\top}\bm{\eta}}{\sqrt{d}}+\sqrt{\Delta}\frac{\bm{\mu}_{1}^{\top}\bm{\xi}}{\sqrt{d}}}{\sigma^{2}(1-\Delta)+\Delta}}\right)$
$\displaystyle\qquad\times\tanh\left({\scriptstyle\frac{\sqrt{1-\Delta}\frac{s\bm{\mu}_{2}^{\top}\bm{\mu}^{\star}}{d}+\sqrt{1-\Delta}\sigma\frac{\bm{\mu}_{2}^{\top}\bm{\eta}}{\sqrt{d}}+\sqrt{\Delta}\frac{\bm{\mu}_{2}^{\top}\bm{\xi}}{\sqrt{d}}}{\sigma^{2}(1-\Delta)+\Delta}}\right)\Bigg{\rangle}_{\bm{\mu}_{1},\bm{\mu}_{2}\sim\mathbb{P}_{b}}$
$\displaystyle=\mathrm{mse}_{\circ}-\mathbb{E}_{\mathcal{D}}\sum\limits_{s=\pm
1}\mathbb{E}_{z\sim\mathcal{N}(0,1)}\frac{\Delta}{\sigma^{2}(1-\Delta)+\Delta}\Bigg{\langle}\left({\scriptstyle\frac{(1-b\sqrt{1-\Delta})s\bm{\mu}^{\top}\bm{\mu}^{\star}}{d}}\right)\tanh\left({\scriptstyle\frac{\sqrt{1-\Delta}\frac{s\bm{\mu}^{\top}\bm{\mu}^{\star}}{d}+\sqrt{(1-\Delta)\sigma^{2}+\Delta}\frac{\lVert\bm{\mu}\lVert}{\sqrt{d}}z}{\sigma^{2}(1-\Delta)+\Delta}}\right)\Bigg{\rangle}_{\bm{\mu}\sim\mathbb{P}_{b}}$
$\displaystyle\qquad+\mathbb{E}_{\mathcal{D}}\frac{1}{2}\sum\limits_{s=\pm
1}\mathbb{E}_{u,v\sim\mathcal{N}(0,\Omega)}\Bigg{\langle}\left(\frac{\Delta}{\sigma^{2}(1-\Delta)+\Delta}\right)^{2}\frac{\bm{\mu}_{1}^{\top}\bm{\mu}_{2}}{d}\tanh\left({\scriptstyle\frac{\sqrt{1-\Delta}\frac{s\bm{\mu}_{1}^{\top}\bm{\mu}^{\star}}{d}+\sqrt{(1-\Delta)\sigma^{2}+\Delta}u}{\sigma^{2}(1-\Delta)+\Delta}}\right)$
$\displaystyle\qquad\times\tanh\left({\scriptstyle\frac{\sqrt{1-\Delta}\frac{s\bm{\mu}_{2}^{\top}\bm{\mu}^{\star}}{d}+\sqrt{(1-\Delta)\sigma^{2}+\Delta}v}{\sigma^{2}(1-\Delta)+\Delta}}\right)\Bigg{\rangle}_{\bm{\mu}_{1},\bm{\mu}_{2}\sim\mathbb{P}_{b}}.$
(91)

We have adopted the shortcuts

$\displaystyle
b=\frac{\sqrt{1-\Delta}\sigma^{2}}{\sigma^{2}(1-\Delta)+\Delta},$
$\displaystyle\Omega=\begin{bmatrix}\frac{\lVert\bm{\mu}_{1}\lVert^{2}}{d}&\frac{\bm{\mu}_{1}^{\top}\bm{\mu}_{2}}{d}\\\
\frac{\bm{\mu}_{1}^{\top}\bm{\mu}_{2}}{d}&\frac{\lVert\bm{\mu}_{2}\lVert^{2}}{d}.\end{bmatrix}.$
(92)

(C.2) shows that the inner averages involved in the Bayes MSE
$\mathrm{mse}_{\mathfrak{b}}$ only depend on the overlaps
$\nicefrac{{\bm{\mu}_{1,2}^{\top}\bm{\mu}_{1,2}}}{{d}}$ for
$\bm{\mu}_{1},\bm{\mu}_{2}$ two independently sampled vectors from
$\mathbb{P}_{b}$. Motivated by similar high-dimensional studies, it is
reasonable to expect these quantities to concentrate as $n,d\to\infty$ in the
measure $\mathbb{E}_{\mathcal{D}}\langle\cdot\rangle_{\mathbb{P}_{b}}$. A more
principle – but much more painstaking– way to derive this assumption is to
introduce a Fourier representation, and carry out the
$\mathbb{E}_{\mathcal{D}}$ average using the replica method. We refer the
interested reader to, for instance, Appendix H of [36], where such a route is
taken.

###  C.3 Derivation

To compute the summary statistics of $\mathbb{P}_{b}$, we again resort to the
replica method to compute the moment-generating function (free entropy), see
also Appendix A. The replicated partition function reads

$\displaystyle\mathbb{E}_{\mathcal{D}}Z^{s}=\int\prod\limits_{a=1}^{s}d\bm{\mu}_{a}e^{-\frac{1}{2\hat{\sigma}^{2}}\lVert\bm{\mu}_{a}\lVert^{2}}\prod\limits_{\mu=1}^{n}\sum\limits_{s=\pm
1}\frac{1}{2}\mathbb{E}_{\bm{\eta}\sim\mathcal{N}(0,\mathbb{I}_{d})}\left[\prod\limits_{a=1}^{s}\cosh\left(\frac{s\frac{\bm{\mu}_{a}^{\top}\bm{\mu}^{\star}}{d}+\sigma\frac{\bm{\mu}_{a}^{\top}\bm{\eta}}{\sqrt{d}}}{\sigma^{2}}\right)\right].$
(93)

Introducing the order parameters

$\displaystyle q_{ab}\equiv\frac{\bm{\mu}_{a}^{\top}\bm{\mu}_{b}}{d},$
$\displaystyle m_{a}\equiv\frac{\bm{\mu}_{a}^{\top}\bm{\mu}^{\star}}{d},$ (94)

one reaches

$\displaystyle\mathbb{E}_{\mathcal{D}}Z^{s}=$
$\displaystyle\int\prod\limits_{a\leq
b}dq_{ab}d\hat{q}_{ab}\prod\limits_{a=1}^{s}dm_{a}d\hat{m}_{a}\underbrace{e^{-d\sum\limits_{a\leq
b}q_{ab}\hat{q}_{ab}-d\sum\limits_{a}m_{a}\hat{m}_{a}}}_{e^{sd\Psi_{t}}}\underbrace{\int\prod\limits_{a=1}^{s}d\bm{\mu}_{a}e^{-\frac{1}{2\hat{\sigma}^{2}}\lVert\bm{\mu}_{a}\lVert^{2}}e^{\sum\limits_{a\leq
b}\hat{q}_{ab}\bm{\mu}_{a}^{\top}\bm{\mu}_{b}+\sum\limits_{a}\bm{\mu}_{a}^{\top}\bm{\mu}^{\star}}}_{e^{sd\Psi_{w}}}$
$\displaystyle\underbrace{\left[n\sum\limits_{s=\pm
1}\frac{1}{2}\int\mathcal{N}\left(\\{\lambda_{a}\\}_{a=1}^{s};0,(q_{ab})_{1\leq
a,b\leq
s}\right)\left[\prod\limits_{a=1}^{s}\cosh\left(\frac{sm_{a}+\sigma\lambda_{a}}{\sigma^{2}}\right)\right]\right]^{n}}_{e^{s\alpha
d\Psi_{y}}}.$ (95)

The computation of $\Psi_{t},\Psi_{w},\Psi_{y}$ is rather standard and follows
the main steps presented in Appendix A. One again assumes the replica
symmetric ansatz [21, 22]

$\displaystyle q_{ab}=(r-q)\delta_{ab}+q,$ $\displaystyle m_{a}=m,$ (96)
$\displaystyle\hat{q}_{ab}=\left(-\frac{\hat{r}}{2}+-\hat{q}\right)\delta_{ab}+\hat{q},$
$\displaystyle\hat{m}_{a}=\hat{m}.$ (97)

Introducing as in Appendix A the variances

$\displaystyle V\equiv r-q,$ $\displaystyle\hat{V}\equiv\hat{r}+\hat{q},$ (98)

one reaches

$\displaystyle\Psi_{t}=\frac{\hat{V}q-V\hat{q}}{2}-m\hat{m},$ (99)
$\displaystyle\Psi_{w}=-\frac{1}{2}\ln\left(1+\hat{V}\hat{\sigma}^{2}\right)+\frac{\hat{\sigma}^{2}}{2}\frac{\hat{q}+\hat{m}^{2}}{1+\hat{V}\hat{\sigma}^{2}},$
(100)
$\displaystyle\Psi_{y}=\mathbb{E}_{\eta\sim\mathcal{N}(0,1)}\ln\left[\int
d\lambda\mathcal{N}(\lambda;\sqrt{q}\eta,V)\cosh\left(\frac{\sigma\lambda+m}{\sigma^{2}}\right)\right].$
(101)

Therefore the free entropy reads

$\displaystyle\Phi=\underset{q,m,V,\hat{q},\hat{m},\hat{V}}{\mathrm{extr}}$
$\displaystyle\frac{\hat{V}q-V\hat{q}}{2}-m\hat{m}-\frac{1}{2}\ln\left(1+\hat{V}\hat{\sigma}^{2}\right)+\frac{\hat{\sigma}^{2}}{2}\frac{\hat{q}+\hat{m}^{2}}{1+\hat{V}\hat{\sigma}^{2}}$
$\displaystyle+\alpha\mathbb{E}_{\eta\sim\mathcal{N}(0,1)}\ln\left[\int
d\lambda\mathcal{N}(\lambda;\sqrt{q}\eta,V)\cosh\left(\frac{\sigma\lambda+m}{\sigma^{2}}\right)\right].$
(102)

The solution of this extremization problem is given by the set of self-
consistent equations, imposing that gradients with respect to all parameters
be zero:

$\displaystyle\begin{cases}V=\frac{\hat{\sigma}^{2}}{1+\hat{V}\hat{\sigma}^{2}}\\\
q=\frac{\hat{\sigma}^{4}(\hat{q}+\hat{m}^{2})}{(1+\hat{V}\hat{\sigma}^{2})^{2}}+\frac{\hat{\sigma}^{2}}{1+\hat{V}\hat{\sigma}^{2}}\\\
m=\frac{\hat{\sigma}^{2}\hat{m}}{1+\hat{V}\hat{\sigma}^{2}}\end{cases}$
$\displaystyle\begin{cases}\hat{V}=-\alpha\frac{1}{\sigma\sqrt{q}}\mathbb{E}_{\eta\sim\mathcal{N}(0,1)}\left[\tanh\left(\frac{\sigma\sqrt{q}\eta+m}{\sigma^{2}}\right)\eta\right]\\\
\hat{q}=\frac{\alpha V}{\sigma^{2}}\\\
\hat{m}=\frac{\alpha}{\sigma^{2}}\mathbb{E}_{\eta\sim\mathcal{N}(0,1)}\left[\tanh\left(\frac{\sigma\sqrt{q}\eta+m}{\sigma^{2}}\right)\right]\end{cases}$
(103)

This is a simple optimization problem over $6$ variables, which can be solved
numerically. One can finally use the obtained summary statistics to evaluate
the Bayes optimal MSE:

$\displaystyle\mathrm{mse}_{\mathfrak{b}}$
$\displaystyle=\mathrm{mse}_{\circ}-2\mathbb{E}_{\mathcal{D}}\mathbb{E}_{z\sim\mathcal{N}(0,1)}\frac{\Delta}{\sigma^{2}(1-\Delta)+\Delta}\left({\scriptstyle(1-b\sqrt{1-\Delta})m}\right)\tanh\left({\scriptstyle\frac{\sqrt{1-\Delta}m+\sqrt{(1-\Delta)\sigma^{2}+\Delta}\sqrt{q+V}z}{\sigma^{2}(1-\Delta)+\Delta}}\right)$
$\displaystyle+\mathbb{E}_{\mathcal{D}}\mathbb{E}_{u,v\sim\mathcal{N}(0,\Omega)}\left({\scriptstyle\frac{\Delta}{\sigma^{2}(1-\Delta)+\Delta}}\right)^{2}q\tanh\left({\scriptstyle\frac{\sqrt{1-\Delta}m+\sqrt{(1-\Delta)\sigma^{2}+\Delta}u}{\sigma^{2}(1-\Delta)+\Delta}}\right)\tanh\left({\scriptstyle\frac{\sqrt{1-\Delta}m+\sqrt{(1-\Delta)\sigma^{2}+\Delta}v}{\sigma^{2}(1-\Delta)+\Delta}}\right)$
(104)

where

$\displaystyle\Omega=\begin{bmatrix}q+V&q\\\ q&q+V\end{bmatrix}$ (105)

Figure 9: (Binary isotropic mixture, $\bm{\Sigma}_{1,2}=0.09\mathbb{I}_{d}$,
$\rho_{1,2}=\nicefrac{{1}}{{2}}$, and $\bm{\mu}_{1}=-\bm{\mu}_{2}$, see also
Fig. 1). (right) MSE (minus the rescaling baseline $\mathrm{mse}_{\circ}$) for
the DAE ($p=1,\sigma=\tanh$) (3) (blue), the oracle Tweedie denoiser (79)
(green), the Bayes denoiser (C.3) (red), and the plug-in oracle-PCA denoiser
(106) (black). Solid lines correspond to the theoretical formulae II.3, (B.2)
and (C.3); simulations correspond to a DAE optimized using the Pytorch
implementation of Adam (weight decay $\lambda=0.1$, learning rate $\eta=0.05$,
full batch, $2000$ epochs, sample complexity $\alpha=1$), averaged over $N=10$
instances, with error bars representing one standard deviation. (right) Cosine
similarity (7) for the same algorithms, with identical color code.

This completes the derivation of the MSE of the Bayes estimator agnostic of
the cluster means.

###  C.4 A simple plug-in denoiser

Another simple denoiser assuming only perfect knowledge of the variance
$\sigma$, but imperfect knowledge of $\bm{\mu}$, is simply to plug the PCA
estimate $\hat{\bm{\mu}}_{\mathrm{PCA}}$ of $\bm{\mu}$ into the Tweedie oracle
(79), i.e.

$\displaystyle
t(\hat{\bm{\mu}}_{\mathrm{PCA}},\tilde{\bm{x}})=\frac{\sqrt{1-\Delta}\sigma^{2}}{\sigma^{2}(1-\Delta)+\Delta}\tilde{\bm{x}}+\frac{\Delta}{\sigma^{2}(1-\Delta)+\Delta}\tanh\left(\frac{\tilde{\bm{x}}^{\top}\hat{\bm{\mu}}_{\mathrm{PCA}}\sqrt{1-\Delta}}{\sigma^{2}(1-\Delta)+\Delta}\right)\times\hat{\bm{\mu}}_{\mathrm{PCA}}.$
(106)

The performance of this plug-in denoiser is plotted in Fig. 9, and contrasted
to the one of the Bayes denoiser $\mathfrak{b}(\cdot)$ (C.3), the oracle
$t(\cdot)$ (B.2) and the DAE II.3. Strikingly, the performance of the PCA
plug-in denoiser is sizeably identical to the one of the Bayes-optimal
denoiser, both in terms of MSE and cosine similarity. It is important to
remind at this point that both the plug-in (106) and the Bayes denoiser (C.3)
still rely on perfect knowledge of the cluster variance $\sigma$ and the noise
level $\Delta$, while the DAE (3) is agnostic to these parameters. Yet, these
two denoisers make for a fairer comparison than the oracle (B.2), as they
importantly take into account the finite training set. The DAE is relatively
close to the Bayes baseline for small noise levels $\Delta$, while a gap can
be seen to open up for larger $\Delta$.

##  Appendix D Details on real data simulations

In this Appendix, we provide further details on the real data experiments
presented in Fig. 2 and Fig. 4.

#### Preprocessing

The original MNIST [15] and FashionMNIST [16] data-sets were flattened
(vectorized), centered, and rescaled by $400$ (MNIST) and $600$
(FashionMNIST). For simplicity and ease of discussion, for each data set, only
two labels were kept, namely $1$s and $7$s for MNIST and boots and shoes for
FashionMNIST. Note that the visual similarity between the two selected classes
should make the denoising problem more challenging.

#### Means and covariance estimation

For each data-set, data from the same class were considered to belong to the
same Gaussian mixture cluster. For simplicity, we kept the same number of
training points in each cluster, so as to obtain balanced clusters
($\rho_{1}=\rho_{2}=\nicefrac{{1}}{{2}}$) for definiteness. For each cluster,
the corresponding mean $\mu$ and covariance $\Sigma$ were numerically
evaluated from the empirical mean and covariance over the $6000$ boots (shoes)
in the FashionMNIST training set, and the $6265$ $1$s ($7$s) in the MNIST
training set. The solid lines in Fig. 2 correspond to using those estimates in
the asymptotic formulae of Result II.3, in their generic form (61) (see
Appendix A) for sample complexity $\alpha=1$ and $\ell_{2}$ regularization
$\lambda=0.1$.

####  Pytorch simulations

The red points in Fig. 2 were obtained from numerical simulations, by
optimizing a $p=1$ DAE (3) with $\sigma=\tanh$ activation using the Pytorch
implementation of the full-batch Adam [43] optimizer, with weight decay
$\lambda=0.1$, learning rate $\eta=0.05$, over $T=2000$ epochs. For each value
of $\Delta$, $n=784$ images were drawn without replacement from the data-set
training set, corrupted by noise, and fed to the DAE along with the clean
image for training. The denoising test MSE, cosine similarity and summary
statistics of the trained DAE were estimated after optimization using
$n_{\mathrm{test}}=1000$ images randomly drawn from the test set (from which
also only the relevant classes were kept, and equally represented– i.e.
balanced). Each value was averaged over $N=10$ instances of the train set,
noise realization and test set. The obtained values were furthermore found to
be robust with respect to the choice of the optimizer, learning rate and
number of epochs.

###  D.1 Complementary simulations

Figure 10: Cosine similarity $\theta$ (7) (green), weight norm
$\nicefrac{{\lVert\hat{\bm{w}}\lVert^{2}}}{{d}}$ and skip connection strength
$\hat{b}$ for the MNIST dataset (left), of which for simplicity only $1$s and
$7$s were kept, and FashionMNIST (right), of which only boots and shoes were
kept. In solid lines, the theoretical predictions resulting from using Result
II.3 in its generic formulation (61) with the empirically estimated
covariances and means. Dots represent numerical simulations of a DAE
($p=1,\sigma=\tanh$) trained with $n=784$ training points, using the Pytorch
implementation of full-batch Adam, with learning rate $\eta=0.05$ over $2000$
epochs, weight decay $\lambda=0.1$, averaged over $N=5$ instances. Error bars
represent one standard deviation. See Appendix D for full details on the
preprocessing.

Fig.2 in the main text addresses the comparison of the theoretical prediction
of Result II.3 with simulations on the real data-set, at the level of the
denoising test MSE (6). For completeness, we show here the same comparison for
the other learning metrics and summary statistics, namely the cosine
similarity $\theta$ (7), weights norm
$\nicefrac{{\lVert\hat{\bm{w}}\lVert^{2}}}{{d}}$ and trained skip connection
strength $\hat{b}$. Theoretical asymptotic characterization are again provided
by Result II.3 (11), (12) and (II.3). They are plotted as solid lines in Fig.
10, and contrasted to numerical simulations (dots) on the real data-sets, by
training a $p=1$ DAE (3) using the Pytorch implementation of the Adam