††thanks: All authors contributed equally and names are ordered
alphabetically. Correspondence should be sent to<EMAIL_ADDRESS>
# Information limits and Thouless-Anderson-Palmer equations
for spiked matrix models with structured noise

Jean Barbier  Francesco Camilli  Yizhou Xu  The Abdus Salam International
Centre for Theoretical Physics
Strada Costiera 11, Trieste 34151, Italy  Marco Mondelli  Institute of Science
and Technology Austria,
Am Campus 1, 3400 Klosterneuburg, Austria

###### Abstract

We consider a prototypical problem of Bayesian inference for a structured
spiked model: a low-rank signal is corrupted by additive noise. While both
information-theoretic and algorithmic limits are well understood when the
noise is i.i.d. Gaussian, the more realistic case of structured noise still
proves to be challenging. To capture the structure while maintaining
mathematical tractability, a line of work has focused on rotationally
invariant noise. However, existing studies either provide sub-optimal
algorithms or they are limited to a special class of noise ensembles. In this
paper, we establish the first characterization of the information-theoretic
limits for a noise matrix drawn from a general trace ensemble. These limits
are then achieved by an efficient algorithm inspired by the theory of adaptive
Thouless-Anderson-Palmer (TAP) equations. Our approach leverages tools from
statistical physics (replica method) and random matrix theory (generalized
spherical integrals), and it unveils the equivalence between the rotationally
invariant model and a surrogate Gaussian model.

††preprint: APS/123-QED

##  I Introduction

Recovering a low-rank signal from a high-dimensional observation corrupted by
noise is a ubiquitous problem, appearing e.g. in sparse principal component
analysis (PCA) [1, 2], community detection [3, 4], group synchronization [5]
and sub-matrix localization or clustering [6]. In this paper, we consider the
prototypical task of estimating the rank-1 signal
${\mathbf{X}}^{*}{\mathbf{X}}^{*\intercal}\in\mathbb{R}^{N\times N}$ from a
symmetric matrix ${\mathbf{Y}}$ of noisy observations given by

$\displaystyle{\mathbf{Y}}=\frac{\lambda}{N}{\mathbf{X}}^{*}{\mathbf{X}}^{*\intercal}+{\mathbf{Z}},$
(1)

where $\lambda\geq 0$ represents the signal-to-noise ratio (SNR) and
${\mathbf{Z}}\in\mathbb{R}^{N\times N}$ is additive noise. This is often
referred to as the Johnstone spiked covariance model [7], and it was
originally formulated as a probabilistic model for PCA. Starting with the
seminal result of [8], the behavior of eigenvalues and eigenvectors of (1) has
been extensively studied in statistics and random matrix theory, see e.g. [9,
10, 11, 12, 13, 14]. Specifically, the authors of [8] identified a phase
transition phenomenon – named BBP after their initials – tuned by the SNR
$\lambda$: above the transition, the largest eigenvalue of ${\mathbf{Y}}$
detaches from the bulk of the spectrum containing the noise eigenvalues, and
the top eigenvector of ${\mathbf{Y}}$ is correlated with ${\mathbf{X}}$; below
the transition, the largest eigenvalue sits at the edge of the bulk, and the
top eigenvector exhibits vanishing correlation with ${\mathbf{X}}$.

Going beyond the estimator obtained from the top eigenvector of
${\mathbf{Y}}$, a line of work has focused on Approximate Message Passing
(AMP) algorithms. Originally proposed in the context of compressed sensing
[15] and CDMA [16], AMP methods have been developed for numerous high-
dimensional inference problems, including the estimation of low-rank matrices
[17, 18] as in (1), generalized linear regression [19, 20] and inference in
multi-layer models [21]. The popularity of the AMP paradigm stems from its
attractive features: _(i)_ AMP can be tailored to take advantage of structural
information about the signal, in the form of a Bayesian prior; _(ii)_ the AMP
performance in the high-dimensional limit is precisely characterized by a low-
dimensional deterministic recursion known as state evolution [22, 23]; and
_(iii)_ using state evolution, it has been proved that AMP achieves Bayes-
optimal performance in a number of settings [24, 18] and, even when
information-theoretic limits are not met, AMP remains optimal among a vast
class of efficient algorithms [25, 26].

However, most theoretical studies on low-rank matrix estimation are limited by
an _i.i.d._ (independently and identically distributed) hypothesis on the
noise matrix components. In this setting, the fundamental limits of inference
are well understood [27], and they are achieved by an AMP algorithm, unless
there is a statistical-to-computational gap. While some of the results on AMP
can be generalized to the broader class of i.i.d. sub-Gaussian matrices via
universality arguments [28, 29], the i.i.d. assumption is rather limiting:
these models remain structureless, and no concrete correlations can appear in
the data matrices. A way to relax the identicality assumption was proposed in
the mathematical physics literature of spin glasses, in the context of the
Sherrington-Kirkpatrick model. Specifically the authors of [30, 31], and later
[32, 33, 34, 35], consider random couplings whose variances depend on the
index labeling the coupled sites. This idea also appeared earlier in the
context of inference, under the name of spatial coupling [36, 37, 38]. Yet, in
the mentioned works the independence hypothesis still remains crucial.

In the seminal papers [39, 40, 41], the authors considered instead a class of
_rotationally invariant matrices_ , which break the independence between the
elements of the coupling matrices, leaving a model which is still tractable.
The amount of works in this setting, or similar ones, see for instance [42,
43, 44] for spin glasses, and [45, 46, 47, 48, 49, 50] in inference, shows a
growing interest towards the topic. Even if the performance of spectral PCA
can be predicted with a fairly generic additive rotationally invariant noise
(see e.g. [11]), establishing how to factor in also the prior information in
the inference procedure, as well as characterizing information-theoretic
limits has proven to be significantly more challenging.

The recent paper [51] takes a step forward by revealing that, in order to
achieve the information-theoretic limits, it is necessary to apply a peculiar
pre-processing to the data that depends on the type of correlations in the
noise. Despite the new mechanism pinpointed by [51], the analysis has remained
limited only to certain classes of noise distributions, until now. The goal of
the present paper is precisely to elaborate a concise theory, and to formulate
implementable algorithms, that can be tailored to treat _any_ kind of
rotationally invariant noise coming from a trace ensemble, i.e., whose
distribution is encoded in a matrix potential.

We finally note that, while this manuscript was in the final stage of
preparation, the concurrent paper [52] focusing on the same setting was posted
on arXiv. Specifically, [52] develops a new class of AMP algorithms, together
with a rigorous state evolution result for them. A fixed point of that AMP is
shown to match the replica predictions of [51] and, in fact, we verified that
such a fixed point matches the replica predictions we make in the present
paper as well. Therefore, the algorithm of [52], just as the one we propose
here, is conjectured to be Bayes-optimal when no statistical-to-computational
gap is present.

###  I.1 Our contributions

  * _i)_

Using the celebrated _replica method_ [53] and the inhomogeneous spherical
integral of [51], we compute the information-theoretic limits for low-rank
estimation in the presence of rotationally invariant noise. Specifically, we
consider the teacher-student scenario in (1), i.e., the teacher plants a
rank-1 spike matrix inside an extensive rank noise bulk, and we compute the
mutual information between the observation ${\mathbf{Y}}$, and the ground
truth signal ${\mathbf{X}}^{*}{\mathbf{X}}^{*\intercal}$.

    * _ia)_

We simplify and solve numerically the fixed point equations coming from the
replica symmetric variational potential for the mutual information.
Remarkably, thanks to some inherent symmetries of the model, called Nishimori
identities, the fixed point equations, which have functions among their
unknowns, reduce to one simple scalar equation.

    * _ib)_

The channel is not Gaussian, so one cannot use the usual _I-MMSE_ relation
[54] to compute the minimum mean square error (MMSE) from the mutual
information. Nevertheless, our final formula for the mutual information
consists in a variational principle, whose order parameters, at their
stationary values, yield the MMSE when properly combined.

  * _ii)_

We express the mutual information between data and ground truth using the
AdaTAP formalism of [43], relying on the validity of the latter in the
presence of a spike. This approach outputs the pre-processing function as a
functional of the matrix potential of the noise ensemble.

  * _iii)_

We run numerical experiments supporting the consistency between the fixed
point of the TAP equations of point _ii)_ and our replica prediction for the
MMSE.

##  II Setting

Let the data be constructed, conditionally on the unknown spike
${\mathbf{X}}^{*}{\mathbf{X}}^{*\intercal}$, according to (1), with noise
matrix

$\displaystyle{\mathbf{Z}}$
$\displaystyle={\mathbf{O}}^{\intercal}{\mathbf{D}}{\mathbf{O}}\sim
C_{V}\exp\Big{(}-\frac{N}{2}\text{Tr}V({\mathbf{Z}})\Big{)}d{\mathbf{Z}},$ (2)

where $d{\mathbf{Z}}=\prod_{i\leq j}^{N}dZ_{ij}$ and $Z_{ij}=Z_{ji}$ for any
$1\leq i,j\leq N$, and $C_{V}$ is a normalizing constant depending on the
matrix potential $V$. The matrices have $O(1)$ spectral norm (i.e.,
independent on $N$) and the signal vector ${\mathbf{X}}^{*}$ has i.i.d.
entries drawn from the prior distribution $P_{X}$ with second moment equal to
one. The prior thus implies $\|{\mathbf{X}}^{*}\|^{2}=N+o_{N}(1)$ with high
probability. With an abuse of notation, we use $P_{X}({\mathbf{x}})$ instead
of $P_{X}^{\otimes N}({\mathbf{x}})$. Here, and throughout the paper, any
function $f$ applied to a symmetric matrix ${\mathbf{M}}$ with
eigendecomposition
${\mathbf{M}}={\mathbf{O}}{\mathbf{D}}{\mathbf{O}}^{\intercal}$ is actually
applied only to its eigenvalues $(D_{i})_{i\leq N}$:
$f({\mathbf{M}}):={\mathbf{O}}f({\mathbf{D}}){\mathbf{O}}^{\intercal}$,
$f({\mathbf{D}})=\text{diag}(f(D_{i}))_{i\leq N}$.

We stress that, compared to the previous work [51] where the matrix potential
was restricted to be a low degree polynomial, here $V$ can be _any_ analytic
function.

The posterior measure of the problem is given by

$dP_{X|Y}({\mathbf{x}}|{\mathbf{Y}})=\frac{C_{V}}{P_{Y}({\mathbf{Y}})}dP_{X}({\mathbf{x}})e^{-\frac{N}{2}\text{Tr}V({\mathbf{Y}}-\frac{\lambda}{N}{\mathbf{x}}{\mathbf{x}}^{\intercal})},$
(3)

where the evidence, i.e. the probability of the observations is

