# Optimal errors and phase transitions in high-dimensional generalized linear models

Jean Barbier, Florent Krzakala, Nicolas Macris, Léo Miolane, Lenka Zdeborová

► To cite this version:

Jean Barbier, Florent Krzakala, Nicolas Macris, Léo Miolane, 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, 2019, 116 (12), pp.5451-5460. 10.1073/pnas.1802705116. cea-01614258

HAL Id: cea-01614258

https://cea.hal.science/cea-01614258

Submitted on 21 Nov 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<table>
  <tbody>
    <tr>
      <td><strong>Appendix B Some properties of the scalar channels</strong></td>
      <td>77</td>
    </tr>
    <tr>
      <td>B.1 The additive Gaussian scalar channel</td>
      <td>77</td>
    </tr>
    <tr>
      <td>B.2 The non-linear scalar channel</td>
      <td>77</td>
    </tr>
    <tr>
      <td><strong>Appendix C Approximation</strong></td>
      <td>84</td>
    </tr>
    <tr>
      <td>C.1 Relaxing the hypotheses on P<sub>0</sub> and Φ</td>
      <td>84</td>
    </tr>
    <tr>
      <td>C.2 Relaxing the hypotheses on φ</td>
      <td>85</td>
    </tr>
    <tr>
      <td>C.3 The case of discrete channels: Removing the Gaussian noise</td>
      <td>87</td>
    </tr>
    <tr>
      <td><strong>Appendix D Some sup-inf formulas</strong></td>
      <td>87</td>
    </tr>
    <tr>
      <td><strong>Appendix E Concentration of free entropy and overlaps</strong></td>
      <td>90</td>
    </tr>
    <tr>
      <td>E.1 Concentration of the free entropy</td>
      <td>90</td>
    </tr>
    <tr>
      <td>E.2 Concentration of the overlap</td>
      <td>95</td>
    </tr>
    <tr>
      <td><strong>Appendix F Details on numerics</strong></td>
      <td>99</td>
    </tr>
    <tr>
      <td>F.1 General purpose algorithms</td>
      <td>99</td>
    </tr>
    <tr>
      <td>F.2 Evaluating the replica formula</td>
      <td>100</td>
    </tr>
    <tr>
      <td>F.3 Breaking the symmetry in GAMP</td>
      <td>100</td>
    </tr>
  </tbody>
</table>

# Part I

# Main part

## 1 Introduction

As datasets grow larger and more complex, modern data analysis requires solving high-dimensional estimation problems with very many parameters. Developing algorithms up to the task and understanding their limitations has become a major challenge in computer science, machine learning, statistics, signal processing, communications and related fields.

In the present contribution, we address this challenge in the case of generalized linear estimation models (GLMs) [1,2] where data are generated as follows: Given a *n*-dimensional vector **X**\*, hidden to the statistician, he/she observes instead a *m*-dimensional vector **Y** where each component reads

$$Y_{\mu} = \varphi\left(\frac{1}{\sqrt{n}}[\mathbf{\Phi}\mathbf{X}^{*}]_{\mu}, A_{\mu}\right), \quad 1 \le \mu \le m, \qquad (1)$$

where **Φ** is a *m* × *n* “measurement” or “data” matrix, the random variables ($A_\mu$) $\stackrel{\text{iid}}{\sim} P_A$ account for noise (or randomness) of the model. The model is “linear” because the output $Y_\mu$ depends on a *linear* combination of the data $z_\mu = \frac{1}{\sqrt{n}}[\mathbf{\Phi}\mathbf{X}^*]_\mu = \frac{1}{\sqrt{n}}\sum_{i=1}^n \Phi_{\mu i} X_i^*$. The GLM generalizes the ordinary linear regression by allowing the output function $\varphi(z, A)$ to be non-linear and/or stochastic; in the case of a deterministic model we simply write $\varphi(z)$. Explicit examples will be given below.

GLMs belong to the realm of supervised learning and arise in a wide variety of scientific fields. In signal processing one usually observes $Y_\mu$ given as a linear combination of the signal-elements **X**\*. In a range of applications these observations are obtained via a non-linear function $\varphi$. In optics or X-ray crystallography one often measures only the amplitude of $[\mathbf{\Phi}\mathbf{X}^*]_\mu$, leading to the phase retrieval problem [3]. A real-valued analog is the problem of sign-retrieval when we only observe $|[\mathbf{\Phi}\mathbf{X}^*]_\mu|$ [4, 5]. Observations are sometimes quantized in order to reduce the storage, leading for instance to the problem of 1-bit compressed sensing [6]. In statistics and machine learning, classification is often described via a GLM where the output function $\varphi$
is discrete and corresponds to the labels that classify the data-points $\mathbf{\Phi}_{\mu}$ [1, 2, 7]. GLMs with non-linear output functions are also the basic building blocks of each layer of neural networks [8]: $\varphi$ corresponds to the activation, the rows of the matrix $\mathbf{\Phi}$ are different data samples while $\mathbf{X}^*$ are the set of synaptic weights to be learned.

There are two main learning problems in GLMs: i) The estimation task requires, knowing the measured vector $\mathbf{Y}$ and the matrix $\mathbf{\Phi}$, to infer the unknown vector $\mathbf{X}^*$; ii) the prediction or generalization task instead requires, again knowing $\mathbf{Y}$ and $\mathbf{\Phi}$, to predict accurately new values $Y_{\text{new}}$ when new rows (i.e. data-points) are added to the matrix $\mathbf{\Phi}$.

In the present paper we build a rigorous theory for both these tasks for random instances of the GLM. In this setting each element $\mathbf{\Phi}_{\mu i}$ of the matrix is sampled independently from a probability distribution of zero mean and unit variance, and the unknown vector $\mathbf{X}^*$ has been also created randomly from a probability distribution $P_0$, with each of its components $X_1^*, \dots, X_n^* \stackrel{\text{iid}}{\sim} P_0$. Since our main aim is to study the intrinsic information-theoretic and algorithmic limitations caused by the lack of samples and/or the amplitude of the noise, we assume throughout this paper that $P_0$ and $\varphi$ are known to the statistician (if they are not the task can only be harder). Our results are derived in the challenging and interesting high-dimensional limit where $m, n \to \infty$ while $m/n \to \alpha$ a constant. Random instances of GLMs are both practically and theoretically relevant in many different contexts:

a) In signal processing, GLM estimation with a random matrix $\mathbf{\Phi}$ has been studied with considerable attention in the context of compressed sensing [9–11] where a $n$-dimensional sparse signal is recovered from $m < n$ noisy measurements. While standard compressed sensing focused on the linear case –where $\varphi(z, A) = z + A$ with a Gaussian noise $A$– the generalized case was also widely studied [12, 13], especially for quantized output [14] and 1-bit compressed sensing [6, 15] where $\varphi(z, A) = \text{sign}(z + A)$, as well as for compressive phase retrieval when $\varphi(z, A) = |z + A|$ [16].

b) In statistical learning, important activity is dedicated to understand the limitation of learning with data generated by GLMs, both in the linear case, e.g. in the context of ridge regression or LASSO [17], or with non-linear probabilistic output, e.g. logistic regression. Random instances were studied in particular in the context of so-called M-estimators [18–21].

c) In studies of artificial neural networks there has been a large amount of works using random instances of GLMs, with $\varphi$ playing the role of a non-linear activation function. In this context the random GLM was introduced as the teacher-student setting for the perceptron in the pioneering work of Gardner and Derrida [22]. Large volume of work followed and is reviewed, e.g., in [23–25]. While initial works concentrated on a simple activation functions $\varphi(z) = \text{sign}(z - K)$ ($K$ is the threshold constant), many other functions were considered, e.g. in [26–28]. Recently, the study of random instances of neural networks have emerged as a key ingredient in understanding the performance of deep learning algorithms [29, 30]. Computing mutual informations in GLMs is also a critical issue in confirming the information bottleneck scenario of [31, 32].

d) In communications, error-correcting codes that use random constructions are particularly efficient, as discussed by Shannon in his seminal paper [33]. Random instances of GLMs describe both the setting of code-division multiple access –a multi-user access method used in communication technologies [34, 35]– as well as an error correction scheme called sparse superposition codes, that have been shown to achieve the Shannon capacity for any type of noisy channel [36–40].

Interestingly there is an important gap in the above volume of work. On the one hand there are studies that rely on the algorithmic performance of the so-called generalized approximate message-passing algorithm (GAMP) [11, 12, 41]. GAMP is remarkable in that its asymptotic $(n, m \to \infty, m/n \to \alpha)$ performance can be analyzed rigorously using the so-called state evolution [42–45]. However, GAMP is not expected to be always information-theoretically optimal. On the other hand, other results are concerned with the linear case of the GLM with additive Gaussian noise for which the information-theoretically optimal performance was established in [46–48] (the methodology of these works unfortunately does not generalize straightforwardly
to the important non-linear case or to other types of additive noise). All the other works, giving information-theoretic results for the non-linear case, are based on powerful and sophisticated but non-rigorous techniques originating in statistical physics of disordered systems, such as the cavity and replica methods [49]. Historically, the first of these non-rigorous, yet correct, results on information-theoretic limitations of learning was for the perceptron with binary weights and was established using the replica method in [22, 50, 51], including a discontinuous phase transition to perfect learning that appears as the ratio between number of samples and the dimension exceeds $\alpha \approx 1.249$.

In the present paper we close the above gap between mathematically rigorous work and conjectures (some of them several decades old) from statistical mechanics. In particular, we prove that the results for GLMs stemming from the replica method are indeed correct and imply the optimal value of both the estimation and generalization error. These results are summarized in section “Main results”. The proof is based on the adaptive interpolation method recently developed in [52] and is of independent interest as it is applicable to a range of other models, see section “Methods and proofs” and the supplementary informations (SI). We compare our information-theoretic results to the performance of the GAMP algorithm and its state evolution (as reviewed briefly in section “Main results”). We determine regions of parameters where this algorithm is or is not information-theoretically optimal. Up to technical assumptions (as specified below), our results apply to all activation functions $\phi$ and priors $P_0$, thus unifying a large volume of previous work where many particular functions have been analyzed on a case by case basis. This generality allows us to provide a unifying understanding of the types of phase transitions and phase diagrams that we can encounter in GLMs, which is as well of independent interest and we devote section “Application to learning and inference” to its presentation.

## 2 Main results

This section summarizes our main results. Their formal statement together with all technical assumptions and full proofs are provided in section “Methods and proofs” and in the SI.

For the random GLM problem as defined in the introduction, the optimal way to estimate the ground-truth signal/weights $\mathbf{X}^*$ relies on its posterior probability distribution

$$P(\mathbf{x}|\mathbf{Y}, \mathbf{\Phi}) = \frac{1}{\mathcal{Z}(\mathbf{Y}, \mathbf{\Phi})} \prod_{i=1}^{n} P_0(x_i) \prod_{\mu=1}^{m} P_{\text{out}}\left(Y_{\mu} \middle| \frac{[\mathbf{\Phi}\mathbf{x}]_{\mu}}{\sqrt{n}}\right) \quad (2)$$

where we used the prior $P_0$ of $\mathbf{X}^*$, and introduced the likelihood $P_{\text{out}}$ that an output $Y_\mu$ is observed given $\frac{1}{\sqrt{n}}[\mathbf{\Phi}\mathbf{x}]_\mu$. $P_{\text{out}}(\cdot | z)$ is the probability density function of $\phi(z, A)$ (where again the r.v. $A \sim P_A$ accounts for noise). This paper is concerned with the so-called Bayes-optimal setting where the prior $P_0$ and the likelihood $P_{\text{out}}$ that appear in the posterior (2) were also used to generate the ground-truth signal $\mathbf{X}^*$ and the labels $\mathbf{Y}$, using a known random matrix $\mathbf{\Phi}$.

A first quantity of interest is the free entropy (which is the free energy up to a sign) defined as $f_n(\mathbf{Y}, \mathbf{\Phi}) = \frac{1}{n} \ln \mathcal{Z}(\mathbf{Y}, \mathbf{\Phi})$. The expectation of the free entropy is equal to minus the conditional entropy density of the observation $-\frac{1}{n}H(\mathbf{Y}|\mathbf{\Phi})$, as well as (up to an additive constant) to the mutual information density between the signal and the observations $\frac{1}{n}I(\mathbf{X}^*; \mathbf{Y}|\mathbf{\Phi})$.

### 2.1 The free entropy

Our first result is the rigorous determination of the free entropy, in the high-dimensional asymptotic regime $n, m \to \infty$, $m/n \to \alpha$. For a random matrix $\mathbf{\Phi}$ with independent entries of zero mean and unit variance, for output $\mathbf{Y}$ that was generated using (1), and under appropriate technical assumptions stated precisely in
section “Methods and proofs”, the free entropy converges in probability to:

$$f_n(\mathbf{Y}, \mathbf{\Phi}) \equiv \frac{1}{n} \ln \mathcal{Z}(\mathbf{Y}, \mathbf{\Phi}) \xrightarrow[n \to \infty]{\mathbb{P}} \sup_{q \in [0, \rho]} \inf_{r \ge 0} f_{\text{RS}}(q, r; \rho) \quad (3)$$

where $\rho \equiv \mathbb{E}_{P_0}[(X^*)^2]$ and where the potential $f_{\text{RS}}(q, r; \rho)$ is

$$f_{\text{RS}}(q, r; \rho) \equiv \psi_{P_0}(r) + \alpha \Psi_{\text{Out}}(q; \rho) - rq/2, \quad (4)$$

$$\psi_{P_0}(r) \equiv \mathbb{E}_{Z_0, X_0} \ln \int dP_0(x) e^{rxX_0 + \sqrt{r}xZ_0 - rx^2/2}, \quad (5)$$

$$\Psi_{P_{\text{Out}}}(q; \rho) \equiv \mathbb{E}_{V,W,\tilde{Y}_0} \ln \int \mathcal{D}w P_{\text{Out}}(\tilde{Y}_0 | \sqrt{q}V + \sqrt{\rho-q}w), \quad (6)$$

where $\mathcal{D}w = dw \exp(-w^2/2)/\sqrt{2\pi}$ is a standard Gaussian measure and the scalar r.v. are independently sampled from $X_0 \sim P_0$, then $V, W, Z_0 \stackrel{\text{iid}}{\sim} \mathcal{N}(0, 1)$ and $\tilde{Y}_0 \sim P_{\text{Out}}(\cdot|\sqrt{q}V + \sqrt{\rho-q}W)$. Only the special linear case with Gaussian $P_{\text{Out}}$ was known rigorously so far [46–48]. Convergence of the *averaged* free entropy is precisely stated in Theorem 1; the one in probability follows from concentration results in the SI.

One can check by explicit comparison that for specific choices of $P_0$ and $P_{\text{Out}}$ the expression (4) is the replica-symmetric free entropy derived in numerous statistical physics papers (thus the RS in $f_{\text{RS}}$), and in particular in [22, 41, 50, 51] for $\varphi(z) = \text{sign}(z)$. The formula for general $P_0$ and $P_{\text{Out}}$ was conjectured based on the statistical physics derivation in [13]. Establishing (3) closes these old conjectures and yields an important step towards vindication of the cavity and replica methods for inference, alongside with e.g. [43, 53]. We now discuss the main consequences of this formula.

## 2.2 Overlap and optimal estimation error

Our second result concerns the overlap between a sample $\mathbf{x}$ from the posterior (2) and the ground-truth. We obtain that as $n, m \to \infty$, $n/m \to \alpha$,

$$\frac{1}{n} |\mathbf{x} \cdot \mathbf{X}^*| \xrightarrow[n \to \infty]{\mathbb{P}} q^* \quad (7)$$

whenever $q^* = q^*(\alpha)$ the maximizer in formula (3) is unique. This is the case for almost every $\alpha$ (see the SI).

It is a simple fact of Bayesian inference that, given the measurements $\mathbf{Y}$ and the measurement matrix $\mathbf{\Phi}$, the estimator $\hat{\mathbf{X}}$ that minimizes the mean-square error with the ground-truth $\mathbf{X}^*$ is the mean of the posterior distribution (2), i.e. $\hat{\mathbf{X}} = \mathbb{E}_{P(\mathbf{x}|\mathbf{Y}, \mathbf{\Phi})}[\mathbf{x}]$. The minimum mean-square error (MMSE) that is achieved by such “Bayes-optimal” estimator is deduced, again in the limit $n \to \infty$, $m/n \to \alpha$, as follows:

$$\text{MMSE} = \frac{1}{n} \mathbb{E} \left[ \| \mathbf{X}^* - \hat{\mathbf{X}} \|^2 \right] \to \rho - q^*. \quad (8)$$

We refer to Theorem 2 in section “Main theorems” for rigorous statements. Again the value of the MMSE was known rigorously so far only for the linear case with Gaussian noise [46–48] (and conjectured for the non-linear case e.g. in [13]).

## 2.3 Optimal generalization error

Our third result concerns the prediction error, also called generalization error. Consider again the statistical model (1). To define the *Bayes-optimal generalization error*, one is given a new row of the matrix/data point, denoted $\mathbf{\Phi}_{\text{new}} \in \mathbb{R}^n$ (in addition to the data $\mathbf{\Phi}$ and associated outputs $\mathbf{Y}$ used for the learning), and is asked to estimate the corresponding output value $Y_{\text{new}}$. We seek for an estimator $\hat{Y}_{\text{new}} = \hat{Y}_{\text{new}}(\mathbf{Y}, \mathbf{\Phi}, \mathbf{\Phi}_{\text{new}})$ that achieves
$\mathcal{E}_{\text{gen}} \equiv \min_{\hat{Y}_{\text{new}}} \mathbb{E}[(Y_{\text{new}} - \hat{Y}_{\text{new}})^2]$, i.e. that minimizes the MSE with the true $Y_{\text{new}}$ obtained using the ground-truth weights $\mathbf{X}^*$. Such estimator is again obtained from the posterior: $\hat{Y}_{\text{new}} = \mathbb{E}_{P_A(a)} \mathbb{E}_{P(\mathbf{x}|\mathbf{Y}, \mathbf{\Phi})} \varphi(\frac{1}{\sqrt{n}} \mathbf{\Phi}_{\text{new}} \cdot \mathbf{x}, a)$. Note that this is different than the plug-in estimator $\tilde{Y}_{\text{new}} = \varphi(\frac{1}{\sqrt{n}} \mathbf{\Phi}_{\text{new}} \cdot \hat{\mathbf{X}})$, which leads to a worse MSE than $\hat{Y}_{\text{new}}$. Yet it is often used in practice for deterministic models since most algorithms for generalized linear regression do not provide the full posterior distribution.

Our result states that the optimal generalization error follows from the I-MMSE theorem [54] applied to the free entropy (3) (see the SI for the details). The optimal generalization error reads as $n \to \infty$, $m/n \to \alpha$ ($q^*$ is the maximizer in (3))

$$ \mathcal{E}_{\text{gen}} \to \mathbb{E}_{V,a} [\varphi(\sqrt{\rho} V, a)^2] - \mathbb{E}_{V} [\mathbb{E}_{w,a} [\varphi(\sqrt{q^*} V + \sqrt{\rho-q^*} w, a)]^2], \quad (9) $$

where $V, w \stackrel{\text{iid}}{\sim} \mathcal{N}(0, 1)$ and $a \sim P_A$. See again Theorem 2 in the section “Main theorems” for the precise statement (and Theorems 3 and 4 in the SI).

Note that for labels $\mathbf{Y}$ belonging to a discrete set the MSE might not be a suitable loss and we are more often interested in maximizing the so-called overlap, i.e. the probability of obtaining the correct label. In that case the Bayes-optimal estimator is computed as the argmax of the posterior marginals, rather than as its mean, i.e. for discrete labels $\tilde{Y}_{\text{new}} = \text{argmax}_{y} \mathbb{P}(y = \varphi(\frac{1}{\sqrt{n}} \mathbf{\Phi}_{\text{new}} \cdot \mathbf{x}, a))$ where again $\mathbf{x}$ is distributed according to (2), $a \sim P_A$. The replica method has been used to compute the optimal generalization error for the perceptron where $\varphi(x) = \text{sign}(z)$ in the pioneering works of [23, 50, 55]. We note that in this special case the plug-in estimator $\tilde{Y}_{\text{new}}$ is actually equal to the optimal one $\tilde{Y}_{\text{new}}$.

A final note concerns the issue of overfitting. In optimization-based approaches to learning overfitting may lead to a generalization error which is too large as compared to the training error. In the Bayes-optimal setting the estimators are constructed in order not to overfit. This is related to general properties of Bayes-optimal inference and learning that are called “Nishimori conditions” in the physics literature [13] and that turn out to be crucial in our proofs.

## 2.4 Optimality of approximate message-passing

While the three results stated above are of an information-theoretic nature, our fourth one concerns the performance of an algorithm to solve random instances of GLMs called generalized approximate message-passing (GAMP) [11–13], which is closely related to the TAP equations developed in statistical physics [41, 56, 57].

The GAMP algorithm can be summarized as follows [11–13]: Given initial estimates $\hat{\mathbf{x}}^0, \mathbf{v}^0$ for the marginal posterior means and variances of the unknown signal vector $\mathbf{X}^*$ entries, GAMP iterates the following equations, with $g_{\mu}^0 = 0$:

$$ \left\{ \begin{array}{l} V^t = \overline{\mathbf{v}^{t-1}} \\ \omega^t = \mathbf{\Phi} \hat{\mathbf{x}}^{t-1} / \sqrt{n} - V^t \mathbf{g}^{t-1} \\ g_{\mu}^t = g_{P_{\text{out}}} (Y_{\mu}, \omega_{\mu}^t, V^t) \quad \forall \mu = 1, \dots, m \\ \lambda^t = \alpha g_{P_{\text{out}}}^2 (\mathbf{Y}, \omega^t, V^t) \\ \mathbf{R}^t = \hat{\mathbf{x}}^{t-1} + (\lambda^t)^{-1} \mathbf{\Phi}^{\top} \mathbf{g}^t / \sqrt{n} \\ \hat{x}_i^t = g_{P_0} (R_i^t, \lambda^t) \quad \forall i = 1, \dots, n \\ v_i^t = (\lambda^t)^{-1} \partial_R g_{P_0} (R, \lambda^t) \Big|_{R=R_i^t} \quad \forall i = 1, \dots, n \end{array} \right. $$

(here we denote by $\bar{\mathbf{u}}$ the average over all the components of a vector $\mathbf{u}$). The so-called thresholding function $g_{P_0}(R, \lambda)$ is defined as the mean of the normalized distribution $\propto P_0(x) \exp(-\lambda(R - x)^2/2)$ and the output function $g_{P_{\text{out}}}(Y, \omega, V)$ is similarly the mean of the normalized distribution (of $x) \propto P_{\text{out}}(Y|\omega +
$$ \sqrt{V}x) \exp(-x^2/2). $$

The heuristic derivation of GAMP in statistical physics [13] suggests via the definition of the function $g_{P_{out}}$ that $\omega$ and $V$ are the estimates of the means and average variance of the components of the variable $\mathbf{z} = \mathbf{\Phi}\mathbf{x}$. This, in turn, suggests a GAMP prediction of labels of new data points:

$$ \hat{Y}_{\text{new}}^{\text{GAMP},t} = \int y P_{\text{out}}(y|\omega_{\text{new}}^t + z\sqrt{V^t}) dy \mathcal{D}z $$

where $\omega_{\text{new}}^t \equiv \frac{1}{\sqrt{n}}\mathbf{\Phi}_{\text{new}} \cdot \hat{\mathbf{x}}^{t-1}$. Comparing it with the test-set labels, this serves to compute GAMP's generalization error.

One of the strongest assets of GAMP is that its performance can be tracked via a closed form procedure known as state evolution (SE), again in the asymptotic limit when $n, m \to \infty$, $m/n \to \alpha$. For proofs of SE see [43, 44] for the linear case, and [45] for the generalized one. In our notations, SE tracks the correlation (or “overlap”) between the true weights $\mathbf{X}^*$ and their estimate $\hat{\mathbf{x}}^t$ defined as $q^t \equiv \lim_{n \to \infty} \frac{1}{n} \mathbf{X}^* \cdot \hat{\mathbf{x}}^t$ via:

$$ q^t = 2\psi'_{P_0}(r^t), \quad r^t = 2\alpha\Psi'_{P_{\text{out}}}(q^{t-1}; \rho). \qquad (10) $$

The derivatives are w.r.t. the first argument. Similarly for the evolution of GAMP's generalization error $\mathcal{E}_{\text{gen}}^{\text{GAMP},t}$ (see SI) we get that it is asymptotically, and with high probability, given by the r.h.s. of formula (9) but with $q^*$ replaced by $q^t$.

It is a simple algebraic fact that the fixed points of the SE equations (10) correspond to the critical points of the potential (4). The question of GAMP achieving asymptotically optimal MMSE or generalization error therefore reduces to the study of the extrema of the two-scalar-variables potential (4). If the SE (10) converges to the same couple $(q, r)$ as the extremizer $(q^*, r^*)$ of (3) then GAMP is optimal, and if it does not then GAMP is sub-optimal. In the next section we illustrate this result on several examples, delimiting regions where GAMP reaches optimality. We note that optimality of AMP-based algorithms in terms of the MMSE on the ground-truth vector $\mathbf{X}^*$ was proven for several cases where the extremizer $q^*$ in (3) is unique, see e.g. [58], or in the linear case of GLM in [47]. Our results allow to complete the characterization of regions of parameters where the algorithm reaches optimal performance in terms of the estimation and generalization errors. While the asymptotic value of the Bayes-optimal generalization error was predicted for some cases of $P_{\text{out}}$ and $P_0$ [55], and TAP-based algorithms were argued to reach this performance in [59, 60], it was not known whether this error can be achieved provably nor for what exact regions of parameters the algorithm is sub-optimal. Our present work settles this question thanks to the state evolution of the GAMP algorithm. Interestingly, heuristic arguments based on the glassy nature of the corresponding probability measure were used to argue that direct sampling or optimization-based approaches will not be able to match this performance [51]. Whether this statement is correct goes beyond the scope of the present paper.

## 3 Application to learning and inference

In this section, we report what our results imply for the information-theoretically optimal errors, and those reached by the GAMP algorithm for several interesting cases of output functions $\varphi$ and prior distributions $P_0$. We do not seek to be exhaustive in any way, we simply aim to illustrate the kind of insights about the GLM that can be obtained from our results. We focus on determination of phase transitions in performance as we vary parameters of the model, e.g. the number of samples or the sparsity of the signal. We use careful numerical procedures to compute the expectations required in the formula (4), and check that the reported results are stable towards the choice of various precision-parameters. In this section we, however, do not seek rigor in bounding formally the corresponding numerical errors. Many of the codes used in this section are
<figure><img src="image_2.png" /><figcaption>Figure 1: Phase diagrams showing boundaries of the region where almost exact recovery is possible (in absence of noise). **Left:** The case of sign-less sparse recovery, $\varphi(x) = |x|$ with a Gauss-Bernoulli signal, as a function of the ratio between number of samples/measurements and the dimension $\alpha = m/n$, and the fraction of non-zero components $\rho$. Evaluating (4) for this case, we find that a recovery of the signal is information-theoretically impossible for $\alpha < \alpha_{IT} = \rho$. Recovery becomes possible starting from $\alpha > \rho$, just as in the canonical compressed sensing. Algorithmically the sign-less case is much harder. Evaluating (11) we conclude that GAMP is not able to perform better than a random guess as long as $\alpha < \alpha_c = 1/2$, and the same is true for spectral algorithms, see [61]. For larger values of $\alpha$, the inference using GAMP leads to better results than a purely random guess. GAMP can recover the signal and generalize perfectly only for values of $\alpha$ larger than $\alpha_{AMP}$ (full red line). The dotted red line shows for comparison the algorithmic phase transition of the canonical compressed sensing. **Center:** Analogous to the left panel, for the ReLU output function, $\varphi(x) = \max(0, x)$. Here it is always possible to perform better than random guessing using GAMP. The dotted red line shows the algorithmic phase transition when using information only about the non-zero observations. **Right:** Phase diagram for the symmetric door output function $\varphi(z) = \text{sign}(|z| - K)$ for a Rademacher signal, as a function of $\alpha$ and $K$. The stability line $\alpha_c$ is depicted in dashed blue, the information-theoretic phase transition to almost exact recovery $\alpha_{IT}$ in black, and the algorithmic one $\alpha_{AMP}$ in red.</figcaption></figure>

given online in a github repository [62].

## 3.1 General observations about fixed points and terminology

* **Non-informative fixed point and its stability:** It is instrumental to analyze under what conditions $q^* = 0$ is the optimizer in (3). Our result (8) about the MMSE implies that if $q^* = 0$ then the MMSE is as large as if we had no samples/measurements at our disposition. A necessary condition for $q^* = 0$ is that it is a fixed point of the state evolution. In turn, a sufficient condition for the state evolution (10) to have such a fixed point is that $i)$ the output density $P_{\text{out}}(y|z)$ is even in the argument $z$, and $ii)$ that the prior $P_0$ has zero mean. A proof of this is given in the SI. In order for $q^* = 0$ to be a fixed point to which the state evolution (10) converges, it needs to be stable. We detail in the SI that under properties $i)$ and $ii)$ this fixed point is stable when