$P_{Y}({\mathbf{Y}})=C_{V}\int_{\mathbb{R}^{N}}dP_{X}({\mathbf{x}})e^{-\frac{N}{2}\text{Tr}V({\mathbf{Y}}-\frac{\lambda}{N}{\mathbf{x}}{\mathbf{x}}^{\intercal})}.$
(4)

Our main object of interest is the free entropy (or minus the free energy)

$\displaystyle
f_{N}:=\frac{1}{N}\mathbb{E}_{\mathbf{Y}}\log\mathcal{Z}({\mathbf{Y}}),$ (5)

where

$\displaystyle\mathcal{Z}({\mathbf{Y}}):=\int_{\mathbb{R}^{N}}dP_{X}({\mathbf{x}})e^{-\frac{N}{2}\big{(}\text{Tr}V({\mathbf{Y}}-\frac{\lambda}{N}{\mathbf{x}}{\mathbf{x}}^{\intercal})-{\rm
Tr}V({\mathbf{Z}})\big{)}}$ (6)

and in particular its high-dimensional limit $f=\lim_{N\to\infty}f_{N}$. From
the above, we can define a Hamiltonian function
$H_{N}({\mathbf{x}},\lambda,{\mathbf{X}}^{*},{\mathbf{Z}})=H_{N}({\mathbf{x}})$
equal to

$\displaystyle\mathcal{H}_{N}({\mathbf{x}})=\frac{N}{2}\big{(}{\rm
Tr}V({\mathbf{Y}}-\frac{\lambda}{N}{\mathbf{x}}{\mathbf{x}}^{\intercal})-{\rm
Tr}V({\mathbf{Z}})\big{)}$ (7)

and interpret the problem as a spin glass model where ${\mathbf{Z}}$ and
${\mathbf{X}}^{*}$ play the role of quenched disorder. The subtraction of the
term ${\rm Tr}V({\mathbf{Z}})$ is needed to ensure that the Hamiltonian
remains of $O(N)$, and thus the free entropy stays $O(1)$ in the thermodynamic
limit $N\to\infty$. In fact,

$\displaystyle\mathcal{H}_{N}({\mathbf{x}})$
$\displaystyle=\frac{\lambda}{2}\int_{0}^{1}dt\,{\rm
Tr}V^{\prime}\big{(}{\mathbf{Z}}+\frac{t\lambda}{N}\big{(}{\mathbf{X}}^{*}{\mathbf{X}}^{*\intercal}-{\mathbf{x}}{\mathbf{x}}^{\intercal}\big{)}\big{)}$
$\displaystyle\qquad\times\big{(}{\mathbf{X}}^{*}{\mathbf{X}}^{*\intercal}-{\mathbf{x}}{\mathbf{x}}^{\intercal}\big{)}=O(N),$
(8)

since the difference between the two projectors in the last line is at most
rank two and the related eigenvalues remain $O(1)$.

The free entropy is intimately connected with the mutual information between
data and signal:

$\displaystyle I({\mathbf{X}}^{*}{\mathbf{X}}^{*\intercal};{\mathbf{Y}})$
$\displaystyle=-\mathbb{E}_{\mathbf{Y}}\log
P_{Y}({\mathbf{Y}})+\mathbb{E}_{\mathbf{Z}}\log C_{V}e^{-\frac{N}{2}{\rm
Tr}V({\mathbf{Z}})}$ $\displaystyle=-Nf_{N}.$ (9)

Finally, let us recall some basic concepts from random matrix theory. For a
square symmetric random matrix ${\mathbf{M}}\in\mathbb{R}^{N\times N}$, denote
the associated resolvent matrix

${\mathbf{G}}_{M}(z):=(zI_{N}-{\mathbf{M}})^{-1}.$ (10)

Notice that the resolvent matrix shares the same eigenvectors with
${\mathbf{M}}$. Assume that the empirical spectral density (ESD)
$\hat{\rho}_{M}^{(N)}$ of ${\mathbf{M}}$ converges weakly almost surely to a
distribution $\rho_{M}$, i.e.

$\displaystyle\hat{\rho}_{M}^{(N)}=\frac{1}{N}\sum_{i=1}^{N}\delta_{\lambda_{i}({\mathbf{M}})}\xrightarrow[w.a.s.]{N\to\infty}\rho_{M}.$
(11)

We can then define the Stieltjes transform associated with the random matrix
ensemble of ${\mathbf{M}}$:

$g_{M}(z):=\mathbb{E}_{D}\frac{1}{z-D},\quad D\sim\rho_{M}.$ (12)

$g_{M}(z)$ is well defined for $z\in\mathbb{C}$ outside the support of
$\rho_{M}$. Under the aforementioned hypothesis on the ESD, the Stieltjes
transform is closely related to the resolvent matrix through

$\lim_{N\to\infty}\frac{1}{N}\text{Tr}{\mathbf{G}}_{M}(z)=g_{M}(z).$ (13)

Denote the inverse function of $g_{M}(z)$ as $\zeta_{M}(g)$, and then the
R-transform of ${\mathbf{M}}$ is given by

$R_{M}(g):=\zeta(g)-\frac{1}{g}.$ (14)

The resolvent and the R-transform play a crucial role in our analysis since
they encode all the relevant combinatorics of the random matrix ensemble. As
we shall see later, the resolvent allows us to define a new family of order
parameters, which are quadratic forms of replicas drawn from the posterior
with the resolvent mediating their product. In one shot, this new family of
overlaps encompasses all the new order parameters that were introduced in
[51], as well as additional ones.

##  III Information limits through
the replica method

In this section we analyse the spiked model with generic rotationally
invariant noise using the powerful replica method from statistical physics of
disordered systems [53, 55]. This method is non-rigorous, but it is believed
to be exact for the asymptotic analysis of a broad class of spin glass,
inference and learning models. Historically, one of the first proofs of its
exactness was given for the Sherrington-Kirkpatrick model by Guerra [56] and
Talagrand [57], and later remarkably refined by Panchenko leveraging
ultrametricity [58]. Moreover, the replica symmetry assumption we are going to
employ during the analysis is intimately connected to concentration-of-measure
phenomena proven in broad generality in optimal Bayesian inference. We
therefore conjecture the analysis below leads to asymptotically exact
formulas. For further dicussions on the topic, we refer the reader to [59, 60,
61, 62, 63, 64, 65].

We now state our main result from the information-theoretic perspective. This
comes in the form of a variational formula for the free entropy. The physical
meaning of some of the order parameters entering these formulas is given in
the replica analysis of the next section.

######  Result (Information-theoretic limits: Replica free entropy and minimum
mean-square error).

Let $\Gamma$ be an arbitrary contour in the complex plane $\mathbb{C}$ that
encircles all eigenvalues of the matrix
${\mathbf{Y}}-{\mathbf{x}}{\mathbf{x}}^{\intercal}/N$ for any choice of
${\mathbf{x}}$ with positive measure according to the prior $P_{X}$. Let
$Z\sim\mathcal{N}(0,1)$ and $X\sim P_{X}$ as well as
$D,D_{1},D_{2}\sim\rho_{Z}$ i.i.d. from the noise asymptotic spectral density.
The replica free entropy at the replica symmetric level (which is exact in
Bayes-optimal inference) is given by

$\displaystyle f$ $\displaystyle={\rm extr}\Big{\\{}\frac{1}{4\pi
i}\oint_{\Gamma}dz\Big{[}V^{\prime}(z)\log(1+\lambda
B(z))-2\frac{\lambda\hat{B}(z)M(z)^{2}}{1-\lambda
g_{Z}(z)}+2\hat{M}(z)M(z)+2\hat{B}(z)B(z)\Big{]}$
$\displaystyle\qquad+\frac{v\hat{v}}{2}-\hat{m}m+\frac{q\hat{q}}{2}+\mathbb{E}_{Z,X}\log\int_{\mathbb{R}}dP_{X}(x)\exp\Big{(}\sqrt{\hat{q}}Zx-\frac{\hat{q}+\hat{v}}{2}x^{2}+\hat{m}Xx\Big{)}+m\bar{m}+\frac{v\bar{v}}{2}-\frac{q\bar{q}}{2}-\frac{1}{2}$
$\displaystyle\qquad-\frac{1}{2}\mathbb{E}_{D}\log(\bar{v}-\bar{q}+2\tilde{B}(D))-\frac{1}{2}\mathbb{E}_{D}\frac{\bar{q}-(\bar{q}+\tilde{M}(D))^{2}}{\bar{v}-\bar{q}+2\tilde{B}(D)}-\frac{1}{2}\log(v-q)-\frac{q-m^{2}}{2(v-q)}\Big{\\}}$
$\displaystyle\qquad+\frac{1}{4\pi
i}\oint_{\Gamma}dzV^{\prime}(z)\log(1-\lambda g_{Z}(z)),$ (15)

with an extremization w.r.t. nine scalar order parameters
$(m,q,v,\hat{m},\hat{q},\hat{v},\bar{m},\bar{q},\bar{v})$ and four functions
$(M,B,\hat{M},\hat{B})$ from $\mathbb{C}$ to $\mathbb{C}$. The extremization
selects the solution of the saddle point equations, obtained by equating to
zero the gradient of the replica potential $\\{\cdots\\}$, which maximizes it.
After simplifications of the replica saddle point equations, this can also be
written as

$\displaystyle f=\max_{\mathcal{M}^{\rm RS}}f^{\rm
RS}(m,\hat{m})-\frac{\lambda}{2}\mathbb{E}_{D}V^{\prime}(D),$ (16)
$\displaystyle f^{\rm
RS}(m,\hat{m}):=-\frac{\lambda^{2}}{2}\mathbb{E}_{D_{1},D_{2}}Q(D_{1})Q(D_{2})H(D_{1})H(D_{2})\frac{V^{\prime}(D_{1})-V^{\prime}(D_{2})}{D_{1}-D_{2}}-\frac{m^{2}}{2(1-m)}-\frac{1}{2}\log(1-m)-\frac{m}{2}$
$\displaystyle\quad+\mathbb{E}_{Z,X}\log\int_{\mathbb{R}}dP_{X}(x)\exp\Big{(}\sqrt{\hat{m}}Zx-\frac{\hat{m}}{2}x^{2}+\hat{m}Xx\Big{)}+\frac{1}{2}\mathbb{E}_{D}\log
H(D)-\frac{1}{2}\mathbb{E}_{D}\Big{(}\hat{m}-\frac{1}{1-m}-Q(D)^{2}\Big{)}H(D),$
(17)

with

$\displaystyle H(x)$
$\displaystyle:=\Big{(}\frac{1}{1-m}-\hat{m}-J(x)\Big{)}^{-1},$ (18)
$\displaystyle J(x)$ $\displaystyle:=\lambda
V^{\prime}(x)-\lambda^{2}\mathbb{E}_{D\sim\rho_{Z}}\frac{V^{\prime}(x)-V^{\prime}(D)}{x-D},$
(19)