$$ \alpha \int dy \frac{\left( \int \mathcal{D}z (z^2 - 1) P_{\text{out}}(y|\sqrt{\rho}z) \right)^2}{\int \mathcal{D}z P_{\text{out}}(y|\sqrt{\rho}z)} < 1. \qquad (11) $$

In what follows we will denote $\alpha_c$ the largest value of $\alpha$ for which the above condition holds. Consequently the error reachable by the GAMP algorithm is as bad as random guessing for both the estimation and generalization errors as long as $\alpha < \alpha_c$. For $\alpha > \alpha_c$, starting with infinitesimal positive $q$ the state evolution will move towards larger $q$ as in [63]. Note that the condition (11) also appears in a recent work [61] as a barrier for performance of spectral algorithms.

Concerning the information-theoretically optimal error, we will call the phase where MMSE = $\rho$, i.e. $q^* =$
0 is the extremizer of (4), the non-informative phase. Existing literature sometimes refers to such behavior as retarded learning phase [64], in the sense that in that case a critical number of samples is required for the generalization error to be better than random guessing. Below we will evaluate condition (11) explicitly for several examples.

• Almost exact recovery fixed point: Another fixed point of (10) that is worth our particular attention is the one corresponding to almost exact recovery, meaning with average error per coordinate going to 0 as $n \to \infty$, where $q^* = \rho$. A sufficient and necessary condition for this to be a fixed point is that $\lim_{q \to \rho} \Psi'_{\text{out}}(q; \rho) = +\infty$. This means that the integral of the Fisher information of the output channel diverges:

$$ \int dyd\omega \frac{e^{-\frac{\omega^2}{2\rho}} P'_{\text{out}}(y|\omega)^2}{\sqrt{2\pi\rho} P_{\text{out}}(y|\omega)} = +\infty, $$

where $P'_{\text{out}}(y|\omega)$ denotes the partial derivative w.r.t. $\omega$. This typically means that the output channel should be noiseless. For example, for the Gaussian channel with noise variance $\Delta$, the above expression equals $1/\Delta$. For the probit channel where $P_{\text{out}}(y|z) = \text{erfc}(-yz/\sqrt{2\Delta})/2$ the above expression at small $\Delta$ is proportional to $1/\sqrt{\Delta}$.

Stability of the almost exact recovery fixed point depends non-trivially on both the properties of the output channel, and of the prior. Below we give several examples where almost exact recovery either is or is not possible. In what follow we call the region of parameters for which MMSE = 0, i.e. $q^* = \rho$ is the extremizer in (3), the almost exact recovery phase.

• Hard phase: As can be anticipated from the statement of our main algorithmic result, there are regions of parameters for which the error reached by GAMP is asymptotically equal to the optimal error, and regions where it is not. We will call hard phase the region of parameters where MMSE < MSE$_{AMP}$ with a strict inequality. Focusing on the ratio $\alpha$ between the number of samples and the dimensionality, we will denote $\alpha_{IT}$ the ratio for which the hard phase appears, and $\alpha_{AMP} > \alpha_{IT}$ the ratio for which it disappears. In other words, the hard phase is an interval $(\alpha_{IT}, \alpha_{AMP})$, and is associated to a first order phase transition in the Bayes-optimal posterior probability distribution.

It remains a formidable open question of average computational complexity whether in the setting of this paper (and for problems that are NP-complete in the worst case) there exists an efficient algorithm that achieves better performance than GAMP in the hard phase. The authors are not aware of any, and tend to conjecture that there is none.

## 3.2 Sensing compressively with non-linear outputs

Existing literature covers in detail the case of noiseless compressed sensing, i.e. when the output function $\varphi(z) = z$. Representative sparse prior distribution is the Gauss-Bernoulli (GB) distribution $P_0 = \rho N(0, 1) + (1 - \rho)\delta_0$, where $\rho$ is the average fraction of non-zeros, which are in this case standard Gaussians. The phase diagram of this case is well known see e.g. [67, 68]. In noiseless compressed sensing with random i.i.d. matrices and GB prior, almost exact recovery of the signal is possible for $\alpha > \alpha_{IT} = \rho$ and GAMP recovers the signal for $\alpha > \alpha_{AMP,CS}$ where $\alpha_{AMP,CS}$ is plotted in Fig. 1 (left) with a dotted red line, thus delimiting the hard phase of compressed sensing. We note that the Donoho-Tanner phase transition [9] known as the performance limit of the LASSO $\ell_1$ regularization is slightly higher than $\alpha_{AMP,CS}$.

• Sign-less output channel: The phase diagram of noiseless compressed sensing changes intriguingly when only the absolute value of the output is measured, i.e. when $\varphi(z) = |z|$ instead of $\varphi(z) = z$. Such an output channel is reminiscent of the widely studied phase retrieval problem [3] where the signal is complex valued and only the amplitude is observed. Generalization of our results to the complex case would require extensions, as done for the algorithmic aspects in [69]. The real-valued case was studied under the name sparse recovery from quadratic measurements in the literature, e.g. [70] and references therein, when the number of
<figure><img src="image_3.png" /><figcaption>Figure 2: Optimal generalization error in three classification problems versus the sample complexity $\alpha$, the size of the training set being $\alpha n$. The red line is the Bayes-optimal generalization error (9) while the green line shows the (asymptotic) performances of GAMP as predicted by the state evolution (10). For comparison, we also show the results of GAMP (black dots) and, in blue, the performance of a standard out-of-the-box solver. **Left:** Perceptron, with $\varphi(x) = \text{sign}(x)$ and a binary Rademacher signal. While a perfect generalization is information-theoretically possible starting from $\alpha_{\text{IT}} \approx 1.249$, the state evolution predicts that GAMP will achieve such perfect prediction only above $\alpha_{\text{AMP}} \approx 1.493$. The results of a logistic regression with fine-tuned regularizations with the software scikit-learn [65] are shown for comparison. **Middle:** Perceptron with Gauss-Bernoulli distribution of the weights. No phase transition is observed in this case, but a smooth decrease of the error with $\alpha$. The results of a logistic regression are very close to optimal. **Right:** The symmetric door activation rule with parameter $K$ chosen in order to observe the same number of occurrence of the two classes. In this case there is a sharp phase transition from as bad as random to perfect generalization at $\alpha_{\text{IT}} = 1$. GAMP identifies the rule perfectly only starting from $\alpha_{\text{AMP}} \approx 1.566$. The non-informative fixed point is stable up to $\alpha_c = 1.36$ (dashed line). Interestingly, this non linear rule seems very hard to learn for standardly used solvers. Using Keras [66], a neural network with 2 hidden layers was able to learn only approximately the rule, only for considerably larger training set sizes and much larger number of iterations than GAMP.</figcaption></figure>

non-zero variables grows slower than linearly with the dimension $n$. Our results give access to the phase diagram of sparse recovery from quadratic (or equivalently sign-less) measurements, that is presented in Fig. 1 (left) for the GB prior.

We observe that the information-theoretical phase transition $\alpha_{\text{IT}}$ is the same in the sign-less sparse recovery as in the canonical linear case, i.e. almost exact recovery is possible whenever $\alpha > \rho$. However, the algorithmic phase transition $\alpha_{\text{AMP}}$ above which GAMP is able to find the sparse signal<fnref n="1" /> is strikingly larger for the sign-less case (red line in left panel of Fig. 1). We note that even for a dense signal $\rho = 1$ almost exact recovery is algorithmically possible only for $\alpha > \alpha_{\text{AMP}}(\rho = 1) \approx 1.128$. For very sparse signals, small $\rho$, the situation is even more striking because measurement rate of at least $\alpha > \alpha_c = 1/2$ is needed for algorithmically tractable almost exact recovery for every $\rho$. This is in sharp contrast with the canonical compressed sensing where $\alpha_{\text{AMP,CS}} \to 0$ as $\rho \to 0$. The nature of this algorithmic difficulty of GAMP is related to the symmetry of the output channel thanks to which the non-informative fixed point is stable for $\alpha < \alpha_c = 1/2$. Summarizing this result in one sentence, tractable compressive sensing is impossible (for $\alpha < 1/2$) if we have lost the signs. We remind that this result holds in the setting of the present paper, i.e. in particular when the sparsity $\rho$ is of constant order. For signals where $\rho = o(1)$ the situation is expected to be different [70].

• **ReLU output channel:** Another case of output channel that attracted our interest is the rectified linear unit (ReLU), $\varphi(z) = \max(0, z)$, as widely used in multi-layer feed-forward neural networks. In the present single-layer case reconstruction with the ReLU output is interesting mathematically. With GB signals, roughly half of the measurements are given without noise, but the only information we have about the other half is its sign. A straightforward upper-bound for both information-theoretic and tractable almost exact recovery is simply twice as many measurements than needed in the canonical noiseless compressed sensing. It is interesting to ask whether this bound is tight. Results in the present paper imply that for the information-theoretic performance this bound indeed is tight. However, the phase transition $\alpha_{\text{AMP}}$ above which almost exact re-

<footnote n="1">¹We note that in order to break the symmetry that prevents GAMP to find the signal in constant number of iteration steps, we mismatch infinitesimally the output function $\varphi$ used in the algorithm from the symmetric one used to generate the data. Another way to deal with this issue is related to a spectral initialization as discussed recently in [61].</footnote>
covery is possible with the GAMP algorithm is strictly lower than twice the phase transition of compressed sensing; both are depicted in the central panel on Fig. 1. This implies that while the negative outputs are not useful information-theoretically, they do help to achieve better performance algorithmically.

## 3.3 Perceptron and alike

* **Binary and Gauss-Bernoulli perceptron:** One of the most studied problems that fits in the setting of the present paper is the problem of perceptron [71], where $\varphi(z) = \text{sign}(z)$, that has been analyzed for random patterns $\Phi$ in the statistical physics literature, see [23–25] for reviews. We plot in Fig. 2 the optimal generalization error (9) as follows from our results for the binary perceptron, i.e. weights taken from the Rademacher distribution $P_0 = \frac{1}{2}\delta_{+1} + \frac{1}{2}\delta_{-1}$ (left panel), and for the GB perceptron where $P_0 = \rho N(0, 1) + (1-\rho)\delta_0$ (central panel). The information-theoretically optimal value of the generalization error that we report and prove agrees with existing predictions obtained by the non-rigorous replica method from [50, 51, 55]. Notably, we see that for the GB case the optimal generalization error decreases smoothly as $\alpha$ increases, while for the binary case the generalization error has a first order (i.e. discontinuous) phase transition towards perfect generalization at $\alpha_{IT} \approx 1.249$ as predicted already in [50]. Our results provide rigorous validation for these old conjectures.

Furthermore, our results together with recent literature on GAMP provide a refreshing clarification of the algorithmic questions. It is natural to ask for what region of parameters the optimal generalization error can be provably achieved with efficient algorithms. This question remained unanswered until now. Indeed, for the spherical perceptron the optimal generalization error was computed in [55], and argued empirically on small instances to be achievable with a TAP-like algorithm [59]. The state evolution of GAMP together with our formulas for the generalization error (9) for the average optimal one and with $q^t$ replacing $q^*$ in this formula for GAMP) imply that the optimal generalization error is indeed achievable asymptotically for all $\alpha$ in the GB perceptron.

For the binary perceptron the optimal generalization error was computed in [50, 51]. By comparison with the state evolution of GAMP we obtain that it can also be asymptotically achieved by GAMP, but this time only outside of the hard phase ($\alpha_{IT}, \alpha_{AMP}$) with $\alpha_{AMP} \approx 1.493$. The past literature was unclear on the algorithmic question, Ref. [50] identified the spinodal of the replica-symmetric solution to be at $\alpha \approx 1.493$, but did not attribute it any algorithmic nor physical meaning. Ref. [51] argues that metastable states exist at least up to $\alpha_{RSB} \approx 1.628$ and speculates that Gibbs sampling based algorithms will not be able to reach perfect generalization before that point [23]. Taking our results into account, the main algorithmic question that remains open is whether efficient algorithms can reach perfect generalization for $\alpha_{IT} < \alpha < \alpha_{AMP}$.

* **Symmetric door:** Out of interest we explored an example of binary output channel for which $P_{\text{out}}(y|z)$ is even in the argument $z$, so that the non-informative fixed point $q^* = 0$ exists. Specifically we analyzed the *symmetric door* channel with $\varphi(z) = \text{sign}(|z| - K)$ and Rademacher prior $P_0$. In existing literature such a perceptron was studied with the replica method in the context of lossy data compression [28]. In Fig. 1 (right panel) we report the phase diagram in terms of the stability line of the non-informative fixed point $\alpha_c$ (below which GAMP is not better than random guesses), the information-theoretic phase transition towards perfect generalization $\alpha_{IT}$, and the phase transition of GAMP to perfect generalization $\alpha_{AMP}$.

A simple counting lower bound states that for binary outputs and weights $X_i^*$ perfect generalization is not possible for $\alpha < 1$. Thus it is interesting to notice that the symmetric door channel is able to saturate this lower-bound for $K \approx 0.6745$ for which the probability of $y_\mu = 1$ is $1/2$. This saturation was already remarked in [28]. Our results also, however, imply that in that case the perfect generalization will not be achievable with GAMP (and we conjecture no other efficient) algorithm unless $\alpha > \alpha_{AMP} \approx 1.566$. The generalization error that GAMP provides for this case is depicted in Fig. 2 (right).
## 3.4 Empirical comparison with general purpose algorithms

In this section we argue that many cases that fit into the setting of the present paper could serve as useful benchmarks for existing machine learning algorithms. We believe that the situation is perhaps similar to Shannon coding theorems that have driven algorithmic developments in error correcting codes, achieving the Shannon bound being the primary goal in many works in communications. In machine learning, classification is a natural task and algorithms are usually benchmarked using open access databases. In current state-of-the-art applications of machine learning we usually have very little insight about what is the sample complexity, i.e. how many samples are truly needed so that a given generalization error can be achieved. In our setting the situation is different: We can present samples $(y_\mu, \Phi_\mu)$ to generic out-of-the-box classification algorithms and see how their performances compare to the information-theoretic optimal performance and to the one of the GAMP algorithm that is fine-tuned to the problem.

In Fig. 2 we present examples of state-of-the-art classification algorithms compared to our results. In the left and center panels we compare the optimal and GAMP performances to a simple logistic regression, fine-tuned by manually optimizing the ridge penalty (for $\ell_2$ regularization) and LASSO penalty (for a sparsity-enhancing $\ell_1$ regularization) with the software scikit-learn [65]. We observe that while for the GB case the logistic regression is comparable to the performance of GAMP, for binary weights perfect generalization is not achieved close to the GAMP phase transition.

In the right panel of Fig. 2 we study classification for labels generated by the symmetric door channel. A general purpose algorithm would not know about the form of the channel. A neural network with only two hidden units is in principle able to represent the corresponding function (each of the hidden neurons can learn one of the two planes that separate data in the symmetric door function). A more intriguing question is whether a more generic multi-layer neural network is indeed able to learn this rule and how many samples does it need? In the example used in Fig. 2, using the software Keras [66] with a tensorflow backend, we show the performance of a network with two hidden layers, ReLU activation and dropout (the details for this particular run can be found on the github repository [62]). The symmetric door function thus provides a challenging benchmark that could be used to study how to improve performance of the general purpose multi-layer neural network classifiers. In the SI we provide additional examples comparing the optimal performance to general-purpose algorithms for regression.

# 4 Methods and proofs

In this section we give the main theorem for the free entropy and main ideas of the proof. An essential tool is the adaptive interpolation method recently introduced in [52] which is a powerful evolution of the Guerra and Toninelli interpolation method developed for spin glasses [72]. Reference [52] analyzed simpler inference problems. In particular the proof for the upper bound in [52] does not apply to GLMs and requires non-trivial new ingredients. One such new ingredient is to work with a potential $f_{\text{RS}}(q, r; \rho)$ depending on two parameters $(q, r)$ instead of a single one as in [52]. This allows us to use convexity arguments that are crucial in order to finish the proof, see the last section "Matching bounds and end of proof". We stress that the present analysis heavily relies on properties of Bayes-optimal inference that translate into remarkable identities between correlation functions (called Nishimori identities by physicists; see SI for their formulation) valid for all values of parameters. These identities are used in the derivation of (17) and (18) below, which are two essential steps of our proof. The formula from Theorem 1 relies on the Nishimori identities and does not hold out of the Bayes-optimal setting.
## 4.1 Main theorems

For the proof it is necessary to work with a slightly different model with an additive regularizing Gaussian noise with variance $\Delta \ge 0$:

$$Y_{\mu} = \varphi\left(\frac{1}{\sqrt{n}}[\mathbf{\Phi}\mathbf{X}^{*}]_{\mu}, A_{\mu}\right) + \sqrt{\Delta}Z_{\mu}, \quad 1 \le \mu \le m, \qquad (12)$$

where $(Z_{\mu}) \stackrel{\text{iid}}{\sim} \mathcal{N}(0, 1)$, and $(A_{\mu}) \stackrel{\text{iid}}{\sim} P_A$ are r.v. that represent the stochastic part of $\varphi$. It is also instrumental to think of the measurements as the outputs of a “channel” $Y_{\mu} \sim P_{\text{out}}(\cdot|\frac{1}{\sqrt{n}}[\mathbf{\Phi}\mathbf{X}^{*}]_{\mu})$ with transition density $P_{\text{out}}(y|z) = (2\pi\Delta)^{-1/2} \int dP_A(a) \exp\{-\frac{1}{2\Delta}(y - \varphi(z, a))^2\}$ if $\Delta > 0$, or $P_{\text{out}}(y|z) = \int dP_A(a)\mathbf{1}(y = \varphi(z, a))$ else, where $\mathbf{1}(\cdot)$ is the indicator function. Our main theorem holds under the following rather general hypotheses:

(h1) The prior distribution $P_0$ admits a finite third moment and has at least two points in its support.

(h2) The sequence $(\mathbb{E}[|\varphi(\frac{1}{\sqrt{n}}[\mathbf{\Phi}\mathbf{X}^{*}]_1, \mathbf{A}_1)|^{2+\gamma}])_{n \ge 1}$ is bounded for some $\gamma > 0$.

(h3) The r.v. $(\Phi_{\mu i})$ are independent with zero mean, unit variance and finite third moment bounded with $n$.

(h4) For almost-all values of $a$ (w.r.t. the distribution $P_A$), the function $x \mapsto \varphi(x, a)$ is continuous almost everywhere.

(h5) $(\Delta > 0)$ or $(\Delta = 0$ and $\varphi$ takes values in $\mathbb{N})$.

In general, when $\varphi$ is continuous the condition $\Delta > 0$ (but arbitrarily small) is necessary for the existence of a finite limit of the free entropy (for particular choices of $(\varphi, P_A)$ this might not be needed, e.g. $\varphi(z, A) = z + A$ with $A \sim \mathcal{N}(0, \sigma^2)$). We also assume that the kernel $P_{\text{out}}$ is informative, i.e. there exists $y$ such that $P_{\text{out}}(y|\cdot)$ is not equal almost everywhere to a constant. If $P_{\text{out}}$ is not informative, it is not difficult to show that estimation is then impossible.

We define the set of the critical points of $f_{\text{RS}}$, (4), also called “state evolution fixed points” (as it is clear from (10)):

$$\Gamma = \{ (q, r) \in [0, \rho] \times (\mathbb{R}_+ \cup \{+\infty\}) \mid \begin{cases} q = 2\psi'_{P_0}(r) \\ r = 2\alpha\Psi'_{P_{\text{out}}}(q; \rho) \end{cases} \}.$$

Define $f_n = \mathbb{E}f_n(\mathbf{Y}, \mathbf{\Phi}) = \frac{1}{n}\mathbb{E}\ln Z(\mathbf{Y}, \mathbf{\Phi})$. Then the main theorem of this paper is stated as follows:

**Theorem 1 (Replica-symmetric free entropy).** Suppose that (h1)-(h2)-(h3)-(h4)-(h5.a) hold. Then, for the GLM (12),

$$\lim_{n \to \infty} f_n = \sup_{q \in [0, \rho]} \inf_{r \ge 0} f_{\text{RS}}(q, r) = \sup_{(q,r) \in \Gamma} f_{\text{RS}}(q, r).$$

Moreover, as one can see in the SI, the “sup inf” and the supremum over $\Gamma$ above are achieved over the same couples. Under stronger assumptions on $P_0$ and $P_{\text{out}}$, one can show (see Theorem 6 in the SI) that $f_n(\mathbf{Y}, \mathbf{\Phi})$ concentrates around its mean $f_n$ and thus obtains convergence in probability (3).

An immediate corollary of Theorem 1 is the limiting expression of the mutual information $I(\mathbf{X}^*; \mathbf{Y}|\mathbf{\Phi}) = \mathbb{E}\ln P(\mathbf{Y}, \mathbf{X}^*|\mathbf{\Phi}) - \mathbb{E}\ln(P(\mathbf{Y}|\mathbf{\Phi})P(\mathbf{X}^*))$ between the observations and the unknown vector:

**Corollary 1 (Mutual information).** Under the same hypotheses as in Theorem 1, the mut. info. for the GLM (12) verifies

$$\lim_{n \to \infty} \frac{1}{n} I(\mathbf{X}^*; \mathbf{Y}|\mathbf{\Phi}) = \inf_{q \in [0, \rho]} \sup_{r \ge 0} i_{\text{RS}}(q, r) = \inf_{(q,r) \in \Gamma} i_{\text{RS}}(q, r), \\ i_{\text{RS}}(q, r) = \alpha \Psi_{P_{\text{out}}}(\rho; \rho) - \alpha \Psi_{P_{\text{out}}}(q; \rho) - \psi_{P_0}(r) + rq/2.$$

Finally, we gather our main results related to the optimal errors in a single theorem, see the SI for more details, including results on the optimality of the GAMP algorithm:
**Theorem 2 (Optimal errors).** Assume the same hypotheses as in Theorem 1. Then formula (9) for the generalization error is true as $n, m \to \infty$, $m/n \to \alpha$ whenever the maximizer $q^*(\alpha)$ of (3) is unique, which is the case for almost every $\alpha$. If moreover all the moments of $P_0$ are finite, then formula (7) for the overlap as well as the following matrix-MMSE formula

$$ \frac{1}{n^2} \mathbb{E} \left[ \| \mathbf{X}^* \mathbf{X}^{*\top} - \mathbb{E}_{P(\mathbf{x}|\mathbf{Y}, \mathbf{\Phi})} [\mathbf{x}\mathbf{x}^\top] \|_{\text{F}}^2 \right] \to \rho^2 - q^*(\alpha)^2 \quad (13) $$

are true, where $\| - \|_{\text{F}}$ is the Frobenius norm.

There are cases of GLMs (e.g. the sign-less output channel $\mathbf{Y} = |\mathbf{\Phi}\mathbf{X}^*|/\sqrt{n} + \mathbf{Z}$) where the sign of $\mathbf{X}^*$ simply cannot be estimated (thus the absolute value in (7)). This is why our general theorem is related to an error metric (13) insensitive to this $\pm$ symmetry. Nevertheless formula (8) for the signal MSE is formally valid when there is no such sign symmetry.

## 4.2 Proof by the adaptive interpolation method

We now give the main ideas behind the proof of Theorem 1. We defer to the SI the details, as well as those of Corollary 1 and Theorem 2.

A word about notation. The r.v. $\mathbf{Y}$ (and also $\mathbf{\Phi}$, $\mathbf{X}^*$, $\mathbf{A}$, $\mathbf{Z}$) are called quenched because once the measurements are acquired they are fixed. The expectation w.r.t. all quenched r.v. will be denoted by $\mathbb{E}$ without subscript. In contrast, expectation of annealed variables w.r.t. a posterior distribution at fixed quenched variables is denoted by Gibbs brackets $\langle - \rangle$.

### 4.2.1 Two scalar inference channels

An important role in the proof is played by two simple scalar inference channels. The free entropy is expressed in terms of the free entropies of these channels. This “decoupling property” stands at the root of the replica approach in statistical physics.

The first scalar channel is an additive Gaussian channel. Suppose that we observe $Y_0 = \sqrt{r} X_0 + Z_0$ where $X_0 \sim P_0$ and $Z_0 \sim \mathcal{N}(0, 1)$ are independent. Consider the inference problem consisting of retrieving $X_0$ from the observation $Y_0$. The free entropy associated with this channel is the expectation of the logarithm of the normalization factor of the associated posterior $dP(x|Y_0)$, that is given by (5) (up to a constant).

The second scalar channel that appears naturally in the problem is linked to the channel $P_{\text{out}}$ through the following inference model. Suppose that $V, W^* \stackrel{\text{iid}}{=} \mathcal{N}(0, 1)$ where $V$ is known while the inference problem is to recover the unknown $W^*$ from the observation $\tilde{Y}_0 \sim P_{\text{out}}(\cdot | \sqrt{q}V + \sqrt{\rho-q}W^*)$ where $\rho > 0$ and $q \in [0, \rho]$. The free entropy for this model, again related to the average logarithm of the normalization factor of the posterior of $w$ given $\tilde{Y}_0$ and $V$, is exactly (6).

### 4.2.2 Interpolating estimation problem

To carry on the proof, we introduce an “interpolating estimation problem” that interpolates between the original problem $Y_\mu \sim P_{\text{out}}(\cdot | \frac{1}{\sqrt{n}}[\mathbf{\Phi}\mathbf{X}^*]_\mu)$ at $t = 0, t \in [0, 1]$ being the interpolation parameter, and the two scalar problems described above at $t = 1$. For $t \in (0, 1)$ the interpolating estimation problem is a mixture of the original and the scalar problems. This interpolation scheme is inspired by the interpolation paths used by Talagrand to study the perceptron, see [73]. Thanks to a novel ingredient specific to the adaptive interpolation method [52], it allows to obtain in a unified manner a complete proof of the replica formula for the free entropy and this in the whole phase diagram.

Let $q_\epsilon : [0, 1] \to [0, \rho]$, $r_\epsilon : [0, 1] \to [0, r_{\text{max}}]$, $r_{\text{max}} \equiv 2\alpha\Psi'_{\text{out}}(\rho; \rho)$, be two continuous “interpolating functions” parametrized by $\epsilon = (\epsilon_1, \epsilon_2) \in \mathcal{B}_n \equiv [s_n, 2s_n]^2$, with $(s_n)_{n \ge 1} \in (0, 1/2]^N$ a sequence
that tends to zero slowly enough. Set $R_1(t, \epsilon) \equiv \epsilon_1 + \int_0^t r_\epsilon(v)dv$, $R_2(t, \epsilon) \equiv \epsilon_2 + \int_0^t q_\epsilon(v)dv$ and define $S_{t,\mu} = S_{t,\mu}(\mathbf{X}^*, W_\mu^*, V_\mu, \mathbf{\Phi})$ as