and $Q(x)$ is the solution of

$Q(x)=\hat{m}-\frac{1}{1-m}+\lambda^{2}\mathbb{E}_{D}\frac{V^{\prime}(x)-V^{\prime}(D)}{x-D}Q(D)H(D).$
(20)

Finally, $\mathcal{M}^{\rm RS}$ represents the set of solution(s) of the
following fixed point equations:

$\displaystyle\begin{cases}\hat{m}=-R_{J({\mathbf{Z}})}(1-m),\\\
m=\mathbb{E}X\langle x\rangle_{\hat{m}},\end{cases}$ (21)

where $\langle\cdot\rangle_{\hat{m}}$ denotes the posterior of a scalar
Gaussian channel with signal-to-noise ratio $\hat{m}$:

$\langle f(x)\rangle_{\hat{m}}:=\frac{\int
dP_{X}(x)e^{\sqrt{\hat{m}}Zx+\hat{m}xX-\frac{\hat{m}}{2}x^{2}}f(x)}{\int
dP_{X}(x)e^{\sqrt{\hat{m}}Zx+\hat{m}xX-\frac{\hat{m}}{2}x^{2}}}.$ (22)

Recall that $\int P_{X}(dx)x^{2}=1$. Let $m_{*}$ be the value of the order
parameter $m$ picked by the above extremization, i.e., the solution of (21)
which maximizes $f^{\rm RS}(m,\hat{m})$. The asymptotic minimum mean-square
error corresponding to the Bayes-optimal estimator (i.e., the posterior mean)
$\mathbb{E}[{\mathbf{X}}^{*}{\mathbf{X}}^{*\intercal}\mid{\mathbf{Y}}]$ is
given by

$\displaystyle\lim_{N\to\infty}\frac{1}{N^{2}}\mathbb{E}\big{\|}{\mathbf{X}}^{*}{\mathbf{X}}^{*\intercal}-\mathbb{E}[{\mathbf{X}}^{*}{\mathbf{X}}^{*\intercal}\mid{\mathbf{Y}}]\big{\|}_{\rm
F}^{2}=1-m_{*}^{2}.$ (23)

###  III.1 Analysis by the replica method

We now provide the derivation of the previous result. Before replicating the
partition function we are going to re-express it in a more amenable form. We
start by extracting the matrix entering the potential in the log-likelihood
term of the partition function using Cauchy’s formula. We will then repeatedly
use Sherman-Morrison’s formula to deal with inverses of rank-one perturbations
of matrices:

$\displaystyle\mathcal{Z}$
$\displaystyle=\int_{\mathbb{R}^{N}}dP_{X}({\mathbf{x}})e^{-\frac{N}{2}\text{Tr}\big{(}V\big{(}{\mathbf{Y}}-\frac{\lambda}{N}{\mathbf{x}}{\mathbf{x}}^{\intercal}\big{)}-{\rm
Tr}V({\mathbf{Z}})\big{)}}$ (24)
$\displaystyle=\mathbb{E}_{\mathbf{x}}e^{\frac{N}{2}{\rm
Tr}V({\mathbf{Z}})-\frac{N}{4\pi
i}\text{Tr}\oint_{\Gamma}dzV(z)\big{(}zI_{N}-{\mathbf{Y}}+\frac{\lambda}{N}{\mathbf{x}}{\mathbf{x}}^{\intercal}\big{)}^{-1}}$
$\displaystyle=\mathbb{E}_{\mathbf{x}}\exp\Big{(}\frac{N}{2}{\rm
Tr}V({\mathbf{Z}})-\frac{N}{4\pi i}\text{Tr}\oint_{\Gamma}dzV(z)$
$\displaystyle\qquad\times\Big{[}{\mathbf{G}}_{Y}(z)-\frac{\lambda}{N}\frac{{\mathbf{G}}_{Y}(z){\mathbf{x}}{\mathbf{x}}^{\intercal}{\mathbf{G}}_{Y}(z)}{1+\frac{\lambda}{N}{\mathbf{x}}^{\intercal}{\mathbf{G}}_{Y}(z){\mathbf{x}}}\Big{]}\Big{)}$
$\displaystyle=C\,\mathbb{E}_{\mathbf{x}}\exp\Big{(}\frac{N}{4\pi
i}\oint_{\Gamma}dzV(z)\frac{\lambda}{N}\frac{{\mathbf{x}}^{\intercal}{\mathbf{G}}_{Y}(z)^{2}{\mathbf{x}}}{1+\frac{\lambda}{N}{\mathbf{x}}^{\intercal}{\mathbf{G}}_{Y}(z){\mathbf{x}}}\Big{)}$
$\displaystyle=C\,\mathbb{E}_{\mathbf{x}}\exp\Big{(}\frac{N}{4\pi
i}\oint_{\Gamma}dzV^{\prime}(z)\log\Big{(}1+\frac{\lambda}{N}{\mathbf{x}}^{\intercal}{\mathbf{G}}_{Y}(z){\mathbf{x}}\Big{)}\Big{)},$

where we used $\partial_{z}{\mathbf{G}}_{Y}(z)=-{\mathbf{G}}_{Y}(z)^{2}$ and
an integration by part in the last equality. Here

$C:=\exp\Big{(}-\frac{N}{2}\big{(}\text{Tr}V({\mathbf{Y}})-{\rm
Tr}V({\mathbf{Z}})\big{)}\Big{)}$

is a multiplicative constant which we will compute separately.

We are now ready to replicate the partition function. We denote with the
replica index 0 the signal ${\mathbf{X}}^{*}={\mathbf{x}}_{0}$. We then get
that

$\displaystyle\mathbb{E}\mathcal{Z}^{n}$
$\displaystyle=C^{n}\int\prod_{a=0}^{n}dP_{X}({\mathbf{x}}_{a})$
$\displaystyle\ \ \times\exp\Big{(}\frac{N}{4\pi
i}\sum_{a=1}^{n}\oint_{\Gamma}dzV^{\prime}(z)\log(1+\lambda
B^{aa}(z))\Big{)},$

where we introduce the following order parameters for $1\leq a\leq n$, which
are generalized data-dependent self-overlap functions,

$B^{aa}(z):=\frac{1}{N}{\mathbf{x}}_{a}^{\intercal}{\mathbf{G}}_{Y}(z){\mathbf{x}}_{a}\in\mathbb{C}.$
(25)

We now introduce delta functions in their Fourier form, together with the
conjugate order parameters. Jointly denote
$\mathcal{D}[\hat{B},B]:=\prod_{a=1}^{n}\mathcal{D}[\hat{B}^{aa},B^{aa}]$ for
the differential element in functional (path) integrals. The replicated
partition function $\mathbb{E}\mathcal{Z}^{n}$ then becomes

$\displaystyle
C^{n}\int\mathcal{D}[\hat{B},B]\prod_{a=0}^{n}dP_{X}({\mathbf{x}}_{a})$
$\displaystyle\ \times\exp\Big{(}\frac{N}{4\pi i}\sum_{a\leq
n}\oint_{\Gamma}dzV^{\prime}(z)\log(1+\lambda B^{aa}(z))\Big{)}$
$\displaystyle\ \times\mathbb{E}_{\mathbf{O}}\exp\Big{(}\sum_{a\leq
n}\oint_{\Gamma}dz\hat{B}^{aa}(z)(NB^{aa}(z)-{\mathbf{x}}_{a}^{\intercal}{\mathbf{G}}_{Y}(z){\mathbf{x}}_{a})\Big{)}.$

The $z$-integral is on the contour $\Gamma$. In order to perform the quenched
average over the Haar distributed eigenbasis ${\mathbf{O}}$ of the noise, we
need to decompose explicitly the data into signal plus noise. The last term
can then be simplified using Sherman-Morrison again:

$\displaystyle\frac{{\mathbf{x}}_{a}^{\intercal}{\mathbf{G}}_{Y}(z){\mathbf{x}}_{a}}{N}\\!$
$\displaystyle=\\!\frac{1}{N}{\mathbf{x}}_{a}^{\intercal}\Big{(}{\mathbf{G}}_{Z}(z)+\frac{\lambda{\mathbf{G}}_{Z}(z){\mathbf{x}}_{0}{\mathbf{x}}_{0}^{\intercal}{\mathbf{G}}_{Z}(z)}{N(1-\frac{\lambda}{N}{\mathbf{x}}_{0}^{\intercal}{\mathbf{G}}_{Z}(z){\mathbf{x}}_{0})}\Big{)}{\mathbf{x}}_{a}$
$\displaystyle\\!=\\!\frac{{\mathbf{x}}_{a}^{\intercal}{\mathbf{G}}_{Z}(z){\mathbf{x}}_{a}}{N}+\frac{\lambda
M^{aa}(z)^{2}}{1-\lambda g_{Z}(z)},$ (26)

where we introduce the main order parameters, i.e., generalized overlaps
between replicas and the ground-truth (which, by Bayes-optimality, also
corresponds to the generalized overlap between different replicas): for $1\leq
a\leq n$,

$M^{aa}(z):=\frac{1}{N}{\mathbf{x}}_{a}^{\intercal}{\mathbf{G}}_{Z}(z){\mathbf{x}}_{0}\in\mathbb{C}.$
(27)

Note that by definition we have

$\displaystyle\lim_{N\to\infty}\frac{1}{N}{\mathbf{x}}_{0}^{\intercal}{\mathbf{G}}_{Z}(z){\mathbf{x}}_{0}=g_{Z}(z).$
(28)

Let us make a remark concerning the generalized overlap function $M^{aa}(z)$.
By expanding in series the resolvent around $z\to+\infty$, we realize that it
corresponds to the generating function for an infinite family of scalar
overlaps
$(\frac{1}{N}{\mathbf{x}}_{a}^{\intercal}{\mathbf{Z}}^{k}{\mathbf{x}}_{0})_{k\geq
0}$. A similar observation can be made for $B^{aa}(z)$ which encodes
$(\frac{1}{N}{\mathbf{x}}_{a}^{\intercal}{\mathbf{Z}}^{k}{\mathbf{x}}_{a})_{k\geq
0}$. The first few of these overlaps are the order parameters identified in
[51]. We note that the analysis of [51] is restricted to low-degree
polynomials for the potential $V$. This is because the number of order
parameters – and, therefore, the number of replica saddle point equations –
grows with the degree of the polynomial, which quickly leads to intractable,
non interpretable, formulas. However, by identifying these generating
functions order parameters, we can easily encode such infinite families of
scalars and write down compact equations. This is one key mechanism that
allows us to treat generic potential functions $V$. Similar ideas have been
used to study gradient-flow dynamics in [66, 67, 68, 69].

We consider a replica symmetric ansatz: for all $1\leq a\leq n$, we set