$$S_{t,\mu} \equiv \sqrt{\frac{1-t}{n}} [\mathbf{\Phi}\mathbf{X}^*]_\mu + \sqrt{R_2(t, \epsilon)} V_\mu + \sqrt{\rho t - R_2(t, \epsilon) + 2s_n} W_\mu^*$$

where $(V_\mu), (W_\mu^*) \stackrel{\text{iid}}{\sim} \mathcal{N}(0, 1)$. Consider the following observation channels, with two types of observations obtained through

$$\left\{ \begin{array}{ll} Y_{t,\mu} & \sim P_{\text{out}}(\cdot | S_{t,\mu}), & \text{for } 1 \le \mu \le m, \\ Y'_{t,i} & = \sqrt{R_1(t, \epsilon)} X_i^* + Z'_i, & \text{for } 1 \le i \le n, \end{array} \right. \qquad (14)$$

where $(Z'_i) \stackrel{\text{iid}}{\sim} \mathcal{N}(0, 1)$. We assume that $\mathbf{V} = (V_\mu)_{\mu=1}^m$ is known and that the inference problem is to recover both $\mathbf{W}^* = (W_\mu^*)_{\mu=1}^m$ and $\mathbf{X}^* = (X_i^*)_{i=1}^n$ from the “t-dependent” observations $\mathbf{Y}_t = (Y_{t,\mu})_{\mu=1}^m$ and $\mathbf{Y}'_t = (Y'_{t,i})_{i=1}^n$.

We now understand that $(R_1, R_2)$ and $1-t$ appearing in the first and second set of measurements in (14) play the role of signal-to-noise ratios (snr) in the interpolating problem, with $t$ giving more and more “power” to the scalar inference channels when increasing. Here is the first crucial ingredient of our interpolation scheme. In classical interpolations, these snr would all take a trivial form, i.e. be linear in $t$, but here, the non-trivial integral dependency in $t$ of $(R_1, R_2)$ allows for much more flexibility when choosing the interpolation path. This allows us to actually choose the “optimal interpolation path”. This will become clear below as well as the role of the “small perturbation” parameters $(\epsilon_1, \epsilon_2)$.

Define $u_y(x) \equiv \ln P_{\text{out}}(y|x)$ and, with a slight abuse of notations, we also define the quantity $s_{t,\mu} = s_{t,\mu}(\mathbf{x}, w_\mu, V_\mu, \mathbf{\Phi}) \equiv S_{t,\mu}(\mathbf{x}, w_\mu, V_\mu, \mathbf{\Phi})$, the expression above with $\mathbf{X}^*, W_\mu^*$ replaced by $\mathbf{x}, w_\mu$. We introduce the interpolating Hamiltonian $\mathcal{H}_t = \mathcal{H}_t(\mathbf{x}, \mathbf{w}; \mathbf{Y}_t, \mathbf{Y}'_t, \mathbf{\Phi}, \mathbf{V})$

$$\mathcal{H}_t \equiv -\sum_{\mu=1}^{m} u_{Y_{t,\mu}}(s_{t,\mu}) + \frac{1}{2} \sum_{i=1}^{n} (Y'_{t,i} - \sqrt{R_1(t, \epsilon)} x_i)^2$$

and the corresponding $(t$-dependent) Gibbs bracket $\langle - \rangle_t$ which is the expectation w.r.t. the joint posterior distribution of $(\mathbf{x}, \mathbf{w})$ given the observations $\mathbf{Y}_t, \mathbf{Y}'_t$ (and $\mathbf{\Phi}, \mathbf{V}$), defined as

$$\langle L(\mathbf{x}, \mathbf{w}) \rangle_t \equiv \frac{1}{\mathcal{Z}_t(\mathbf{Y}_t, \mathbf{Y}'_t, \mathbf{\Phi}, \mathbf{V})} \int dP_0(\mathbf{x}) \mathcal{D}\mathbf{w} L(\mathbf{x}, \mathbf{w}) e^{-\mathcal{H}_t},$$

for every continuous bounded test function $L$. Here $\mathcal{Z}_t \equiv \int dP_0(\mathbf{x}) \mathcal{D}\mathbf{w} \exp\{-\mathcal{H}_t(\mathbf{x}, \mathbf{w}; \mathbf{Y}_t, \mathbf{Y}'_t, \mathbf{\Phi}, \mathbf{V})\}$ is the appropriate normalization, $\mathcal{D}\mathbf{w}$ is the standard Gaussian measure. Finally we introduce

$$f_{n,\epsilon}(t) \equiv \frac{1}{n} \mathbb{E} \ln \mathcal{Z}_t(\mathbf{Y}_t, \mathbf{Y}'_t, \mathbf{\Phi}, \mathbf{V})$$

which is the interpolating free entropy. One verifies that

$$\left\{ \begin{array}{l} f_{n,\epsilon}(0) = f_n - \frac{1}{2} + \mathcal{O}(s_n), \\ f_{n,\epsilon}(1) = \psi_{P_0}(\int_0^1 r_\epsilon(t)dt) - \frac{1}{2}(1 + \rho \int_0^1 r_\epsilon(t)dt) + \frac{m}{n} \Psi_{P_{\text{out}}}(\int_0^1 q_\epsilon(t)dt; \rho) + \mathcal{O}(s_n), \end{array} \right. \qquad (15)$$

where $|\mathcal{O}(s_n)| \le C s_n$ for a constant $C > 0$. Now comes another crucial property of the interpolating model: It is such that at $t = 0$ we recover the original problem and $f_{n,\epsilon}(0) = f_n - 1/2 + \mathcal{O}(s_n)$ (the constant 1/2 comes from the purely noisy measurements of the second channel in (14)), while at $t = 1$ we have two scalar inference channels and thus the associated terms $\psi_{P_0}$ and $\Psi_{P_{\text{out}}}$ appear in $f_{n,\epsilon}(1)$. These are precisely the
terms appearing in the free entropy potential (4).

### 4.2.3 Entropy variation along the interpolation

From the understanding of the previous section, it is natural to evaluate the variation of entropy along the interpolation, which allows to “compare” the original and purely scalar models thanks to the identity

$$f_n = f_{n,\epsilon}(1) - \int_0^1 \frac{df_{n,\epsilon}(t)}{dt} dt + \frac{1}{2} + \mathcal{O}(s_n). \qquad (16)$$

Then by choosing the optimal interpolation path thanks to the non-trivial snr dependencies in $t$, we are able to show the equality between the replica formula and the free entropy $f_n$.

We thus compute the $t$-derivative of the free entropy (see the SI for the details of this calculation). It is given by

$$\frac{df_{n,\epsilon}(t)}{dt} = \frac{r_{\epsilon}(t)}{2}(q_{\epsilon}(t) - \rho) - \frac{1}{2}\mathbb{E}\left\langle \left(\frac{1}{n}\sum_{\mu=1}^{m} u'_{Y_{t,\mu}}(S_{t,\mu})u'_{Y_{t,\mu}}(s_{t,\mu}) - r_{\epsilon}(t)\right)(Q - q_{\epsilon}(t)) \right\rangle_t + \mathcal{O}_n(1), \quad (17)$$

where $\mathcal{O}_n(1)$ is a quantity that goes to 0 in the $n, m \to \infty$ limit, uniformly in $t, \epsilon$, and in the interpolating functions $q_\epsilon, r_\epsilon$. The overlap is $Q = Q_n \equiv \sum_{i=1}^n X_i^* x_i/n$.

We now state a crucial result in an informal way and refer to the SI for precise statements. Formally, the overlap concentrates around its mean (for all $t \in [0, 1]$), a behaviour called “replica-symmetric” in statistical physics. In order to make this statement rigorous, one has to include the $\epsilon$-dependent small perturbation in (14) which effectively adds “side-information” about $X^*$ (e.g. think of $t = 0$) without affecting the asymptotic free entropy density. This perturbation forces the overlap to concentrate. We prove that: If for each $t$ the map $R^t: (\epsilon_1, \epsilon_2) \in \mathcal{B}_n \mapsto (R_1(t, \epsilon), R_2(t, \epsilon)) \in R^t(\mathcal{B}_n)$ is a $C^1$ diffeomorphism whose Jacobian has determinant greater or equal to 1, then we have for $s_n = \frac{1}{2}n^{-1/16}$ (see Proposition 4 of Sec. 4.3 in the SI for the precise statement)

$$\frac{1}{s_n^2} \int_{\mathcal{B}_n} d\epsilon \int_0^1 dt \mathbb{E}\left\langle (Q - \mathbb{E}\langle Q \rangle_t)^2 \right\rangle_t = \mathcal{O}(n^{-1/8}). \qquad (18)$$

As will be seen below it is possible to choose interpolating functions that satisfy the required condition.

### 4.2.4 Canceling the remainder

Note from (15) and (4) that the first term appearing in (17) is precisely the missing one to obtain the expression of the potential on the r.h.s. of (16). Thus we would like to “cancel” the Gibbs bracket in (17). This term is called *remainder*. In order to prove the replica formula, we have to show that this remainder vanishes: Thanks to the freedom of choice of interpolation paths $(r_\epsilon, q_\epsilon)$, we are able to do so by “adapting” the interpolation. Thus we would like to choose $q_\epsilon(t) = \mathbb{E}\langle Q \rangle_t \approx Q$ because of (18). However, $\mathbb{E}\langle Q \rangle_t$ is a function of $\int_0^t q_\epsilon(v) dv$ (and of $t, \epsilon$ and $\int_0^t r_\epsilon(v) dv$). The equation $q_\epsilon(t) = \mathbb{E}\langle Q \rangle_t$ is therefore a first order differential equation over $t \mapsto \int_0^t q_\epsilon(v) dv$. Assuming for the moment that this differential equation has a solution, the Cauchy-Schwarz inequality applied to the remainder together with (18) allows to show that the absolute value of this remainder integrated over $(\epsilon, t) \in \mathcal{B}_n \times [0, 1]$ is $\mathcal{O}(s_n^2 n^{-1/16})$. Combining this result with (15) and (16) leads to the following fundamental sum rule (Proposition 5 of Sec. 4.3 in SI):

$$f_n = \frac{1}{s_n^2} \int_{\mathcal{B}_n} d\epsilon \{ \psi_{P_0}(\int_0^1 r_\epsilon(t) dt) + \alpha \Psi_{P_{\text{out}}}(\int_0^1 q_\epsilon(t) dt; \rho) - \frac{1}{2} \int_0^1 r_\epsilon(t) q_\epsilon(t) dt \} + \mathcal{O}_n(1). \qquad (19)$$
### 4.2.5 Matching bounds and end of proof

We now possess all the necessary tools to finish the proof of Theorem 1. We first prove that $\lim_{n \to \infty} f_n = \sup_{r \ge 0} \inf_{q \in [0, \rho]} f_{RS}(q, r)$. Then in the SI, we show that $i$) this is also equal to $\sup_{q \in [0, \rho]} \inf_{r \ge 0} f_{RS}(q, r)$ which gives the first equality of the theorem; $ii$) that this sup inf is attained at the supremum of the state evolution fixed points, which gives the second equality.

* **Lower bound:** We choose the constant function $r_\epsilon(t) = r$ for $t \in [0, 1]$. In the SI we show, using the Cauchy-Lipschitz theorem and the Liouville formula, that the differential equation $q_\epsilon(t) = \mathbb{E}\langle Q \rangle_t$ possesses a unique solution and that the map $R^t: (\epsilon_1, \epsilon_2) \mapsto (\epsilon_1 + rt, \epsilon_2 + \int_0^t q_\epsilon(v)dv)$ is a $C^1$ diffeomorphism with Jacobian greater than 1 (so (18) is valid). Identity (19) then implies $\liminf_{n \to \infty} f_n \ge \inf_{q \in [0, \rho]} f_{RS}(q, r)$ for all $r \in [0, r_{\max}]$. Thus $\liminf_{n \to \infty} f_n \ge \sup_{r \in [0, r_{\max}]} \inf_{q \in [0, \rho]} f_{RS}(q, r)$. In the SI an easy argument shows the r.h.s. is in fact equal to $\sup_{r \ge 0} \inf_{q \in [0, \rho]} f_{RS}(q, r)$.

* **Upper bound:** We choose the interpolating functions as solutions of the following system of 1st order differential equations: $r_\epsilon(t) = 2\alpha \Psi'_{P_{\text{out}}}(\mathbb{E}\langle Q \rangle_t)$, $q_\epsilon(t) = \mathbb{E}\langle Q \rangle_t$. Again, applying the Cauchy-Lipschitz theorem and the Liouville formula we show in the SI that this system admits a unique solution and the map $R^t: (\epsilon_1, \epsilon_2) \mapsto (R_1(t, \epsilon), R_2(t, \epsilon))$ is a $C^1$ diffeomorphism with determinant greater or equal to 1. So with this choice of interpolating functions (18) is valid and we have (19). We show in the SI (Proposition 18) that $\Psi_{P_{\text{out}}}(q; \rho)$ is convex in $q$ and thus $g: q \in [0, \rho] \mapsto 2\alpha \Psi_{P_{\text{out}}}(q; \rho) - r_\epsilon(t)q$ is convex too. Since by the differential equations $r_\epsilon(t) = 2\alpha \Psi'_{P_{\text{out}}}(q_\epsilon(t))$, the function $g$ must attain its minimum at $q = q_\epsilon(t)$. By Proposition 17 in the SI $\psi_{P_0}(r)$ is convex, thus from Jensen and the last remark, the integrand $\{\dots\}$ in (19) is bounded as

$$ \begin{aligned} & \psi_{P_0}(\int_0^1 r_\epsilon(t)dt) + \alpha \Psi_{P_{\text{out}}}(\int_0^1 q_\epsilon(t)dt; \rho) - \frac{1}{2} \int_0^1 r_\epsilon(t)q_\epsilon(t)dt \\ & \le \int_0^1 dt \{\psi_{P_0}(r_\epsilon(t)) + \alpha \Psi_{P_{\text{out}}}(q_\epsilon(t); \rho) - \frac{1}{2}r_\epsilon(t)q_\epsilon(t)\} \\ & = \int_0^1 dt \{\psi_{P_0}(r_\epsilon(t)) + \inf_{q \in [0, \rho]} (\alpha \Psi_{P_{\text{out}}}(q; \rho) - \frac{1}{2}r_\epsilon(t)q)\} \\ & \le \sup_{r \ge 0} \inf_{q \in [0, \rho]} \{\psi_{P_0}(r) + \alpha \Psi_{P_{\text{out}}}(q; \rho) - \frac{1}{2}rq\} \end{aligned} $$

which implies $\limsup_{n \to \infty} f_n \le \sup_{r \ge 0} \inf_{q \in [0, \rho]} f_{RS}(q, r)$.

## Acknowledgments

J.B. acknowledges funding from the SNSF (grant 200021-156672). F.K. and L.Z. acknowledge funding from the ERC under the European Union's 7th Framework Programme Grant Agreement 307087-SPARCS and under Horizon 2020 Research and Innovation Programme Grant Agreement 714608-SMiLe. F.K. and N.M. acknowledge support from the ANR-PAIL. Part of this work was done while L.M. was visiting EPFL.

## References

[1] J. Nelder and R. Wedderburn. Generalized linear models. *Journal of the Royal Statistical Society. Series A (General)*, 135(3):370–384, 1972.

[2] P. McCullagh. Generalized linear models. *European Journal of Operational Research*, 16(3):285–292, 1984.

[3] J. R. Fienup. Phase retrieval algorithms: a comparison. *Applied optics*, 21(15), 1982.
[4] L. Demanet and P. Hand. Stable optimizationless recovery from phaseless linear measurements. *Journal of Fourier Analysis and Applications*, 20(1):199-221, 2014.

[5] E. J. Candes, T. Strohmer, and V. Voroninski. Phaselift: Exact and stable signal recovery from magnitude measurements via convex programming. *Communications on Pure and Applied Mathematics*, 66(8):1241-1274, 2013.

[6] P. T. Boufounos and R. G. Baraniuk. 1-bit compressive sensing. In *42nd Annual Conference on Information Sciences and Systems (CISS)*, pages 16–21. IEEE, 2008.

[7] P. Bühlmann and S. Van De Geer. *Statistics for high-dimensional data: methods, theory and applications*. Springer Science & Business Media, 2011.

[8] Y. LeCun, Y. Bengio, and G. Hinton. Deep learning. *Nature*, 521(7553):436-444, 2015.

[9] D. L. Donoho and J. Tanner. Sparse nonnegative solution of underdetermined linear equations by linear programming. *Proceedings of the National Academy of Sciences of the United States of America*, 102(27):9446-9451, 2005.

[10] E. J. Candes and T. Tao. Near-optimal signal recovery from random projections: Universal encoding strategies? *IEEE Transactions on Information Theory*, 52(12):5406-5425, Dec 2006.

[11] D. L. Donoho, A. Maleki, and A. Montanari. Message-passing algorithms for compressed sensing. *Proceedings of the National Academy of Sciences*, 106(45):18914-18919, Nov 2009.

[12] S. Rangan. Generalized approximate message passing for estimation with random linear mixing. In *IEEE International Symposium on Information Theory Proceedings (ISIT)*, pages 2168-2172, July 2011.

[13] L. Zdeborová and F. Krzakala. Statistical physics of inference: thresholds and algorithms. *Advances in Physics*, 65(5):453-552, 2016.

[14] U. Kamilov, V. K. Goyal, and S. Rangan. Optimal quantization for compressive sensing under message passing reconstruction. In *IEEE International Symposium on Information Theory Proceedings (ISIT)*, pages 459-463. IEEE, 2011.

[15] Y. Xu, Y. Kabashima, and L. Zdeborová. Bayesian signal reconstruction for 1-bit compressed sensing. *Journal of Statistical Mechanics: Theory and Experiment*, 2014(11):P11015, 2014.

[16] P. Schniter and S. Rangan. Compressive phase retrieval via generalized approximate message passing. *IEEE Transactions on Signal Processing*, 63(4):1043-1055, 2015.

[17] M. Bayati and A. Montanari. The lasso risk for gaussian matrices. *IEEE Transactions on Information Theory*, 58(4):1997-2017, April 2012.

[18] N. El Karoui, D. Bean, P. J. Bickel, C. Lim, and B. Yu. On robust regression with high-dimensional predictors. *Proceedings of the National Academy of Sciences*, 110(36):14557-14562, 2013.

[19] D. Donoho and A. Montanari. High dimensional robust m-estimation: asymptotic variance via approximate message passing. *Probability Theory and Related Fields*, 166(3):935-969, Dec 2016.

[20] R. Gribonval and P. Machart. Reconciling" priors" & "priors" without prejudice? In *Advances in Neural Information Processing Systems*, pages 2193-2201, 2013.

[21] M. Advani and S. Ganguli. An equivalence between high dimensional bayes optimal inference and m-estimation. In *Advances in Neural Information Processing Systems*, 2016.

[22] E. Gardner and B. Derrida. Three unfinished works on the optimal storage capacity of networks. *Journal of Physics A: Mathematical and General*, 22(12):1983, 1989.
[23] H. S. Seung, H. Sompolinsky, and N. Tishby. Statistical mechanics of learning from examples. *Phys. Rev. A*, 45:6056–6091, Apr 1992.

[24] T. L. H. Watkin, A. Rau, and M. Biehl. The statistical mechanics of learning a rule. *Rev. Mod. Phys.*, 65:499–556, Apr 1993.

[25] A. Engel and C. Van den Broeck. *Statistical mechanics of learning*. Cambridge University Press, 2001.

[26] A. Engel and L. Reimers. Reliability of replica symmetry for the generalization problem in a toy multilayer neural network. *Europhys. Lett.*, 28:531, 1994.

[27] G. J. Bex, R. Serneels, and C. V. den Broeck. Storage capacity and generalization error for the reversed-wedge ising perceptron. *Phys. Rev. E*, 51:6309, 1995.

[28] T. Hosaka, Y. Kabashima, and H. Nishimori. Statistical mechanics of lossy data compression using a nonmonotonic perceptron. *Physical Review E*, 66(6):066126, 2002.

[29] C. Baldassi, C. Borgs, J. T. Chayes, A. Ingrosso, C. Lucibello, L. Saglietti, and R. Zecchina. Unreasonable effectiveness of learning neural networks: From accessible states and robust ensembles to basic algorithmic schemes. *Proceedings of the National Academy of Sciences*, 113(48):E7655–E7662, 2016.

[30] C. H. Martin and M. W. Mahoney. Rethinking generalization requires revisiting old ideas: statistical mechanics approaches and complex learning behavior. 2017.

[31] N. Tishby, F. C. Pereira, and W. Bialek. The information bottleneck method. In *Proc. of the 37-th Annual Allerton Conference on Communication, Control and Computing*, pages 368–377, 1999.

[32] R. Shwartz-Ziv and N. Tishby. Opening the black box of deep neural networks via information. CoRR, abs/1703.00810, 2017.

[33] C. E. Shannon. A mathematical theory of communication, part i, part ii. *Bell Syst. Tech. J.*, 27:623–656, 1948.

[34] T. Tanaka. A statistical-mechanics approach to large-system analysis of cdma multiuser detectors. *IEEE Transactions on Information Theory*, 48(11):2888–2910, Nov 2002.

[35] D. Guo and S. Verdú. Randomly spread cdma: Asymptotics via statistical physics. *IEEE Transactions on Information Theory*, 51(6):1983–2010, June 2005.

[36] A. R. Barron and A. Joseph. Toward fast reliable communication at rates near capacity with gaussian noise. In *IEEE International Symposium on Information Theory (ISIT)*, pages 315–319, June 2010.

[37] J. Barbier and F. Krzakala. Approximate message-passing decoder and capacity-achieving sparse superposition codes. *IEEE Transactions on Information Theory*, 2017.

[38] C. Rush, A. Greig, and R. Venkataramanan. Capacity-achieving sparse superposition codes via approximate message passing decoding. *IEEE Transactions on Information Theory*, 63(3):1476–1500, 2017.

[39] J. Barbier, M. Dia, and N. Macris. Threshold saturation of spatially coupled sparse superposition codes for all memoryless channels. In *IEEE Information Theory Workshop (ITW)*, Sept 2016.

[40] J. Barbier, M. Dia, and N. Macris. Universal sparse superposition codes with spatial coupling and GAMP decoding. CoRR, abs/1707.04203, 2017.

[41] M. Mézard. The space of interactions in neural networks: Gardner’s computation with the cavity method. *Journal of Physics A: Mathematical and General*, 22(12):2181–2190, 1989.

[42] E. Bolthausen. An iterative construction of solutions of the tap equations for the sherrington-kirkpatrick model. *Communications in Mathematical Physics*, 325(1):333–366, 2014.
[43] M. Bayati and A. Montanari. The dynamics of message passing on dense graphs, with applications to compressed sensing. *IEEE Transactions on Information Theory*, 57(2):764–785, Feb 2011.

[44] M. Bayati, M. Lelarge, and A. Montanari. Universality in polytope phase transitions and message passing algorithms. *The Annals of Applied Probability*, 25(2):753–822, 2015.

[45] A. Javanmard and A. Montanari. State evolution for general approximate message passing algorithms, with applications to spatial coupling. *Information and Inference: A Journal of the IMA*, 2(2):115–144, 2013.

[46] J. Barbier, M. Dia, N. Macris, and F. Krzakala. The mutual information in random linear estimation. In *54th Annual Allerton Conference on Communication, Control, and Computing*, pages 625–632. IEEE, 2016.

[47] J. Barbier, N. Macris, M. Dia, and F. Krzakala. Mutual information and optimality of approximate message-passing in random linear estimation. arXiv:1701.05823, 2017.

[48] G. Reeves and H. D. Pfister. The replica-symmetric prediction for compressed sensing with gaussian matrices is exact. In *IEEE International Symposium on Information Theory (ISIT)*, pages 665–669, 2016.

[49] M. Mézard, G. Parisi, and M.-A. Virasoro. *Spin glass theory and beyond*. World Scientific Publishing, 1987.

[50] G. Györgyi. First-order transition to perfect generalization in a neural network with binary synapses. *Physical Review A*, 41(12):7097, 1990.

[51] H. Sompolinsky, N. Tishby, and H. S. Seung. Learning from examples in large neural networks. *Physical Review Letters*, 65(13):1683, 1990.

[52] J. Barbier and N. Macris. The adaptive interpolation method: a simple scheme to prove replica formulas in bayesian inference. *Probability Theory and Related Fields*, Oct 2018.

[53] A. Coja-Oghlan, F. Krzakala, W. Perkins, and L. Zdeborova. Information-theoretic thresholds from the cavity method. In *Proceedings of the 49th Annual ACM SIGACT Symposium on Theory of Computing (STOC)*, pages 146–157, 2017.

[54] D. Guo, S. Shamai, and S. Verdú. Mutual information and minimum mean-square error in gaussian channels. *IEEE Transactions on Information Theory*, 51(4):1261–1282, April 2005.

[55] M. Opper and D. Haussler. Generalization performance of bayes optimal classification algorithm for learning a perceptron. *Physical Review Letters*, 66(20):2677, 1991.

[56] D. J. Thouless, P. W. Anderson, and R. G. Palmer. Solution of ‘solvable model of a spin glass’. *Philosophical Magazine*, 35(3):593–601, 1977.

[57] Y. Kabashima. Inference from correlated patterns: a unified theory for perceptron learning and linear vector channels. *Journal of Physics: Conference Series*, 95(1):012001, 2008.

[58] D. L. Donoho, A. Javanmard, and A. Montanari. Information-theoretically optimal compressed sensing via spatial coupling and approximate message passing. *IEEE Transactions on Information Theory*, 59(11):7434–7464, Nov 2013.

[59] M. Opper and O. Winther. Mean field approach to bayes learning in feed-forward neural networks. *Physical review letters*, 76(11):1964, 1996.

[60] M. Opper and O. Winther. Tractable approximations for probabilistic models: The adaptive thouless-anderson-palmer mean field approach. *Physical Review Letters*, 86(17):3695, 2001.
[61] M. Mondelli and A. Montanari. Fundamental limits of weak recovery with applications to phase retrieval. *arXiv preprint arXiv:1708.05932*, 2017.

[62] Github repository of numerical experiments for the current paper. https://github.com/sphinxteam/GeneralizedLinearModel2017, 2017.

[63] A. K. Fletcher and S. Rangan. Iterative reconstruction of rank-one matrices in noise. *Information and Inference: A Journal of the IMA*, 7(3):531–562, 2018.

[64] D. Hansel, G. Mato, and C. Meunier. Memorization without generalization in a multilayered neural network. *EPL (Europhysics Letters)*, 20(5):471, 1992.

[65] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python. *Journal of Machine Learning Research*, 12:2825–2830, 2011.

[66] F. Chollet. *keras*. https://github.com/fchollet/keras, 2015.

[67] Y. Wu and S. Verdú. Rényi information dimension: Fundamental limits of almost lossless analog compression. *IEEE Transactions on Information Theory*, 56(8):3721–3748, 2010.

[68] F. Krzakala, M. Mézard, F. Sausset, Y. Sun, and L. Zdeborová. Statistical-physics-based reconstruction in compressed sensing. *Phys. Rev. X*, 2:021005(18), May 2012.

[69] A. Maleki, L. Anitori, Z. Yang, and R. G. Baraniuk. Asymptotic analysis of complex lasso via complex approximate message passing (camp). *IEEE Transactions on Information Theory*, 59(7):4290–4308, 2013.

[70] M. Soltanolkotabi. Structured signal recovery from quadratic measurements: Breaking sample complexity barriers via nonconvex optimization. *arXiv preprint arXiv:1702.06175*, 2017.

[71] F. Rosenblatt. *The perceptron, a perceiving and recognizing automaton Project Para*. Cornell Aeronautical Laboratory, 1957.

[72] F. Guerra and F. L. Toninelli. The thermodynamic limit in mean field spin glass models. *Communications in Mathematical Physics*, 230(1):71–79, 2002.

[73] M. Talagrand. *Mean field models for spin glasses: Volume I: Basic examples*, volume 54. Springer Science & Business Media, 2010.

[74] M. Mezard and A. Montanari. *Information, physics, and computation*. Oxford University Press, 2009.

[75] C. Baldassi, A. Braunstein, N. Brunel, and R. Zecchina. Efficient supervised learning in networks with binary synapses. *Proceedings of the National Academy of Sciences*, 104(26):11079–11084, 2007.

[76] T. Richardson and R. Urbanke. *Modern coding theory*. Cambridge university press, 2008.

[77] M. J. Wainwright, M. I. Jordan, et al. Graphical models, exponential families, and variational inference. *Foundations and Trends® in Machine Learning*, 1(1-2):1–305, 2008.

[78] J. P. Vila and P. Schniter. Expectation-maximization gaussian-mixture approximate message passing. *IEEE Transactions on Signal Processing*, 61(19):4658–4672, 2013.

[79] J. Ziniel, P. Schniter, and P. Sederberg. Binary linear classification and feature selection via generalized approximate message passing. In *Information Sciences and Systems (CISS), 2014 48th Annual Conference on*, pages 1–6. IEEE, 2014.

[80] F. Krzakala, M. Mézard, F. Sausset, Y. Sun, and L. Zdeborová. Probabilistic reconstruction in compressed sensing: algorithms, phase diagrams, and threshold achieving matrices. *Journal of Statistical Mechanics: Theory and Experiment*, 2012(08):P08009(57), 2012.
[81] D. Donoho and J. Tanner. Observed universality of phase transitions in high-dimensional geometry, with implications for modern data analysis and signal processing. *Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences*, 367(1906):4273-4293, 2009.

[82] S. Oymak, A. Jalali, M. Fazel, Y. C. Eldar, and B. Hassibi. Simultaneously structured models with application to sparse and low-rank matrices. *IEEE Transactions on Information Theory*, 61(5):2886-2908, 2015.

[83] E. B. Baum and Y.-D. Lyuu. The transition to perfect generalization in perceptrons. *Neural computation*, 3(3):386-401, 1991.

[84] S. Diamond and S. Boyd. CVXPY: A Python-embedded modeling language for convex optimization. *Journal of Machine Learning Research*, 2016. To appear.

[85] T. Goldstein and C. Studer. Phasemax: Convex phase retrieval via basis pursuit. *arXiv preprint arXiv:1610.07531*, 2016.

[86] J. Barbier, N. Macris, and L. Miolane. The Layered Structure of Tensor Estimation and its Mutual Information. In *47th Annual Allerton Conference on Communication, Control, and Computing (Allerton)*, September 2017.

[87] A. Montanari. Estimating random variables from random sparse observations. *European Transactions on Telecommunications*, 19(4):385-403, 2008.

[88] S. B. Korada and N. Macris. Exact solution of the gauge symmetric p-spin glass model on a complete graph. *Journal of Statistical Physics*, 136(2):205-230, 2009.

[89] Y. Wu and S. Verdú. Optimal phase transitions in compressed sensing. *IEEE Transactions on Information Theory*, 58(10):6241-6263, Oct 2012.

[90] N. Macris. Griffith-kelly-sherman correlation inequalities: A useful tool in the theory of error correcting codes. *IEEE Transactions on Information Theory*, 53(2):664-683, Feb 2007.

[91] S. B. Korada and N. Macris. Tight bounds on the capacity of binary input random cdma systems. *IEEE Transactions on Information Theory*, 56(11):5590-5613, Nov 2010.

[92] P. Hartman. *Ordinary Differential Equations*. Society for Industrial and Applied Mathematics, 2002, 1964.

[93] P. Milgrom and I. Segal. Envelope theorems for arbitrary choice sets. *Econometrica*, 70(2):583-601, 2002.

[94] D. Guo, Y. Wu, S. S. Shitz, and S. Verdú. Estimation in gaussian noise: Properties of the minimum mean-square error. *IEEE Transactions on Information Theory*, 57(4):2371-2385, 2011.

[95] S. B. Korada and A. Montanari. Applications of the Lindeberg principle in communications and statistical learning. *IEEE Transactions on Information Theory*, 57(4):2440-2450, 2011.

[96] Y. Wu and S. Verdú. Functional properties of minimum mean-square error and mutual information. *IEEE Transactions on Information Theory*, 58(3):1289-1301, 2012.

[97] L. Miolane. Phase transitions in spiked matrix estimation: information-theoretic analysis. *arXiv preprint arXiv:1806.04343*, 2018.

[98] V. I. Bogachev. *Measure theory*, volume 1. Springer Science & Business Media, 2007.

[99] M. Lelarge and L. Miolane. Fundamental limits of symmetric low-rank matrix estimation. *ArXiv e-prints*, November 2016.
[100] S. Chatterjee et al. A generalization of the lindeberg principle. *The Annals of Probability*, 34(6):2061–2076, 2006.

[101] R. T. Rockafellar. *Convex analysis*. Princeton university press, 2015.

[102] S. Boucheron, G. Lugosi, and O. Bousquet. Concentration inequalities. In *Advanced Lectures on Machine Learning*, pages 208–240. Springer, 2004.

[103] J. Bezanzon, S. Karpinski, V. Shah, and A. Edelman. Julia: A fast dynamic language for technical computing. In *Lang.NEXT*, April 2012.

[104] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Watenberg, M. Wicke, Y. Yu, and X. Zheng. TensorFlow: Large-scale machine learning on heterogeneous systems, 2015. Software available from tensorflow.org.
# Part II

# Supplementary informations

## 1 Setting

### 1.1 Generalized linear estimation: Problem statement

We give a formal description of the observation model to which our results apply. The generalized linear model covers both the estimation (or inference) problem and the supervised learning problems (see Sec. 3).

Let $n, m \in \mathbb{N}^*$. Let $P_0$ be a probability distribution over $\mathbb{R}$ and let $(X_i^*)_{i=1}^n \stackrel{\text{iid}}{\sim} P_0$ be the components of a signal vector $\mathbf{X}^*$ (this is also denoted $\mathbf{X}^* \stackrel{\text{iid}}{\sim} P_0$). We fix a function $\varphi: \mathbb{R} \times \mathbb{R}^{k_A} \to \mathbb{R}$ and consider $(\mathbf{A}_\mu)^m_{\mu=1} \stackrel{\text{iid}}{\sim} P_A$, where $P_A$ is a probability distribution over $\mathbb{R}^{k_A}$, $k_A \in \mathbb{N}$. We acquire $m$ measurements through

$$Y_{\mu} = \varphi\left(\frac{1}{\sqrt{n}}[\mathbf{\Phi}\mathbf{X}^*]_{\mu}, \mathbf{A}_{\mu}\right) + \sqrt{\Delta}Z_{\mu}, \quad 1 \le \mu \le m, \qquad (20)$$

where $(Z_\mu)_{\mu=1}^m \stackrel{\text{iid}}{\sim} \mathcal{N}(0, 1)$ is an additive Gaussian noise, $\Delta \ge 0$, and $\mathbf{\Phi}$ is a $m \times n$ measurement matrix with independent entries that have zero mean and unit variance. The estimation problem is to recover $\mathbf{X}^*$ from the knowledge of $\mathbf{Y} = (Y_\mu)_{\mu=1}^m$, $\varphi$, $\mathbf{\Phi}$, $\Delta$, $P_0$ and $P_A$ (the realization of the random stream $\mathbf{A}$ itself, if present in the model, is unknown). We use the notation $[\mathbf{\Phi}\mathbf{X}^*]_\mu = \sum_{i=1}^n \Phi_{\mu i} X_i^*$. When $\varphi(x, \mathbf{A}) = \varphi(x) = x$ we have a random linear estimation problem, whereas if, say, $\varphi(x) = \text{sgn}(x)$ we have a noisy single layer perceptron. Sec. 3 discusses various examples related to non-linear estimation and supervised learning.

Denote the prior over the signal as $dP_0(\mathbf{x}) = \prod_{i=1}^n dP_0(x_i)$, and similarly $dP_A(\mathbf{a}) = \prod_{\mu=1}^m dP_A(\mathbf{a}_\mu)$. It is also fruitful to think of the measurements as the outputs of a “channel”,

$$Y_{\mu} \sim P_{\text{out}}\left(\cdot \mid \frac{1}{\sqrt{n}}[\mathbf{\Phi}\mathbf{X}^*]_{\mu}\right). \qquad (21)$$

When $\Delta > 0$ the transition kernel $P_{\text{out}}$ admits a transition density with respect to (w.r.t.) Lebesgue's measure, given by

$$P_{\text{out}}(y|x) = \frac{1}{\sqrt{2\pi\Delta}} \int dP_A(\mathbf{a}) e^{-\frac{1}{2\Delta}(y-\varphi(x,\mathbf{a}))^2}. \qquad (22)$$

When $\Delta = 0$, we will only consider discrete channels where $\varphi$ takes values in $\mathbb{N}^2$. In that case $P_{\text{out}}$ admits a transition density with respect the counting measure on $\mathbb{N}$ given by (here $\mathbf{1}(\cdot)$ is the indicator function)

$$P_{\text{out}}(y|x) = \int dP_A(\mathbf{a}) \mathbf{1}(y = \varphi(x, \mathbf{a})). \qquad (23)$$

Note that for deterministic models, $\mathbf{A}$ in (20) is absent and thus the associated $\int dP_A(\mathbf{a})$ integral in (22)-(23) simply disappears. In fact (20) is sometimes called a “random function representation” of a transition kernel $P_{\text{out}}$. Our analysis uses both representations (20) and (21).

Throughout this paper we often adopt the language of statistical mechanics. In particular the random variables $\mathbf{Y}$ (and also $\mathbf{\Phi}$, $\mathbf{X}^*$, $\mathbf{A}$, $\mathbf{Z}$) are called quenched variables because once the measurements are acquired they have a “fixed realization”. An expectation taken w.r.t. all quenched random variables appearing in an

<footnote n="2">²Notice that this allows to study any channel whose outputs belong to a countable set $S$ by applying a injection $u: S \to \mathbb{N}$ to the outputs.</footnote>
expression will simply be denoted by $\mathbb{E}$ without subscript. Subscripts are only used when the expectation carries over a subset of random variables appearing in an expression or when some confusion could arise.

A fundamental role is played by the posterior distribution of (the signal) $\mathbf{x}$ given the quenched measurements $\mathbf{Y}$ (recall that $\mathbf{X}^*$, $\mathbf{A}$ and $\mathbf{Z}$ are unknown). According to the Bayes formula this posterior is

$$dP(\mathbf{x}|\mathbf{Y}, \mathbf{\Phi}) = \frac{1}{\mathcal{Z}(\mathbf{Y}, \mathbf{\Phi})} dP_0(\mathbf{x}) \prod_{\mu=1}^{m} P_{\text{out}}\left(Y_{\mu} \middle| \frac{1}{\sqrt{n}}[\mathbf{\Phi}\mathbf{x}]_{\mu}\right) \quad (24)$$

$$= \frac{1}{\mathcal{Z}(\mathbf{Y}, \mathbf{\Phi})} dP_0(\mathbf{x}) e^{-\mathcal{H}(\mathbf{x}; \mathbf{Y}, \mathbf{\Phi})} \quad (25)$$

where the Hamiltonian is defined as

$$\mathcal{H}(\mathbf{x}; \mathbf{Y}, \mathbf{\Phi}) := - \sum_{\mu=1}^{m} \ln P_{\text{out}}\left(Y_{\mu} \middle| \frac{1}{\sqrt{n}}[\mathbf{\Phi}\mathbf{x}]_{\mu}\right) \quad (26)$$

and the partition function (the normalization factor) is defined as

$$\mathcal{Z}(\mathbf{Y}, \mathbf{\Phi}) := \int dP_0(\mathbf{x}) e^{-\mathcal{H}(\mathbf{x}; \mathbf{Y}, \mathbf{\Phi})}. \quad (27)$$

From the point of view of statistical mechanics (25) is a Gibbs distribution and the integration over $dP_0(\mathbf{x})$ in the partition function is best thought as a “sum over annealed or fluctuating degrees of freedom”. Let us introduce a standard statistical mechanics notation for the expectation w.r.t. the posterior (24), the so called Gibbs bracket $\langle - \rangle$ defined as

$$\langle g(\mathbf{x}) \rangle := \mathbb{E}[g(\mathbf{X})|\mathbf{Y}, \mathbf{\Phi}] = \int dP(\mathbf{x}|\mathbf{Y}, \mathbf{\Phi}) g(\mathbf{x}) \quad (28)$$

for any continuous bounded function $g$. The main quantity of interest here is the associated averaged free entropy (or minus the averaged free energy)

$$f_n := \frac{1}{n} \mathbb{E} \ln \mathcal{Z}(\mathbf{Y}, \mathbf{\Phi}). \quad (29)$$

It is perhaps useful to stress that $\mathcal{Z}(\mathbf{Y}, \mathbf{\Phi})$ is nothing else than the density of $\mathbf{Y}$ conditioned on $\mathbf{\Phi}$ so we have the explicit representation (used later on)

$$\begin{aligned} f_n &= \frac{1}{n} \mathbb{E}_{\mathbf{\Phi}} \int d\mathbf{Y} \mathcal{Z}(\mathbf{Y}, \mathbf{\Phi}) \ln \mathcal{Z}(\mathbf{Y}, \mathbf{\Phi}) \\ &= \frac{1}{n} \mathbb{E}_{\mathbf{\Phi}} \int d\mathbf{Y} dP_0(\mathbf{X}^*) e^{-\mathcal{H}(\mathbf{X}^*; \mathbf{Y}, \mathbf{\Phi})} \ln \int dP_0(\mathbf{x}) e^{-\mathcal{H}(\mathbf{x}; \mathbf{Y}, \mathbf{\Phi})}, \end{aligned} \quad (30)$$

where $d\mathbf{Y} = \prod_{\mu=1}^{m} dY_{\mu}$. Thus $f_n$ is minus the conditional entropy $-H(\mathbf{Y}|\mathbf{\Phi})/n$ of the measurements. One of the main contributions of this paper is the derivation, thanks to the adaptive interpolation method, of the thermodynamic limit $\lim_{n \to \infty} f_n$ in the “high-dimensional regime”, namely when $n, m \to \infty$ while $m/n \to \alpha > 0$ ($\alpha$ is sometimes referred to as the “measurement rate” or “sampling rate”).

## 1.2 The teacher-student scenario

We now describe an important conceptual setting, the teacher-student scenario (also called planted model), that allows to then define the optimal generalization error. We voluntarily employ terms coming from machine learning instead of the signal processing terminology used until here.
First the teacher randomly generates a classifier $\mathbf{X}^* \in \mathbb{R}^n$ (the signal in the estimation problem) with $\mathbf{X}^* \stackrel{\text{iid}}{\sim} P_0$ and an ensemble of $m$ patterns (row-vectors) $\mathbf{\Phi}_\mu \in \mathbb{R}^n$ for $\mu = 1, \dots, m$ such that $\mathbf{\Phi}_\mu \stackrel{\text{iid}}{\sim} \mathcal{N}(\mathbf{0}, \mathbf{I}_n)$. The teacher then chooses a model $(\varphi, P_A, \Delta)$ or equivalently $P_{\text{out}}$, which are linked through (22)-(23). The teacher then output labels $Y_\mu \in \mathbb{R}$ through (20) or (21) for $\mu = 1, \dots, m$.

The student is given the distribution $P_0$, the model $(\varphi, P_A, \Delta)$ or equivalently $P_{\text{out}}$ and the training data composed of the pattern-label pairs $\{(Y_\mu, \mathbf{\Phi}_\mu)\}_{\mu=1}^m$ generated by the teacher. His (supervised) learning task is then to predict the labels associated with new, yet unseen, patterns from all this knowledge.

How does the teacher may evaluate the student's prediction capabilities? The teacher starts by randomly generating a new line of the matrix, or pattern, $\mathbf{\Phi}_{\text{new}}$. Then, still using the same $\mathbf{X}^*$, he generates the associated new label $Y_{\text{new}} \sim P_{\text{out}}(\cdot | \mathbf{\Phi}_{\text{new}} \cdot \mathbf{X}^*/\sqrt{n})$. He is now ready to evaluate the student generalization performance. For that purpose, an important quantity is the *generalization error* (or prediction error). If we denote $\hat{Y}(\mathbf{\Phi}_{\text{new}}, \mathbf{\Phi}, \mathbf{Y})$ the estimator used by the student (which is thus a measurable function of the observations), the generalization error is defined as

$$ \mathcal{E}_{\text{gen}}(\hat{Y}) := \mathbb{E}[(Y_{\text{new}} - \hat{Y}(\mathbf{\Phi}_{\text{new}}, \mathbf{\Phi}, \mathbf{Y}))^2]. \quad (31) $$

The optimal generalization error is then defined as the minimum of $\mathcal{E}_{\text{gen}}$ over all estimators $\hat{Y}(\mathbf{\Phi}_{\text{new}}, \mathbf{\Phi}, \mathbf{Y})$:

$$ \mathcal{E}_{\text{gen}}^{\text{opt}} := \min_{\hat{Y}} \mathcal{E}_{\text{gen}}(\hat{Y}) = \text{MMSE}(Y_{\text{new}} | \mathbf{\Phi}_{\text{new}}, \mathbf{\Phi}, \mathbf{Y}) = \mathbb{E}[(Y_{\text{new}} - \mathbb{E}[Y_{\text{new}} | \mathbf{\Phi}_{\text{new}}, \mathbf{\Phi}, \mathbf{Y}])^2]. \quad (32) $$

Here, and for the rest of the paper, we define the minimum mean-square error (MMSE) function as follows: Given two random variables $\mathbf{A}, \mathbf{B}$, the MMSE in estimating $\mathbf{A}$ given $\mathbf{B}$ is defined as

$$ \text{MMSE}(\mathbf{A}|\mathbf{B}) := \mathbb{E}[\|\mathbf{A} - \mathbb{E}[\mathbf{A}|\mathbf{B}]\|^2], \quad (33) $$

where $\mathbb{E}[\mathbf{A}|\mathbf{B}]$ is the expectation of $\mathbf{A}$ with respect to its posterior given $\mathbf{B}$.

*A word about notations:* Let us emphasize on the link between the different notations that we use in the present supplementary material and in the main text. E.g., the expectation w.r.t. to the posterior of $Y_{\text{new}}$ appearing in (32) can be written equivalently as:

$$ \mathbb{E}[Y_{\text{new}} | \mathbf{\Phi}_{\text{new}}, \mathbf{\Phi}, \mathbf{Y}] = \mathbb{E}_{P_A(\mathbf{a})} \mathbb{E}_{P(\mathbf{x}|\mathbf{\Phi}, \mathbf{Y})} \varphi\left(\frac{\mathbf{\Phi}_{\text{new}} \cdot \mathbf{x}}{\sqrt{n}}, \mathbf{a}\right) = \mathbb{E}_{P_A(\mathbf{a})} \left\langle \varphi\left(\frac{\mathbf{\Phi}_{\text{new}} \cdot \mathbf{x}}{\sqrt{n}}, \mathbf{a}\right) \right\rangle. \quad (34) $$

To see that, just write:

$$ \begin{aligned} \mathbb{E}[Y_{\text{new}} | \mathbf{\Phi}_{\text{new}}, \mathbf{\Phi}, \mathbf{Y}] &:= \int dY_{\text{new}} Y_{\text{new}} P(Y_{\text{new}} | \mathbf{\Phi}_{\text{new}}, \mathbf{\Phi}, \mathbf{Y}) \\ &= \int dY_{\text{new}} Y_{\text{new}} P_{\text{out}}\left(Y_{\text{new}} \middle| \frac{1}{\sqrt{n}} \mathbf{\Phi}_{\text{new}} \cdot \mathbf{x}\right) dP(\mathbf{x}|\mathbf{Y}, \mathbf{\Phi}) \\ &= \left\langle \int dY_{\text{new}} Y_{\text{new}} P_{\text{out}}\left(Y_{\text{new}} \middle| \frac{1}{\sqrt{n}} \mathbf{\Phi}_{\text{new}} \cdot \mathbf{x}\right) \right\rangle \\ &= \left\langle \int dP_A(\mathbf{a}) dY_{\text{new}} Y_{\text{new}} \frac{1}{\sqrt{2\pi\Delta}} e^{-\frac{1}{2\Delta}\{Y_{\text{new}} - \varphi\left(\frac{1}{\sqrt{n}}\mathbf{\Phi}_{\text{new}} \cdot \mathbf{x}, \mathbf{a}\right)\}^2} \right\rangle \\ &= \mathbb{E}_{P_A(\mathbf{a})} \left\langle \varphi\left(\frac{1}{\sqrt{n}}\mathbf{\Phi}_{\text{new}} \cdot \mathbf{x}, \mathbf{a}\right) \right\rangle. \end{aligned} \quad (35) $$

Here we used definition (22) for the transition kernel, but using instead (23) would lead to the same identity.
## 1.3 Two scalar inference channels

An important role in our proof of the asymptotic expression of the free entropy is played by simple scalar inference channels. As we will see, the free entropy is expressed in terms of the free entropy of these channels. This “decoupling property” results from the mean-field approach in statistical physics, used through in the replica method to perform a formal calculation of the free entropy of the model [49,74]. Let us now introduce these two scalar denoising models.

The first one is an additive Gaussian channel. Let $r \ge 0$, which plays the role of a signal-to-noise ratio (snr). Suppose that $X_0 \sim P_0$ and that we observe

$$Y_0 = \sqrt{r} X_0 + Z_0, \qquad (36)$$

where $Z_0 \sim \mathcal{N}(0, 1)$ independently of $X_0$. Consider the inference problem consisting of retrieving $X_0$ from the observations $Y_0$. The associated posterior distribution is

$$dP(x|Y_0) = \frac{dP_0(x)e^{\sqrt{r}Y_0x - rx^2/2}}{\int dP_0(x)e^{\sqrt{r}Y_0x - rx^2/2}}. \qquad (37)$$

In this expression all the $x$-independent terms have been simplified between the numerator and the normalization. The free entropy associated with this channel is just the expectation of the logarithm of the normalization factor

$$\psi_{P_0}(r) := \mathbb{E} \ln \int dP_0(x) e^{\sqrt{r} Y_0 x - rx^2/2}. \qquad (38)$$

The basic properties of $\psi_{P_0}$ are presented in Appendix B.1.

The second scalar channel that appears naturally in the problem is linked to the transition kernel $P_{\text{out}}$ through the following inference model. Suppose that $V, W^* \stackrel{\text{iid}}{\sim} \mathcal{N}(0, 1)$ where $V$ is known while the inference problem is to recover the unknown $W^*$ from the following observation

$$\tilde{Y}_0 \sim P_{\text{out}}(\cdot | \sqrt{q}V + \sqrt{\rho-q}W^*), \qquad (39)$$

where $\rho > 0, q \in [0, \rho]$. Notice that the channel (39) is equivalent to $\tilde{Y}_0 = \varphi(\sqrt{q}V + \sqrt{\rho-q}W^*, \mathbf{A}) + \sqrt{\Delta}Z$ with $\Delta \ge 0$ and where $(\mathbf{A}, Z) \sim P_A \otimes \mathcal{N}(0, 1)$, independently of $V, W^*$. The free entropy for this model, again related to the normalization of the posterior $dP(w|\tilde{Y}_0, V)$, is

$$\Psi_{P_{\text{out}}}(q; \rho) = \Psi_{P_{\text{out}}}(q) := \mathbb{E} \ln \int \mathcal{D}w P_{\text{out}}(\tilde{Y}_0 | \sqrt{q}V + \sqrt{\rho-q}w), \qquad (40)$$

where $\mathcal{D}w := dw(2\pi)^{-1/2}e^{-w^2/2}$ is the standard Gaussian measure. In (40) above, $P_{\text{out}}$ denotes either the transition density with respect to Lebesgue's measure (given by (22)) in the case $\Delta > 0$, or the density with respect to the counting measure over $\mathbb{N}$ (given by (23)), in the case of a “discrete” channel ($\varphi$ takes values in $\mathbb{N}$ and $\Delta = 0$). We prove in Appendix B.2 that this function is convex, differentiable and non-decreasing w.r.t. its first argument.
# 2 Main results

## 2.1 Replica-symmetric formula and mutual information

Let us now introduce our first main result, a single-letter replica-symmetric formula for the asymptotic free entropy of model (20), (21). The result holds under the following rather general hypotheses. We will consider two cases, that is when there is some Gaussian noise ($\Delta > 0$, see (h5.a) below) and the case without Gaussian noise ($\Delta = 0$, see (h5.b) below):

(h1) The prior distribution $P_0$ admits a finite third moment and has at least two points in its support.

(h2) There exists $\gamma > 0$ such that the sequence $(\mathbb{E}[|\varphi(\frac{1}{\sqrt{n}}[\mathbf{\Phi}\mathbf{X}^*], \mathbf{A}_1)|^{2+\gamma}])_{n\ge 1}$ is bounded.

(h3) The random variables $(\Phi_{\mu i})$ are independent with zero mean, unit variance and finite third moment that is bounded with $n$.

(h4) For almost-all values of $\mathbf{a} \in \mathbb{R}^{k_A}$ (w.r.t. $P_A$), the function $x \mapsto \varphi(x, \mathbf{a})$ is continuous almost everywhere.

We will also assume that one of the two following hypotheses hold:

(h5.a) $\Delta > 0$.

(h5.b) $\Delta = 0$ and $\varphi$ takes values in $\mathbb{N}$.

**Remark 1.** The above hypotheses are here stated using the “random function” representation of (20). In many cases, it can be useful to state them using the “transition kernel” representation of (21). The hypotheses (h2) and (h4) are respectively equivalent<fnref n="3" /> to:

(h2') There exists $\gamma > 0$ such that $\mathbb{E}[|Y_1|^{2+\gamma}]$ remains bounded with $n$.

(h4') $x \in \mathbb{R} \mapsto P_{\text{out}}(\cdot|x)$ is continuous almost everywhere for the weak convergence.

Under the above hypothesis (h5.a) (respectively (h5.b)), the transition kernel $P_{\text{out}}$ admits a density with respect to Lebesgue's measure on $\mathbb{R}$ (resp. the counting measure on $\mathbb{N}$) that will be denoted by $P_{\text{out}}(\cdot|x)$. We will call the kernel $P_{\text{out}}$ informative if there exists $y \in \mathbb{R}$ (resp. $y \in \mathbb{N}$) such that $P_{\text{out}}(y|\cdot)$ is not equal almost everywhere to a constant. If $P_{\text{out}}$ is not informative, it is not difficult to show that estimation is then impossible.

Let us define the replica-symmetric potential (or just potential). Call $\rho := \mathbb{E}[(X^*)^2]$ where $X^* \sim P_0$. Then the potential is

$$f_{\text{RS}}(q, r; \rho) = f_{\text{RS}}(q, r) := \psi_{P_0}(r) + \alpha \Psi_{P_{\text{out}}}(q; \rho) - \frac{rq}{2}. \quad (41)$$

We define also $f_{\text{RS}}(\rho, +\infty) = \lim_{r \to \infty} f_{\text{RS}}(\rho, r)$. From now on denote $\psi'_{P_0}(r)$ and $\Psi'_{P_{\text{out}}}(q) = \Psi'_{P_{\text{out}}}(q; \rho)$ the derivatives of $\psi_{P_0}(r)$ and $\Psi_{P_{\text{out}}}(q; \rho)$ w.r.t. their first argument. We need also to define the set of the critical points of $f_{\text{RS}}$:

$$\Gamma := \left\{ (q, r) \in [0, \rho] \times (\mathbb{R}_+ \cup \{+\infty\}) \mid \begin{array}{l} q = 2\psi'_{P_0}(r) \\ r = 2\alpha\Psi'_{P_{\text{out}}}(q; \rho) \end{array} \right\}, \quad (42)$$

where, with a slight abuse of notation, we define $\psi'_{P_0}(+\infty) = \lim_{r \to \infty} \psi'_{P_0}(r)$ and $\Psi'_{P_{\text{out}}}(\rho) = \lim_{q \to \rho} \Psi'_{P_{\text{out}}}(q)$. These limits are well defined by convexity of $\psi_{P_0}$ and $\Psi_{P_{\text{out}}}$. The elements of $\Gamma$ are called “fixed points of the state evolution”. Our first main result is

<footnote n="3">³The implications (h2) ⇔ (h2') and (h4) ⇒ (h4') are obvious. If (h4') holds one can show, by inverting cumulative distribution functions, that there exists a function $\varphi: \mathbb{R} \times [0, 1] \to \mathbb{R}$ such that (20) holds for $A_\mu \stackrel{\text{iid}}{\sim} P_A = \text{Unif}([0, 1])$ and that (h4) is verified.</footnote>
**Theorem 1** (Replica-symmetric formula for the free entropy). Suppose that hypotheses (h1)-(h2)-(h3)-(h4) hold. Suppose that either hypothesis (h5.a) or (h5.b) holds. Then, for the generalized linear estimation model (20), (21) the thermodynamic limit of the free entropy (29) verifies

$$f_{\infty} := \lim_{n \to \infty} f_n = \sup_{q \in [0, \rho]} \inf_{r \ge 0} f_{\text{RS}}(q, r) = \sup_{(q,r) \in \Gamma} f_{\text{RS}}(q, r). \quad (43)$$

Moreover, if $P_{\text{out}}$ is informative, then the “sup inf” and the supremum over $\Gamma$ in (43) are achieved over the same couples $(q, r)$.

An immediate corollary of Theorem 1 is the limiting expression of the mutual information between the signal and the observations. To state the result, we need to introduce two mutual informations associated to the two scalar channels presented in Sec. 1.3, namely

$$I_{P_0}(r) := I(X_0; \sqrt{r} X_0 + Z_0) = \frac{r\rho}{2} - \psi_{P_0}(r) \quad (44)$$

for the channel (36) and

$$\mathcal{I}_{P_{\text{out}}}(q) := I(W^*; \tilde{Y}_0|V) = \Psi_{P_{\text{out}}}(\rho) - \Psi_{P_{\text{out}}}(q) \quad (45)$$

for the channel (39).

**Corollary 2** (Single-letter formula for the mutual information). The thermodynamic limit of the mutual information for model (20), (21) between the observations and the hidden variables verifies

$$i_{\infty} := \lim_{n \to \infty} \frac{1}{n} I(\mathbf{X}^*; \mathbf{Y} | \mathbf{\Phi}) = \inf_{q \in [0, \rho]} \sup_{r \ge 0} i_{\text{RS}}(q, r) = \inf_{(q,r) \in \Gamma} i_{\text{RS}}(q, r), \quad (46)$$

where

$$i_{\text{RS}}(q, r) := I_{P_0}(r) + \alpha \mathcal{I}_{P_0}(q) - \frac{r}{2}(\rho - q). \quad (47)$$

*Proof.* This follows from a simple calculation:

$$\begin{aligned} \frac{1}{n} I(\mathbf{X}^*; \mathbf{Y} | \mathbf{\Phi}) &= \frac{1}{n} H(\mathbf{Y} | \mathbf{\Phi}) - \frac{1}{n} H(\mathbf{Y} | \mathbf{X}^*, \mathbf{\Phi}) = -f_n + \frac{1}{n} \mathbb{E} \ln P(\mathbf{Y} | \mathbf{X}^*, \mathbf{\Phi}) \\ &= -f_n + \frac{m}{n} \mathbb{E} \ln P_{\text{out}}(Y_1 | \mathbf{\Phi}_1 \cdot \mathbf{X}^*/\sqrt{n}). \end{aligned} \quad (48)$$

By the central limit theorem (that we can apply under hypotheses (h1)-(h3)) we have

$$S_n := \frac{1}{\sqrt{n}} \mathbf{\Phi}_1 \cdot \mathbf{X}^* = \frac{1}{\sqrt{n}} \sum_{i=1}^{n} \Phi_{1,i} X_i^* \xrightarrow{n \to \infty} \mathcal{N}(0, \rho).$$

Now, under the hypotheses (h2)-(h4) and either (h5.a) or (h5.b) it is not difficult to verify that

$$\begin{aligned} \mathbb{E} \ln P_{\text{out}}(Y_1 | \mathbf{\Phi}_1 \cdot \mathbf{X}^*/\sqrt{n}) &= \mathbb{E} \int dY P_{\text{out}}(Y | S_n) \ln P_{\text{out}}(Y | S_n) \\ &\xrightarrow{n \to \infty} \mathbb{E} \int dY P_{\text{out}}(Y | \sqrt{\rho} V) \ln P_{\text{out}}(Y | \sqrt{\rho} V) = \Psi_{P_{\text{out}}}(\rho) \end{aligned}$$

where $V \sim \mathcal{N}(0, 1)$. We conclude, using (48):

$$\frac{1}{n} I(\mathbf{X}^*; \mathbf{Y} | \mathbf{\Phi}) = -f_n + \alpha \Psi_{P_{\text{out}}}(\rho) + o_n(1) \quad (49)$$

where $\lim_{n \to \infty} o_n(1) = 0$.

<figure><img src="image_4.png" /></figure>
The next proposition, proved in Appendix A.2, states that for almost every $\alpha > 0$ there is one unique optimizer $q^*$ in (43) (or equivalently in (46)):

**Proposition 1.** Assume that the assumptions of Theorem 1 hold and that $P_{\text{out}}$ is informative. Define

$$D^* := \{\alpha > 0 \mid (43) \text{ (or equivalently (46)) admits a unique optimizer } q^*(\alpha)\}. \quad (50)$$

the set $D^*$ is equal to $\mathbb{R}_+^*$ minus some countable set. Moreover $\alpha \mapsto q^*(\alpha)$ is continuous on $D^*$.

As an application of Theorem 1 we can compute the free entropy of the “planted perceptron” on the hypercube and the sphere. This perceptron model has already been studied in physics [22] and more recently in statistics, where it is known as “one-bit compressed sensing” [6, 15]. The limit of the free entropy follows from an application of Theorem 1 with $\varphi(x) = \text{sgn}(x)$ and $P_0 = \frac{1}{2}\delta_{-1} + \frac{1}{2}\delta_1$ (for the hypercube) or $P_0 = \mathcal{N}(0, 1)$ (for the sphere). For $\mu \in \{1, \dots, m\}$ we define

$$S_{\mu} := \{ \mathbf{x} \in \mathbb{R}^n \mid \text{sgn}(\mathbf{x} \cdot \mathbf{\Phi}_{\mu}) = \text{sgn}(\mathbf{X}^* \cdot \mathbf{\Phi}_{\mu}) \}. \quad (51)$$

We will use the notation $\mathcal{N}(x) = \mathbb{P}(Z \le x)$ for $Z \sim \mathcal{N}(0, 1)$. Let $\mathbb{S}_n$ be the unit sphere in $\mathbb{R}^n$ and $\mu_n$ the uniform probability measure on $\mathbb{S}_n$.

**Corollary 3 (Free entropy of the planted perceptron).** Let $Z, V \stackrel{\text{id}}{\sim} \mathcal{N}(0, 1)$. We have

$$\frac{1}{n} \mathbb{E} \ln \left( \# \bigcap_{\mu=1}^{m} S_{\mu} \cap \{-1, 1\}^n \right) \xrightarrow{n \to \infty} \ln(2) + \sup_{q \in [0,1]} \inf_{r \ge 0} \left\{ \mathbb{E} \ln \cosh(\sqrt{r}Z + r) + 2\alpha \mathbb{E} \left[ \mathcal{N}\left(\frac{\sqrt{q}V}{\sqrt{1-q}}\right) \ln \mathcal{N}\left(\frac{\sqrt{q}V}{\sqrt{1-q}}\right) \right] - \frac{r(q+1)}{2} \right\}, \quad (52)$$

$$\frac{1}{n} \mathbb{E} \ln \mu_n \left( \bigcap_{\mu=1}^{m} S_{\mu} \cap \mathbb{S}_n \right) \xrightarrow{n \to \infty} \sup_{q \in [0,1)} \left\{ \frac{1}{2} \ln(1-q) + 2\alpha \mathbb{E} \left[ \mathcal{N}\left(\frac{\sqrt{q}V}{\sqrt{1-q}}\right) \ln \mathcal{N}\left(\frac{\sqrt{q}V}{\sqrt{1-q}}\right) \right] + \frac{q}{2} \right\}. \quad (53)$$

## 2.2 Optimal reconstruction (or estimation) error

We first consider the problem of estimating $\mathbf{X}^*$ given $\mathbf{Y}$ and $\mathbf{\Phi}$. The following theorem states that the optimizer $q^*(\alpha)$ of the replica-symmetric formula (43) gives the asymptotic correlation between the planted solution $\mathbf{X}^*$ and a typical sample from the posterior distribution $P(\cdot | \mathbf{Y}, \mathbf{\Phi})$:

**Theorem 2 (Limit of the overlap).** Assume that all the moments of $P_0$ are finite and that $P_{\text{out}}$ is informative. Assume that (h1)-(h2)-(h3)-(h4) hold and that either (h5.a) or (h5.b) holds. Then for all $\alpha \in D^*$,

$$\frac{1}{n} |\mathbf{x} \cdot \mathbf{X}^*| = \frac{1}{n} \left| \sum_{i=1}^{n} x_i X_i^* \right| \xrightarrow{n \to \infty} q^*(\alpha), \quad \text{in probability,} \quad (54)$$

where $\mathbf{x} = (x_1, \dots, x_n)$ is sampled from the posterior distribution of the signal $P(\cdot | \mathbf{Y}, \mathbf{\Phi})$ given by (24), independently of everything else.

Theorem 2 is proved in Sec. 5.3. Notice that in all generality it is only possible to estimate $\mathbf{X}^*$ up to its sign (think for instance to $\mathbf{Y} = |\mathbf{\Phi}\mathbf{X}^*|/\sqrt{n} + \sqrt{\Delta Z}$), this is why the absolute values in (54) are needed. For this reason, the usual MSE on $\mathbf{X}^*$

$$\text{mse}(\hat{\mathbf{X}}) := \frac{1}{n} \mathbb{E} \left[ \| \mathbf{X}^* - \hat{\mathbf{X}}(\mathbf{Y}, \mathbf{\Phi}) \|^2 \right]$$
is not (in all generality) an appropriate error metric. Indeed, in the case where $\mathbf{Y} = |\mathbf{\Phi}\mathbf{X}^*|/\sqrt{n} + \sqrt{\Delta}\mathbf{Z}$, where $\mathbf{\Phi}, \mathbf{X}^*, \mathbf{Z}$ have all independent $\mathcal{N}(0, 1)$ entries, then $\mathbb{E}[\mathbf{X}^*|\mathbf{Y}, \mathbf{\Phi}] = 0$ and $\min_{\hat{\mathbf{X}}} \text{mse}(\hat{\mathbf{X}}) = 1$. This means that the minimum mean-square error is always equal to the variance and thus, in this sense, it is never possible to estimate the signal better than trivial estimators. For this reason, the appropriate error metric for the reconstruction problem is the MSE on $\mathbf{X}^*\mathbf{X}^*^\top$. From Theorem 2 one deduces the limit of the MMSE in estimating $\mathbf{X}^*\mathbf{X}^*^\top$:

**Corollary 4 (Matrix minimum mean-square error).** Under the same conditions as in Theorem 2, for all $\alpha \in D^*$ we have

$$ \text{MMSE}_n := \frac{1}{n^2} \mathbb{E} \left[ \| \mathbf{X}^* \mathbf{X}^*^\top - \mathbb{E}[\mathbf{X}^* \mathbf{X}^*^\top | \mathbf{Y}, \mathbf{\Phi}] \|_{\text{F}}^2 \right] \underset{n \to \infty}{\longrightarrow} \rho^2 - q^*(\alpha)^2, \quad (55) $$

where $\| - \|_{\text{F}}$ denotes the Frobenius norm.

## 2.3 Optimal generalization (or prediction) error

In order to express the optimal generalization error we introduce the following function (recall that $\tilde{Y}_0$ is drawn from the channel (39)):

$$ \mathcal{E}(q) := \text{MMSE}(\tilde{Y}_0|V) = \mathbb{E}[(\tilde{Y}_0 - \mathbb{E}[\tilde{Y}_0|V])^2] \quad (56) $$

$$ = \mathbb{E}_V \int dY Y^2 P_{\text{out}}(Y|\sqrt{\rho}V) - \mathbb{E}_V \left[ \mathbb{E}_W \left[ \int dY Y P_{\text{out}}(Y|\sqrt{q}V + \sqrt{\rho-q}W) \right]^2 \right] \quad (57) $$

$$ = \mathbb{E}[\varphi(\sqrt{\rho}V, \mathbf{A})^2] - \mathbb{E}_V [\mathbb{E}_{W,\mathbf{A}}[\varphi(\sqrt{q}V + \sqrt{\rho-q}W, \mathbf{A})]^2] + \Delta \quad (58) $$

where $V, W \stackrel{\text{iid}}{\sim} \mathcal{N}(0, 1)$, $\mathbf{A} \sim P_A$ are independent random variables, $\mathbb{E}_{W,\mathbf{A}}$ denotes the expectation w.r.t. $W$ and $\mathbf{A}$ only and $\mathbb{E}_W[-]^2 = (\mathbb{E}_W[-])^2$. We recall $\rho := \mathbb{E}[(X^*)^2]$ with $X^* \sim P_0$. With a slight abuse of notation $\int dY$ denotes in (57) either the integration w.r.t. Lebesgue's measure on $\mathbb{R}$ in the case $\Delta > 0$ or the integration w.r.t. the counting measure on $\mathbb{N}$ (in the case $\Delta = 0$).

Recall the teacher-student setting of Sec. 1.2: The generalization error is related to the estimation of a new output $Y_{\text{new}} \sim P_{\text{out}}(\cdot | \mathbf{\Phi}_{\text{new}} \cdot \mathbf{X}^*/\sqrt{n})$ where $\mathbf{\Phi}_{\text{new}}$ is a new row of the matrix, and is defined by (32).

**Theorem 3 (Optimal generalization error).** Assume that $P_{\text{out}}$ is informative, that (h1)-(h2)-(h3)-(h4) hold and that either (h5.a) or (h5.b) hold. Then for all $\alpha \in D^*$ we have

$$ \mathcal{E}_{\text{gen}}^{\text{opt}}(\alpha) \underset{n \to \infty}{\longrightarrow} \mathcal{E}(q^*(\alpha)) \quad (59) $$

where $q^*(\alpha)$ is the optimizer of the replica-symmetric formula (43), see Proposition 1.

Theorem 3 follows from a more general result, that we state now. Let $f: \mathbb{R} \to \mathbb{R}$ and consider the generalized optimal generalization error

$$ \mathcal{E}_{f,n}(\alpha) := \text{MMSE}(f(Y_{\text{new}})|\mathbf{\Phi}_{\text{new}}, \mathbf{Y}, \mathbf{\Phi}) = \mathbb{E} \left[ (f(Y_{\text{new}}) - \mathbb{E}[f(Y_{\text{new}})|\mathbf{\Phi}_{\text{new}}, \mathbf{Y}, \mathbf{\Phi}])^2 \right] \quad (60) $$

which is the minimum mean-square error on $f(Y_{\text{new}})$. In particular $\mathcal{E}_{\text{gen}}^{\text{opt}}(\alpha) = \mathcal{E}_{f,n}(\alpha)$ for $f: x \mapsto x$. We define also

$$ \mathcal{E}_f(q) := \text{MMSE}(f(\tilde{Y}_0)|V) = \mathbb{E}[(f(\tilde{Y}_0) - \mathbb{E}[f(\tilde{Y}_0)|V])^2] \quad (61) $$

$$ = \mathbb{E}[f(\varphi(\sqrt{\rho}V, \mathbf{A}) + \sqrt{\Delta}Z)^2] - \mathbb{E}_V [\mathbb{E}_{W,Z,\mathbf{A}}[f(\varphi(\sqrt{q}V + \sqrt{\rho-q}W, \mathbf{A}) + \sqrt{\Delta}Z)]^2], \quad (62) $$

where $\tilde{Y}_0$ is the output of the second scalar channel (39).
**Theorem 4** (Generalized optimal generalization error). Let $f: \mathbb{R} \to \mathbb{R}$ be a measurable function such that $\mathbb{E}[|f(Y_1)|^{2+\gamma}]$ remains bounded as $n$ grows, for some $\gamma > 0$. Assume that $P_{\text{out}}$ is informative and that (h1)-(h2)-(h3)-(h4) hold and that either (h5.a) or (h5.b) holds. Then for all $\alpha \in D^*$ we have

$$ \mathcal{E}_{f,n}(\alpha) \xrightarrow{n \to \infty} \mathcal{E}_f(q^*(\alpha)) \qquad (63) $$

where $q^*(\alpha)$ is the optimizer of the replica-symmetric formula (43), see Proposition 1.

Theorem 4 is proved in Sec. 5.1.

## 2.4 Optimality of the generalized approximate message-passing algorithm

### 2.4.1 The generalized approximate message-passing algorithm

While the main results presented until now are information-theoretic, our next one concerns the performance of a popular algorithm to solve random instances of generalized linear problems, called generalized approximate message-passing (GAMP). We shall not re-derive its properties here, and instead refer to the original papers for details. This approach has a long history, especially in statistical physics [41, 56, 57, 75], error correcting codes [76], and graphical models [77]. For a modern derivation in the context of linear models, see [11, 68, 78]. The case of generalized linear models was discussed by Rangan in [12], and has been used for classification purpose in [79].

We first need to define two so-called threshold functions that are associated to the two scalar channels (36) and (39). The first one is the posterior mean of the signal in channel (36) with signal-to-noise ratio $r$:

$$ g_{P_0}(y, r) := \mathbb{E}[X_0|Y_0=y]. \qquad (64) $$

The second one is the posterior mean of $W^*$ in channel (39) with “noise level” $\eta = \rho - q$:

$$ g_{P_{\text{out}}}(\tilde{y}, v, \eta) := \mathbb{E}[W^*|\tilde{Y}_0 = \tilde{y}, \sqrt{q}V = v]. \qquad (65) $$

These functions act componentwise when applied to vectors.

Given initial estimates $(\hat{\mathbf{x}}^0, \mathbf{v}^0)$ for the means and variances of the elements of the signal vector $\mathbf{X}^*$, GAMP takes as input the observation vector $\mathbf{Y}$ and then iterates the following equations with initialization $g_{\mu}^0 = 0$ for all $\mu = 1, \dots, m$ (we denote by $\bar{\mathbf{u}}$ the average over all the components of the vector $\mathbf{u}$ and $\mathbf{\Phi}^\top$ is the transpose of the matrix $\mathbf{\Phi}$): From $t=1$ until convergence,

$$ \left\{ \begin{array}{l} V^t = \overline{\mathbf{v}^{t-1}} \\ \boldsymbol{\omega}^t = \mathbf{\Phi} \hat{\mathbf{x}}^{t-1} / \sqrt{n} - V^t \mathbf{g}^{t-1} \\ g_{\mu}^t = g_{P_{\text{out}}} (Y_{\mu}, \boldsymbol{\omega}_{\mu}^t, V^t) \qquad \forall \mu = 1, \dots, m \\ \lambda^t = \alpha g_{P_{\text{out}}}^2 (\mathbf{Y}, \boldsymbol{\omega}^t, V^t) \\ \mathbf{R}^t = \hat{\mathbf{x}}^{t-1} + (\lambda^t)^{-1} \mathbf{\Phi}^\top \mathbf{g}^t / \sqrt{n} \\ \hat{\mathbf{x}}_i^t = g_{P_0} (R_i^t, \lambda^t) \qquad \forall i = 1, \dots, n \\ \mathbf{v}_i^t = (\lambda^t)^{-1} \partial_R g_{P_0} (R, \lambda^t) \Big|_{R=R_i^t} \qquad \forall i = 1, \dots, n \end{array} \right. \qquad (66) $$

One of the strongest asset of GAMP is that its performance can be tracked rigorously in the limit $n, m \to \infty$ while $m/n \to \alpha$ via a procedure known as state evolution (SE), see [43, 44] for the linear case, and [12, 45] for the generalized one. In our notations, state evolution tracks the asymptotic value of the overlap between the true hidden value $\mathbf{X}^*$ and its estimate by GAMP $\hat{\mathbf{x}}^t$ defined as $q^t := \lim_{n \to \infty} \mathbf{X}^* \cdot \hat{\mathbf{x}}^t/n$ (that is related to the asymptotic mean-square error (MSE) $E^t$ between $\mathbf{X}^*$ and its estimate $\hat{\mathbf{x}}^t$ by $E^t = \rho - q^t$, where recall that
$\rho := \mathbb{E}[(X^*)^2]$ with $X^* \sim P_0$) via:

$$ \left\{ \begin{array}{l} q^{t+1} = 2\psi'_{P_0}(r^t), \\ r^t = 2\alpha\Psi'_{P_{\text{out}}}(q^t; \rho). \end{array} \right. \qquad (67) $$

From Theorem 1 we realize that the fixed points of these equations correspond to the critical points of the asymptotic free entropy in (43). In fact, in the replica heuristic, the optimizer $q^*$ of the potential is conjectured to give the optimal value of the overlap, a fact that was proven for the linear channel [46-48]. We will see in Sec. 3 that $q^* \xrightarrow{t \to \infty} q^*$ for a large set of parameters.

### 2.4.2 Estimation and generalization error of GAMP

Perhaps more surprisingly, one can use GAMP in the teacher-student scenario described in Sec. 1.2 in order to provide an estimation of a new output $Y_{\text{new}} \sim P_{\text{out}}(\cdot | \mathbf{\Phi}_{\text{new}} \cdot \mathbf{X}^*/\sqrt{n})$ where $\mathbf{\Phi}_{\text{new}}$ is a new row of the matrix. As $\hat{\mathbf{x}}^t$ is the GAMP estimate of the posterior expectation of $\mathbf{X}^*$, with estimated variance $\mathbf{v}^t$, the natural heuristic is to consider for the posterior probability distribution of the random variable $\mathbf{\Phi}_{\text{new}} \cdot \mathbf{X}^*/\sqrt{n}$ a Gaussian with mean $\mathbf{\Phi}_{\text{new}} \cdot \hat{\mathbf{x}}^{t-1}/\sqrt{n}$ and variance $V^t = E^t = \rho - q^t$ (the fact that the variance and MSE are equal follows from the Nishimori identity of Proposition 12 but applied to GAMP instead of the Gibbs measure, see e.g. [13] where this is shown). This allows to estimate the posterior mean of the output, which leads to the GAMP prediction (recall the $P_{\text{out}}$ definition (22)-(23)):

$$ \hat{Y}^{\text{GAMP},t} := \int y P_{\text{out}}\left(y \mid \frac{1}{\sqrt{n}} \mathbf{\Phi}_{\text{new}} \cdot \hat{\mathbf{x}}^{t-1} + \sqrt{\rho - q^t} w\right) \mathcal{D}w dy, \qquad (68) $$

where $\mathcal{D}w$ denotes the standard Gaussian measure. The following claim, from [12], gives the precise estimation error of GAMP. It is stated there as a claim because some steps of the proof are missing. The paper [45] affirms in its abstract to prove the claim of [12], but without further details. For these reasons, we believe that the claim holds, however we prefer to state it here as a claim (instead of a theorem).

**Claim 1 (GAMP estimation error, [12]).** We have almost surely for all $t \in \mathbb{N}$,

$$ \lim_{n \to \infty} \frac{1}{n} \hat{\mathbf{x}}^t \cdot \mathbf{X}^* = \lim_{n \to \infty} \frac{1}{n} \|\hat{\mathbf{x}}^t\|^2 = q^t, \qquad (69) $$

as well as

$$ \lim_{n \to \infty} \frac{1}{n^2} \mathbb{E} \left[ \| \mathbf{X}^* \mathbf{X}^{*\top} - \hat{\mathbf{x}}^t (\hat{\mathbf{x}}^t)^\top \|^2 \right] = \rho^2 - (q^t)^2. \qquad (70) $$

Comparing (70) with the MMSE given by Corollary 5, we see that if $\lim_{t \to \infty} q^t = q^*(\alpha)$, then GAMP achieves the MMSE. Provided that Claim 1 holds we can deduce the generalization error of GAMP:

**Proposition 2 (GAMP generalization error).** Suppose that hypotheses (h1)-(h2)-(h4) hold. Moreover suppose that either (h5.a) or (h5.b) holds. Assume that $(\mathbf{\Phi}_{\mu i}) \stackrel{\text{iid}}{\sim} \mathcal{N}(0, 1)$, and that $x \mapsto P_{\text{out}}(\cdot | x)$ is continuous almost everywhere for the Wasserstein distance of order 2. Let $t \in \mathbb{N}$. Assume that the limit (69) holds in probability and that there exists $\eta > 0$ such that $\mathbb{E}[|\hat{Y}^{\text{GAMP},t}|^{2+\eta}]$ remains bounded (as $n$ grows). Then we have

$$ \lim_{n \to \infty} \mathcal{E}_{\text{gen}}^{\text{GAMP},t} := \lim_{n \to \infty} \mathbb{E} \left[ (Y_{\text{new}} - \hat{Y}^{\text{GAMP},t})^2 \right] = \mathcal{E}(q^t). \qquad (71) $$

**Remark 2.** If we modify slightly the GAMP estimator of (68) by changing the first $y$ into $f(y)$, it is not difficult to show (following the steps of Proposition 2) that this new estimator achieves an asymptotic error of $\mathcal{E}_f(q^t)$, given by (61), for estimating $f(Y_{\text{new}})$.
Proposition 2 is proved in Sec. 5.2. We see that this formula matches the one for the Bayes-optimal generalization error, see Theorem 3, up to the fact that instead of $q^*(\alpha)$ (the optimizer of the replica formula (43)) appearing in the optimal error formula, here it is $q^t$ which appears. Thus clearly, when $q^t$ converges to $q^*(\alpha)$ (we shall see that this is the case in many situations in the examples of Sec. 3) this yields a very interesting and non trivial result: GAMP achieves the Bayes-optimal generalization error in a plethora of models (a task again often believed to be intractable) and this for large sets of parameters.

## 2.5 Optimal denoising error

Another interesting error measure to study is the following “denoising error”. Assume that the observations are noisy, i.e. $\Delta > 0$ in (20). The goal here is to denoise the observations $Y_\mu$ and estimate the signal which in this case is $\varphi(\frac{1}{\sqrt{n}}[\mathbf{\Phi}\mathbf{X}^*]_\mu, \mathbf{A}_\mu)$.

The minimum denoising error (in $L^2$ sense) is actually a simple corollary from the replica-symmetric formula of Theorem 1 and follows from a so-called “I-MMSE relation”, see Proposition 13. We will need the joint posterior distribution of $(W^*, \mathbf{A})$ given $(V, \tilde{Y}_0)$ for the scalar channel (39). So we define the Gibbs bracket for the scalar channel by (here $\mathbf{a} \in \mathbb{R}^{k_A}$):

$$ \langle g(w, \mathbf{a}) \rangle_{\text{sc}} := \mathbb{E}[g(W^*, \mathbf{A})|\tilde{Y}_0, V] = \frac{\int \mathcal{D}wdP_A(\mathbf{a})g(w, \mathbf{a})e^{-\frac{1}{2\Delta}\{\tilde{Y}_0 - \varphi(\sqrt{q}V + \sqrt{\rho-q}w, \mathbf{a})\}^2}}{\int \mathcal{D}wdP_A(\mathbf{a})e^{-\frac{1}{2\Delta}\{\tilde{Y}_0 - \varphi(\sqrt{q}V + \sqrt{\rho-q}w, \mathbf{a})\}^2}}, \quad (72) $$

for any continuous bounded function $g$. When the function depends only on $w$ it may be re-written as

$$ \langle g(w) \rangle_{\text{sc}} = \frac{\int \mathcal{D}wg(w)P_{\text{out}}(\tilde{Y}_0|\sqrt{q}V + \sqrt{\rho-q}w)}{\int \mathcal{D}wP_{\text{out}}(\tilde{Y}_0|\sqrt{q}V + \sqrt{\rho-q}w)}. \quad (73) $$

**Corollary 5 (Optimal denoising error).** Suppose that hypotheses (h1)-(h2)-(h3)-(h4) hold. Suppose that either hypothesis (h5.a) or (h5.b) holds. Then for almost every $\Delta > 0$, for any optimal couple $(q^*, r^*)$ of (43),

$$ \lim_{n \to \infty} \frac{1}{m} \text{MMSE} \left( \varphi \left( \frac{1}{\sqrt{n}} \mathbf{\Phi} \mathbf{X}^*, \mathbf{A} \right) \middle| \mathbf{\Phi}, \mathbf{Y} \right) = \text{MMSE} \left( \varphi(\sqrt{q^*} V + \sqrt{\rho-q^*} W^*, \mathbf{A}) \middle| \tilde{Y}_0, V \right) \\ = \mathbb{E} \left[ \varphi(\sqrt{\rho} V, \mathbf{A})^2 \right] - \mathbb{E} \left[ \langle \varphi(\sqrt{q^*} V + \sqrt{\rho-q^*} w, \mathbf{a}) \rangle_{\text{sc}}^2 \right], \quad (74) $$

where $\langle - \rangle_{\text{sc}}$ acts jointly on $(w, \mathbf{a})$ and is defined by (72), and $V, W^* \stackrel{\text{iid}}{=} \mathcal{N}(0, 1)$.

Note that the joint posterior over both the signal $\mathbf{X}^*$ and the random stream $\mathbf{A}$ is simply expressed as

$$ dP(\mathbf{x}, \mathbf{a} | \mathbf{Y}, \mathbf{\Phi}) \propto dP_0(\mathbf{x})dP_A(\mathbf{a}) \prod_{\mu=1}^{m} e^{-\frac{1}{2\Delta} \{Y_\mu - \varphi(\frac{1}{\sqrt{n}}[\mathbf{\Phi}\mathbf{x}]_\mu, \mathbf{a}_\mu)\}^2}. \quad (75) $$

The proof of Corollary 5 is presented in Sec. 5.4.

## 3 Application to concrete situations

In this section, we show how our main results can be applied to several models of interest in fields ranging from machine learning to signal processing, and unveil several interesting new phenomena in learning of generalized linear models. For various specific cases of prior $P_0$ and output $P_{\text{out}}$, we evaluate numerically the free entropy potential (41), its stationary points $\Gamma$ and identify which of them gives the information-theoretic results, i.e. is the optimizer in (43). We also identify which of the stationary points corresponds to the result
obtained asymptotically by the GAMP algorithm, i.e. the fixed point of the state evolution (67). Finally we compute the corresponding generalization error (71). We stress that in this section the results are based on numerical investigation of the resulting formulas: We do not aim at rigor that would involve precise bounds and more detailed analytical control for the corresponding integrals.

## 3.1 Generic observations

Using the functions $g_{P_{\text{out}}}$ and $g_{P_0}$ introduced in Sec. 2.4 we can rewrite the fixed point equations (42) as

$$q = 2\psi'_{P_0}(r) = \mathbb{E}[g_{P_0}(Y_0, r)^2], \qquad (76)$$

$$r = 2\alpha\Psi'_{P_{\text{out}}}(q) = \frac{\alpha}{\rho - q}\mathbb{E}[g_{P_{\text{out}}}(\tilde{Y}_0, \sqrt{q}V, \rho - q)^2], \qquad (77)$$

where the expectation in (76) corresponds to the scalar channel (36) and the expectation in (77) corresponds to the second scalar channel (39).

**Non-informative fixed point and its stability:** It is interesting to analyze under what conditions $q^* = 0$ is the optimizer of (43). Notice that $q^* = 0$ corresponds to the error on the recovery of the signal as large as it would be if we had no observations at our disposition. Theorem 1 gives that any optimal couple $(q^*, r^*)$ of (43) should be a fixed point of the state evolution equations (76)–(77). A sufficient condition for $(q, r) = (0, 0)$ to be a fixed point of (76)–(77) is that:

(a) The transition density $P_{\text{out}}(y|z)$ is even in the argument $z$.

(b) The prior $P_0$ has zero mean.

In order to see this, notice that if $P_{\text{out}}(y|z)$ is even in $z$ then from the definition (65) of the function $g_{P_{\text{out}}}$ we have $g_{P_{\text{out}}}(y, 0, \rho) = 0$ and consequently from (77) we have $\Psi'_{P_{\text{out}}}(0) = 0$. From the second point, notice that we have $\psi'_{P_0}(0) = \frac{1}{2}\mathbb{E}_{P_0}[X_0]^2 = 0$.

We assume now that the transition density $P_{\text{out}}(y|z)$ is even in the argument $z$ and that the prior $P_0$ has zero mean. In order for $q = 0$ to be the global maximizer $q^*$ of (43) or to be a relevant fixed point of the state evolution (67) (relevant in the sense that GAMP might indeed converge to it in a practical setting) we need $q = 0$ to be a stable fixed point of the above equations (76)–(77). We therefore need to expand (76)–(77) around $q = 0$, and doing so, we obtain that $q = 0$ is stable if

$$2\alpha\Psi''_{P_{\text{out}}}(0) \times 2\psi''_{P_0}(0) = \alpha\mathbb{E}\left[(\langle w^2 \rangle_{\text{sc}} - \langle w \rangle_{\text{sc}}^2 - 1)^2\right] < 1, \qquad (78)$$

where the expectation corresponds to the scalar channel (39) with $q = 0$ and the Gibbs bracket $\langle - \rangle_{\text{sc}}$ is given by (72). The expectation quantifies how the observation of $\tilde{Y}_0$ in the scalar channel (39) modifies the variance of $W^*$ (which is 1 without any observation). Rewriting this condition more explicitly into a form that is convenient for numerical evaluation we get (recalling (73), $q = 0$ and condition (a))

$$\alpha \int dy \frac{\left( \int \mathcal{D}z (z^2 - 1) P_{\text{out}}(y | \sqrt{\rho}z) \right)^2}{\int \mathcal{D}z P_{\text{out}}(y | \sqrt{\rho}z)} < 1, \qquad (79)$$

where recall that $\mathcal{D}z$ is a standard Gaussian measure. We conjecture that the condition (79) delimits precisely the region where polynomial-time algorithms do not perform better than “random guessing” (see the discussion below, where we will make this stability condition explicit for several examples of symmetric output channels). Note that the condition (79) also appears in a recent work [61] as a barrier for performance of spectral algorithms.
**Exact recovery fixed point:** Another particular fixed point of (76)–(77) that we observe is the one corresponding to exact recovery $q^* = \rho$. A sufficient and necessary condition for this to be a fixed point is that $\lim_{q \to \rho} \Psi'_{\text{out}}(q) = +\infty$. Heuristically, this means that the integral of the Fisher information of the output channel should diverge:

$$ \int dyd\omega \frac{e^{-\frac{\omega^2}{2\rho}} P'_{\text{out}}(y|\omega)^2}{\sqrt{2\pi\rho} P_{\text{out}}(y|\omega)} = +\infty, \qquad (80) $$

where $P'_{\text{out}}(y|\omega)$ denotes the partial derivative w.r.t. $\omega$. This typically means that the channel should be noiseless. For the Gaussian channel with noise variance $\Delta$, the above expression equals $1/\Delta$. For the probit channel where $P_{\text{out}}(y|z) = \text{erfc}(-yz/\sqrt{2\Delta})/2$ the above expression at small $\Delta$ is proportional to $1/\sqrt{\Delta}$.

Stability of the exact recovery fixed point was also investigated, but we did not obtain any unified expression. The stability depends non-trivially on both the properties of the output channel, but also on the properties of the prior. Below we give several examples where exact recovery either is or is not possible, or where there is a phase transition between the two regimes.

## 3.2 Phase diagram of perfect learning

In this section we consider deterministic (noiseless) output channels and ask: How many measurements are needed in order to perfectly recover the signal?

Our crucial point is to compare with the well explored phase diagram of Bayesian (noiseless) compressed sensing in the case of the linear channel [68, 80]. As the number of samples (measurements) varies we encounter five different regimes of parameters:

* The tractable recovery phase: This is the region in the parameter space where GAMP achieves perfect reconstruction.

* The non-informative phase: Region where perfect reconstruction is information-theoretically impossible and moreover even the Bayes-optimal estimator is as bad as a random guess based on the prior information and on the knowledge of the output function.

* The no recovery phase: Region where perfect reconstruction is information-theoretically impossible, but an estimator positively correlated with the ground truth exists.

* The hard phase: Region where the perfect reconstruction is information-theoretically possible, but where GAMP is unable to achieve it. At the same time, in this region GAMP leads to a better generalization error than the one corresponding to the non-informative fixed point. It remains a challenging open question whether polynomial-time algorithms can achieve perfect reconstruction in this regime.

* The hard non-informative phase: This phase corresponds to the region where perfect reconstruction is information-theoretically possible but where GAMP only achieves an error as bad as randomly guessing, given by the trivial fixed point. In this phase as well, the existence of polynomial-time exact recovery algorithms is an open question. This phase does not exist for the linear channel.

Some of the codes used in this section can be consulted online on the github repository [62].

### 3.2.1 The linear channel

The case of exact recovery of a sparse signal after it passed through a noiseless linear channel, i.e. $\phi(x) = x$, is studied in the literature in great details, especially in the context of compressed sensing [10, 11]. For a signal with a fraction $\rho$ of non-zero entries it is found that as soon as $\alpha > \rho$, perfect reconstruction is theoretically possible, although it may remain computationally difficult. The whole field of compressed sensing builds on
<figure><img src="image_5.png" /><figcaption>Figure 3: Phase diagrams showing boundaries of the region where exact recovery is possible (in absence of noise). **Left:** The case of sign-less sparse recovery, $\varphi(x) = |x|$ with a Gauss-Bernoulli signal, as a function of the ratio between number of samples/measurements and the dimension $\alpha = m/n$, and the fraction of non-zero components $\rho$. Evaluating the free entropy for this case, we find that a recovery of the signal is information-theoretically impossible for $\alpha < \alpha_{\text{IT}} = \rho$. Recovery becomes possible starting from $\alpha > \rho$, just as in the canonical compressed sensing. Algorithmically the sign-less case is much harder. Evaluating (79) we conclude that GAMP is not able to perform better than a random guess as long as $\alpha < \alpha_c = 1/2$. For larger values of $\alpha$, the inference using GAMP leads to better results than a purely random guess. GAMP can exactly recover the signal and generalize perfectly only for values of $\alpha$ larger than $\alpha_{\text{AMP}}$ (full red line). The dotted red line shows for comparison the algorithmic phase transition of the canonical compressed sensing. **Center:** Analogous to the left panel, for the ReLU output function, $\varphi(x) = \max(0, x)$. Here it is always possible to perform better than random guessing using GAMP. The dotted red line shows the algorithmic phase transition when using information only about the non-zero observations. **Right:** Phase diagram for the symmetric door output function $\varphi(z) = \operatorname{sign}(|z| - K)$ for a Rademacher signal, as a function of $\alpha$ and $K$. The stability line $\alpha_c$ is depicted in dashed blue, the information-theoretic phase transition to exact recovery $\alpha_{\text{IT}}$ in black, and the algorithmic one $\alpha_{\text{AMP}}$ in red.</figcaption></figure>

the realization that, using the $\ell_1$ norm minimization technique, one can efficiently recover the signal for larger $\alpha$ values than the so-called Donoho-Tanner transition [11, 81].

In the context of the present paper, when the empirical distribution of the signal is known, one can fairly easily beat the $\ell_1$ transition and reconstruct the signal up to lower values of $\alpha$ using the Bayesian GAMP algorithm [11, 12, 68, 80]. In this case, three different phases are present [68, 80]: i) For $\alpha < \rho$, perfect reconstruction is impossible; ii) for $\rho < \alpha < \alpha_s$ reconstruction is possible, but not with any known polynomial-complexity algorithm; iii) for $\alpha > \alpha_s$, the so-called *spinodal transition* computed with state evolution, GAMP provides a polynomial-complexity algorithm able to reach perfect reconstruction. The line $\alpha_s(\rho)$ depends on the distribution of the signal. For a Gauss-Bernoulli signal with a fraction $\rho$ of non-zero (Gaussian) values we compare the GAMP performance to the optimal one in Fig. 3 (left and right). This is the same figure as in the main text. We copy it here so that the SI is self-contained.

### 3.2.2 The rectified linear unit (ReLU) channel

Let us start by discussing the case of a generalized linear model with the ReLU output channel, i.e. $\varphi(x) = \max(0, x)$, with a signal coming from a Gauss-Bernoulli distribution $P_0 = \rho N(0, 1) + (1 - \rho)\delta_0$ with a fraction $\rho$ of non-zero (Gaussian) values. We are motivated by the omnipresent use of the ReLU activation function in deep learning, and explore its properties for GLMs that can be seen as a simple single layer neural network.

Our analysis shows that a perfect generalization (and thus a perfect reconstruction of the signal as well) is possible whenever the number of samples per dimension (measurement rate) $\alpha > 2\rho$, and impossible when $\alpha < 2\rho$. This is very intuitive, since half of the measurements (those non-zero) are giving as much information as in the linear case, thus the factor 2.
How hard is it to actually solve the problem with an efficient algorithm? The answer is given by applying the state evolution analysis to GAMP, which tells us that only for even larger values of $\alpha$, beyond the spinodal transition, does GAMP reach a perfect recovery. Notice, however, that this spinodal transition occurs at a significantly lower measurement rate $\alpha$ than one would reach just keeping the non-zero measurements. This shows that, actually, these zero measurements contain a useful information for the algorithm. The situation is shown in the center panel of Fig. 3: The zero measurements do not help information-theoretically but they, however, do help algorithmically.

### 3.2.3 The sign-less channel

We now discuss the sign-less channel where only the absolute value of the linear mixture is observed, i.e. $\varphi(x) = |x|$. This case can be seen as the real-valued analog of the famous phase retrieval problem. We again consider the signal to come from a Gauss-Bernoulli distribution with a fraction $\rho$ of non-zero (Gaussian) values.

Sparse phase retrieval has been well explored in the literature in the regime where the number $s$ of non-zeros is sub-leading in the dimension, $s = O(n)$. This case is known to present a large algorithmic gap. While analogously to compressed sensing exact recovery is information-theoretically possible for a number of measurement $\Omega(s \ln(n/s))$, best known algorithms achieve it only with $\Omega(s^2/\ln n)$ measurements [82], see also [70] and references therein for a good discussion of other related literature. This is sometimes referred to as the $s^2$ barrier. We are not aware of a study where, as in our setting, the sparsity is $s = \rho n$ and the number of measurements is $\alpha n$ with $\alpha = \Omega(1)$. Our analysis in this regime hence sheds a new light on the hardness of the problem of recovering a sparse signal from sign-less measurements.

Our analysis of the mutual information shows that a perfect reconstruction is information-theoretically possible as soon as $\alpha > \rho$: In other words, the problem is –information-theoretically– as easy, or as hard as the compressed sensing one. This is maybe less surprising when one thinks of the following algorithm: Try all $2^m$ choices of the possible signs for the $m$ outputs, and solve a compressed sensing problem for each of them. Clearly, this should yields a perfect solution only in the case of the actual combination of signs.

Algorithmically, however, the problem is much harder than for the linear output channel. As shown in the left side of Fig. 3, for small $\rho$ one requires a much larger fraction $\alpha$ of measurements in order for GAMP to recover the signal. For the linear channel the algorithmic transition $\alpha_s(\rho) \to 0$ as $\rho \to 0$, while for the sign-less channel we get $\alpha_s(\rho) \to 1/2$ as $\rho \to 0$. In other words if one looses the signs one cannot perform recovery in compressed sensing with less than $n/2$ measurements.

What we observe in this example for $\alpha < 1/2$ is in the statistical physics literature on neural networks known as retarded learning [64]. This appears in problems where the $\varphi(x)$ function is symmetric, as seen at the beginning of this section: There is always a critical point of the mutual information with an overlap value $q = 0$. For this problem, this critical point is actually “stable” (meaning that it is actually a local minimum in $q$ in the mutual information (46)) for all $\alpha < 1/2$ independently of $\rho$. To see that, we have to go back to (78). For the absolute value channel, the posterior distribution of $W^*$ given $\tilde{Y}_0 = \sqrt{\rho}|W^*|$ in the second scalar channel (39) is a mean of two Dirac masses on $-W^*$ and $W^*$. Thus, the posterior variance is $\langle w^2 \rangle_{sc} - \langle w \rangle_{sc}^2 = (W^*)^2$. Consequently, (78) leads to the stability condition $\alpha < 1/2$.

This has the two following implications: i) In the non-informative phase, when $\alpha < 1/2$ and $\rho > \alpha$, the minimum at $q = 0$ is actually the global one. In this case, the MMSE on $X^*$ and the generalization error are the ones given by using 0 as a guess for each element of $X^*$; in other words, there is no useful information that one can exploit and no algorithmic approach can be better than a random guess. ii) In the hard non-informative phase when $\alpha < 1/2$, GAMP initialized at random, i.e. close to the $q = 0$ fixed point, will remain there. This suggests that in this region, even if a perfect reconstruction is information-theoretically possible, it will still be very hard to beat a random guess with a tractable algorithm.
<figure><img src="image_6.png" /><figcaption>Figure 4: Generalization error in three classification problems as a function of the number of data-samples per dimension $\alpha$. The red line is the Bayes-optimal generalization error, while the green one shows the (asymptotic) performances of GAMP as predicted by the state evolution (SE), when different. For comparison, we also show the result of GAMP (black dots) and, in blue, the performance of a standard out-of-the-box solver, both tested on a single randomly generated instance. **Left:** Perceptron, with $\varphi(x) = \sgn(x)$ and a Rademacher ($\pm 1$) signal. While a perfect generalization is information-theoretically possible starting from $\alpha = 1.249(1)$, the state evolution predicts that GAMP will allow such perfect prediction only from $\alpha = 1.493(1)$. The results of a logistic regression with fine-tuned ridge penalty with the software scikit-learn [65] are shown for comparison. **Middle:** Perceptron with Gauss-Bernoulli coefficients for the signal. No phase transition is observed in this case, but a smooth decrease of the error with $\alpha$. The results of a logistic regression with fine-tuned $\ell_1$ sparsity-enhancing penalty (again with [65]) are very close to optimal. **Right:** The symmetric door activation rule with parameter $K = 0.67449$ chosen in order to observe the same number of occurrence of the two classes. In this case there is a sharp phase transition at $\alpha = 1$ from a situation where it is impossible to learn the rule, so that the generalization is not better than a random guess, to a situation where the optimal generalization error drops to zero. However, GAMP identifies the rule perfectly only starting from $\alpha_s = 1.566(1)$ (GAMP error stays 1 up to $\alpha_{stab} = 1.36$, see the black dashed curve). Interestingly, this non-linear rule seems very hard to learn for other existing algorithms. Using Keras [66], a neural network with two hidden layers was able to learn approximately the rule, but only for much larger training set sizes (shown in inset, the Keras/tensorflow code for this particular run can be found on the github repository [62]).</figcaption></figure>

### 3.2.4 The symmetric door channel

The third output channel we study in detail is the symmetric door channel, where $\varphi(x) = \sgn(|x| - K)$. In case of channels with discrete set of outputs exact recovery is only possible when the prior is also discrete. In the present case we consider the signal to be Rademacher, where each element is chosen at random between 1 and $-1$, i.e. $P_0 = \frac{1}{2}\delta_{+1} + \frac{1}{2}\delta_{-1}$. This channel was studied previously using the replica method in the context of optimal data compression [28].

This output channel is in the class of symmetric channels for which overlap $q = 0$ is a fixed point. This fixed point is stable for $\alpha < \alpha_c(K)$. Exact recovery is information-theoretically possible above $\alpha_{IT}(K)$ and tractable with the GAMP algorithm above the spinodal transition $\alpha_s(K)$. The values of these three transition lines are depicted in the right panel of Fig. 3.

We note that $\alpha_{IT} \ge 1$ is a generic bound on exact recovery for every $K$, required by a simple counting argument. While a-priori it is not clear whether this bound is saturated for some $K$, we observe that it is for $K = 0.67449$ such that half of the observed measurements are negative and the rest positive. This is, however, not efficiently achievable with GAMP. The saturation of the $\alpha_{IT} \ge 1$ bound was remarked previously in the context of the work [28] on optimal data compression. Our work predicts that this information-theoretic result will not be achievable with known efficient algorithms.

## 3.3 Examples of optimal generalization error

Besides the formula for the mutual information, the main result of this paper is the Theorem 3 for the optimal generalization error, and formula (71) for the generalization error achieved by the GAMP algorithm. In this section we evaluate both these generalization errors for several cases of priors and output functions. We study both regression problems, where the output is real-valued, and classification problems, where the output is
discrete.

While in realistic regression and classification problems the matrix $\Phi$ corresponds to the data, and is thus not i.i.d. random, we view the practical interest of our theory as a benchmark for state-of-the-art algorithms. Our work provides an exact asymptotic analysis of optimal generalization error and sample complexity for a range of simple rules where a teacher uses random data to generate labels. The challenge for state-of-the-art multi-purpose algorithms is to try to match as closely as possible the performance that can be obtained with GAMP that is fine-tuned to the specific form of the output and prior.

### 3.3.1 Threshold output: The perceptron

The example of non-linear output that is the most widely explored in the literature is the threshold output, where the deterministic output (or “activation”) function is $\varphi(x) = \text{sgn}(x)$. This output in the teacher-student setting of the present paper is known as the perceptron problem [22], or equivalently, the one-bit compressed sensing in signal processing [6]. Its solution has been discussed in details within the replica formalism (see for instance [15, 55, 57, 83]) and we confirm all of these heuristic computations within our approach. Let $V \sim \mathcal{N}(0, 1)$. The formula (57), (63) for the generalization error then reduces to (recall $q^* = q^*(\alpha)$ is an optimizer of (43))

$$ \lim_{n \to \infty} \mathcal{E}_{\text{gen}}^{\text{opt}} = 1 - \int \mathcal{D}V \left( \frac{2}{\sqrt{\pi}} \int_{0}^{V} \sqrt[2]{\frac{q^*}{2(\rho - q^*)}} dt e^{-t^2/2} \right)^2 = 1 - \mathbb{E} \left[ \text{erf} \left( V \sqrt{\frac{q^*}{2(\rho - q^*)}} \right)^2 \right]. \quad (81) $$

In Fig. 4 (left) we plot the optimal generalization error of the perceptron with a Rademacher signal, the state evolution prediction of the generalization error of the GAMP algorithm, together with the error actually achieved by GAMP on one randomly generated instance of the problem. We also compare these to the performance of a standard logistic regression. As expected from existing literature [22, 23] we confirm that in this case the information-theoretic transition appears at a number of samples per dimension $\alpha_{IT} = 1.249(1)$, while the algorithmic transition is at $\alpha_s = 1.493(1)$. Logistic regression does not seem to be able to match the performance on GAMP in this case.

In Fig. 4 (center) we plot the generalization error for a Gauss-Bernoulli signal with density $\rho = 0.2$. Cases as this one were studied in detail in the context of one-bit compressed sensing [15] and GAMP was found to match the optimal generalization performance with no phase transitions observed, which is confirmed by our analysis. In this case the logistic regression is rather close to the performance of GAMP.

### 3.3.2 Symmetric Door

The next classification problem, i.e. discrete output rule, we study is the symmetric door function $\varphi(x) = \text{sgn}(|x| - K)$. In this case the generalization error (57) becomes (here again $V \sim \mathcal{N}(0, 1)$)

$$ \lim_{n \to \infty} \mathcal{E}_{\text{gen}}^{\text{opt}} = 1 - \mathbb{E}_V \left[ \left\{ \text{erf}\left(\frac{K - \sqrt{q^*}V}{\sqrt{2(\rho - q^*)}}\right) - \text{erf}\left(-\frac{K + \sqrt{q^*}V}{\sqrt{2(\rho - q^*)}}\right) - 1 \right\}^2 \right]. \quad (82) $$

In Fig. 4 (right) we plot the generalization error for $K = 0, 67449$ such that $1/2$ of the outputs are 1 and $1/2$ are $-1$. The symmetric door output is an example of function for which the optimal generalization error for $\alpha < \alpha_{IT} = 1$ (for that specific value of $K$, see phase diagram in the right panel of Fig. 3) is as bad as if we were guessing randomly. The GAMP algorithm still achieves such a bad generalization until $\alpha_{stab} = 1.36$, and achieves perfect generalization only for $\alpha > \alpha_s = 1.566(1)$.

Interestingly, labels created from this very simple symmetric door rule seem to be very challenging to learn for general purpose algorithms. We tried to optimize parameters of a two-layers neural network and
<figure><img src="image_7.png" /><figcaption>Figure 5: The generalization error for three regression problems is plotted as a function of the number of samples per dimension $\alpha$. The red line is again the Bayes-optimal generalization error, while the green one shows the (asymptotic) performances of GAMP as predicted by the state evolution (SE), when different. Again, we also show the result of GAMP on a particular instance (black dots) and, in blue, the performance of an out-of-the-box solver. **Left:** White Gaussian noise output and a Gauss-Bernoulli signal. For this choice of noise, there is no sharp transition (as opposed to what happens at smaller noises). The results of a LASSO with fine-tuned $\ell_1$ sparsity-enhancing penalty (with [65]) are very close to optimal. **Middle:** Here we analyze a ReLU output function $\varphi(x) = \max(0, x)$, still with a Gauss-Bernoulli signal. Now there is an information-theoretic phase transition at $\alpha = 2\rho = 0.4$, but GAMP requires $\alpha_s = 0.589(1)$ to reach perfect recovery. We show for comparison the results of maximum likelihood estimation performed with CVXPY—a powerful python-embedded language for convex optimization [84]— using two methods that are both amenable to convex optimization: In CVX-1 we use only the non-zero values of $\mathbf{Y}$, and perform a minimization of the $\ell_1$ norm of $\mathbf{x}$ subject to $Y_\mu = \mathbf{\Phi}_\mu \cdot \mathbf{x}$ for $\mu \in \{1, \dots, m\}$ such that $Y_\mu \neq 0$, while in CVX-2, we use all the dataset, with the constraint that $Y_\mu = \mathbf{\Phi}_\mu \cdot \mathbf{x}$ for $\mu \in \{1, \dots, m\}$ such that $Y_\mu \neq 0$ (as before) and the additional restriction $\mathbf{\Phi}_\mu \cdot \mathbf{x} \le 0$ for $\mu \in \{1, \dots, m\}$ such that $Y_\mu = 0$. In both case, a perfect generalization is obtained only for $\alpha \gtrsim 1$. **Right:** The sign-less output function $\varphi(x) = |x|$. The information-theoretic perfect recovery starts at $\alpha = \rho = 0.5$, but the problem is again harder algorithmically for GAMP that succeeds only above $\alpha_s = 0.90(1)$. Again, the problem appears to be hard for other solvers. In inset, we show the performance for the estimation problem using PhaseMax [85], which is able to learn the rule only using about four times as many measurements than needed information-theoretically.</figcaption></figure>

only managed to get the performances shown in the inset of Fig. 4 (right). It is an interesting theoretical challenge whether a deeper neural network can learn this simple rule from fewer samples.

### 3.3.3 Linear regression

The additive white Gaussian noise (AWGN) channel, or linear regression, is defined by $\varphi(x, A) = x + \sigma A$ with $A \sim \mathcal{N}(0, 1)$. This models the (noisy) linear regression problem, as well as noisy random linear estimation and compressed sensing. In this case (57) leads to

$$ \lim_{n \to \infty} \mathcal{E}_{\text{gen}}^{\text{opt}} = \rho - q^* + \sigma^2. \qquad (83) $$

This result agrees with the generalization error analyzed heuristically in [23] in the limit $\sigma \to 0$. Fig. 5 (left) depicts the generalization error for this example. The performance of GAMP in this case is very close to the one of LASSO.

### 3.3.4 Rectified linear unit (ReLU)

In Fig. 5 (center) we analyze the generalization error for the ReLU output function, $\varphi(x) = \max(0, x)$. This channel models the behavior of a single neuron with the rectified linear unit activation [8] widely used in multilayer neural networks. In this case (57) becomes after simple algebra and Gaussian integration by parts
(again $V \sim \mathcal{N}(0, 1))$,

$$ \lim_{n \to \infty} \mathcal{E}_{\text{gen}}^{\text{opt}} = \frac{\rho}{2} - \frac{q^*}{4} \left( 1 + \mathbb{E}_V \left[ V^2 \text{erf} \left( V \sqrt{\frac{q^*}{2(\rho - q^*)}} \right)^2 \right] \right) - \frac{(\rho - q^*)^{3/2}}{\sqrt{\rho + q^*}} \left( \frac{1}{2\pi} + \frac{q^*}{\rho\pi} \sqrt{\frac{\rho + q^*}{\rho - q^*}} \right). \quad (84) $$

For sparse Gauss-Bernoulli signals in Fig. 5 (center) we observe again the information-theoretic transition to perfect generalization to be distinct from the algorithmic one. At the same time our test with existing algorithms were not able to closely match the performance of GAMP. This hence also remains an interesting benchmark.

### 3.3.5 Sign-less channel

In Fig. 5 (right) we analyze the generalization error for the sign-less output function where $\varphi(x) = |x|$. This models a situation similar to compressed sensing, except that the sign of the output has been lost. This is a real-valued analog of the phase retrieval problem as discussed in Sec. 3.2.3. In this case the generalization error (57) becomes (again $V \sim \mathcal{N}(0, 1)$)

$$ \lim_{n \to \infty} \mathcal{E}_{\text{gen}}^{\text{opt}} = \rho - \mathbb{E}_V \left[ b\left(V\sqrt{q^*}, \rho - q^*\right)^2 \right], \quad (85) $$

where

$$ b(x, y) = \sqrt{\frac{2y}{\pi}} e^{-\frac{x^2}{2y}} + \frac{x}{2} \text{erfc} \left( -\frac{x}{\sqrt{2y}} \right) - \frac{x}{2} \left\{ 1 + \text{erf} \left( -\frac{x}{\sqrt{2y}} \right) \right\}. \quad (86) $$

Our comparison with the performance of a state-of-the-art algorithm PhaseMax [85] suggests that also for this simple benchmark there is room for improvement in term of matching the performance of GAMP.

### 3.3.6 Sigmoid, or logistic regression

Let us also consider an output function with auxiliary randomization. After having generated the classifier $\mathbf{X}^*$, the teacher randomly associates the label +1 to the pattern $\mathbf{\Phi}_\mu$ with probability $f_\lambda(n^{-1/2}\mathbf{\Phi}_\mu \cdot \mathbf{X}^*)$, where $f_\lambda(x) = (1 + \exp(-\lambda x))^{-1} \in [0, 1]$ is the sigmoid of parameter $\lambda > 0$, and the label -1 with probability $1 - f_\lambda(n^{-1/2}\mathbf{\Phi}_\mu \cdot \mathbf{X}^*)$. One of the (many) possible ways for the teacher to do so is by selecting $\varphi(x, A) = \mathbf{1}(A \le f_\lambda(x)) - \mathbf{1}(A > f_\lambda(x))$, where $\mathbf{1}(E)$ is the indicator function of the event $E$. He then generates a stream of uniform random numbers $A \stackrel{\text{iid}}{\sim} \mathcal{U}_{[0,1]}$ and obtains the labels through (20) (with $\Delta = 0$). Let $V, w \stackrel{\text{iid}}{\sim} \mathcal{N}(0, 1)$. In this setting the error (57) becomes

$$ \lim_{n \to \infty} \mathcal{E}_{\text{gen}}^{\text{opt}} = 2 - 4 \mathbb{E}_V \left[ \left\{ \mathbb{E}_w f_\lambda \left( \sqrt{q^*} V + \sqrt{\rho - q^*} w \right) \right\}^2 \right]. \quad (87) $$

This formula reduces to (81) when $\lambda \to \infty$ as it should.

## 4 Proof of the replica formula by the adaptive interpolation method

We now prove Theorem 1. Our main tool will be an interpolation method recently introduced in [52] and called “adaptive interpolation method”. Here we formulate the method as a direct evolution of the Guerra and Toninelli interpolation method developed in the context of spin glasses [72]. In contrast with the discrete and more pedestrian version of the adaptive interpolation method presented in [52], here we employ a continuous
approach which is more straightforward (see [52] for the links between the discrete and continuous versions of the method) and that has also been recently used in [86] for studying non-symmetric tensor estimation.

We will prove Theorem 1 under the following hypotheses:

(H1) The support of the prior distribution $P_0$ is included in $[-S, S]$, for some $S > 0$.

(H2) $\varphi$ is a bounded $C^2$ function with bounded first and second derivatives w.r.t. its first argument.

(H3) $(\Phi_{\mu i})_{\text{id}} \sim \mathcal{N}(0, 1)$.

These stronger assumptions will then be relaxed in Appendix C to the weaker assumptions (h1)-(h2)-(h3)-(h4) and (h5.a) or (h5.b). Since the observations (20) are equivalent to the rescaled observations

$$ \tilde{Y}_{\mu} := \Delta^{-1/2} Y_{\mu} = \Delta^{-1/2} \varphi \left( \frac{1}{\sqrt{n}} [\mathbf{\Phi} \mathbf{X}^*]_{\mu}, \mathbf{A}_{\mu} \right) + Z_{\mu}, \quad 1 \le \mu \le m, \qquad (88) $$

the variance $\Delta$ of the Gaussian noise can be “incorporated” inside the function $\varphi$. Thus, it suffices to prove Theorem 1 for $\Delta = 1$ and we suppose, for the rest of the proof, that we are in this equivalent case.

## 4.1 Interpolating estimation problem

We introduce an “interpolating estimation problem” that interpolates between the original problem (21) at $t = 0, t \in [0, 1]$ being the interpolation parameter, and the two scalar problems described in Sec. 1.3 at $t = 1$ which are analytically tractable. For $t \in (0, 1)$ the interpolating estimation problem is a mixture of the original and scalar problems. This interpolation scheme is inspired from the interpolation paths used by Talagrand to study the perceptron, see [73]. There are two major differences between the “non-planted perceptron” studied by Talagrand, and the “planted perceptron” that we are investigating:

* In the planted case, the presence of a planted solution forces (under small perturbations) the correlations to vanish for all values of the parameters, see [87, 88]. In the non-planted case, proving such decorrelation is much more involved, and is proved only in a limited region of the parameter space, see [73].

* However, in the planted case, there can be arbitrarily many solutions to the state evolution equations (42) (see Remark 21 in [89]), whereas in the region studied by [73], there is only one solution. For this reason, our interpolation method needs to be more sophisticated in order to interpolate with the “right fixed point”.

We fix a sequence $(s_n)_{n \ge 1} \in (0, 1/2]^N$ that converges to 0 as $n$ goes to infinity ($s_n$ will be chosen in Sec. 4.3 below to be equal to $\frac{1}{2}n^{-1/16}$). We define $\mathcal{B}_n := [s_n, 2s_n]^2$. For all $\epsilon = (\epsilon_1, \epsilon_2) \in \mathcal{B}_n$, we consider two continuous “interpolation functions” $q_{\epsilon} : [0, 1] \to [0, \rho]$ and $r_{\epsilon} : [0, 1] \to [0, r_{\text{max}}]$, where $r_{\text{max}} := 2\alpha \sup_{q \in [0, \rho]} \Psi'_{\text{Pout}}(q; \rho) = 2\alpha \Psi'_{\text{Pout}}(\rho; \rho)$ (recall that by Proposition 18, $\Psi'_{\text{Pout}}$ is non-decreasing). We define also for all $t \in [0, 1]$ and all $\epsilon \in \mathcal{B}_n$

$$ R_1(t, \epsilon) := \epsilon_1 + \int_0^t r_{\epsilon}(v)dv, \quad R_2(t, \epsilon) := \epsilon_2 + \int_0^t q_{\epsilon}(v)dv. \qquad (89) $$

We will be mainly interested in functions $r_{\epsilon}, q_{\epsilon}$ that satisfy some regularity properties. We will use the following definition:

**Definition 1 (Regularity).** We say that the families of functions $(q_{\epsilon})_{\epsilon \in \mathcal{B}_n}$ and $(r_{\epsilon})_{\epsilon \in \mathcal{B}_n}$, taking values respectively in $[0, \rho]$ and $[0, r_{\text{max}}]$, are regular if for all $t \in [0, 1]$ the mapping

$$ R^t: \begin{array}{rcl} (s_n, 2s_n)^2 & \rightarrow & R^t((s_n, 2s_n)^2) \\ \epsilon & \mapsto & (R_1(t, \epsilon), R_2(t, \epsilon)) \end{array} \qquad (90) $$
is a $C^1$ diffeomorphism, whose Jacobian is greater or equal to 1.

Define

$$S_{t,\mu} := \sqrt{\frac{1-t}{n}} [\mathbf{\Phi}\mathbf{X}^*]_{\mu} + \sqrt{R_2(t, \epsilon)} V_{\mu} + \sqrt{\rho t - R_2(t, \epsilon) + 2s_n} W_{\mu}^* \quad (91)$$

where $V_{\mu}, W_{\mu}^*$ are iid $N(0, 1)$. Consider the following observation channels, with two types of observations obtained through

$$\left\{ \begin{array}{ll} Y_{t,\mu} & \sim P_{\text{out}}(\cdot | S_{t,\mu}), \quad 1 \le \mu \le m, \\ Y'_{t,i} & = \sqrt{R_1(t, \epsilon)} X_i^* + Z'_i, \quad 1 \le i \le n, \end{array} \right. \quad (92)$$

where $(Z'_i)_{i=1}^n \stackrel{\text{iid}}{\sim} N(0, 1)$. We assume that $V = (V_{\mu})_{\mu=1}^m$ is known. Then the inference problem is to recover both unknowns $W^* = (W_{\mu}^*)_{\mu=1}^m$ and $X^* = (X_i^*)_{i=1}^n$ from the knowledge of $V, \Phi$ and the “time-dependent” observations $Y_t = (Y_{t,\mu})_{\mu=1}^m$ and $Y'_t = (Y'_{t,i})_{i=1}^n$.

We now understand that $R_1(t, \epsilon)$ appearing in the second set of measurements in (92), and the terms $1-t, R_2(t, \epsilon)$ and $\rho t - R_2(t, \epsilon) + 2s_n$ appearing in the first set all play the role of signal-to-noise ratios in the interpolating model, with $t$ giving more and more “power” (or weight) to the scalar inference channels when increasing. Here is the first crucial and novel ingredient of our interpolation scheme. In the classical interpolation method, these signal intensities would all take a trivial form (i.e. would be linear in $t$) but here, the non-trivial (integral) dependency in $t$ of the intensities through the use of the interpolation functions $q$ and $r$ allows for much more flexibility when choosing the interpolation path. This will allow us to actually choose the “optimal interpolation path” (this will become clear soon).

Define $u_y(x) := \ln P_{\text{out}}(y|x)$ and, with a slight abuse of notations,

$$s_{t,\mu} = s_{t,\mu}(\mathbf{x}, w_{\mu}) := \sqrt{\frac{1-t}{n}} [\mathbf{\Phi}\mathbf{x}]_{\mu} + \sqrt{R_2(t, \epsilon)} V_{\mu} + \sqrt{\rho t - R_2(t, \epsilon) + 2s_n} w_{\mu}. \quad (93)$$

We introduce the interpolating Hamiltonian

$$\mathcal{H}_{t,\epsilon}(\mathbf{x}, \mathbf{w}; \mathbf{Y}_t, \mathbf{Y}'_t, \mathbf{\Phi}, \mathbf{V}) := -\sum_{\mu=1}^{m} \ln P_{\text{out}}(Y_{t,\mu}|s_{t,\mu}) + \frac{1}{2} \sum_{i=1}^{n} (Y'_{t,i} - \sqrt{R_1(t, \epsilon)} x_i)^2. \quad (94)$$

The dependence in $\mathbf{\Phi}$ and $\mathbf{V}$ of the Hamiltonian is through the $(s_{t,\mu})_{\mu=1}^m$. It becomes, when the observations are replaced by their expression (92),

$$\mathcal{H}_{t,\epsilon}(\mathbf{x}, \mathbf{w}; \mathbf{Y}_t, \mathbf{Y}'_t, \mathbf{\Phi}, \mathbf{V}) = -\sum_{\mu=1}^{m} u_{Y_{t,\mu}}(s_{t,\mu}) + \frac{1}{2} \sum_{i=1}^{n} (\sqrt{R_1(t, \epsilon)}(X_i^* - x_i) + Z'_i)^2. \quad (95)$$

We also introduce the corresponding Gibbs bracket $\langle - \rangle_{n,t,\epsilon}$ which is the expectation operator w.r.t. the $(t, \epsilon)$-dependent posterior distribution of $(\mathbf{X}^*, \mathbf{W}^*)$ given $(\mathbf{Y}_t, \mathbf{Y}'_t, \mathbf{\Phi}, \mathbf{V})$. It is defined as

$$\langle g(\mathbf{x}, \mathbf{w}) \rangle_{n,t,\epsilon} := \frac{1}{\mathcal{Z}_{t,\epsilon}(\mathbf{Y}_t, \mathbf{Y}'_t, \mathbf{\Phi}, \mathbf{V})} \int dP_0(\mathbf{x}) \mathcal{D}\mathbf{w} g(\mathbf{x}, \mathbf{w}) e^{-\mathcal{H}_{t,\epsilon}(\mathbf{x}, \mathbf{w}; \mathbf{Y}_t, \mathbf{Y}'_t, \mathbf{\Phi}, \mathbf{V})}, \quad (96)$$

for every continuous bounded function $g$ on $\mathbb{R}^n \times \mathbb{R}^m$. In (96) $\mathcal{D}\mathbf{w} = (2\pi)^{-m/2} \prod_{\mu=1}^m dw_{\mu}e^{-w_{\mu}^2/2}$ is the $m$-dimensional standard Gaussian distribution and $\mathcal{Z}_{t,\epsilon}(\mathbf{Y}_t, \mathbf{Y}'_t, \mathbf{\Phi}, \mathbf{V})$ is the appropriate normalization (or
partition function):

$$Z_{t,\epsilon}(\mathbf{Y}_t, \mathbf{Y}'_t, \mathbf{\Phi}, \mathbf{V}) := \int dP_0(\mathbf{x}) D\mathbf{w} e^{-\mathcal{H}_{t,\epsilon}(\mathbf{x},\mathbf{w};\mathbf{Y}_t,\mathbf{Y}'_t,\mathbf{\Phi},\mathbf{V})}. \quad (97)$$

Finally the interpolating free entropy is

$$f_{n,\epsilon}(t) := \frac{1}{n} \mathbb{E} \ln Z_{t,\epsilon}(\mathbf{Y}_t, \mathbf{Y}'_t, \mathbf{\Phi}, \mathbf{V}). \quad (98)$$

Note that the presence of the perturbation $\epsilon = (\epsilon_1, \epsilon_2)$ induces only a small change in the free entropy, namely of the order of $s_n$:

**Lemma 1** (Small free entropy variation under perturbation). For all $\epsilon_1, \epsilon_2 \in [s_n, 2s_n]$,

$$|f_{n,\epsilon}(0) - f_{n,\epsilon=(0,0)}(0)| \le C s_n \quad (99)$$

for some constant $C$ that only depends on $S$, $\alpha$ and $\varphi$.

*Proof.* Let us compute

$$\left| \frac{df_{n,\epsilon}(0)}{d\epsilon_1} \right| = \frac{1}{2} \left| \mathbb{E} \langle Q \rangle_{n,0,\epsilon} \right| \le \frac{S^2}{2}, \quad (100)$$

by hypothesis (H1). Next we compute

$$\left| \frac{df_{n,\epsilon}(0)}{d\epsilon_2} \right| = \frac{1}{2n} \sum_{\mu=1}^{m} \left| \mathbb{E} \left[ u'_{Y_{0,\mu}}(S_{0,\mu}) \langle u'_{Y_{0,\mu}}(s_{0,\mu}) \rangle_{n,0,\epsilon} \right] \right|. \quad (101)$$

This identity is obtained using very similar steps as in Sec. A.5 to which we refer. Under hypothesis (H2) this quantity is bounded by a constant that only depends on $\alpha$ and $\varphi$. Then by the mean value theorem we obtain $|f_{n,\epsilon}(0) - f_{n,(0,0)}(0)| \le C \|\epsilon\| \le 2\sqrt{2}s_n C$ for some constant $C$ that only depends on $S$, $\alpha$ and $\varphi$. $\square$

One verifies easily, using the Lemma 1, that for all $\epsilon \in \mathcal{B}_n$

$$\left\{
\begin{aligned}
f_{n,\epsilon}(0) &= f_{n,(0,0)}(0) + \mathcal{O}(s_n) \\
&= f_n - \frac{1}{2} + \mathcal{O}(s_n), \\
f_{n,\epsilon}(1) &= \psi_{P_0}(R_1(1, \epsilon)) - \frac{1}{2}(1 + \rho R_1(1, \epsilon)) + \frac{m}{n}\Psi_{P_{\text{out}}}(R_2(1, \epsilon); \rho + 2s_n) \\
&= \psi_{P_0}(\int_0^1 r_\epsilon(t)dt) - \frac{1}{2}(1 + \rho \int_0^1 r_\epsilon(t)dt) + \frac{m}{n}\Psi_{P_{\text{out}}}(\int_0^1 q_\epsilon(t)dt; \rho) + \mathcal{O}(s_n).
\end{aligned}
\right.
\quad (102)$$

where $f_n$ is given by (30) and where $\mathcal{O}(s_n)$ denotes a quantity that is bounded by $Cs_n$ for some constant $C > 0$ that only depends on $S$, $\varphi$ and $\alpha$. For the last equality we used Proposition 17 in Appendix B.1 which says that $\psi_{P_0}$ is $\frac{\rho}{2}$-Lipschitz and, similarly to Proposition 18, it is not difficult to verify that $(q_1, q_2) \mapsto \Psi_{P_{\text{out}}}(q_1; q_2)$ is $C^1$ on the compact set $\{(q_1, q_2) | 0 \le q_1 \le q_2 \le \rho+1\}$ and is thus Lipschitz. We emphasize a crucial property of the interpolating model: It is such that at $t=0$ we recover the original model and thus $f_{n,\epsilon}(0) \approx f_n - 1/2$ (the trivial constant comes from the purely noisy measurements of the second channel in (92)), while at $t=1$ we have the two scalar inference channels and thus the associated terms $\psi_{P_0}$ and $\Psi_{P_{\text{out}}}$ discussed in Sec. 1.3 appear in $f_{n,\epsilon}(1)$. These are precisely the terms appearing in the potential (41).

## 4.2 Free entropy variation along the interpolation path

From the understanding of the previous section, it is at this stage very natural to evaluate the variation of free entropy along the interpolation path, which allows to “compare” the original and purely scalar models
thanks to the identity

$$f_n = f_{n,\epsilon}(0) + \frac{1}{2} + \mathcal{O}(s_n) = f_{n,\epsilon}(1) - \int_0^1 \frac{df_{n,\epsilon}(t)}{dt} dt + \frac{1}{2} + \mathcal{O}(s_n), \quad (103)$$

where the first equality follows from (102). As discussed above, part of the potential (41) appears in $f_{n,\epsilon}(1)$. If the interpolation is properly done, the missing terms required to obtain the potential on the r.h.s. of (103) should naturally appear. Then by choosing the optimal interpolation path thanks to the non-trivial snr dependencies in $t$ (i.e. by selecting the proper interpolating functions $q$ and $r$), we will be able to show the equality between the replica formula and the free entropy $\lim_{n \to \infty} f_n$.

We thus now compute the $t$-derivative of the free entropy along the interpolation path (see Appendix A.5 for the proof). Let $u'_y(x)$ be the derivative (w.r.t. $x$) of $u_y(x)$. Then we have the following.

**Proposition 3 (Free entropy variation).** The derivative of the free entropy (98) verifies, for all $\epsilon \in \mathcal{B}_n$ and all $t \in (0, 1)$

$$\frac{df_{n,\epsilon}(t)}{dt} = -\frac{1}{2} \mathbb{E} \left\langle \left( \frac{1}{n} \sum_{\mu=1}^{m} u'_{Y_{t,\mu}}(S_{t,\mu}) u'_{Y_{t,\mu}}(s_{t,\mu}) - r_{\epsilon}(t) \right) (Q - q_{\epsilon}(t)) \right\rangle_{n,t,\epsilon} + \frac{r_{\epsilon}(t)}{2} (q_{\epsilon}(t) - \rho) + \mathcal{O}_n(1), \quad (104)$$

where $\mathcal{O}_n(1)$ is a quantity that goes to 0 in the $n, m \to \infty$ limit, uniformly in $t \in (0, 1)$, $\epsilon \in \mathcal{B}_n$ and uniformly in the choice of the functions $q_\epsilon$ and $r_\epsilon$. The overlap is

$$Q_n = Q := \frac{1}{n} \mathbf{X}^* \cdot \mathbf{x} = \frac{1}{n} \sum_{i=1}^{n} X_i^* x_i \quad (105)$$

where $\mathbf{x}$ is a sample from the posterior of model (92) associated with the Gibbs bracket $\langle - \rangle_{n,t,\epsilon}$, see (96).

## 4.3 Overlap concentration and fundamental sum rule

The next lemma plays a key role in our proof. Essentially it states that the overlap concentrates around its mean, a behavior called “replica symmetric” in statistical physics. Similar results have been obtained in the context of the analysis of spin glasses [53, 73]. Here we use a formulation taylored to Bayesian inference problems as developed in the context of LDPC codes, random linear estimation [47] and Nishimori symmetric spin glasses [88, 90, 91].

**Proposition 4 (Overlap concentration).** Assume that the interpolation functions $(q_\epsilon)$, $(r_\epsilon)$ are regular, see Definition 1. Let $s_n = \frac{1}{2}n^{-1/16}$ for all $n \ge 1$. Under assumptions (H1), (H2) and (H3) there exists a constant $C(\varphi, S, \alpha)$ that depends only on $S, \varphi$ and $\alpha$ such that

$$\frac{1}{s_n^2} \int_{\mathcal{B}_n} d\epsilon \int_0^1 dt \mathbb{E} \langle (Q - \mathbb{E}\langle Q \rangle_{n,t,\epsilon})^2 \rangle_{n,t,\epsilon} \le \frac{C(\varphi, S, \alpha)}{n^{1/8}}. \quad (106)$$

Proposition 4 follows from Proposition 29 proved in Appendix E.2, combined with (240) and Fubini's theorem. Note from (102) and (41) that the second term appearing in (104) is precisely the missing one that is required in order to obtain the expression of the potential on the r.h.s. of (103). Thus in order to prove Theorem 1 we would like to “cancel” the Gibbs bracket in (104), which is the so called remainder (once integrated over $t$). This is made possible thanks to the adaptive interpolating functions.

One possible way to cancel the remainder is to choose $q_\epsilon(t) = \mathbb{E}\langle Q \rangle_{n,t,\epsilon}$, which is approximately equal to $Q$ because it concentrates by Proposition 4. However, $\mathbb{E}\langle Q \rangle_{n,t,\epsilon}$ depends on $\int_0^t q_\epsilon(v)dv$ (and on $t$, $\int_0^t r_\epsilon(v)dv$ and $\epsilon$ too). The equation $q_\epsilon(t) = \mathbb{E}\langle Q \rangle_{n,t,\epsilon}$ is therefore a first order differential equation over $t \mapsto \int_0^t q_\epsilon(v)dv$.
We will see in details in Sec. 4.4 that it possesses a solution, but for the moment we just assume it exists in order to derive the following fundamental sum rule, which is a core identity in the proof scheme:

**Proposition 5 (Fundamental sum rule).** Assume that the interpolation functions $(q_\epsilon)$ and $(r_\epsilon)$ are regular (see Definition 1). Assume that for all $t \in [0, 1]$ and $\epsilon \in \mathcal{B}_n$ we have $q_\epsilon(t) = \mathbb{E}\langle Q \rangle_{n,t,\epsilon}$. Then

$$f_n = \frac{1}{s_n^2} \int_{\mathcal{B}_n} \left\{ \psi_{P_0} \left( \int_0^1 r_\epsilon(t) dt \right) + \alpha \Psi_{P_{\text{out}}} \left( \int_0^1 q_\epsilon(t) dt; \rho \right) - \frac{1}{2} \int_0^1 q_\epsilon(t) r_\epsilon(t) dt \right\} d\epsilon + \mathcal{O}_n(1), \quad (107)$$

where $\mathcal{O}_n(1)$ denotes a quantity that goes to 0 as $n \to \infty$ uniformly w.r.t. the choice of the interpolation functions.

*Proof.* By the Cauchy-Schwarz inequality

$$\left( \frac{1}{s_n^2} \int_{\mathcal{B}_n} d\epsilon \int_0^1 dt \mathbb{E} \left\langle \left( \frac{1}{n} \sum_{\mu=1}^m u'_{Y_{t,\mu}}(S_{t,\mu}) u'_{Y_{t,\mu}}(s_{t,\mu}) - r_\epsilon(t) \right) (Q - q_\epsilon(t)) \right\rangle_{n,t,\epsilon} \right)^2 \\ \le \frac{1}{s_n^2} \int_{\mathcal{B}_n} d\epsilon \int_0^1 dt \mathbb{E} \left\langle \left( \frac{1}{n} \sum_{\mu=1}^m u'_{Y_{t,\mu}}(S_{t,\mu}) u'_{Y_{t,\mu}}(s_{t,\mu}) - r_\epsilon(t) \right)^2 \right\rangle_{n,t,\epsilon} \times \frac{1}{s_n^2} \int_{\mathcal{B}_n} d\epsilon \int_0^1 dt \mathbb{E} \left\langle (Q - q_\epsilon(t))^2 \right\rangle_{n,t,\epsilon}.$$

The first term of this product is bounded by some constant $C(\varphi, \alpha)$ that only depend on $\varphi$ and $\alpha$, see Appendix A.6. The second term is bounded by $C(\varphi, S, \alpha)n^{-1/8}$ by Proposition 4, since we assumed that for all $\epsilon \in \mathcal{B}_n$ and all $t \in [0, 1]$ we have $q_\epsilon(t) = \mathbb{E}\langle Q \rangle_{n,t,\epsilon}$. We have therefore

$$\left| \frac{1}{s_n^2} \int_{\mathcal{B}_n} d\epsilon \int_0^1 dt \mathbb{E} \left\langle \left( \frac{1}{n} \sum_{\mu=1}^m u'_{Y_{t,\mu}}(S_{t,\mu}) u'_{Y_{t,\mu}}(s_{t,\mu}) - r_\epsilon(t) \right) (Q - q_\epsilon(t)) \right\rangle_{n,t,\epsilon} \right| \le \frac{C(\varphi, S, \alpha)}{n^{1/16}}.$$

Therefore from (104)

$$\frac{1}{s_n^2} \int_{\mathcal{B}_n} d\epsilon \int_0^1 dt \frac{df_{n,\epsilon}(t)}{dt} = \frac{1}{2s_n^2} \int_{\mathcal{B}_n} d\epsilon \int_0^1 dt \{q_\epsilon(t)r_\epsilon(t) - r_\epsilon(t)\rho\} + \mathcal{O}_n(1) + \mathcal{O}(n^{-1/16}). \quad (108)$$

Here the small terms are going to 0 both uniformly w.r.t. to the choice of $q_\epsilon$ and $r_\epsilon$. When replacing (108) in (103) and combining it with (102) we reach the claimed identity (107), but up to the fact that $\Psi_{P_{\text{out}}} \left( \int_0^1 q_\epsilon(t) dt; \rho \right)$ is multiplied by $m/n$ instead of $\alpha$. Recalling that $m/n \to \alpha$ as $m, n \to \infty$ allows to finish the argument (notice that $\Psi_{P_{\text{out}}}$ is continuous and hence bounded on $[0, \rho]$, see Proposition 18). $\square$

We are now ready to prove matching bounds.

## 4.4 Lower and upper matching bounds

We now possess all the necessary tools to prove Theorem 1 in three steps.

(i) We prove that, under assumptions (H1), (H2) and (H3), $\lim_{n \to \infty} f_n = \sup_{r \ge 0} \inf_{q \in [0, \rho]} f_{RS}(q, r)$.

(ii) Under hypothesis (H2), the function $\Psi_{P_{\text{out}}}$ is convex, Lipschitz and non-decreasing (Proposition 19). We thus apply Corollary 8 of Appendix D to get $\sup_{r \ge 0} \inf_{q \in [0, \rho]} f_{RS}(q, r) = \sup_{q \in [0, \rho]} \inf_{r \ge 0} f_{RS}(q, r)$. We then deduce from (i) that $\lim_{n \to \infty} f_n = \sup_{q \in [0, \rho]} \inf_{r \ge 0} f_{RS}(q, r)$ under (H1)-(H2)-(H3).

(iii) Finally, the approximation arguments given in Appendix C permit to relax (H1)-(H2) to the weaker hypotheses (h1)-(h2) and allow to replace the Gaussian assumption (H3) on $\Phi$ by (h1)-(h3)-(h4). The fact that for discrete channels the Gaussian noise can then be removed, allowing to replace (h5.a) (i.e. $\Delta > 0$ treated until here) to (h5.b) (i.e. $\Delta = 0$ and $\varphi$ takes values in $\mathbb{N}$), is proven in Sec. C.3. This
proves the first equality of Theorem 1. The last equality in (43) and the remaining part of Theorem 1 follow then from Lemma 23.

It thus remains to tackle (i), but before that we need a definition. For $t \in [0, 1]$ and $\epsilon \in \mathcal{B}_n$, we write $R^t(\epsilon) = (R_1(t, \epsilon), R_2(t, \epsilon))$. The quantity $\mathbb{E}\langle Q \rangle_{n,t,\epsilon}$ is a function of $n, t, R^t(\epsilon)$ that we write $\mathbb{E}\langle Q \rangle_{n,t,\epsilon} = F_n(t, R^t(\epsilon))$, where $F_n$ is a function defined on

$$D_n := \{(t, r_1, r_2) \in [0, 1] \times \mathbb{R}_+ \times \mathbb{R}_+ \mid r_2 \le \rho t + 2s_n\}. \quad (109)$$

The following proposition, proven in Appendix A.7, will be useful.

**Proposition 6.** $F_n$ is a continuous function from $D_n$ to $[0, \rho]$. Let $D_n^\circ$ denotes the interior of $D_n$. $F_n$ admits partial derivatives with respect to its second and third argument on $D_n^\circ$. These partial derivatives are both continuous and non-negative on $D_n^\circ$.

Let us now start with the lower bound.

### 4.4.1 Lower bound

**Proposition 7 (Lower bound).** *The free entropy (29) verifies*

$$\liminf_{n \to \infty} f_n \ge \sup_{r \ge 0} \inf_{q \in [0, \rho]} f_{RS}(q, r). \quad (110)$$

*Proof.* We consider, for $(\epsilon_1, \epsilon_2) \in \mathcal{B}_n$ and a fixed value $r \in [0, r_{\max}]$, the following 1st order differential equation:

$$y(0) = (\epsilon_1, \epsilon_2) \quad \text{and} \quad \forall t \in [0, 1], \quad y'(t) = (r, F_n(t, y(t))). \quad (111)$$

By the Cauchy-Lipschitz Theorem (see for instance Theorem 3.1 in Chapter V from [92]) this equation admits a (unique) solution that we write $y(\cdot, \epsilon) = (y_1(\cdot, \epsilon), y_2(\cdot, \epsilon))$. The hypotheses of the Cauchy-Lipschitz Theorem are verified, because of Proposition 6. We define then, for all $t \in [0, 1]$,

$$r_{\epsilon}(t) = y'_{1}(t, \epsilon) = r \quad \text{and} \quad q_{\epsilon}(t) = y'_{2}(t, \epsilon) = F_{n}(t, y(t, \epsilon)) \in [0, \rho].$$

We have therefore $R_1(t, \epsilon) = \epsilon_1 + \int_0^t y'_1(s, \epsilon)ds = y_1(t, \epsilon)$ and similarly $R_2(t, \epsilon) = y_2(t, \epsilon)$. We obtain that for all $t \in [0, 1]$,

$$q_{\epsilon}(t) = F_{n}(t, y(t, \epsilon)) = F_{n}(t, (R_{1}(t, \epsilon), R_{2}(t, \epsilon))) = \mathbb{E}\langle Q \rangle_{n,t,\epsilon}.$$

Let us show now that the functions $(q_\epsilon)$ and $(r_\epsilon)$ are regular (see Definition 1). Let $t \in [0, 1]$. The function $R^t: \epsilon \mapsto (R_1(t, \epsilon), R_2(t, \epsilon)) = y(t, \epsilon)$ is the flow of (111) and is thus injective (by unicity of the solution) and $C^1$ because of the regularity properties (see Proposition 6) of $F_n$. The Jacobian of the flow is given by the Liouville formula (see Corollary 3.1 in Chapter V from [92]):

$$\det\left(\frac{\partial R^t}{\partial \epsilon}(\epsilon)\right) = \exp\left(\int_0^t dv \frac{\partial F_n}{\partial y_2}(v, y(v, \epsilon))\right) \ge 1,$$

because by Proposition 6 we have $\partial_{y_2} F_n \ge 0$. We obtain (by the local inversion Theorem) that $R^t$ is a $C^1$ diffeomorphism, and since its Jacobian is greater or equal to 1 the functions $(q_\epsilon)$ and $(r_\epsilon)$ are regular.
We have seen that for all $\epsilon \in \mathcal{B}_n$ and all $t \in [0, 1]$, $q_\epsilon(t) = \mathbb{E}\langle Q \rangle_{t,n,\epsilon}$, so we can apply Proposition 5 to get

$$
\begin{align*}
f_n &= \frac{1}{s_n^2} \int_{\mathcal{B}_n} \left\{ \psi_{P_0}(r) + \alpha \Psi_{P_{\text{out}}}(\int_0^1 q_\epsilon(t) dt; \rho) - \frac{r}{2} \int_0^1 q_\epsilon(t) dt \right\} d\epsilon + \mathcal{O}_n(1) \\
&= \frac{1}{s_n^2} \int_{\mathcal{B}_n} f_{\text{RS}}(\int_0^1 q_\epsilon(t) dt, r) d\epsilon + \mathcal{O}_n(1) \\
&\geq \inf_{q \in [0, \rho]} f_{\text{RS}}(q, r) + \mathcal{O}_n(1)
\end{align*}
$$

and thus $\liminf_{n \to \infty} f_n \ge \inf_{q \in [0, \rho]} f_{RS}(q, r)$. This is true for all $r \in [0, r_{\max}]$ so we get

$$
\liminf_{n \to \infty} f_n \ge \sup_{r \in [0, r_{\max}]} \inf_{q \in [0, \rho]} f_{RS}(q, r). \quad (112)
$$

Let $r \ge r_{\max}$. We have for all $q \in [0, \rho]$, $\partial_q f_{RS}(q, r) = \alpha \Psi'_{P_{\text{out}}}(q) - \frac{r}{2} \le 0$, because $r \ge r_{\max} \ge 2\alpha \Psi'_{P_{\text{out}}}(q)$.
Therefore for all $r \ge r_{\max}$, $\inf_{q \in [0, \rho]} f_{RS}(q, r) = f_{RS}(\rho, r)$ and

$$
\frac{\partial}{\partial r} \inf_{q \in [0, \rho]} f_{\text{RS}}(q, r) = \frac{\partial}{\partial r} f_{\text{RS}}(\rho, r) = \psi'_{P_0}(r) - \frac{\rho}{2} \le 0,
$$

because by Proposition 17, $\psi_{P_0}$ is $\frac{\rho}{2}$-Lipschitz. The function $r \mapsto \inf_{q \in [0, \rho]} f_{RS}(q, r)$ is therefore non-increasing on $[r_{\max}, +\infty)$. Going back to (112), we conclude

$$
\liminf_{n \to \infty} f_n \ge \sup_{r \in [0, r_{\max}]} \inf_{q \in [0, \rho]} f_{RS}(q, r) = \sup_{r \ge 0} \inf_{q \in [0, \rho]} f_{RS}(q, r). \quad (113)
$$

□

4.4.2 Upper bound

**Proposition 8 (Upper bound).** The free entropy (29) verifies

$$
\limsup_{n \to \infty} f_n \le \sup_{r \ge 0} \inf_{q \in [0, \rho]} f_{RS}(q, r). \quad (114)
$$

*Proof.* We consider, for $(\epsilon_1, \epsilon_2) \in \mathcal{B}_n$, the following order-1 system of differential equations:

$$
y(0) = (\epsilon_1, \epsilon_2) \quad \text{and} \quad \forall t \in [0, 1], \quad y'(t) = \begin{pmatrix} 2\alpha \Psi'_{P_{\text{out}}}(F_n(t, y(t))) \\ F_n(t, y(t)) \end{pmatrix}. \quad (115)
$$

By Proposition 18 the function $\Psi'_{P_{\text{out}}}$ is $C^1$ and takes values in $[0, r_{\max}]$. By Proposition 6, the function $F_n$ is continuous, bounded and admits partial derivatives w.r.t. its second and third arguments, that are continuous. We can therefore apply the Cauchy-Lipschitz Theorem as in the proof of Proposition 7: The equation (115) admits a (unique) solution that we write $y(\cdot, \epsilon) = (y_1(\cdot, \epsilon), y_2(\cdot, \epsilon))$. We define then, for all $t \in [0, 1]$,

$$
r_{\epsilon}(t) = y'_{1}(t, \epsilon) = 2\alpha \Psi'_{P_{\text{out}}}(F_{n}(t, y(t, \epsilon))) \in [0, r_{\text{max}}] \quad \text{and} \quad q_{\epsilon}(t) = y'_{2}(t, \epsilon) = F_{n}(t, y(t, \epsilon)) \in [0, \rho].
$$

We have therefore $R_1(t, \epsilon) = \epsilon_1 + \int_0^t y'_1(s, \epsilon)ds = y_1(t, \epsilon)$ and similarly $R_2(t, \epsilon) = y_2(t, \epsilon)$. We obtain that for all $t \in [0, 1]$,

$$
q_{\epsilon}(t) = F_{n}(t, y(t, \epsilon)) = F_{n}(t, (R_{1}(t, \epsilon), R_{2}(t, \epsilon))) = \mathbb{E}\langle Q \rangle_{n,t,\epsilon}.
$$

Let us show now that the functions ($q_\epsilon$) and ($r_\epsilon$) are regular (see Definition 1). Let $t \in [0, 1]$. The function $R^t : \epsilon \mapsto (R_1(t, \epsilon), R_2(t, \epsilon)) = y(t, \epsilon)$ is the flow of (115) and is thus injective and $C^1$ because of the regularity
properties (see Proposition 6) of $F_n$. The Jacobian of the flow is again given by the Liouville formula:

$$ \det\left(\frac{\partial R^t}{\partial \epsilon}(\epsilon)\right) = \exp\left(\int_0^t dv 2\alpha \frac{\partial F_n}{\partial y_1}(v, y(v, \epsilon)) \Psi_{P_{\text{out}}}''(F_n(v, y(v, \epsilon))) + \int_0^t dv \frac{\partial F_n}{\partial y_2}(v, y(v, \epsilon))\right) \geq 1, $$

because by Proposition 6, $\frac{\partial F_n}{\partial y_1}$ and $\frac{\partial F_n}{\partial y_2}$ are both non negative and since $\Psi_{P_{\text{out}}}$ is convex (see Proposition 18), we have also $\Psi_{P_{\text{out}}}'' \geq 0$. We obtain (by the local inversion Theorem) that $R^t$ is a $C^1$ diffeomorphism. Its Jacobian is greater or equal to 1, and the functions ($q_\epsilon$) and ($r_\epsilon$) are therefore regular.

We have seen that for all $\epsilon \in \mathcal{B}_n$ and all $t \in [0, 1]$, $q_\epsilon(t) = \mathbb{E}\langle Q \rangle_{t,n,\epsilon}$, so we can apply Proposition 5 to get

$$ \begin{aligned} f_n &= \frac{1}{s_n^2} \int_{\mathcal{B}_n} \left\{ \psi_{P_0}(\int_0^1 r_\epsilon(t) dt) + \alpha \Psi_{P_{\text{out}}}(\int_0^1 q_\epsilon(t) dt; \rho) - \frac{1}{2} \int_0^1 q_\epsilon(t) r_\epsilon(t) dt \right\} d\epsilon + \mathcal{O}_n(1) \\ &\le \frac{1}{s_n^2} \int_{\mathcal{B}_n} \int_0^1 \left\{ \psi_{P_0}(r_\epsilon(t)) + \alpha \Psi_{P_{\text{out}}}(q_\epsilon(t); \rho) - \frac{1}{2} q_\epsilon(t) r_\epsilon(t) \right\} dt d\epsilon + \mathcal{O}_n(1) \end{aligned} \quad (116) $$

by Jensen's inequality, because by Propositions 17 and 18 the functions $\psi_{P_0}$ and $\Psi_{P_{\text{out}}}$ are convex.

Let us fix $\epsilon \in \mathcal{B}_n$ and $t \in [0, 1]$. By definition of $r_\epsilon$ and $q_\epsilon$, we have

$$ r_{\epsilon}(t) = 2\alpha \Psi'_{P_{\text{out}}}(F_n(t, y(t, \epsilon))) = 2\alpha \Psi'_{P_{\text{out}}}(q_{\epsilon}(t)). \quad (117) $$

The function $g: q \in [0, \rho] \mapsto 2\alpha\Psi_{P_{\text{out}}}(q; \rho) - r_\epsilon(t)q$ is convex by Proposition 18. By equation (117) above, we see that $g'(q_\epsilon(t)) = 0$ and therefore:

$$ \alpha \Psi_{P_{\text{out}}}(q_{\epsilon}(t); \rho) - \frac{1}{2} q_{\epsilon}(t) r_{\epsilon}(t) = \inf_{q \in [0, \rho]} \left\{ \alpha \Psi_{P_{\text{out}}}(q; \rho) - \frac{1}{2} q r_{\epsilon}(t) \right\}. $$

This holds for all $\epsilon \in \mathcal{B}_n$ and all $t \in [0, 1]$. Plugging this back in (116), we get:

$$ \begin{aligned} f_n &\le \frac{1}{s_n^2} \int_{\mathcal{B}_n} \int_0^1 \inf_{q \in [0, \rho]} \left\{ \psi_{P_0}(r_\epsilon(t)) + \alpha \Psi_{P_{\text{out}}}(q; \rho) - \frac{1}{2} q r_\epsilon(t) \right\} dt d\epsilon + \mathcal{O}_n(1) \\ &\le \sup_{r \ge 0} \inf_{q \in [0, \rho]} \left\{ \psi_{P_0}(r) + \alpha \Psi_{P_{\text{out}}}(q; \rho) - \frac{1}{2} q r \right\} + \mathcal{O}_n(1). \end{aligned} $$

This proves Proposition 8. $\square$

From the arguments given at the beginning of the section, this ends the proof of Theorem 1.

# 5 Proofs of the limits of optimal errors

## 5.1 Optimal generalization error: Proof of Theorem 4

### 5.1.1 Formal derivation and proof idea: A teacher-student scenario with side information

Before proving Theorem 4 rigorously, we find useful to provide a conceptual framework allowing to formally derive the generalization error (a framework that will actually serve as a basis for the rigorous derivation presented in the next section). In order to obtain the (generalized) optimal generalization error, we need first to assume that the new “test labels” are also observed by the student in the teacher-student scenario of Sec. 1.2 but with a very low signal-to-noise ratio. The presence of this side information will allow us to use the I-MMSE relation (Proposition 13) to obtain the generalization error when small, but non-zero, information about the test labels is known by the student. Then, by formally taking the limit of vanishing side information
on the resulting expression (and assuming that the large $n$ and vanishing side information limits commute), we will recover the generalization error. We thus now introduce the following “train-test” observation model.

The set of patterns and labels are divided into two sets by the teacher: The training set $S^{\text{tr}}$ of size $m$ that will be used as the main source of information by the student in order to then generalize, and the test set $S^{\text{te}}$ of size $m' = \epsilon n$ that will be used by the teacher in order to evaluate the performance of the student, but also by the student as small additional side information. Let us be more precise: The teacher gives to the student both the patterns and associated labels of the training set, namely $S^{\text{tr}} := \{(Y_\mu; \Phi_\mu)\}_{\mu=1}^m$ (recall the labels are given by (20), (21)). For the test set, the test patterns to classify are given to the student but the associated labels are (almost) not: Let $\epsilon, \lambda \ge 0$. Instead of the test labels $\{\tilde{Y}_\mu\}_{\mu=1}^{m'}$ (that should be totally unknown to the student in the ideal setting), what is given to the student is

$$ U_{\mu} = \sqrt{\lambda} Y'_{\mu} + Z'_{\mu}, \quad \text{for} \quad 1 \le \mu \le m' = \epsilon n, \qquad (118) $$

where $Z'_{\mu} \stackrel{\text{iid}}{\sim} \mathcal{N}(0, 1)$, and $Y'_{\mu}$ is given by

$$ Y'_{\mu} = f(\tilde{Y}_{\mu}), \quad \tilde{Y}_{\mu} \sim P_{\text{out}}\left(\cdot \mid \frac{\Phi'_{\mu} \cdot \mathbf{X}^*}{\sqrt{n}}\right), \qquad (119) $$

where $f: \mathbb{R} \to \mathbb{R}$ is a continuous bounded function and $\Phi'_{\mu} \stackrel{\text{iid}}{\sim} \mathcal{N}(0, \mathbf{I}_n)$ independently of everything else. We will first prove Theorem 4 for continuous bounded functions $f$, and then relax this at the end of the proof. The test set given to the student, in addition of the training set, is $S^{\text{te}} = S^{\text{te}}(\lambda, \epsilon) := \{(U_{\mu} = \sqrt{\lambda} Y'_{\mu} + Z'_{\mu}; \Phi'_{\mu})\}_{\mu=1}^{\epsilon n}$ where $\lambda$ is typically very small. Indeed, we are particularly interested in the case $\lambda, \epsilon \to 0$ when the student has no information about the test labels, which is the ideal setting we want to study. But in order to employ the I-MMSE relation we consider instead very small $\lambda > 0$.

The learning of the classifier $\mathbf{X}^*$ given $S^{\text{tr}}$ and $S^{\text{te}}$ is a slight extension of model (20). Define $\mathbf{Y}' = (Y'_{\mu})_{\mu=1}^{m'}$ as the vector of test labels (before they are corrupted by additional noise through (118)). Then the (generalized) optimal generalization error with side information (i.e. at $\lambda, \epsilon > 0$) in this “train-test” observation model is

$$ \mathcal{E}_{f,n}^{\text{side}}(\lambda, \epsilon) := \min_{\hat{\mathbf{Y}}'} \frac{1}{\epsilon n} \mathbb{E}[\|\mathbf{Y}' - \hat{\mathbf{Y}}'(S^{\text{te}}, S^{\text{tr}})\|^2] = \frac{1}{\epsilon n} \mathbb{E}[\|\mathbf{Y}' - \mathbb{E}[\mathbf{Y}'|S^{\text{te}}, S^{\text{tr}}]\|^2]. \qquad (120) $$

The “true” generalization error (60) is recovered by defining instead $S^{\text{te}} = \Phi'$ or equivalently letting $\lambda, \epsilon \to 0$, i.e. when only $S^{\text{tr}}$ and the test patterns are given to the student: $\lim_{\lambda, \epsilon \to 0} \mathcal{E}_{f,n}^{\text{side}}(\lambda, \epsilon) = \mathcal{E}_{f,n}$. Note that the $f$ function only plays a role in the test set, while the labels of the training data are generated through the “pure model” (20), (21).

From there one can use the I-MMSE relation of Proposition 13 in order to formally compute the limiting $n \to \infty$ expression of (120). Indeed,

$$ \frac{\partial}{\partial \lambda} \frac{1}{n} I(\mathbf{Y}'; \sqrt{\lambda} \mathbf{Y}' + \mathbf{Z}' | \mathbf{Y}, \mathbf{\Phi}, \mathbf{\Phi}') = \frac{\epsilon}{2} \mathcal{E}_{f,n}^{\text{side}}(\lambda, \epsilon). \qquad (121) $$

Fortunately, by a straightforward extension of the interpolation method presented in Sec. 4 one can generalize Theorem 1 to take into account this additional side information and access this mutual information (see the end of the section for the proof):

**Lemma 2.** For all $\epsilon, \lambda \ge 0$ we have

$$ \frac{1}{n} I(\mathbf{Y}'; \sqrt{\lambda} \mathbf{Y}' + \mathbf{Z}' | \mathbf{Y}, \mathbf{\Phi}, \mathbf{\Phi}') \xrightarrow{n \to \infty} \inf_{q \in [0, \rho]} \sup_{r \ge 0} \tilde{i}_{\text{RS}}(q, r, \lambda) - i_{\infty}, \qquad (122) $$
where $i_\infty$ is given by Corollary 2 and

$$ \tilde{i}_{\text{RS}}(q, r, \lambda) := i_{\text{RS}}(q, r) + \epsilon I(f(Y^{(q)}); \sqrt{\lambda} f(Y^{(q)}) + Z'|V) \quad (123) $$

$$ = I_{P_0}(r) + \alpha \mathcal{I}_{P_{\text{out}}}(q; \rho) + \epsilon I(f(Y^{(q)}); \sqrt{\lambda} f(Y^{(q)}) + Z'|V) - \frac{r}{2}(\rho - q). \quad (124) $$

Recall that $Y^{(q)}$ is sampled from the “second scalar channel” (39): $Y^{(q)} \sim P_{\text{out}}(\cdot | \sqrt{q}V + \sqrt{\rho-q}W^*)$, where $V, W^* \stackrel{\text{iid}}{\sim} \mathcal{N}(0, 1)$.

Define the following MMSE function:

$$ M_f : (\lambda, q) \mapsto \text{MMSE}(f(Y^{(q)}|\sqrt{\lambda}f(Y^{(q)}) + Z', V). \quad (125) $$

By concavity arguments detailed in the next section, we have almost everywhere

$$ \lim_{n \to \infty} \frac{\partial}{\partial \lambda} \frac{1}{n} I(\mathbf{Y}'; \sqrt{\lambda} \mathbf{Y}' + \mathbf{Z}' | \mathbf{Y}, \mathbf{\Phi}, \mathbf{\Phi}') = \frac{\partial}{\partial \lambda} \lim_{n \to \infty} \frac{1}{n} I(\mathbf{Y}'; \sqrt{\lambda} \mathbf{Y}' + \mathbf{Z}' | \mathbf{Y}, \mathbf{\Phi}, \mathbf{\Phi}') = \frac{\partial}{\partial \lambda} \inf_{q \in [0, \rho]} \sup_{r \ge 0} \tilde{i}_{\text{RS}}(q, r, \lambda) $$

using Lemma 2 for the last equality. Assuming $\partial_\lambda \inf_{q \in [0, \rho]} \sup_{r \ge 0} \tilde{i}_{\text{RS}}(q, r, \lambda) = \partial_\lambda \tilde{i}_{\text{RS}}(q, r, \lambda)|_{(q_\lambda^*, r_\lambda^*)}$, where $(q_\lambda^*, r_\lambda^*)$ is an optimal couple, (121) and the last identity combined lead to

$$ \frac{\epsilon}{2} \lim_{n \to \infty} \mathcal{E}_{f,n}^{\text{side}}(\lambda, \epsilon) = \epsilon \frac{\partial}{\partial \lambda} I(f(Y^{(q)}); \sqrt{\lambda} f(Y^{(q)}) + Z'|V) \Big|_{q_\lambda^*} = \frac{\epsilon}{2} M_f(\lambda, q_\lambda^*), \quad (126) $$

using again the I-MMSE relation for the last equality. Thus $\lim_{n \to \infty} \mathcal{E}_{f,n}^{\text{side}}(\lambda, \epsilon) = M_f(\lambda, q_\lambda^*)$.

A formal calculation of the vanishing side information limit of $M_f(\lambda, q_\lambda^*)$ gives back $\mathcal{E}_f(q^*(\alpha))$ (recall (61) and $q^*(\alpha)$ is the optimizer of the replica-symmetric formula (43)), so that $\lim_{\lambda, \epsilon \to 0} \lim_{n \to \infty} \mathcal{E}_{f,n}^{\text{side}}(\lambda, \epsilon) = \mathcal{E}_f(q^*(\alpha))$. It is very natural to believe that the vanishing side information limit of $\lim_{n \to \infty} \mathcal{E}_{f,n}^{\text{side}}(\lambda, \epsilon)$ should give back the true asymptotic generalization error. So if one could justify the commutation of limits

$$ \lim_{\lambda, \epsilon \to 0} \lim_{n \to \infty} \mathcal{E}_{f,n}^{\text{side}}(\lambda, \epsilon) = \lim_{n \to \infty} \lim_{\lambda, \epsilon \to 0} \mathcal{E}_{f,n}^{\text{side}}(\lambda, \epsilon) = \lim_{n \to \infty} \mathcal{E}_{f,n} $$

this would end the proof. We prove this point in the next section.

**Proof of Lemma 2:** Extending the interpolation method presented in Sec. 4, one can generalize Theorem 1 to take into account this additional side information. This gives directly

$$ \frac{1}{n} I(\mathbf{X}^*; \mathbf{Y}, \sqrt{\lambda}\mathbf{Y}' + \mathbf{Z}' | \mathbf{\Phi}, \mathbf{\Phi}') \underset{n \to \infty}{\longrightarrow} \tilde{I}_{\infty}(\alpha, \epsilon, \lambda) := \inf_{q \in [0, \rho]} \sup_{r \ge 0} \tilde{i}_{\text{RS}}(q, r, \lambda) \quad (127) $$

where $\tilde{I}_{\text{RS}}(q, r, \lambda)$ is given by

$$ \tilde{I}_{\text{RS}}(q, r, \lambda) := I_{P_0}(r) + \alpha \mathcal{I}_{P_{\text{out}}}(q; \rho) + \epsilon I(W^*; \sqrt{\lambda} f(Y^{(q)}) + Z'|V) - \frac{r}{2}(\rho - q). \quad (128) $$

Conditionally on $(V, f(Y^{(q)}))$, the random variables $W^*$ and $\sqrt{\lambda}f(Y^{(q)}) + Z'$ are independent, therefore

$$ I(f(Y^{(q)}); \sqrt{\lambda}f(Y^{(q)}) + Z'|V) = I(W^*, f(Y^{(q)}); \sqrt{\lambda}f(Y^{(q)}) + Z'|V). $$

Now, by the chain rule of the mutual information we have

$$ I(W^*, f(Y^{(q)}); \sqrt{\lambda}f(Y^{(q)}) + Z'|V) = I(W^*; \sqrt{\lambda}f(Y^{(q)}) + Z'|V) + I(f(Y^{(q)}); \sqrt{\lambda}f(Y^{(q)}) + Z'|V, W^*). $$