${\rm RS\
ansatz}\mbox{:}\quad\begin{cases}(B^{aa}(z),\hat{B}^{aa}(z))=(B(z),\hat{B}(z)),\\\
(M^{aa}(z),\hat{M}^{aa}(z))=(M(z),\hat{M}(z)).\end{cases}$

This implies the following simplifications for the replicated partition
function: letting $\mathcal{D}[\cdots]=\mathcal{D}[\hat{B},B,\hat{M},M]$, we
have

$\displaystyle\mathbb{E}\mathcal{Z}^{n}$
$\displaystyle=C^{n}\int\mathcal{D}[\cdots]\exp\Big{(}\frac{Nn}{4\pi
i}\oint_{\Gamma}dzV^{\prime}(z)\log(1+\lambda B(z))$ $\displaystyle\qquad-4\pi
i\frac{\lambda\hat{B}(z)M(z)^{2}}{1-\lambda
g_{Z}(z)}+Nn\oint_{\Gamma}dz\hat{M}(z)M(z)$
$\displaystyle\qquad+Nn\oint_{\Gamma}dz\hat{B}(z)B(z)+NnI^{\rm RS}\Big{)},$
(29)

where we also define

$\displaystyle\exp(nNI^{\rm RS})$
$\displaystyle:=\mathbb{E}_{\mathbf{O}}\exp\Big{(}-\sum_{a\leq
n}\oint_{\Gamma}dz\Big{(}\hat{B}(z){\mathbf{x}}_{a}^{\intercal}{\mathbf{G}}_{Z}(z){\mathbf{x}}_{a}$
$\displaystyle\qquad+\hat{M}(z){\mathbf{x}}_{0}^{\intercal}{\mathbf{G}}_{Z}(z){\mathbf{x}}_{a}\Big{)}\Big{)}.$
(30)

The above integral $I^{\rm RS}$ is an instance of the _inhomogeneous spherical
integral_ defined and analysed in [51]. A key property of this integral is
that it depends on the replicas only through their overlap structure which,
under a replica symmetric ansatz, reads

$\displaystyle{\rm RS\
ansatz}\mbox{:}\quad\begin{cases}{\mathbf{x}}_{a}^{\intercal}{\mathbf{x}}_{b}/N=q,\quad
1\leq a<b\leq n,\\\ {\mathbf{x}}_{0}^{\intercal}{\mathbf{x}}_{a}/N=m,\quad
1\leq a\leq n,\\\ {\mathbf{x}}_{a}^{\intercal}{\mathbf{x}}_{a}/N=v,\quad 1\leq
a\leq n.\end{cases}$ (31)

Let us define

$\displaystyle\tilde{M}(x)$ $\displaystyle:=\frac{1}{2\pi
i}\oint_{\Gamma}\frac{\hat{M}(z)dz}{z-x},\quad\tilde{B}(x):=\frac{1}{2\pi
i}\oint_{\Gamma}\frac{\hat{B}(z)dz}{z-x}.$ (32)

It is important to note that in general $\tilde{B}(x)\neq\hat{B}(x)$ and
$\tilde{M}(x)\neq\hat{M}(x)$ because they might not be holomorphic. Then, the
result of the inhomogeneous spherical integral reads

$\displaystyle I^{\rm
RS}(q,m,v,\tilde{B},\tilde{M})=\text{extr}_{(\bar{m},\bar{v},\bar{q})}\Big{\\{}m\bar{m}+\frac{v\bar{v}}{2}-\frac{q\bar{q}}{2}$
$\displaystyle\quad-\frac{1}{2}\mathbb{E}_{D}\log(\bar{v}-\bar{q}+2\tilde{B}(D))-\frac{1}{2}\mathbb{E}_{D}\frac{\bar{q}-(\bar{q}+\tilde{M}(D))^{2}}{\bar{v}-\bar{q}+2\tilde{B}(D)}\Big{\\}}$
$\displaystyle\quad-\frac{1}{2}-\frac{1}{2}\log(v-q)-\frac{q-m^{2}}{2(v-q)}+O(n).$
(33)

Now, fixing the overlap definitions using additional delta functions in
Fourier form, under the same replica ansatz for the Fourier conjugates, we
reach

$\displaystyle\mathbb{E}\mathcal{Z}^{n}$
$\displaystyle=C^{n}\int\mathcal{D}[\cdots]\exp\Big{(}\frac{Nn}{4\pi
i}\oint_{\Gamma}dzV^{\prime}(z)\log(1+\lambda B(z))-4\pi
i\frac{\lambda\hat{B}(z)M(z)^{2}}{1-\lambda g_{Z}(z)}$
$\displaystyle\qquad+Nn\oint_{\Gamma}dz\hat{M}(z)M(z)+Nn\oint_{\Gamma}dz\hat{B}(z)B(z)+NnI^{RS}(q,m,v,\tilde{B},\tilde{M})+\frac{Nn}{2}v\hat{v}-Nnm\hat{m}$
$\displaystyle\qquad-\frac{Nn(n-1)}{2}q\hat{q}+N\log\int_{\mathbb{R}^{n+1}}\prod_{a=0}^{n}dP_{X}(x_{a})\exp\Big{(}-\frac{\hat{v}}{2}\sum_{a=1}^{n}x_{a}^{2}+\hat{m}\sum_{a=1}^{n}x_{0}x_{a}+\hat{q}\sum_{1\leq
a<b\leq n}x_{a}x_{b}\Big{)}\Big{)}.$ (34)

In order to decouple the replicas in the last term we use an Hubbard-
Stratonovitch transform: with $Z\sim\mathcal{N}(0,1)$,

$\displaystyle\mathbb{E}_{(x_{a})}e^{-\frac{\hat{v}}{2}\sum_{a=1}^{n}x_{a}^{2}+\hat{m}\sum_{a=1}^{n}x_{0}x_{a}+\hat{q}\sum_{1\leq
a<b\leq n}x_{a}x_{b}}$
$\displaystyle\quad=\mathbb{E}_{Z,x_{0}}\big{(}\mathbb{E}_{x}e^{-\frac{\hat{v}+\hat{q}}{2}x^{2}+\hat{m}x_{0}x+\sqrt{\hat{q}}Zx}\big{)}^{n}.$

The final steps are an integration with respect to the order parameters by
saddle point, followed by an application of the replica trick (assuming
commutation of thermodynamic and replica limits for the saddle point
integration):

$\displaystyle\lim_{N\to\infty}\frac{1}{N}\mathbb{E}\ln\mathcal{Z}$
$\displaystyle=\lim_{N\to\infty}\lim_{n\to
0}\frac{1}{Nn}\ln\mathbb{E}\mathcal{Z}^{n}$ $\displaystyle=\lim_{n\to
0}\lim_{N\to\infty}\frac{1}{Nn}\ln\mathbb{E}\mathcal{Z}^{n}.$ (35)

This finally gives the average limiting free entropy expression (15) after
having changed variables as follows $(2\pi i\hat{M},2\pi
i\hat{B})\to(\hat{M},\hat{B})$.

To conclude, let us compute the contribution related to the constant $C$. One
may think it is irrelevant, however it is a function of the SNR and concurs to
the value of the mutual information and some of its fundamental properties
(for instance, monotonicity and concavity). We can again use Cauchy’s integral
representation and the Shermann-Morrison inversion formula for
$V({\mathbf{Y}})$:

$\displaystyle\frac{1}{N}\mathbb{E}\log C$
$\displaystyle=-\frac{1}{2}\mathbb{E}{\rm
Tr}\big{(}V({\mathbf{Y}})-V({\mathbf{Z}})\big{)}$
$\displaystyle=-\frac{1}{4\pi i}\mathbb{E}{\rm
Tr}\oint_{\Gamma}dzV(z)\big{(}{\mathbf{G}}_{Y}(z)-{\mathbf{G}}_{Z}(z)\big{)}$
$\displaystyle=-\frac{1}{4\pi
i}\mathbb{E}\oint_{\Gamma}dzV(z)\frac{\lambda}{N}\frac{{\mathbf{x}}_{0}^{\intercal}{\mathbf{G}}^{2}_{Z}(z){\mathbf{x}}_{0}}{(1-\frac{\lambda}{N}{\mathbf{x}}_{0}^{\intercal}{\mathbf{G}}_{Z}(z){\mathbf{x}}_{0})},$

which in the limit converges to

$\displaystyle\frac{1}{N}\mathbb{E}\log C$ $\displaystyle\to\frac{1}{4\pi
i}\oint_{\Gamma}dzV(z)\frac{\lambda\partial_{z}g_{Z}(z)}{(1-\lambda
g_{Z}(z))}$ $\displaystyle=\frac{1}{4\pi
i}\oint_{\Gamma}dzV^{\prime}(z)\log(1-\lambda g_{Z}(z))\,.$ (36)

###  III.2 Simplifying the saddle point equations

In this section we show how to go from the variational formulation (15) for
the free entropy, to a simpler formula (16), (17) with only two order
parameters. Before giving the complete set of saddle point equations derived
from (15), we stress that the physical meaning of some order parameters makes
it possible to fix directly their values to their expectation (assuming
concentration), obtainable using the Nishimori identities, see [70,
Proposition 15] for a proof.

Nishimori identity. For any bounded function $f$ of the signal
${\mathbf{X}}^{*}$, the data ${\mathbf{Y}}$ and of conditionally i.i.d.
samples from the posterior ${\mathbf{x}}^{j}\sim P_{X\mid
Y}(\,\cdot\mid{\mathbf{Y}})$, $j=1,2,\ldots,n$, we have that

$\displaystyle\mathbb{E}\langle
f({\mathbf{Y}},{\mathbf{X}}^{*},{\mathbf{x}}^{2},\ldots,{\mathbf{x}}^{n})\rangle=\mathbb{E}\langle
f({\mathbf{Y}},{\mathbf{x}}^{1},{\mathbf{x}}^{2},\ldots,{\mathbf{x}}^{n})\rangle,$
(37)

where the bracket notation $\langle\,\cdot\,\rangle$ is used for the joint
expectation over the posterior samples $({\mathbf{x}}^{j})_{j\leq n}$, and
$\mathbb{E}$ is over the signal ${\mathbf{X}}^{*}$ and data ${\mathbf{Y}}$.

To begin with, recall that we fixed $v$ to be the squared norm of a sample
from the posterior re-scaled by the number of components. Assume that
concentration effects take place, i.e. that the order parameters of the
problem are limiting values of self-averaging quantities, as they should in
this optimal setting [62], and denote

$\langle f(x)\rangle:=\frac{\int
dP_{X}(x)e^{\sqrt{\hat{q}}Zx+\hat{m}xX-\frac{\hat{q}+\hat{v}}{2}x^{2}}f(x)}{\int
dP_{X}(x)e^{\sqrt{\hat{q}}Zx+\hat{m}xX-\frac{\hat{q}+\hat{v}}{2}x^{2}}}.$ (38)

Using the Nishimori identity, we have that

$v=\lim_{N\to\infty}\frac{1}{N}\mathbb{E}\langle\|{\mathbf{x}}\|^{2}\rangle=\lim_{N\to\infty}\frac{1}{N}\mathbb{E}\|{\mathbf{X}}^{*}\|^{2}=1.$
(39)

We have $\hat{v}=0$ because by Bayes-optimality the constraint $v=1$ is
already enforced by the prior without the need of a delta constraint. The
Nishimori identity also imposes

$m=q.$ (40)

Moreover, $B(z)$ is also fixed by the Nishimori identity (below $N$ is large
and equalities are understood up to a vanishing correction as $N\to\infty$):

$\displaystyle B(z)$
$\displaystyle=\frac{1}{N}\mathbb{E}\langle{\mathbf{x}}^{\intercal}{\mathbf{G}}_{Y}(z){\mathbf{x}}\rangle$
$\displaystyle=\frac{1}{N}\mathbb{E}{\mathbf{X}}^{*\intercal}{\mathbf{G}}_{Y}(z){\mathbf{X}}^{*}$
$\displaystyle=\frac{1}{N}\mathbb{E}{\mathbf{X}}^{*\intercal}\Big{[}{\mathbf{G}}_{Z}(z)+\frac{\lambda}{N}\frac{{\mathbf{G}}_{Z}(z){\mathbf{X}}^{*}{\mathbf{X}}^{*\intercal}{\mathbf{G}}_{Z}(z)}{1-\frac{\lambda}{N}{\mathbf{X}}^{*\intercal}{\mathbf{G}}_{Z}(z){\mathbf{X}}^{*}}\Big{]}{\mathbf{X}}^{*}$
$\displaystyle=g_{Z}(z)+\frac{\lambda g_{Z}(z)^{2}}{1-\lambda
g_{Z}(z)}=\frac{g_{Z}(z)}{1-\lambda g_{Z}(z)},$ (41)

where we used the Nishimori identity in the second equality, Sherman-Morrison
in the third one, and

$\frac{1}{N}\mathbb{E}{\mathbf{X}}^{*\intercal}{\mathbf{G}}_{Z}(z){\mathbf{X}}^{*}=g_{Z}(z)$

in the fourth (by independence of the signal and noise).

We now state the complete set of saddle point equations obtained by cancelling
the gradient of the replica free entropy potential $\\{\cdots\\}$ in (15)
w.r.t. the order parameters. The parameter w.r.t. which the gradient is
computed in order to obtain a certain saddle point equation is reported in the
round parenthesis. Let

$\displaystyle H(x)$ $\displaystyle:=(\bar{v}-\bar{q}+2\tilde{B}(x))^{-1},$
(42) $\displaystyle R(x)$
$\displaystyle:=\bar{q}-v(\bar{m}+\tilde{M}(x))^{2}.$ (43)

Let $D\sim\rho_{Z}$ be drawn from the spectral distribution of the noise. Then
the saddle point equations read

$\displaystyle(\hat{m}):\ m=\mathbb{E}X\langle x\rangle$
$\displaystyle(\hat{q}):\ q=\mathbb{E}\langle x\rangle^{2}$
$\displaystyle(\hat{v}):\ v=\mathbb{E}\langle x^{2}\rangle=1$
$\displaystyle(\bar{m}):\ m=-\mathbb{E}_{D}(\bar{m}+\hat{M}(D))H(D)$
$\displaystyle(\bar{q}):\ q=\mathbb{E}_{D}H(D)^{2}R(D)$
$\displaystyle(\bar{v}):\ v=\mathbb{E}_{D}H(D)(1-H(D)R(D))$
$\displaystyle(m):\ -\hat{m}+\bar{m}+\frac{m}{1-q}=0$ $\displaystyle(q):\
\hat{q}-\bar{q}=\frac{q}{1-q}$ $\displaystyle(v):\ \bar{v}=1$
$\displaystyle(\hat{B}):\ B(z)-\frac{\lambda M(z)^{2}}{1-\lambda g_{Z}(z)}$
$\displaystyle\qquad\qquad=-\mathbb{E}_{D}\Big{[}\frac{1}{D-z}(H(D)-R(D)H(D)^{2})\Big{]}$
$\displaystyle(\hat{M}):\
M(z)=\mathbb{E}_{D}\Big{[}\frac{1}{D-z}(\bar{m}+\tilde{M}(D))H(D)\Big{]}$
$\displaystyle(B):\ \frac{\lambda V^{\prime}(z)}{1+\lambda
B(z)}+2\hat{B}(z)=0$ $\displaystyle(M):\
-\frac{2\lambda\hat{B}(z)M(z)}{1-\lambda g_{Z}(z)}+\hat{M}(z)=0,$

where we used $v=1$, $\hat{v}=0$.

We can now simplify. Firstly, from $(\hat{m})$, $(\hat{q})$ and (40), we have

$m=q=\mathbb{E}X\langle x\rangle$ (44)

and $\hat{m}=\hat{q}$. Then, from $(m)$ and $(q)$, we have

$\bar{q}=\bar{m}=\hat{m}-\frac{m}{1-m}.$ (45)

From $(\bar{q})$ and $(\bar{v})$, we have

$m=1-\mathbb{E}_{D}H(D).$ (46)

From $(B)$ and $(M)$, we have

$\displaystyle\hat{B}(z)$
$\displaystyle=-\frac{\lambda}{2}V^{\prime}(z)(1-\lambda g_{Z}(z)),$ (47)
$\displaystyle\hat{M}(z)$ $\displaystyle=-\lambda^{2}V^{\prime}(z)M(z).$ (48)

Let us keep in mind that $\hat{B}(z)$ is not holomorphic, and $\hat{M}(z)$ is
in general not holomorphic.

The complex integrals in $\tilde{B}$ and $\tilde{M}$ can be performed by the
residue theorem.

$\displaystyle\tilde{B}(x)$ $\displaystyle=\frac{1}{2\pi
i}\oint_{\Gamma}\frac{dz}{z-x}\Big{[}-\frac{\lambda}{2}V^{\prime}(z)(1-\lambda
g_{Z}(z))\Big{]}$ $\displaystyle=-\frac{\lambda}{2}V^{\prime}(x)+\frac{1}{4\pi
i}\oint_{\Gamma}\frac{dz}{z-x}\Big{[}\lambda^{2}\mathbb{E}_{D}\frac{V^{\prime}(z)}{D-z}\Big{]}$
$\displaystyle=-\frac{\lambda}{2}V^{\prime}(x)+\frac{\lambda^{2}}{2}\mathbb{E}_{D}\frac{V^{\prime}(x)-V^{\prime}(D)}{x-D},$
(49)

and similarly for the other function

$\displaystyle\tilde{M}(x)$ $\displaystyle=-\frac{1}{2\pi
i}\oint_{\Gamma}\frac{dz}{z-x}\lambda^{2}V^{\prime}(z)M(z)$
$\displaystyle=\frac{1}{2\pi
i}\oint_{\Gamma}\frac{dzV^{\prime}(z)}{z-x}\mathbb{E}\Big{[}\frac{\lambda^{2}}{z-D}(\bar{m}+\tilde{M}(D))H(D)\Big{]}$
$\displaystyle=\lambda^{2}\mathbb{E}_{D}\frac{V^{\prime}(x)-V^{\prime}(D)}{x-D}(\bar{m}+\tilde{M}(D))H(D),$
(50)

where we used $(\hat{M})$, (47) and (48).

Finally, let us denote

$J(x):=-2\tilde{B}(x),$ (51)

which leads to (19) according to (49). Recall that $\bar{v}=1$ according to
$(v)$. Then, the definition of $H(x)$ and (45) lead to (18). Combining (18)
with (46) gives

$\hat{m}=-R_{J({\mathbf{Z}})}(1-m).$ (52)

This forms the fixed point equations together with

$m=\mathbb{E}X\langle x\rangle_{\hat{m}},$ (53)

after noticing that the posterior mean (38) simplifies to (22) when
$\hat{q}=\hat{m}$ and $\hat{v}=0$. The above analysis gives the simplified
saddle point equations (21).

To obtain (17), we simply represent all order parameters through $m$,
$\hat{m}$ and $Q(x):=\bar{m}+\tilde{M}(x)$, while (20) is obtained from (50).
We also simplify two contour integrals as follows. The first integral is

$\displaystyle-\frac{1}{2\pi i}\oint_{\Gamma}dz\hat{M}(z)M(z)$
$\displaystyle=\frac{1}{2\pi
i}\oint_{\Gamma}dz\lambda^{2}V^{\prime}(z)M(z)^{2}$
$\displaystyle=\frac{1}{2\pi
i}\oint_{\Gamma}dz\lambda^{2}V^{\prime}(z)\mathbb{E}_{D1,D2}\frac{Q(D_{1})Q(D_{2})H(D_{1})H(D_{2})}{(D_{1}-z)(D_{2}-z)}$
$\displaystyle=\lambda^{2}\mathbb{E}Q(D_{1})Q(D_{2})H(D_{1})H(D_{2})\frac{V^{\prime}(D_{1})-V^{\prime}(D_{2})}{D_{1}-D_{2}},$

with $D_{1},D_{2}$ i.i.d. from $\rho_{Z}$ and where we have used (48),
$(\hat{M})$ and the residue theorem. The second integral is

$\displaystyle\frac{1}{2\pi i}\oint_{\Gamma}dz\hat{B}(z)B(z)$
$\displaystyle=-\frac{\lambda}{4\pi i}\oint_{\Gamma}dzV^{\prime}(z)g_{Z}(z)$
$\displaystyle=-\frac{\lambda}{2}\mathbb{E}_{D}V^{\prime}(D),$ (54)

where we have used (47) and the residue theorem.

Finally, notice that at the saddle point, specifically using (41), the first
term in the free entropy (15) is precisely

$\displaystyle-\frac{1}{4\pi i}\oint_{\Gamma}dzV^{\prime}(z)\log(1-\lambda
g_{Z}(z)),$ (55)

which cancels with the constant evaluated in (III.1).

###  III.3 Relation to a Gaussian surrogate model

We can notice that the replica saddle point equations (as well as the TAP
equations defined in the next section) are closely related to those appearing
in the the Gaussian noise case. In fact, the replica saddle point equations
for Gaussian noise (with SNR $\tilde{\lambda}$) reads as follows:

$\hat{m}=\tilde{\lambda}^{2}m,\qquad m=\mathbb{E}X\langle x\rangle_{\hat{m}}.$
(56)

Therefore, by choosing

$\tilde{\lambda}=\sqrt{-R_{J({\mathbf{Z}})}(1-m_{*})/m_{*}},$

where $m_{*}$ takes the value at the extremizer of (21), the Gaussian model
and the rotational invariant model share the same fixed point (and thus same
minimum mean-square error, but necessarily the same mutual information).

##  IV Thouless-Anderson-Palmer
free entropy and equations

Along the lines of [51], we employ here the adaTAP approach [43] as an
alternative to the replica method. AdaTAP offers the advantage of expressing
the free entropy as a variational principle over an extensive number of
parameters, which can be interpreted as site marginal means and variances for
every variable in the system, namely signal components in our setting.
Strictly speaking, before our work, the validity of this approach was verified
only for spin glass models containing at most two-body interactions, mediated
by rotationally invariant matrices. In contrast, our model is not quadratic,
but the precise point of the first steps of the replica computation is to make
it quadratic. Hence, we take again as a starting point (recall that notation
$\mathbb{E}_{\mathbf{x}}$ means integration against the prior $P_{X}^{\otimes
N}$):

$\displaystyle\mathcal{Z}\propto\mathbb{E}_{\mathbf{x}}\exp\Big{(}\frac{N}{4\pi
i}\oint_{\Gamma}dzV^{\prime}(z)\log(1+\frac{\lambda}{N}{\mathbf{x}}^{\intercal}{\mathbf{G}}_{Y}(z){\mathbf{x}})\Big{)}$
$\displaystyle\propto\mathbb{E}_{\mathbf{x}}\int\mathcal{D}[B,i\hat{B}]\exp\Big{(}\frac{N}{4\pi
i}\oint_{\Gamma}dzV^{\prime}(z)\log(1+\lambda B(z))\Big{)}$
$\displaystyle\qquad\times\exp\Big{(}\oint_{\Gamma}\hat{B}(z)(NB(z)-{\mathbf{x}}^{\intercal}{\mathbf{G}}_{Y}(z){\mathbf{x}})dz\Big{)}$
$\displaystyle=\int\mathcal{D}[B,i\hat{B}]\exp\Big{(}\frac{N}{4\pi
i}\oint_{\Gamma}dzV^{\prime}(z)\log(1+\lambda B(z))$
$\displaystyle\qquad+N\oint_{\Gamma}\hat{B}(z)B(z)dz\Big{)}\mathbb{E}_{\mathbf{x}}\exp\Big{(}\frac{1}{2}{\mathbf{x}}^{\intercal}J({\mathbf{Y}}){\mathbf{x}}\Big{)},$
(57)

where $J({\mathbf{Y}})=-2\oint_{\Gamma}\hat{B}(z){\mathbf{G}}_{Y}(z)dz$ ends
up being equal to (19). The last factor in (57) is precisely a two-body
(quadratic) model, whose partition function is therefore computable via the
adaTAP approach. Define

$\displaystyle\varphi({\mathbf{Y}}):=\log\mathbb{E}_{\mathbf{x}}\exp\Big{(}\frac{1}{2}{\mathbf{x}}^{\intercal}J({\mathbf{Y}}){\mathbf{x}}\Big{)}.$
(58)

Then, following Opper and Winther’s prescription, the TAP representation of
the free entropy $\varphi({\mathbf{Y}})$ reads

$\displaystyle\varphi_{\rm TAP}({\mathbf{Y}})$
$\displaystyle=\sum_{i=1}^{N}\Big{[}\lambda_{i}m_{i}+\frac{\gamma_{i}}{2}(\sigma_{i}+m_{i}^{2})+\frac{c_{i}\sigma_{i}-\log\sigma_{i}-1}{2}$
$\displaystyle+\frac{1}{2}{\mathbf{m}}^{\intercal}J({\mathbf{Y}}){\mathbf{m}}+\log\int
dP_{X}(x)e^{-\lambda_{i}x-\frac{\gamma_{i}}{2}x^{2}}\Big{]}$
$\displaystyle-\frac{1}{2}\log\det(\text{diag}(\mathbf{c})-J({\mathbf{Y}}))$
(59)

where ${\mathbf{m}}=(m_{i})_{i\leq N},\mathbf{c}=(c_{i})_{i\leq N}$ and an
extremization w.r.t. the parameters
$\lambda_{i},m_{i},\gamma_{i},\sigma_{i},c_{i}$ is intended. Since we are
interested only in leading terms, we can carry out some simplifications of the
above.

First of all, a common assumption in the thermodynamic limit (see [43, 71]) is
that of homogeneous variances $\sigma_{i}=\sigma$ together with
$\gamma_{i}=\gamma$, which in turn yields $c_{i}=c$. Let us now focus on the
determinant term in $\varphi_{\rm TAP}$, which is supposed to reconstruct the
_Onsager reaction term_ in the TAP equations. We argue that at leading order
it does not depend on the spike, nor on the specific realization of
${\mathbf{Z}}$. The leading contribution is determined just by the spectral
distribution of ${\mathbf{Z}}$. Assume $J({\mathbf{Y}})$ is a regular enough
non-linearity applied to the eigenvalues of a matrix whose spectrum consists
of a bulk of eigenvalues inherited by ${\mathbf{Z}}$, plus possibly one spike
detached from the mentioned bulk. The non-linearity changes the shape of the
spectrum, but it preserves the bulk-plus-spike structure. A spike of one or
few eigenvalues cannot alter the spectral distribution of the overall matrix.
From these considerations we get

$\displaystyle\frac{1}{2}\log\det(\text{diag}(\mathbf{c})-J({\mathbf{Y}}))$
$\displaystyle\simeq\frac{1}{2}\log\det(cI_{N}-J({\mathbf{Z}}))$
$\displaystyle\simeq\frac{N}{2}\mathbb{E}\log(c-J(D)),$ (60)

where $\mathbb{E}$ is intended over $D$, distributed according to the
asymptotic spectral density of the noise. Hence, the TAP representation of the
overall free entropy of the model at leading order in $N$ reads (equalities
are up to a constant and a vanishing correction $o_{N}(1)$):

$\displaystyle Nf_{\rm TAP}$
$\displaystyle=\sum_{i=1}^{N}\Big{[}\lambda_{i}m_{i}+\frac{\gamma}{2}(m_{i}^{2}+\sigma)+\frac{c\sigma-1-\log\sigma}{2}$
$\displaystyle+\log\int
dP_{X}(x)e^{-\lambda_{i}x-\frac{\gamma}{2}x^{2}}\Big{]}-\frac{N}{2}\mathbb{E}\log(c-J(D))$
$\displaystyle+\frac{1}{2}{\mathbf{m}}^{\intercal}J({\mathbf{Y}}){\mathbf{m}}+\frac{N}{4\pi
i}\oint_{\Gamma}dzV^{\prime}(z)\log(1+\lambda B(z))$
$\displaystyle+N\oint_{\Gamma}\hat{B}(z)B(z)dz.$ (61)

Again, extremization is intended w.r.t. $\lambda_{i},m_{i},\gamma,\sigma,c$
and the two functions $\hat{B},B$. As anticipated, extremizing w.r.t. $B$ and
$\hat{B}$ only results in matching the coupling matrix $J({\mathbf{Y}})$ with
the pre-processed matrix using (19).

We can now write the remaining TAP equations. Define for future convenience
the Bayes “denoiser”

$\displaystyle\eta(a,b):=\frac{\int dP_{X}(x)e^{ax-\frac{bx^{2}}{2}}\,x}{\int
dP_{X}(y)e^{ay-\frac{by^{2}}{2}}}.$ (62)

Extremization w.r.t. $c$ yields $\sigma=g_{J({\mathbf{Z}})}(c)$, namely

$\displaystyle c=\frac{1}{\sigma}+R_{J({\mathbf{Z}})}(\sigma).$ (63)

Recall that we are looking for equilibrium configurations that satisfy
Nishimori identities, so in the limit we must have

$\displaystyle\frac{1}{N}\sum_{i=1}^{N}(\sigma+m_{i}^{2})=1\text{, that is
}\sigma=1-\tilde{q},$ (64)

where $\tilde{q}:=\frac{1}{N}\sum_{i=1}^{N}m_{i}^{2}$. Cancelling the
$\sigma$-derivative one then gets

$\displaystyle\gamma=-c+\frac{1}{\sigma}=-R_{J({\mathbf{Z}})}(\sigma)=-R_{J({\mathbf{Z}})}(1-\tilde{q}).$
(65)

Finally, extremizing w.r.t. $\lambda_{i}$ and $m_{i}$ yields the final TAP
equations for our model

$\displaystyle{\mathbf{m}}=\eta(J({\mathbf{Y}}){\mathbf{m}}+\gamma{\mathbf{m}},\gamma),\quad\gamma=-R_{J({\mathbf{Z}})}(1-\bar{q}),$
(66)

where $\eta$ is applied component-wise to the vector in the first entry.

The outcome of this analysis is a fundamental equivalence between the original
model with non-linear likelihood governed by $V$ and a model quadratic in
${\mathbf{x}}$, with effective interaction matrix $J({\mathbf{Y}})$. The
equivalence is information-theoretic: the two models have asymptotically the
same free entropy and, therefore, mutual information and minimum mean-square
error. The main advantage of this equivalence resides in the fact that since
the effective model is quadratic, we are able to employ known analytical and
algorithmic approaches in the next section.

Figure 1: The performance of the TAP iterations (dots) matches well the
replica prediction for the minimum mean-square error (solid lines), for
various eigenvalue distributions for the noise. Top left: Quartic potential.
Top right: Sestic potential. Bottom left: Marchenko–Pastur distribution of
eigenvalues. Bottom right: Normal distribution of eigenvalues. Green dashed
lines (which overlap perfectly the red solid lines) denote the theoretical
performance of spectral PCA as predicted by [11]. We do not include the
performance of spectral PCA for the normal distribution of eigenvalues due to
numerical instabilities.

######  Algorithm (Optimal data pre-processing, and TAP iterations).

Define the optimal pre-processing function $J(x)$ as in (19). Let
${\mathbf{m}}^{0}=\sqrt{N}v_{1}({\mathbf{Y}})$ with $v_{1}({\mathbf{Y}})$ the
unit norm first principal component of ${\mathbf{Y}}$. For $t\geq 1$ the TAP
iterations are defined as

$\displaystyle{\mathbf{m}}^{t+1}=\tau{\mathbf{m}}^{t}+(1-\tau)\eta(J({\mathbf{Y}}){\mathbf{m}}^{t}+\gamma^{t}{\mathbf{m}}^{t-1},\gamma^{t}),\quad\gamma^{t}=-R_{J({\mathbf{Z}})}(1-\tilde{q}^{t}),\quad\tilde{q}^{t}=\frac{\|{\mathbf{m}}^{t}\|^{2}}{N},$
(67)

where $\eta$ is defined in (62), we use a damping parameter $\tau\in[0,1)$ for
improved numerical stability, and $R_{J({\mathbf{Z}})}$ is the $R$-transform
of the asymptotic spectral density of $J({\mathbf{Z}})$.

##  V From TAP equations to an efficient iterative algorithm

Now that the information-theoretic analysis has been performed through the
replica method, we switch towards algorithmic aspects and focus on how to
efficiently match the performance predicted by our theory, based on the
Thouless-Anderson-Palmer formalism [72, 73, 43].

###  V.1 TAP iterations

We can now state our second main result (Algorithm above), which is of an
algorithmic nature. This comes in the form of a Bayes-optimal pre-processing
function to be applied to the data matrix, and an efficient iterative
algorithm exploiting it, in order to reach a solution of the TAP equations.

We draw attention to the time indexing in the algorithm (67). The update rule
is inspired by that of a usual AMP algorithm, and as supported by our
numerical experiments, it proves to be effective to match the results
predicted by the replica analysis. Despite its similarity, with an evident
candidate Onsager reaction term, $\gamma^{t}{\mathbf{m}}^{t-1}$, (67) cannot
be really regarded as an AMP algorithm per se, since we have no theoretical
guarantee that the components of the iterates
$J({\mathbf{Y}}){\mathbf{m}}^{t}+\gamma^{t}{\mathbf{m}}^{t-1}$ have Gaussian
statistics.

###  V.2 Numerical experiments

To verify the match between replica theory and algorithm, we choose four
concrete examples for the noise potential. $(i)$ A quartic potential
$V(x)=\gamma x^{4}/4$ with $\gamma=16/{27}$. Its eigenvalue distribution is
given by

$\rho_{Z}(x)=\frac{1}{2\pi}(2a^{2}\gamma+\gamma x^{2})\sqrt{4a^{2}-x^{2}},$
(68)

where $a=3/4$. $(ii)$ A sestic potential $V(x)=\xi x^{6}/6$ with $\xi=27/80$.
Its eigenvalue distribution is given by

$\rho_{Z}(x)=\frac{1}{2\pi}(6a^{4}\xi+2a^{2}\xi x^{2}+\xi
x^{4})\sqrt{4a^{2}-x^{2}},$ (69)

where $a=\sqrt{2/3}$. In both cases, the constants are chosen in order to
enforce unit variance for the spectral densities. These two cases are (among)
those studied in the previous paper [51], the sestic potential being the
highest degree of a polynomial that the techniques in the reference allowed to
study. With the present contribution we can now analyse arbitrary potentials,
even non-polynomial ones such as $(iii)$ eigenvalues following the
Marchenko–Pastur distribution:

$\rho_{Z}(x)=\frac{1}{2\pi\sigma^{2}}\frac{\sqrt{(\lambda_{+}-x)(x-\lambda_{-})}}{\alpha
x},$ (70)

where $\lambda_{\pm}=\sigma^{2}(1\pm\sqrt{\alpha})^{2}$ and $\alpha=0.2$. The
associated potential is given by $V(x)=[(1-1/\alpha)/x+1/\alpha]/\sigma^{2}$.
Finally, as an example of more exotic distribution of eigenvalues we
considered $(iv)$ eigenvalues following a standard normal distribution
truncated between $[-5,5]$. Its potential has probably no analytical
expression, so we numerically calculated its derivative through [74]

$V^{\prime}(x)=2\,\text{P.V.}\int\frac{\rho_{Z}(d\lambda)}{x-\lambda},$ (71)

where P.V. denotes the Cauchy principal value. Thus, we are able to calculate
$J({\mathbf{Z}})$ and its R-transform. In all cases, the noise is properly
normalized such that $\mathbb{E}(D-\mathbb{E}D)^{2}=1$, which is also how we
determine $\sigma^{2}$ in (70). We consider two choices for the prior $P_{X}$:
the standard Gaussian law $\mathcal{N}(0,1)$ and Rademacher law
$\frac{1}{2}\delta_{-1}+\frac{1}{2}\delta_{1}$.

The algorithm uses a PCA initialization [49, 48], that can be obtained
efficiently via the power method. For the normally distributed eigenvalues,
however, we manually choose an initialization having positive correlation
$\sqrt{0.5}$ with the ground truth for numerical stability. In all
experiments, we use $N=2000$ and show the results averaged over $10$ trials
and the corresponding standard deviation. In some cases, about $20\%$ of the
trials does not converge to the right fixed point, which we exclude when
gathering statistics. Moreover, we fix the Onsager coefficient to its fixed
point value predicted by the replica theory for numerical stability. We use a
damping $\tau=0.9$ in all experiments.

Fig. 1 shows that in all successful cases, our algorithm approaches the Bayes-
optimal performance predicted by our replica-based theory. We therefore
conjecture that the fixed point performance of our algorithm matches that of
the minimum mean-square error estimator, when no statistical-to-computational
gap is present, as in all the above experiments.

Note that, as already predicted in [51], PCA remains Bayes-optimal when the
prior of the signal is Gaussian (or more generically rotationally invariant),
regardless of the noise eigenvalue distribution.

Another remark is that, as can be seen from Fig. 1, a Rademacher prior leads
to better numerical stability, due to a more attractive fixed point for the
dynamics. In fact, the Rademacher prior – being more informative – constrains
the signal estimate more strongly than a Gaussian prior.

Finally, the equivalence between the models with structured and Gaussian
noises does not only hold at the level of static (thermodynamic) properties.
Indeed, Fig. 2 numerically verifies that the gap between the TAP iterations
run for the Gaussian model and those run for the rotational invariant model is
small.

Figure 2: Comparison between the TAP iterations for the quartic model (with
$\lambda=2$) and its information-theoretic equivalent Gaussian surrogate
model. The error bars represent the standard deviation computed over 10
trials. The dashed black line represents the MMSE predicted by the replica
theory.

## Acknowledgments

J.B., F.C. and Y.X. were funded by the European Union (ERC, CHORAL, project
number 101039794). Views and opinions expressed are however those of the
authors only and do not necessarily reflect those of the European Union or the
European Research Council. Neither the European Union nor the granting
authority can be held responsible for them. M.M. was supported by the 2019
Lopez-Loreta Prize. J.B. acknowledges discussions with TianQi Hou at the
initial stage of the project, as well as with Antoine Bodin.

## References

  * Johnstone and Lu [2009a] I. M. Johnstone and A. Y. Lu, On consistency and sparsity for principal components analysis in high dimensions, Journal of the American Statistical Association 104, 682 (2009a).
  * Johnstone and Lu [2009b] I. M. Johnstone and A. Y. Lu, Sparse principal components analysis, arXiv preprint arXiv:0901.4392 (2009b).
  * Abbe [2018] E. Abbe, Community detection and stochastic block models: Recent developments, Journal of Machine Learning Research 18, 1 (2018).
  * Moore [2017] C. Moore, The computer science and physics of community detection: landscapes, phase transitions, and hardness, arXiv preprint, arXiv:1702.00467 (2017).
  * Perry  _et al._ [2018] A. Perry, A. S. Wein, A. S. Bandeira, and A. Moitra, Message-passing algorithms for synchronization problems over compact groups, Communications on Pure and Applied Mathematics 71, 2275 (2018).
  * Lesieur  _et al._ [2015] T. Lesieur, F. Krzakala, and L. Zdeborová, Mmse of probabilistic low-rank matrix estimation: Universality with respect to the output channel, in  _Annual Allerton Conference_ (2015).
  * Johnstone [2001] I. M. Johnstone, On the distribution of the largest eigenvalue in principal components analysis, The Annals of statistics 29, 295 (2001).
  * Baik  _et al._ [2005] J. Baik, G. B. Arous, and S. Péché, Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices, The Annals of Probability 33, 1643 (2005).
  * Bai and Yao [2012] Z. Bai and J. Yao, On sample eigenvalues in a generalized spiked population model, Journal of Multivariate Analysis 106, 167 (2012).
  * Baik and Silverstein [2006] J. Baik and J. W. Silverstein, Eigenvalues of large sample covariance matrices of spiked population models, Journal of multivariate analysis 97, 1382 (2006).
  * Benaych-Georges and Nadakuditi [2011] F. Benaych-Georges and R. R. Nadakuditi, The eigenvalues and eigenvectors of finite, low rank perturbations of large random matrices, Advances in Mathematics 227, 494 (2011).
  * Capitaine  _et al._ [2009] M. Capitaine, C. Donati-Martin, and D. Féral, The largest eigenvalues of finite rank deformation of large Wigner matrices: convergence and nonuniversality of the fluctuations, The Annals of Probability 37, 1 (2009).
  * Féral and Péché [2007] D. Féral and S. Péché, The largest eigenvalue of rank one deformation of large Wigner matrices, Communications in mathematical physics 272, 185 (2007).
  * Knowles and Yin [2013] A. Knowles and J. Yin, The isotropic semicircle law and deformation of Wigner matrices, Communications on Pure and Applied Mathematics 66, 1663 (2013).
  * Donoho  _et al._ [2009] D. L. Donoho, A. Maleki, and A. Montanari, Message-passing algorithms for compressed sensing, Proceedings of the National Academy of Sciences 106, 18914 (2009).
  * Kabashima [2003] Y. Kabashima, A cdma multiuser detection algorithm on the basis of belief propagation, Journal of Physics A: Mathematical and General (2003).
  * Deshpande and Montanari [2014] Y. Deshpande and A. Montanari, Information-theoretically optimal sparse PCA, in  _IEEE International Symposium on Information Theory_ (2014) pp. 2197–2201.
  * Montanari and Venkataramanan [2021] A. Montanari and R. Venkataramanan, Estimation of low-rank matrices via approximate message passing, Annals of Statistics 45, 321 (2021).
  * Rangan [2011] S. Rangan, Generalized approximate message passing for estimation with random linear mixing, in  _International Symposium on Information Theory_ (2011) pp. 2168–2172.
  * Sur and Candès [2019] P. Sur and E. J. Candès, A modern maximum-likelihood theory for high-dimensional logistic regression, Proceedings of the National Academy of Sciences 116, 14516 (2019).
  * Manoel  _et al._ [2017] A. Manoel, F. Krzakala, M. Mézard, and L. Zdeborová, Multi-layer generalized linear estimation, in  _2017 IEEE International Symposium on Information Theory (ISIT)_ (IEEE, 2017) pp. 2098–2102.
  * Bayati and Montanari [2011] M. Bayati and A. Montanari, The dynamics of message passing on dense graphs, with applications to compressed sensing, IEEE Transactions on Information Theory 57, 764 (2011).
  * Bolthausen [2014] E. Bolthausen, An iterative construction of solutions of the TAP equations for the Sherrington–Kirkpatrick model, Communications in Mathematical Physics 325, 333 (2014).
  * Barbier  _et al._ [2019] J. Barbier, F. Krzakala, N. Macris, L. Miolane, and L. Zdeborová, Optimal errors and phase transitions in high-dimensional generalized linear models, Proceedings of the National Academy of Sciences 116, 5451 (2019).
  * Celentano  _et al._ [2020] M. Celentano, A. Montanari, and Y. Wu, The estimation error of general first order methods, in  _Conference on Learning Theory (COLT)_ (2020) pp. 1078–1141.
  * Montanari and Wein [2022] A. Montanari and A. S. Wein, Equivalence of approximate message passing and low-degree polynomials in rank-one matrix estimation, arXiv preprint arXiv:2212.06996 (2022).
  * Lelarge and Miolane [2018] M. Lelarge and L. Miolane, Fundamental limits of symmetric low-rank matrix estimation, Probability Theory and Related Fields 173, 859 (2018).
  * Bayati  _et al._ [2015] M. Bayati, M. Lelarge, and A. Montanari, Universality in polytope phase transitions and message passing algorithms, Annals of Applied Probability 25, 753 (2015).
  * Chen and Lam [2021] W. K. Chen and W.-K. Lam, Universality of approximate message passing algorithms, Electronic Journal of Probability 26, 36 (2021).
  * Barra  _et al._ [2013] A. Barra, P. Contucci, E. Mingione, and D. Tantari, Multi-species mean field spin glasses. rigorous results, Annales Henri Poincaré 16, 691 (2013).
  * Panchenko [2013a] D. Panchenko, The free energy in a multi-species sherrington-kirkpatrick model, The Annals of Probability 43 (2013a).
  * Alberici  _et al._ [2021a] D. Alberici, F. Camilli, P. Contucci, and E. Mingione, The multi-species mean-field spin-glass on the nishimori line, Journal of Statistical Physics 182 (2021a).
  * Alberici  _et al._ [2021b] D. Alberici, F. Camilli, P. Contucci, and E. Mingione, The solution of the deep boltzmann machine on the nishimori line, Communications in Mathematical Physics 387 (2021b).
  * Bates and Sohn [2022] E. Bates and Y. Sohn, Free energy in multi-species mixed p-spin spherical models, Electronic Journal of Probability 27, 1 (2022).
  * Guionnet  _et al._ [2022] A. Guionnet, J. Ko, F. Krzakala, and L. Zdeborová, Low-rank matrix estimation with inhomogeneous noise, arXiv preprint, ariXiv:2208.05918 , arXiv:2208.05918 (2022).
  * Felstrom and Zigangirov [1999] A. J. Felstrom and K. S. Zigangirov, Time-varying periodic convolutional codes with low-density parity-check matrix, IEEE Transactions on Information Theory 45, 2181 (1999).
  * Kudekar  _et al._ [2011] S. Kudekar, T. Richardson, and R. Urbanke, Threshold saturation via spatial coupling: Why convolutional ldpc ensembles perform so well over the bec, IEEE Trans. Info. Th. 57, 803 (2011).
  * Barbier  _et al._ [2016a] J. Barbier, M. Dia, N. Macris, F. Krzakala, T. Lesieur, and L. Zdeborová, Mutual information for symmetric rank-one matrix estimation: A proof of the replica formula, in  _Advances in Neural Information Processing Systems_ (2016).
  * Marinari  _et al._ [1994a] E. Marinari, G. Parisi, and F. Ritort, Replica field theory for deterministic models: I. binary sequences with low autocorrelation, Journal of Physics A 27, 7615 (1994a).
  * Marinari  _et al._ [1994b] E. Marinari, G. Parisi, and F. Ritort, Replica field theory for deterministic models: II. A non-random spin glass with glassy behaviour, Journal of Physics A 27, 7647 (1994b).
  * Parisi and Potters [1999] G. Parisi and M. Potters, Mean-field equations for spin models with orthogonal interaction matrices, Journal of Physics A: Mathematical and General 28, 5267 (1999).
  * Opper  _et al._ [2016] M. Opper, B. Cakmak, and O. Winther, A theory of solving tap equations for ising models with general invariant random matrices, Journal of Physics A: Mathematical and Theoretical 49, 114002 (2016).
  * Opper and Winther [2011] M. Opper and O. Winther, Adaptive and self-averaging thouless-anderson-palmer mean-field theory for probabilistic modeling, Physical Review E 64 (2011).
  * Maillard  _et al._ [2019a] A. Maillard, L. Foini, A. L. Castellanos, F. Krzakala, M. Mézard, and L. Zdeborová, High-temperature expansions and message passing algorithms, Journal of Statistical Mechanics: Theory and Experiment 2019, 113301 (2019a).
  * Rangan  _et al._ [2019] S. Rangan, P. Schniter, and A. K. Fletcher, Vector approximate message passing, IEEE Transactions on Information Theory 65, 6664 (2019).
  * Gerbelot  _et al._ [2020] C. Gerbelot, A. Abbara, and F. Krzakala, Asymptotic errors for teacher-student convex generalized linear models (or: how to prove Kabashima’s replica formula), arXiv preprint, arXiv:2006.06581 (2020).
  * Fan [2022] Z. Fan, Approximate message passing algorithms for rotationally invariant matrices, The Annals of Statistics 50, 197 (2022).
  * Zhong  _et al._ [2021] X. Zhong, T. Wang, and Z. Fan, Approximate message passing for orthogonally invariant ensembles: Multivariate non-linearities and spectral initialization, arXiv:2110.02318 (2021).
  * Mondelli and Venkataramanan [2021] M. Mondelli and R. Venkataramanan, PCA initialization for approximate message passing in rotationally invariant models, in  _Advances in Neural Information Processing Systems_ , Vol. 34 (2021) pp. 29616–29629.
  * Venkataramanan  _et al._ [2022] R. Venkataramanan, K. Kögler, and M. Mondelli, Estimation in rotationally invariant generalized linear models via approximate message passing, in  _International Conference on Machine Learning_ (2022) pp. 22120–22144.
  * Barbier  _et al._ [2023] J. Barbier, F. Camilli, M. Mondelli, and M. Sáenz, Fundamental limits in structured principal component analysis and how to reach them, Proceedings of the National Academy of Sciences 120, e2302028120 (2023).
  * Dudeja  _et al._ [2024] R. Dudeja, S. Liu, and J. Ma, Optimality of approximate message passing algorithms for spiked matrix models with rotationally invariant noise, arXiv preprint arXiv:2405.18081 (2024).
  * Mézard  _et al._ [1990] M. Mézard, G. Parisi, and M.-A. Virasoro,  _Spin glass theory and beyond_ (World Scientific Publishing Co., Inc., Pergamon Press, 1990).
  * Guo  _et al._ [2005] D. Guo, S. Shamai, and S. Verdú, Mutual information and minimum mean-square error in gaussian channels, IEEE Transactions on Information Theory 51, 1261 (2005).
  * Mézard and Montanari [2009] M. Mézard and A. Montanari,  _Information, physics, and computation_ (Oxford Uni. Press, 2009).
  * Guerra [2003] F. Guerra, Broken replica symmetry bounds in the mean field spin glass model, Communications in mathematical physics 233, 1 (2003).
  * Talagrand [2006] M. Talagrand, The parisi formula, Annals of mathematics , 221 (2006).
  * Panchenko [2013b] D. Panchenko, The parisi ultrametricity conjecture, Annals of Mathematics , 383 (2013b).
  * Barbier  _et al._ [2016b] J. Barbier, M. Dia, and N. Macris, Proof of threshold saturation for spatially coupled sparse superposition codes, in  _IEEE International Symposium on Information Theory (ISIT)_ (2016) pp. 1173–1177.
  * El Alaoui and Krzakala [2018] A. El Alaoui and F. Krzakala, Estimation in the spiked wigner model: a short proof of the replica formula, in  _IEEE International Symposium on Information Theory (ISIT)_ (IEEE, 2018) pp. 1874–1878.
  * Barbier and Macris [2019] J. Barbier and N. Macris, The adaptive interpolation method for proving replica formulas. Applications to the Curie–Weiss and Wigner spike models, Journal of Physics A: Mathematical and Theoretical 52, 294002 (2019).
  * Barbier and Panchenko [2022] J. Barbier and D. Panchenko, Strong replica symmetry in high-dimensional optimal bayesian inference, Communications in Mathematical Physics 393, 1 (2022).
  * Barbier [2020] J. Barbier, Overlap matrix concentration in optimal Bayesian inference, Information and Inference: A Journal of the IMA 10, 597 (2020).
  * Alaoui  _et al._ [2020] A. E. Alaoui, F. Krzakala, and M. Jordan, Fundamental limits of detection in the spiked Wigner model, The Annals of Statistics 48, 863 (2020).
  * Barbier  _et al._ [2021] J. Barbier, D. Panchenko, and M. Sáenz, Strong replica symmetry for high-dimensional disordered log-concave Gibbs measures, Information and Inference: A Journal of the IMA 11, 1079 (2021).
  * Bodin and Macris [2021a] A. Bodin and N. Macris, Rank-one matrix estimation: analytic time evolution of gradient descent dynamics, arXiv preprint arXiv:2105.12257 (2021a).
  * Bodin and Macris [2021b] A. Bodin and N. Macris, Model, sample, and epoch-wise descents: exact solution of gradient flow in the random feature model, Advances in Neural Information Processing Systems 34, 21605 (2021b).
  * Bodin [2024] A. P. M. Bodin,  _Random matrix methods for high-dimensional machine learning models_ , Tech. Rep. (EPFL, 2024).
  * Bodin and Macris [2022] A. Bodin and N. Macris, Gradient flow in the gaussian covariate model: exact solution of learning curves and multiple descent structures, arXiv preprint arXiv:2212.06757 (2022).
  * Lelarge and Miolane [2019] M. Lelarge and L. Miolane, Fundamental limits of symmetric low-rank matrix estimation, Probability Theory and Related Fields 173, 859 (2019).
  * Maillard  _et al._ [2019b] A. Maillard, L. Foini, A. L. Castellanos, F. Krzakala, M. Mézard, and L. Zdeborová, High-temperature expansions and message passing algorithms, Journal of Statistical Mechanics: Theory and Experiment 2019, 113301 (2019b).
  * Thouless  _et al._ [1977] D. J. Thouless, P. W. Anderson, and R. G. Palmer, Solution of solvable model of a spin glass’, Philosophical Magazine 35, 593 (1977).
  * Mézard  _et al._ [1987] M. Mézard, G. Parisi, and M. A. Virasoro,  _Spin-Glass Theory and Beyond_ , Lecture Notes in Physics, Vol. 9 (World Scientific, Singapore, 1987).
  * Potters and Bouchaud [2020] M. Potters and J.-P. Bouchaud,  _A First Course in Random Matrix Theory: For Physicists, Engineers and Data Scientists_ (Cambridge University Press, 2020).