# Matrix Denoising with Doubly Heteroscedastic Noise: Fundamental Limits and
Optimal Spectral Methods

Yihan Zhang
Institute of Science and Technology Austria
<EMAIL_ADDRESS>
&Marco Mondelli
Institute of Science and Technology Austria
<EMAIL_ADDRESS>
###### Abstract

We study the matrix denoising problem of estimating the singular vectors of a
rank-$1$ signal corrupted by noise with both column and row correlations.
Existing works are either unable to pinpoint the exact asymptotic estimation
error or, when they do so, the resulting approaches (e.g., based on whitening
or singular value shrinkage) remain vastly suboptimal. On top of this, most of
the literature has focused on the special case of estimating the left singular
vector of the signal when the noise only possesses row correlation (one-sided
heteroscedasticity). In contrast, our work establishes the information-
theoretic and algorithmic limits of matrix denoising with doubly
heteroscedastic noise. We characterize the exact asymptotic minimum mean
square error, and design a novel spectral estimator with rigorous optimality
guarantees: under a technical condition, it attains positive correlation with
the signals whenever information-theoretically possible and, for one-sided
heteroscedasticity, it also achieves the Bayes-optimal error. Numerical
experiments demonstrate the significant advantage of our theoretically
principled method with the state of the art. The proofs draw connections with
statistical physics and approximate message passing, departing drastically
from standard random matrix theory techniques.

##  1 Introduction

Matrix denoising is a central primitive in statistics and machine learning,
and the problem is to recover a signal $X\in{\mathbb{R}}^{n\times d}$ from an
observation $A=X+W$ corrupted by additive noise $W$. This finds applications
across multiple domains of sciences, e.g., imaging [21, 60], biology [13, 42]
and astronomy [67, 5]. When $X$ has low rank and $W$ i.i.d. entries, $A$ is
the standard model for principal component analysis, typically referred to as
the Johnstone spiked covariance model [38]. When $n,d$ are both large and
proportional, which corresponds to the most sample-efficient regime, its
Bayes-optimal limits are well understood [48], and it has been established how
to achieve them efficiently [53]. Minimax/non-asymptotic guarantees are also
available in special cases, such as sparse PCA [17], Gaussian mixtures [69]
and certain joint scalings of $(n,d)$ [54].

However, in most applications, noise is highly structured and correlated,
thereby calling for more realistic assumptions on $W$ than having i.i.d.
entries. A recent line of work addresses this concern by studying matrix
denoising with heteroscedastic noise [1, 66, 29, 40, 23], resting on two basic
ideas: whitening and singular value shrinkage. Whitening refers to multiplying
the data matrix by the square root of the inverse covariance, in order to
reduce the model to one with i.i.d. noise; and singular value shrinkage
retains the singular vectors of the data while deflating the singular values
to correct for the noise. Though the exact asymptotic performance of these
algorithms has been derived [66, 29, 40, 23], their optimality is yet to be
determined from a Bayesian standpoint. In fact, we will prove that whitening
and shrinkage are _not_ the correct way to approach Bayes optimality.

#### Main contributions.

We focus on the prototypical model $A=X+W$, where
$X=\frac{\lambda}{n}u^{*}{v^{*}}^{\top}$ is a rank-$1$ signal, $\lambda$ is
the signal-to-noise ratio (SNR), and $W=\Xi^{1/2}\widetilde{W}\Sigma^{1/2}$ is
doubly heterogeneous noise. Here $u^{*},v^{*}$ follow i.i.d. priors;
$\widetilde{W}$ contains i.i.d. Gaussian entries; the covariance matrices
$\Xi,\Sigma$ capture column and row correlations; and we consider the typical
high-dimensional regime in which $n,d$ are both large and proportional. Our
main results are summarized below.

  1. 1.

We design an efficient spectral estimator to recover $u^{*},v^{*}$, and we
provide a precise asymptotic analysis of its performance, see Theorem 5.1.
This estimator is given by the top singular vectors of a matrix obtained by
carefully pre-processing $A$, see 5.3.

  2. 2.

When the priors of $u^{*},v^{*}$ are standard Gaussian, we show in Corollary
5.2 that the spectral estimator above is optimal in the following sense: _(i)_
under a technical condition, it achieves the _optimal weak recovery threshold_
, namely its mean square error is non-trivial as soon as this is information-
theoretically possible; _(ii)_ it achieves the _Bayes-optimal error_ for
$u^{*}$ (resp. $v^{*}$) when $\Xi$ (resp. $\Sigma$) is the identity. These
optimality guarantees follow from rigorously obtaining the asymptotic minimum
mean square error (MMSE) for the estimation of the whitened signals
$\Xi^{-1/2}u^{*}$ and $\Sigma^{-1/2}v^{*}$, see Theorem 4.2.

Our spectral estimator only involves matrix multiplication and computing
principal singular vectors. Practically, this can be efficiently done using
standard SVD algorithms or power iteration [44]. For both one-sided and double
heteroscedasticity, numerical experiments in Figures 2 and 3 show significant
advantage of our spectral estimator for moderate SNRs over HeteroPCA [72] and
shrinkage-based methods, i.e., Whiten-Shrink-reColor [40, 41], OptShrink [56],
and ScreeNOT [24].

#### Proof techniques.

We take a completely different route from classical approaches in statistics
and random matrix theory (e.g., whitening and shrinkage), and instead exploit
tools from statistical physics and the theory of approximate message passing.
In particular, the MMSE for the whitened signals
$\Xi^{-1/2}u^{*},\Sigma^{-1/2}v^{*}$ is obtained via an interpolation argument
[9, 48, 49]. This result allows us to derive the weak recovery threshold for
estimating the true signals $u^{*},v^{*}$. Moreover, for one-sided
heteroscedasticity, this MMSE coincides with that for estimating the true
signal on the homoscedastic side. Evaluating the Bayes-optimal estimators
requires solving high-dimensional integrals that are computationally
intractable. To circumvent this issue, we propose an efficient spectral method
that still enjoys optimality guarantees. Its design and analysis draw
connections with a family of iterative algorithms called approximate message
passing (AMP) [10, 28].

##  2 Related work

Research on matrix denoising in the homoscedastic case
($\Xi=I_{n},\Sigma=I_{d}$) has a rich history, and in random matrix theory
properties of the spectrum and eigenspaces of $A$ have been studied
exhaustively. Most prominently, the BBP phase transition phenomenon [4] (and
its finite-sample counterpart [57]) unveils a threshold of the SNR $\lambda$
above which a pair of outlier singular value and singular vector emerge. Under
i.i.d. priors, the asymptotic Bayes-optimal estimation error has been derived
[48, 49], rigorously justifying predictions from statistical physics [43]. The
proof uses the interpolation method due to Guerra [31], originally developed
in the context of mean-field spin glasses. Besides low-rank matrix estimation,
this method (including its adaptive variant [9] and the Aizenman–Sims–Starr
scheme [2]) has also been applied to a range of problems, including spiked
tensor estimation [45], generalized linear models [8], stochastic block models
[71] and group synchronization [70].

Moving to the heteroscedastic case, an active line of work concerns optimal
singular value shrinkage methods [40, 29, 41, 66, 56, 23]. These methods can
be regarded as a special family of rotationally invariant estimators, which
apply a univariate function $\eta\colon{\mathbb{R}}_{\geq 0}\to{\mathbb{R}}$
to each empirical singular value. An example widely employed by practitioners
is the thresholding function
$\eta_{\theta}(y)=y\mathds{1}{\left\\{y>\theta\right\\}}$ [24]. In the
presence of noise heteroscedasticity, most of these results are based on
whitening [39]. Another model of noise heterogeneity common in the literature
takes $W=\widetilde{W}\circ\Delta^{\circ 1/2}$, where $\widetilde{W}$ has
i.i.d. Gaussian entries, $\Delta$ is a deterministic block matrix with fixed
(i.e., constant with respect to $n,d$) number of blocks, and $\circ$ denotes
the element-wise product. This means that the entries of the noise are
independent but non-identically distributed, and they follow the variance
profile $\Delta$. The corresponding low-rank perturbation $A$, known as a
spiked inhomogeneous matrix, has attracted attention from both the
information-theoretic [11, 63, 32] and the algorithmic sides [34, 46, 58].
Spiked inhomogeneous matrices have some connections with the model considered
in this paper: if $\Delta$ has rank $1$, such $A$ can be realized by taking
$\Xi,\Sigma$ to be diagonal with suitable block structures. Finally, non-
asymptotic results for the heteroscedastic and the inhomogeneous models have
been derived in varying generality in [72, 78, 20, 1, 16]. We highlight that
our paper is the _first_ to establish information-theoretic and algorithmic
limits for doubly heteroscedastic noise.

Our characterization of the spectral estimator relies on an AMP algorithm that
converges to it by performing power iteration. AMP refers to a family of
iterative procedures, whose performance in the high-dimensional limit is
precisely characterized by a low-dimensional deterministic recursion called
state evolution [10, 14]. Originally introduced for compressed sensing [25],
AMP algorithms have been developed for various settings, including low-rank
estimation [53, 27, 6] and inference in generalized linear models [61, 62,
68]. Beyond statistical estimation, AMP proves its versatility as both an
efficient algorithm and a proof technique for studying e.g. posterior sampling
[55], spectral universality [26], first order methods with random data [19],
mismatched estimation [7], spectral estimators for generalized linear models
[75, 76] and their combination with linear estimators [50].

##  3 Problem setup

Consider the following rank-$1$ rectangular matrix estimation problem with
doubly heteroscedastic noise where we observe

$\displaystyle A$
$\displaystyle=\frac{\lambda}{n}u^{*}{v^{*}}^{\top}+W\in{\mathbb{R}}^{n\times
d},$ (3.1)

and aim to estimate $u^{*},v^{*}$. The following assumptions are imposed
throughout the paper. The dimensions $n,d\to\infty$ obey the proportional
scaling $n/d\to\delta\in(0,\infty)$, where $\delta$ is the aspect ratio. The
SNR $\lambda\in[0,\infty)$ is a known constant (relative to $n,d$). The
signals $(u^{*},v^{*})\sim P^{\otimes n}\otimes Q^{\otimes d}$ have i.i.d.
priors, where $P,Q$ are distributions on ${\mathbb{R}}$ with mean $0$ and
variance $1$. The unknown noise matrix has the form
$W=\Xi^{1/2}\widetilde{W}\Sigma^{1/2}\in{\mathbb{R}}^{n\times d}$, with
$\widetilde{W}_{i,j}\overset{\mathrm{i.i.d.}}{\sim}{\mathcal{N}}(0,1/n)$
independent of $(u^{*},v^{*})$. The covariances $\Xi\in{\mathbb{R}}^{n\times
n},\Sigma\in{\mathbb{R}}^{d\times d}$ are known, deterministic,111All our
results hold verbatim if $\Xi,\Sigma$ are random matrices independent of each
other and of $u^{*},v^{*},\widetilde{W}$.  strictly positive definite and
satisfy

$\displaystyle\lim_{n\to\infty}\frac{1}{n}\operatorname{Tr}(\Xi)$
$\displaystyle=\lim_{d\to\infty}\frac{1}{d}\operatorname{Tr}(\Sigma)=1.$ (3.2)

Their empirical spectral distributions (ESD) converge (as $n,d\to\infty$ s.t.
$n/d\to\infty$) weakly to the laws of the random variables $\overline{\Xi}$
and $\overline{\Sigma}$. Furthermore,
$\left\|\Xi\right\|_{2},\left\|\Sigma\right\|_{2}$ are uniformly bounded over
$d$. The supports of $\overline{\Xi},\overline{\Sigma}$ are compact subsets of
$(0,\infty)$. For all $\varepsilon>0$, there exists $d_{0}\in{\mathbb{N}}$
s.t. for all $d\geq d_{0}$,

$\displaystyle\operatorname{supp}(\operatorname{ESD}(\Xi))\subset\operatorname{supp}(\overline{\Xi})+[-\varepsilon,\varepsilon],\quad\operatorname{supp}(\operatorname{ESD}(\Sigma))\subset\operatorname{supp}(\overline{\Sigma})+[-\varepsilon,\varepsilon].$
(3.3)

The trace assumption 3.2 on the covariances is for normalization purposes
since the values of the traces, if not $1$, can be absorbed into $\lambda$.
The support assumption 3.3 excludes outliers in the spectra of covariances
which may contribute to undesirable spikes in $A$ [66].

##  4 Information-theoretic limits

In this section, we switch to an equivalent rescaled model

$\displaystyle Y$
$\displaystyle\coloneqq\sqrt{n}A=\sqrt{\frac{\gamma}{n}}u^{*}{v^{*}}^{\top}+\Xi^{1/2}Z\Sigma^{1/2}\in{\mathbb{R}}^{n\times
d},$ (4.1)

where $\gamma\coloneqq\lambda^{2}$ and $Z=\sqrt{n}\widetilde{W}$ contains
i.i.d. elements $Z_{i,j}\overset{\mathrm{i.i.d.}}{\sim}{\mathcal{N}}(0,1)$.
Abusing terminology, we refer to $\gamma$ as the SNR of $Y$. Define also
$\alpha\coloneqq 1/\delta\in(0,\infty)$ so that $d/n\to\alpha$. The scaling of
the parameters in 4.1 turns out to be more convenient for presenting the
results in this section. Results for $Y$ can be easily translated to $A$ by a
change of variables.

Let $\widetilde{u}^{*}\coloneqq\Xi^{-1/2}u^{*}$ and
$\widetilde{v}^{*}\coloneqq\Sigma^{-1/2}v^{*}$ denote the whitened signals.
The main result of this section is Theorem 4.2, which characterizes the
performance of the matrix minimum mean square error (MMSE) associated to the
estimation of
$\widetilde{u}^{*}(\widetilde{v}^{*})^{\top},\widetilde{u}^{*}(\widetilde{u}^{*})^{\top}$
and $\widetilde{v}^{*}(\widetilde{v}^{*})^{\top},$ via the corresponding
Bayes-optimal estimators:

$\displaystyle\operatorname{MMSE}_{n}(\gamma)$
$\displaystyle\coloneqq\frac{1}{nd}\mathbb{E}\left[\left\|\widetilde{u}^{*}(\widetilde{v}^{*})^{\top}-\mathbb{E}\left[\widetilde{u}^{*}(\widetilde{v}^{*})^{\top}\,|\,Y\right]\right\|_{\mathrm{F}}^{2}\right],$
(4.2) $\displaystyle\operatorname{MMSE}_{n}^{u}(\gamma)$
$\displaystyle\coloneqq\frac{1}{n^{2}}\mathbb{E}\left[\left\|\widetilde{u}^{*}(\widetilde{u}^{*})^{\top}-\mathbb{E}\left[\widetilde{u}(\widetilde{u}^{*})^{\top}\,|\,Y\right]\right\|_{\mathrm{F}}^{2}\right],$
(4.3) $\displaystyle\operatorname{MMSE}_{n}^{v}(\gamma)$
$\displaystyle\coloneqq\frac{1}{d^{2}}\mathbb{E}\left[\left\|\widetilde{v}^{*}(\widetilde{v}^{*})^{\top}-\mathbb{E}\left[\widetilde{v}^{*}(\widetilde{v}^{*})^{\top}\,|\,Y\right]\right\|_{\mathrm{F}}^{2}\right].$
(4.4)

Our characterization involves a pair of parameters
$(q_{u}^{*},q_{v}^{*})\in{\mathbb{R}}_{\geq 0}^{2}$ defined as the largest
solution to

$\displaystyle q_{u}$ $\displaystyle=\mathbb{E}\left[\frac{\alpha\gamma
q_{v}\overline{\Xi}^{-2}}{1+\alpha\gamma
q_{v}\overline{\Xi}^{-1}}\right],\qquad q_{v}=\mathbb{E}\left[\frac{\gamma
q_{u}\overline{\Sigma}^{-2}}{1+\gamma q_{u}\overline{\Sigma}^{-1}}\right].$
(4.5)

The proposition below, proved in Appendix A, justifies the existence of the
solution to 4.5 and identifies when a non-trivial solution emerges.

######  Proposition 4.1.

The fixed point equation 4.5 always has a trivial solution $(0,0)$. There
exists a non-trivial solution $(q_{u}^{*},q_{v}^{*})\in{\mathbb{R}}_{>0}^{2}$
if and only if

$\displaystyle\alpha\gamma^{2}\mathbb{E}\left[\overline{\Sigma}^{-2}\right]\mathbb{E}\left[\overline{\Xi}^{-2}\right]$
$\displaystyle>1,$ (4.6)

in which case the non-trivial solution is unique.

We are now ready to state our main result on the MMSE.

######  Theorem 4.2.

Assume $P=Q={\mathcal{N}}(0,1)$. For almost every $\gamma>0$,

$\displaystyle\begin{split}\lim_{n\to\infty}\operatorname{MMSE}_{n}(\gamma)&=\mathbb{E}\left[\overline{\Xi}^{-1}\right]\mathbb{E}\left[\overline{\Sigma}^{-1}\right]-q_{u}^{*}q_{v}^{*},\end{split}$
(4.7)
$\displaystyle\begin{split}\lim_{n\to\infty}\operatorname{MMSE}_{n}^{u}(\gamma)&=\mathbb{E}\left[\overline{\Xi}^{-1}\right]^{2}-{q_{u}^{*}}^{2},\qquad\lim_{n\to\infty}\operatorname{MMSE}_{n}^{u}(\gamma)=\mathbb{E}\left[\overline{\Sigma}^{-1}\right]^{2}-{q_{v}^{*}}^{2}.\end{split}$
(4.8)

We note that

$\lim_{n\to\infty}\frac{1}{nd}\mathbb{E}\left[\left\|\widetilde{u}^{*}(\widetilde{v}^{*})^{\top}\right\|_{\mathrm{F}}^{2}\right]=\lim_{n\to\infty}\frac{1}{nd}\mathbb{E}\left[\left\|\widetilde{u}^{*}\right\|_{2}^{2}\right]\mathbb{E}\left[\left\|\widetilde{v}^{*}\right\|_{2}^{2}\right]=\mathbb{E}\left[\overline{\Xi}^{-1}\right]\mathbb{E}\left[\overline{\Sigma}^{-1}\right],$
(4.9)

where the last step follows from Proposition G.2. This quantity represents the
trivial error in the estimation of
$\widetilde{u}^{*}(\widetilde{v}^{*})^{\top}$, which is achieved by the all-0
estimator. Analogous considerations hold for
$\widetilde{u}^{*}(\widetilde{u}^{*})^{\top}$ and
$\widetilde{v}^{*}(\widetilde{v}^{*})^{\top}$, for which the trivial
estimation error is $\mathbb{E}\left[\overline{\Xi}^{-1}\right]^{2}$ and
$\mathbb{E}\left[\overline{\Sigma}^{-1}\right]^{2}$, respectively. Thus,
Proposition 4.1 and Theorem 4.2 identify 4.6 as the condition for non-trivial
estimation, and the smallest $\gamma$ that satisfies 4.6 gives the _weak
recovery threshold_.

We show below that the weak recovery threshold is the same for the estimation
of the true signals $u^{*}{v^{*}}^{\top},u^{*}{u^{*}}^{\top}$ and
$v^{*}{v^{*}}^{\top}$. In this case, since the signal priors are Gaussian,
using the same passages as in 4.9 one has that the trivial estimation error
for $u^{*}{v^{*}}^{\top},u^{*}{u^{*}}^{\top}$ and $v^{*}{v^{*}}^{\top}$ is
always equal to $1$.

######  Corollary 4.3.

Assume $P=Q={\mathcal{N}}(0,1)$. The MMSE associated to the estimation of
$u^{*}{v^{*}}^{\top}$ is non-trivial, i.e,

$\displaystyle\lim_{n\to\infty}\frac{1}{nd}\mathbb{E}\left[\left\|u^{*}{v^{*}}^{\top}-\mathbb{E}\left[u^{*}{v^{*}}^{\top}\,|\,Y\right]\right\|_{\mathrm{F}}^{2}\right]<1$
(4.10)

if and only if 4.6 holds. The same result holds for the MMSE of
$u^{*}{u^{*}}^{\top}$ and $v^{*}{v^{*}}^{\top}$.

#### Proof strategy.

To derive the characterizations in Theorem 4.2, we write the posterior
distribution of $u^{*},v^{*}$ given $Y$ in a Gibbs form, i.e., its density is
the exponential of a Hamiltonian normalized by a partition function. The
interpolation argument relates the log-partition function (also referred to as
the ‘free energy’) of the posterior to that of the posteriors of two Gaussian
location models. Since i.i.d. Gaussianity is key to this approach, the
challenge is to handle noise covariances. Our idea is to incorporate the
covariances into the priors. In terms of the Hamiltonian, the model is
equivalent to the estimation of the whitened signals
$\Xi^{-1/2}u^{*},\Sigma^{-1/2}v^{*}$, whose priors have covariances, in the
presence of i.i.d. Gaussian noise. We then manage to carry out the
interpolation argument for the equivalent model and evaluate the free energy
of the corresponding Gaussian location models.

Specifically, let us starting by writing down the expression of the posterior
distribution after setting up some notation. For
$u\in{\mathbb{R}}^{n},v\in{\mathbb{R}}^{d}$, let
$\widetilde{u}\coloneqq\Xi^{-1/2}u,\widetilde{v}\coloneqq\Sigma^{-1/2}v$.
Define the densities

$\displaystyle\,\mathrm{d}\widetilde{P}(\widetilde{u})\coloneqq\sqrt{\det(\Xi)}\,\mathrm{d}P^{\otimes
n}(\Xi^{1/2}\widetilde{u}),\qquad\,\mathrm{d}\widetilde{Q}(\widetilde{v})\coloneqq\sqrt{\det(\Sigma)}\,\mathrm{d}Q^{\otimes
d}(\Sigma^{1/2}\widetilde{v}),$

where the determinant factors ensure that the integrals equal $1$. With
$P=Q={\mathcal{N}}(0,1)$, we have
$\widetilde{P}={\mathcal{N}}(0_{n},\Xi^{-1}),\widetilde{Q}={\mathcal{N}}(0_{d},\Sigma^{-1})$,
and from Bayes’ rule the posterior of $(u^{*},v^{*})$ given $Y$ is

$\displaystyle\,\mathrm{d}P(u,v\,|\,Y)$
$\displaystyle=\frac{1}{{\mathcal{Z}}_{n}(\gamma)}\exp\left({\mathcal{H}}_{n}(\Xi^{-1/2}u,\Sigma^{-1/2}v)\right)\,\mathrm{d}P^{\otimes
n}(u)\,\mathrm{d}Q^{\otimes d}(v),$ (4.11)

where the Hamiltonian and the partition function are given respectively by

$\displaystyle{\mathcal{H}}_{n}(\widetilde{u},\widetilde{v})\coloneqq\sqrt{\frac{\gamma}{n}}\widetilde{u}^{\top}Z\widetilde{v}+\frac{\gamma}{n}\widetilde{u}^{\top}\widetilde{u}^{*}\widetilde{v}^{\top}\widetilde{v}^{*}-\frac{\gamma}{2n}\left\|\widetilde{u}\right\|_{2}^{2}\left\|\widetilde{v}\right\|_{2}^{2},$
(4.12)
$\displaystyle{\mathcal{Z}}_{n}(\gamma)\coloneqq\int\int\exp\left({\mathcal{H}}_{n}(\Xi^{-1/2}u,\Sigma^{-1/2}v)\right)\,\mathrm{d}P^{\otimes
n}(u)\,\mathrm{d}Q^{\otimes
d}(v)=\int\int\exp\left({\mathcal{H}}_{n}(\widetilde{u},\widetilde{v})\right)\,\mathrm{d}\widetilde{P}(\widetilde{u})\,\mathrm{d}\widetilde{Q}(\widetilde{v}).$
(4.13)

Define the free energy as

$\displaystyle{\mathcal{F}}_{n}(\gamma)$
$\displaystyle\coloneqq\frac{1}{n}\mathbb{E}\left[\log{\mathcal{Z}}_{n}(\gamma)\right].$
(4.14)

The major technical step is to characterize ${\mathcal{F}}_{n}(\gamma)$ in the
large $n$ limit in terms of a bivariate functional ${\mathcal{F}}$ introduced
below. This is the core component to derive the MMSE characterization.

For a positive random variable $\overline{\Sigma}$ subject to the conditions
in Section 3, let

$\displaystyle\psi_{\overline{\Sigma}}(\gamma)$
$\displaystyle\coloneqq\frac{1}{2}\left(\gamma\mathbb{E}\left[\overline{\Sigma}^{-1}\right]-\mathbb{E}\left[\log\left(1+\gamma\overline{\Sigma}^{-1}\right)\right]\right).$
(4.15)

As shown in Appendix B, $\psi_{\overline{\Sigma}}(\gamma)$ is the limiting
free energy of a Gaussian channel, in which one wishes to estimate
$x^{*}\in{\mathbb{R}}^{n}$ from the observation
$Y=\sqrt{\gamma}x^{*}+\Sigma^{1/2}Z$ corrupted by anisotropic Gaussian noise
with covariance $\Sigma$. Using 4.15, let us define the replica symmetric
potential ${\mathcal{F}}$:

$\displaystyle{\mathcal{F}}(q_{u},q_{v})$
$\displaystyle\coloneqq\psi_{\overline{\Xi}}(\alpha\gamma
q_{v})+\alpha\psi_{\overline{\Sigma}}(\gamma
q_{u})-\frac{\alpha\gamma}{2}q_{u}q_{v},$

and the set of critical points of ${\mathcal{F}}$:

$\displaystyle{\mathcal{C}}(\gamma,\alpha)$
$\displaystyle\coloneqq\left\\{(q_{u},q_{v})\in{\mathbb{R}}_{\geq
0}^{2}:\partial_{1}{\mathcal{F}}(q_{u},q_{v})=0,\partial_{2}{\mathcal{F}}(q_{u},q_{v})=0\right\\}$
$\displaystyle=\left\\{(q_{u},q_{v})\in{\mathbb{R}}_{\geq
0}^{2}:q_{u}=2\psi_{\overline{\Xi}}^{\prime}(\alpha\gamma
q_{v}),q_{v}=2\psi_{\overline{\Sigma}}^{\prime}(\gamma q_{u})\right\\}$ (4.16)
$\displaystyle=\left\\{(q_{u},q_{v})\in{\mathbb{R}}_{\geq
0}^{2}:(q_{u},q_{v})\textnormal{ solves
}\lx@cref{creftypecap~refnum}{eqn:fp_it}\right\\},$

where the last equality is a direct calculation of
$\psi_{\overline{\Xi}}^{\prime},\psi_{\overline{\Sigma}}^{\prime}$. The
following result, proved in Appendix C, shows that the limit of
${\mathcal{F}}_{n}(\gamma)$ is given by a dimension-free variational problem
involving ${\mathcal{F}}(q_{u},q_{v})$.

######  Theorem 4.4 (Free energy).

Assume $P=Q={\mathcal{N}}(0,1)$. Then, we have

$\displaystyle\lim_{n\to\infty}{\mathcal{F}}_{n}(\gamma)$
$\displaystyle=\sup_{q_{v}\geq 0}\inf_{q_{u}\geq
0}{\mathcal{F}}(q_{u},q_{v})=\sup_{(q_{u},q_{v})\in{\mathcal{C}}(\gamma,\alpha)}{\mathcal{F}}(q_{u},q_{v}),$

and $\sup_{q_{v}}\inf_{q_{u}}$ and $\sup_{(q_{u},q_{v})}$ are achieved by the
same $(q_{u}^{*},q_{v}^{*})$ in Proposition 4.1.

######  Remark 4.1 (Equivalent models).

Informally, the above result says that the matrix model 4.1 is equivalent at
the level of Hamiltonian to the following two statistically uncorrelated
vector models:

$\displaystyle Y^{u}$ $\displaystyle\coloneqq\sqrt{\alpha\gamma
q_{v}^{*}}u^{*}+\Xi^{1/2}Z^{u}\in{\mathbb{R}}^{n},\quad
Y^{v}\coloneqq\sqrt{\gamma
q_{u}^{*}}v^{*}+\Sigma^{1/2}Z^{v}\in{\mathbb{R}}^{d},$ (4.17)

with $q_{u}^{*},q_{v}^{*}$ the largest solution to 4.5 and
$(u^{*},v^{*},Z^{u},Z^{v})\sim P^{\otimes n}\otimes Q^{\otimes
d}\otimes{\mathcal{N}}(0_{n},I_{n})\otimes{\mathcal{N}}(0_{d},I_{d}).$

######  Remark 4.2 (Gaussian priors).

Theorem 4.4 crucially relies on having Gaussian priors $P,Q$. This assumption
is mainly used to derive single-letter (i.e., dimension-free) expressions of
the free energy of the vector models in 4.17 which, under Gaussian priors, are
nothing but Gaussian integrals. The free energy, and hence the MMSE, are
expected to be sensitive to the priors. Indeed, this is already the case in
the homoscedastic setting $\Xi=I_{n},\Sigma=I_{d}$ [48]. An extension towards
general i.i.d. priors is a challenging open problem and, in fact, without
posing additional assumptions on $\Xi,\Sigma$, it is unclear whether a single-
letter expression for free energy and MMSE is possible.

At this point, the MMSE can be derived from the above characterization of free
energy. Indeed, let

$\displaystyle{\mathcal{D}}(\alpha)$
$\displaystyle\coloneqq\left\\{\gamma>0:{\mathcal{F}}\textnormal{ has a unique
maximizer }(q_{u}^{*},q_{v}^{*})\textnormal{ over
}{\mathcal{C}}(\gamma,\alpha)\right\\}.$

The envelope theorem [47, Corollary 4] ensures that ${\mathcal{D}}(\alpha)$ is
equal to ${\mathbb{R}}_{>0}$ up to a countable set. Using algebraic relations
between free energy and MMSE, we prove 4.7 and 4.8 for all
$\gamma\in{\mathcal{D}}(\alpha)$ (and, thus, for almost every $\gamma>0$).
Then, using the Nishimori identity and the fact that the ESDs of $\Xi,\Sigma$
are upper and lower bounded by constants independent of $n$ and $d$, Corollary
4.3 also follows. The formal arguments are contained in Appendix D.

##  5 Spectral estimator

This section introduces a spectral estimator that meets the weak recovery
threshold and, for one-sided heteroscedasticity, attains the Bayes-optimal
error. Suppose that the following condition holds

$\displaystyle\frac{\lambda^{4}}{\delta}\mathbb{E}\left[\overline{\Sigma}^{-2}\right]\mathbb{E}\left[\overline{\Xi}^{-2}\right]$
$\displaystyle>1,$ (5.1)

which is equivalent to 4.6. Under this condition, the fixed point equations
4.5 have a unique pair of positive solutions $(q_{u}^{*},q_{v}^{*})$. For
convenience, we also define the rescalings $\mu^{*}\coloneqq\lambda
q_{v}^{*}/\delta,\nu^{*}\coloneqq\lambda q_{u}^{*}$, and the auxiliary
quantities

$b^{*}\coloneqq\frac{1}{\delta}\mathbb{E}\left[\frac{\lambda}{\lambda\nu^{*}+\overline{\Sigma}}\right],\qquad
c^{*}\coloneqq\mathbb{E}\left[\frac{\lambda}{\lambda\mu^{*}+\overline{\Xi}}\right].$
(5.2)

Now, we pre-process the data matrix $A$ as

$\displaystyle
A^{*}\coloneqq\lambda(\lambda(\mu^{*}+b^{*})I_{n}+\Xi)^{-1/2}\Xi^{-1/2}A\Sigma^{-1/2}(\lambda(\nu^{*}+c^{*})I_{d}+\Sigma)^{-1/2},$
(5.3)

from which we obtain the spectral estimators

$\displaystyle\begin{split}\widehat{u}&\coloneqq\eta_{u}\sqrt{n}\frac{\Xi^{1/2}(\lambda(\mu^{*}+b^{*})I_{n}+\Xi)^{-1/2}(\lambda\mu^{*}I_{n}+\Xi)u_{1}(A^{*})}{\left\|\Xi^{1/2}(\lambda(\mu^{*}+b^{*})I_{n}+\Xi)^{-1/2}(\lambda\mu^{*}I_{n}+\Xi)u_{1}(A^{*})\right\|_{2}},\end{split}$
(5.4a)
$\displaystyle\begin{split}\widehat{v}&\coloneqq\eta_{v}\sqrt{d}\frac{\Sigma^{1/2}(\lambda(\nu^{*}+c^{*})I_{d}+\Sigma)^{-1/2}(\lambda\nu^{*}I_{d}+\Sigma)v_{1}(A^{*})}{\left\|\Sigma^{1/2}(\lambda(\nu^{*}+c^{*})I_{d}+\Sigma)^{-1/2}(\lambda\nu^{*}I_{d}+\Sigma)v_{1}(A^{*})\right\|_{2}},\end{split}$
(5.4b)

where $u_{1}(\cdot)/v_{1}(\cdot)$ denote the top left/right singular vectors
and

$\displaystyle\begin{split}\eta_{u}&\coloneqq\sqrt{\frac{\lambda\mu^{*}}{\lambda\mu^{*}+1}},\qquad\eta_{v}\coloneqq\sqrt{\frac{\lambda\nu^{*}}{\lambda\nu^{*}+1}}.\end{split}$
(5.5)

Note that $\eta_{u},\eta_{v}>0$, provided that 5.1 holds. The pre-processing
of $A$ in 5.3 and the form of the spectral estimators in 5.4 come from the
derivation of a suitable AMP algorithm, and they are discussed at the end of
the section. We finally defer to Section E.3 the definition of the scalar
quantity $\sigma_{2}^{*}$ obtained via a fixed point equation depending only
on $\overline{\Xi},\overline{\Sigma},\lambda,\delta$, see E.26 for details.

Our main result, Theorem 5.1, shows that, under the criticality condition 5.1,
the matrix $A^{*}$ exhibits a spectral gap between the top two singular
values, and it characterizes the performance of the spectral estimators in
5.4, proving that they achieve weak recovery of $u^{*}$ and $v^{*}$,
respectively.

######  Theorem 5.1.

Suppose that 5.1 holds and that, for any $c>0$,

$\displaystyle\lim_{\beta\downarrow
s}\mathbb{E}\left[\frac{\overline{\Sigma}^{*}}{\beta-c\overline{\Sigma}^{*}}\right]=\lim_{\beta\downarrow
s}\mathbb{E}\left[\left(\frac{\overline{\Sigma}^{*}}{\beta-c\overline{\Sigma}^{*}}\right)^{2}\right]=\infty,\quad\lim_{\alpha\downarrow\sup\operatorname{supp}(\overline{\Xi}^{*})}\mathbb{E}\left[\frac{\overline{\Xi}^{*}}{\alpha-\overline{\Xi}^{*}}\right]=\infty,$
(5.6)

where
$\overline{\Xi}^{*}\coloneqq\frac{\lambda}{\lambda(\mu^{*}+b^{*})+\overline{\Xi}}$,
$\overline{\Sigma}^{*}\coloneqq\frac{\lambda}{\lambda(\nu^{*}+c^{*})+\overline{\Sigma}}$
and $s\coloneqq c\cdot\sup\operatorname{supp}(\overline{\Sigma}^{*})$. Let
$A^{*},\widehat{u},\widehat{v},\sigma_{2}^{*}$ be defined in 5.3, 5.4 and
E.26, and $\sigma_{i}(A^{*})$ denote the $i$-th largest singular value of
$A^{*}$. Then, if $\sigma_{2}^{*}<1$, the following limits hold in
probability:

$\displaystyle\begin{split}\lim_{n\to\infty}\sigma_{1}(A^{*})=1>\sigma_{2}^{*}=\lim_{n\to\infty}\sigma_{2}(A^{*}),\end{split}$
(5.7)
$\displaystyle\begin{split}\lim_{n\to\infty}\frac{\left|\left\langle\widehat{u},u^{*}\right\rangle\right|}{\left\|\widehat{u}\right\|_{2}\left\|u^{*}\right\|_{2}}=\eta_{u},\quad\lim_{d\to\infty}\frac{\left|\left\langle\widehat{v},v^{*}\right\rangle\right|}{\left\|\widehat{v}\right\|_{2}\left\|v^{*}\right\|_{2}}=\eta_{v}\end{split}$
(5.8)
$\displaystyle\begin{split}\lim_{n\to\infty}\frac{1}{n^{2}}\left\|u^{*}{u^{*}}^{\top}-\widehat{u}\widehat{u}^{\top}\right\|_{\mathrm{F}}^{2}=1-\eta_{u}^{4},\quad\lim_{d\to\infty}\frac{1}{d^{2}}\left\|v^{*}{v^{*}}^{\top}-\widehat{v}\widehat{v}^{\top}\right\|_{\mathrm{F}}^{2}=1-\eta_{v}^{4},\end{split}$
(5.9)
$\displaystyle\begin{split}\lim_{n\to\infty}\frac{1}{nd}\left\|u^{*}{v^{*}}^{\top}-\widehat{u}\widehat{v}^{\top}\right\|_{\mathrm{F}}^{2}=1-\eta_{u}^{2}\eta_{v}^{2}.\end{split}$
(5.10)

######  Remark 5.1 (Assumptions).

To guarantee a spectral gap for $A^{*}$ and the weak recoverability of
$u^{*},v^{*}$ via the proposed spectral method, we also require the algebraic
condition $\sigma_{2}^{*}<1$. We conjecture that this condition is implied by
5.1, and we have verified that this is the case in all our numerical
experiments (see Figure 1 for two concrete examples). The additional
assumption 5.6 is a mild regularity condition on the covariances. It ensures
that the densities of $\overline{\Xi}^{*},\overline{\Sigma}^{*}$ decay
sufficiently slowly at the edges of the support, so that $\sigma_{2}^{*}$ is
well-posed [75].

(a) $\Xi=I_{n}$ and $\Sigma$ a Toeplitz matrix with $\rho=0.9$.

(b) $\Xi$ a circulant matrix with $c=0.1,\ell=5$ and $\Sigma$ a Toeplitz
matrix with $\rho=0.5$.

Figure 1: Top two singular values of $A^{*}$ in 5.3, where $d=4000,\delta=4$
and each simulation is averaged over $10$ i.i.d. trials. The singular values
computed experimentally (‘sim’ in the legends and $\times$ in the plots)
closely match our theoretical prediction in 5.7 (‘thy’ in the legends and
solid curves with the same color in the plots). The threshold $\lambda^{*}$ is
such that equality holds in 5.1. We note that the green curve corresponding to
$\sigma_{2}^{*}$ is smaller than $1$ for $\lambda>\lambda^{*}$, i.e., when 5.1
holds.

######  Remark 5.2 (Signal priors).

Theorem 5.1 does not require the prior distributions $P,Q$ to be Gaussian, and
it is valid for any i.i.d. prior with mean $0$ and variance $1$.

On the one hand, Corollary 4.3 shows that, if 5.1 is violated, the problem is
information-theoretically impossible, i.e., no estimator achieves non-trivial
error. On the other hand, Theorem 5.1 exhibits a pair of estimators that
achieves non-trivial error as soon as 5.1 holds – under the additional
assumption $\sigma_{2}^{*}<1$ which we conjecture to be equivalent. Thus, the
spectral method in 5.4 is optimal in terms of weak recovery threshold. Though
such estimators do not attain the optimal error, when both priors are Gaussian
and $\Xi=I_{n}$, $\widehat{u}\widehat{u}^{\top}$ is the Bayes-optimal estimate
for $u^{*}{u^{*}}^{\top}$.

######  Corollary 5.2.

Assume $P=Q={\mathcal{N}}(0,1)$, and consider the setting of Theorem 5.1 with
the additional assumption $\Xi=I_{n}$. Then, $\eta_{u}=\sqrt{q_{u}^{*}}$,
i.e., $\widehat{u}\widehat{u}^{\top}$ achieves the MMSE for
$u^{*}{u^{*}}^{\top}$.

The claim readily follows by noting that, when $\Xi=I_{n}$, the first equation
in 4.5 becomes

$\displaystyle q_{u}^{*}$ $\displaystyle=\frac{\alpha\gamma
q_{v}^{*}}{1+\alpha\gamma
q_{v}^{*}}=\frac{(\lambda^{2}/\delta)(\delta\mu^{*}/\lambda)}{1+(\lambda^{2}/\delta)(\delta\mu^{*}/\lambda)}=\frac{\lambda\mu^{*}}{1+\lambda\mu^{*}}=\eta_{u}^{2},$

where the last equality is by the definition 5.5 of $\eta_{u}$. Let us
highlight that, even if $\Xi=I_{n}$, $\widehat{u}$ still makes non-trivial use
of the other covariance $\Sigma^{1/2}$. At the information-theoretic level,
this is reflected by the fact that $\Sigma^{1/2}$ enters $q_{u}^{*}$ through
the fixed point equations 4.5. Therefore, even though the matrix model in 4.1
is equivalent to a pair of uncorrelated vector models in 4.17 in the sense of
the free energy, the tasks of estimating $u^{*}$ and $v^{*}$ cannot be
decoupled.

(a) Normalized correlation with $u^{*}$

(b) Normalized correlation with $v^{*}$

(c) Matrix MSE for $u^{*}{v^{*}}^{\top}$

Figure 2: Performance comparison when $\Xi=I_{n}$ and $\Sigma$ is a circulant
matrix. The numerical results closely follow the predictions of Theorem 5.1,
and our spectral estimators in 5.4 outperform all other methods (Leeb–Romanov,
OptShrink, ScreeNOT, and HeteroPCA), especially at low SNR.

(a) Normalized correlation with $u^{*}$

(b) Normalized correlation with $v^{*}$

(c) Matrix MSE for $u^{*}{v^{*}}^{\top}$

Figure 3: Performance comparison when $\Xi$ is a Toeplitz matrix and $\Sigma$
is circulant. The numerical results closely follow the predictions of Theorem
5.1, and our spectral estimators in 5.4 outperform all other methods (Leeb,
OptShrink, and ScreeNOT), especially at low SNR.

#### Numerical experiments.

Figures 2 and 3 demonstrate the advantage of our method over existing
approaches, and they display an accurate agreement between simulations (‘sim’
in the legends and $\times$ in the plots) and the theoretical predictions of
Theorem 5.1 (‘thy’ in the legends and solid curves with the same color in the
plots), both plotted as a function of $\lambda$. In both figures,
$n=4000,d=2000$ (so $\delta=2$), and $P=Q={\mathcal{N}}(0,1)$. Each data point
is computed from $20$ i.i.d. trials and error bars are reported at $1$
standard deviation. We let $\Xi$ be either the identity or a Toeplitz matrix
[73, 37, 18], i.e., $\Xi_{i,j}=\rho^{\left|i-j\right|}$ with $\rho=0.9$. We
let $\Sigma$ be a circulant matrix [36, 35]: the first row has $1$ in the
first position, $c=0.0078$ in the second through $(\ell+1)$-st position and in
the last $\ell$ positions ($\ell=300$), with the remaining entries being $0$;
for $2\leq i\leq d$, the $i$-th row is a cyclic shift of the $(i-1)$-st row to
the right by $1$ position. Both matrices satisfy 5.6 and the conditions of
Section 3.

Our spectral estimator outperforms all other approaches: Leeb–Romanov [40],
OptShrink [56], ScreeNOT [24], and HeteroPCA [72] in the one-sided
heteroscedastic case (Figure 2); Leeb [41], OptShrink, and ScreeNOT in the
doubly heteroscedastic case (Figure 3). When computing the normalized
correlation with the signals (left/right overlap), the performance of
Leeb–Romanov and Leeb is the same as the estimators
$\Xi^{1/2}u_{1}(\Xi^{-1/2}A\Sigma^{-1/2}),\Sigma^{1/2}v_{1}(\Xi^{-1/2}A\Sigma^{-1/2})$
referred to as ‘whiten’ in Figures 2(a) and 2(b); the performance of OptShrink
and ScreeNOT is the same as the estimators $u_{1}(A),v_{1}(A)$ referred to as
‘vanilla’ in Figures 3(a) and 3(b). The advantage of our approach (in black)
is especially significant at low SNR; as SNR increases, Leeb-Romanov and Leeb
(in red) achieve similar performance; a much larger SNR ($>2$ and $>3$ in
Figures 2 and 3) is required by HeteroPCA, OptShrink and ScreeNOT (in magenta,
blue and green) to perform comparably.

#### Proof strategy.

The design and analysis of the spectral estimator in 5.4 comprise two steps,
detailed in Appendix E. The _first step_ is to present an AMP algorithm dubbed
Bayes-AMP for matrix denoising with doubly heteroscedastic noise.
Specifically, its iterates are updated as

$\displaystyle u^{t}$
$\displaystyle=\Xi^{-1}A\Sigma^{-1}\widetilde{v}^{t}-b_{t}\Xi^{-1}\widetilde{u}^{t-1},\quad\widetilde{u}^{t}=g_{t}^{*}(u^{t}),\quad
c_{t}=\frac{1}{n}\operatorname{Tr}((\nabla g_{t}^{*}(u^{t}))\Xi^{-1}),$ (5.11)
$\displaystyle v^{t+1}$
$\displaystyle=\Sigma^{-1}A^{\top}\Xi^{-1}\widetilde{u}^{t}-c_{t}\Sigma^{-1}\widetilde{v}^{t},\quad\widetilde{v}^{t+1}=f_{t+1}^{*}(v^{t+1}),\quad
b_{t+1}=\frac{1}{n}\operatorname{Tr}((\nabla
f_{t+1}^{*}(v^{t+1}))\Sigma^{-1}),$

where $\nabla$ denotes the Jacobian matrix, and the functions
$g_{t}^{*},f_{t+1}^{*}$ are specified below in 5.12. As common in AMP
algorithms, the iterates 5.11 are accompanied with a state evolution which
accurately tracks their behavior via a simple deterministic recursion: the
joint empirical distribution of $(u^{*},v^{*},u^{t},v^{t+1})$ converges to the
random variables $(U^{*},V^{*},U_{t},V_{t+1})$, see Proposition E.1 for a
formal statement and the recursive description of the laws of such random
variables. Then, the name ‘Bayes-AMP’ is motivated by the fact that
$g_{t}^{*},f_{t+1}^{*}$ are the posterior-mean denoisers given by

$\displaystyle g_{t}^{*}(u)$
$\displaystyle\coloneqq\mathbb{E}\left[U^{*}\,|\,U_{t}=u\right],\qquad
f_{t+1}^{*}(v)\coloneqq\mathbb{E}\left[V^{*}\,|\,V_{t+1}=v\right].$ (5.12)

Remarkably, Bayes-AMP operates on $\Xi^{-1}A\Sigma^{-1}$, as opposed to the
widely adopted ansatz of considering the whitened matrix
$\Xi^{-1/2}A\Sigma^{-1/2}$. The advantage of operating on
$\Xi^{-1}A\Sigma^{-1}$ is that the fixed point of the corresponding state
evolution matches the extremizers of the free energy in 4.5. This would _not_
be the case if Bayes-AMP used the whitening $\Xi^{-1/2}A\Sigma^{-1/2}$.

The design of Bayes-AMP and the proof of its state evolution follow a two-step
reduction detailed in Appendix F. Using a change of variables, we show in
Section F.2 that Bayes-AMP can be realized by an auxiliary AMP with non-
separable denoising functions (meaning that $g_{t},f_{t+1}$ cannot be written
as univariate functions applied component-wise) operating on
$\Xi^{-1/2}A\Sigma^{-1/2}=\frac{\lambda}{n}\widetilde{u}^{*}(\widetilde{v}^{*})^{\top}+\widetilde{W}$.
Then, in Section F.1 we simulate the auxiliary AMP using a standard AMP
operating on the i.i.d. Gaussian matrix $\widetilde{W}$, whose state evolution
has been established in [12, 30].

However, Bayes-AMP by itself is not a practical algorithm since it needs a
warm start, i.e., an initialization that achieves non-trivial error. Thus, the
_second step_ is to design a spectral estimator that solves the fixed point
equation of Bayes-AMP, which turns out to be an eigen-equation for $A^{*}$.

We now heuristically derive the form 5.3 of $A^{*}$ and the expression 5.4 of
the spectral estimator. To do so, we note that the large-$n$ limits of
$c_{t},b_{t+1}$ coincide with the auxiliary quantities $c^{*},b^{*}$ defined
in 5.2. Furthermore, when the priors of $u^{*},v^{*}$ are Gaussian, 5.12
reduces to

$g_{t}^{*}(u)=\lambda(\lambda\mu^{*}\Xi^{-1}+I_{n})^{-1}u,\qquad
f_{t+1}^{*}(v)=\lambda(\lambda\nu^{*}\Sigma^{-1}+I_{d})^{-1}v,$

where we recall that $\mu^{*}=\lambda q_{v}^{*}/\delta$ and $\nu^{*}=\lambda
q_{u}^{*}$ are rescalings of the non-trivial solution $(q_{u}^{*},q_{v}^{*})$
of 4.5. Denoting by $u,v$ the fixed points of the iteration 5.11, after some
manipulations we have

$\mathfrak{g}(\Xi)u=A^{*}\mathfrak{f}(\Sigma)v,\qquad\mathfrak{f}(\Sigma)v={A^{*}}^{\top}\mathfrak{g}(\Xi)u,$

where $A^{*}$ is given in 5.3 and

$\displaystyle\mathfrak{g}(\Xi)$
$\displaystyle\coloneqq\sqrt{\lambda}(\lambda(\mu^{*}+b^{*})I_{n}+\Xi)^{1/2}(\lambda\mu^{*}I_{n}+\Xi)^{-1}\Xi^{1/2},$
$\displaystyle\mathfrak{f}(\Sigma)$
$\displaystyle\coloneqq\sqrt{\lambda}(\lambda(\nu^{*}+c^{*})I_{d}+\Sigma)^{1/2}(\lambda\nu^{*}I_{d}+\Sigma)^{-1}\Sigma^{1/2}.$

This suggests that $A^{*}$ has top singular value equal to $1$ and
$(\mathfrak{g}(\Xi)u,\mathfrak{f}(\Sigma)v)$ are aligned with the
corresponding singular vectors $(u_{1}(A^{*}),v_{1}(A^{*}))$. Moreover, state
evolution implies that the distribution of the fixed point $(u,v)$ is close to
that of

$(\mu^{*}\Xi^{-1}u^{*}+\sqrt{\mu^{*}/\lambda}w_{u},\nu^{*}\Sigma^{-1}v^{*}+\sqrt{\nu^{*}/\lambda}w_{v}),$

with
$(w_{u},w_{v})\sim{\mathcal{N}}(0_{n},\Xi^{-1})\otimes{\mathcal{N}}(0_{d},\Sigma^{-1})$
independent of $u^{*},v^{*}$. Thus, to obtain estimates of $(u^{*},v^{*})$, we
take
$(\Xi\mathfrak{g}(\Xi)^{-1}u_{1}(A^{*}),\Sigma\mathfrak{f}(\Sigma)^{-1}v_{1}(A^{*}))$
and suitably rescale their norm, which leads to the expressions in 5.4. More
details on the above heuristics are discussed in Section E.2.

The most outstanding step remains to make the heuristics rigorous. This
involves proving that $\Xi u^{t},\Sigma v^{t+1}$ are aligned with the proposed
spectral estimator, which allows for a performance characterization via state
evolution. The formal argument is carried out in Section E.4.

##  6 Concluding remarks

In this work, we establish information-theoretic limits and propose an
efficient spectral method with optimality guarantees, for matrix estimation
with doubly heteroscedastic noise. On the one hand, under Gaussian priors, we
give a rigorous characterization of the MMSE; on the other hand, we present a
spectral estimator that _(i)_ achieves the information-theoretic weak recovery
threshold, and _(ii)_ is Bayes-optimal for the estimation of one of the
signals, when the noise is heteroscedastic only on the other side. While our
analysis focuses on rank-1 estimation, we expect that all results admit proper
extensions to rank-$r$ signals, where $r$ is a constant independent on $n,d$.

The design and analysis of the spectral estimator draws connections with
approximate message passing and, along the way, we introduce a Bayes-AMP
algorithm which could be of independent interest. In this paper, we employ
Bayes-AMP solely as a proof technique. However, one could use the spectral
method designed here as an initialization of Bayes-AMP itself, after suitably
correcting its iterates. This strategy has been successfully carried out for
i.i.d. Gaussian noise in [53] and for rotationally invariant noise in [51,
77]. Bayes-AMP is well equipped to exploit signal priors more informative than
the Gaussian one, and AMP algorithms are known to achieve the information-
theoretically optimal estimation error for low-rank matrix inference [53, 6].
Nevertheless, we point out two obstacles towards doing so in the presence of
doubly heteroscedastic noise. First, for general priors, establishing the
information-theoretic limits remains a challenging open problem, and it is
unclear whether a low-dimensional characterization of the free energy (and,
hence, of the MMSE) is possible. Second, even for Gaussian priors, Bayes-AMP
reduces to the proposed spectral estimator, which is not Bayes-optimal for the
general case of doubly heteroscedastic noise.

Finally, the proposed spectral estimator makes non-trivial use of the
covariances $\Xi,\Sigma$, which are assumed to be known. However, when $n$ and
$d$ grow proportionally, such matrices – if unknown – cannot be consistently
estimated from the data. Thus, a challenging open problem is to construct
estimators that retain comparable performance without knowing the noise
covariances. The paper [52] takes a step in this direction by developing
methods that achieve the minimax risk when the noise is i.i.d. with an unknown
non-Gaussian distribution.

## Acknowledgments and Disclosure of Funding

YZ thanks Shashank Vatedka for discussions at the early stage of this project.
MM thanks Jean Barbier for sharing his insights into the interpolation
argument.

This research is partially supported by the 2019 Lopez-Loreta Prize and by the
Interdisciplinary Projects Committee (IPC) at the Institute of Science and
Technology Austria (ISTA).

## References

  * [1] Joshua Agterberg, Zachary Lubberts, and Carey E. Priebe.  Entrywise estimation of singular vectors of low-rank matrices with heteroskedasticity and dependence.  IEEE Trans. Inform. Theory, 68(7):4618–4650, 2022.
  * [2] Michael Aizenman, Robert Sims, and Shannon L. Starr.  Extended variational principle for the sherrington-kirkpatrick spin-glass model.  Phys. Rev. B, 68:214403, Dec 2003.
  * [3] Z. D. Bai and Y. Q. Yin.  Limit of the smallest eigenvalue of a large-dimensional sample covariance matrix.  Ann. Probab., 21(3):1275–1294, 1993.
  * [4] Jinho Baik, Gérard Ben Arous, and Sandrine Péché.  Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices.  Ann. Probab., 33(5):1643–1697, 2005.
  * [5] Stephen Bailey.  Principal component analysis with noisy and/or missing data.  Publications of the Astronomical Society of the Pacific, 124(919):1015, sep 2012.
  * [6] Jean Barbier, Francesco Camilli, Marco Mondelli, and Manuel Sáenz.  Fundamental limits in structured principal component analysis and how to reach them.  Proc. Natl. Acad. Sci. USA, 120(30):Paper No. e2302028120, 7, 2023\.
  * [7] Jean Barbier, TianQi Hou, Marco Mondelli, and Manuel Saenz.  The price of ignorance: how much does it cost to forget noise structure in low-rank matrix estimation?  In Advances in Neural Information Processing Systems, volume 35, pages 36733–36747, 2022.
  * [8] Jean Barbier, Florent Krzakala, Nicolas Macris, Léo Miolane, and Lenka Zdeborová.  Optimal errors and phase transitions in high-dimensional generalized linear models.  Proc. Natl. Acad. Sci. USA, 116(12):5451–5460, 2019.
  * [9] Jean Barbier and Nicolas Macris.  The adaptive interpolation method: a simple scheme to prove replica formulas in Bayesian inference.  Probab. Theory Related Fields, 174(3-4):1133–1185, 2019.
  * [10] Mohsen Bayati and Andrea Montanari.  The dynamics of message passing on dense graphs, with applications to compressed sensing.  IEEE Trans. Inform. Theory, 57(2):764–785, 2011.
  * [11] Joshua K. Behne and Galen Reeves.  Fundamental limits for rank-one matrix estimation with groupwise heteroskedasticity.  In International Conference on Artificial Intelligence and Statistics, pages 8650–8672, 2022.
  * [12] Raphaël Berthier, Andrea Montanari, and Phan-Minh Nguyen.  State evolution for approximate message passing with non-separable functions.  Inf. Inference, 9(1):33–79, 2020.
  * [13] Tejal Bhamre, Teng Zhang, and Amit Singer.  Denoising and covariance estimation of single particle cryo-em images.  Journal of Structural Biology, 195(1):72–81, 2016.
  * [14] Erwin Bolthausen.  An iterative construction of solutions of the TAP equations for the Sherrington-Kirkpatrick model.  Comm. Math. Phys., 325(1):333–366, 2014.
  * [15] Stéphane Boucheron, Gábor Lugosi, and Pascal Massart.  Concentration inequalities.  Oxford University Press, Oxford, 2013.
  * [16] T. Tony Cai, Rungang Han, and Anru R. Zhang.  On the non-asymptotic concentration of heteroskedastic Wishart-type matrix.  Electron. J. Probab., 27:Paper No. 29, 40, 2022.
  * [17] T. Tony Cai, Zongming Ma, and Yihong Wu.  Sparse PCA: optimal rates and adaptive estimation.  Ann. Statist., 41(6):3074–3110, 2013.
  * [18] T. Tony Cai, Zhao Ren, and Harrison H. Zhou.  Optimal rates of convergence for estimating Toeplitz covariance matrices.  Probab. Theory Related Fields, 156(1-2):101–143, 2013.
  * [19] Michael Celentano, Chen Cheng, and Andrea Montanari.  The high-dimensional asymptotics of first order methods with random data.  arXiv preprint arXiv:2112.07572, 2021.
  * [20] Chen Cheng, Yuting Wei, and Yuxin Chen.  Tackling small eigen-gaps: fine-grained eigenvector estimation and inference under heteroscedastic noise.  IEEE Trans. Inform. Theory, 67(11):7380–7419, 2021.
  * [21] Lucilio Cordero-Grande, Daan Christiaens, Jana Hutter, Anthony N. Price, and Jo V. Hajnal.  Complex diffusion-weighted image estimation via matrix recovery under general noise models.  NeuroImage, 200:391–404, 2019.
  * [22] Romain Couillet and Walid Hachem.  Analysis of the limiting spectral measure of large random matrices of the separable covariance type.  Random Matrices Theory Appl., 3(4):1450016, 23, 2014.
  * [23] Xiucai Ding, Yun Li, and Fan Yang.  Eigenvector distributions and optimal shrinkage estimators for large covariance and precision matrices.  arXiv preprint arXiv:2404.14751, 2024.
  * [24] David Donoho, Matan Gavish, and Elad Romanov.  ScreeNOT: exact MSE-optimal singular value thresholding in correlated noise.  Ann. Statist., 51(1):122–148, 2023.
  * [25] David L. Donoho, Arian Maleki, and Andrea Montanari.  Message passing algorithms for compressed sensing.  Proceedings of the National Academy of Sciences, 106:18914–18919, 2009.
  * [26] Rishabh Dudeja, Subhabrata Sen, and Yue M Lu.  Spectral universality of regularized linear regression with nearly deterministic sensing matrices.  IEEE Transactions on Information Theory, 2024.
  * [27] Zhou Fan.  Approximate message passing algorithms for rotationally invariant matrices.  The Annals of Statistics, 50(1):197–224, 2022.
  * [28] Oliver Y Feng, Ramji Venkataramanan, Cynthia Rush, Richard J Samworth, et al.  A unifying tutorial on approximate message passing.  Foundations and Trends® in Machine Learning, 15(4):335–536, 2022.
  * [29] Matan Gavish, William Leeb, and Elad Romanov.  Matrix denoising with partial noise statistics: optimal singular value shrinkage of spiked F-matrices.  Inf. Inference, 12(3):Paper No. iaad028, 46, 2023.
  * [30] Cédric Gerbelot and Raphaël Berthier.  Graph-based approximate message passing iterations.  Inf. Inference, 12(4):Paper No. iaad020, 67, 2023.
  * [31] Francesco Guerra.  Broken replica symmetry bounds in the mean field spin glass model.  Comm. Math. Phys., 233(1):1–12, 2003.
  * [32] Alice Guionnet, Justin Ko, Florent Krzakala, and Lenka Zdeborová.  Low-rank matrix estimation with inhomogeneous noise.  arXiv preprint arXiv:2208.05918, 2022.
  * [33] Philip Hartman.  Ordinary differential equations, volume 38.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2002.
  * [34] David Hong, Fan Yang, Jeffrey A. Fessler, and Laura Balzano.  Optimally weighted PCA for high-dimensional heteroscedastic data.  SIAM J. Math. Data Sci., 5(1):222–250, 2023.
  * [35] Adel Javanmard and Andrea Montanari.  Confidence intervals and hypothesis testing for high-dimensional regression.  J. Mach. Learn. Res., 15:2869–2909, 2014.
  * [36] Adel Javanmard and Andrea Montanari.  Hypothesis testing in high-dimensional regression under the Gaussian random design model: asymptotic theory.  IEEE Trans. Inform. Theory, 60(10):6522–6554, 2014.
  * [37] Adel Javanmard and Andrea Montanari.  Debiasing the Lasso: optimal sample size for Gaussian designs.  Ann. Statist., 46(6A):2593–2622, 2018.
  * [38] Iain M. Johnstone.  On the distribution of the largest eigenvalue in principal components analysis.  Ann. Statist., 29(2):295–327, 2001.
  * [39] Boris Landa, Thomas T. C. K. Zhang, and Yuval Kluger.  Biwhitening reveals the rank of a count matrix.  SIAM J. Math. Data Sci., 4(4):1420–1446, 2022.
  * [40] William Leeb and Elad Romanov.  Optimal spectral shrinkage and PCA with heteroscedastic noise.  IEEE Trans. Inform. Theory, 67(5):3009–3037, 2021.
  * [41] William E. Leeb.  Matrix denoising for weighted loss functions and heterogeneous signals.  SIAM J. Math. Data Sci., 3(3):987–1012, 2021.
  * [42] Jeffrey T. Leek.  Asymptotic conditional singular value decomposition for high-dimensional genomic data.  Biometrics, 67(2):344–352, 2011.
  * [43] Thibault Lesieur, Florent Krzakala, and Lenka Zdeborová.  Mmse of probabilistic low-rank matrix estimation: Universality with respect to the output channel.  In 2015 53rd Annual Allerton Conference on Communication, Control, and Computing (Allerton), pages 680–687, 2015.
  * [44] Yue M. Lu.  Householder dice: a matrix-free algorithm for simulating dynamics on Gaussian and random orthogonal ensembles.  IEEE Trans. Inform. Theory, 67(12):8264–8272, 2021.
  * [45] Clément Luneau, Jean Barbier, and Nicolas Macris.  Mutual information for low-rank even-order symmetric tensor estimation.  Inf. Inference, 10(4):1167–1207, 2021.
  * [46] Pierre Mergny, Justin Ko, and Florent Krzakala.  Spectral phase transition and optimal pca in block-structured spiked models.  arXiv preprint arXiv:2403.03695, 2024.
  * [47] Paul Milgrom and Ilya Segal.  Envelope theorems for arbitrary choice sets.  Econometrica, 70(2):583–601, 2002.
  * [48] Léo Miolane.  Fundamental limits of low-rank matrix estimation: the non-symmetric case.  arXiv preprint arXiv:1702.00473, 2017.
  * [49] Léo Miolane.  Fundamental limits of inference: A statistical physics approach. Theses, Ecole normale supérieure - ENS PARIS ; Inria Paris, June 2019.
  * [50] Marco Mondelli, Christos Thrampoulidis, and Ramji Venkataramanan.  Optimal combination of linear and spectral estimators for generalized linear models.  Foundations of Computational Mathematics, pages 1–54, 2021.
  * [51] Marco Mondelli and Ramji Venkataramanan.  Pca initialization for approximate message passing in rotationally invariant models.  In Advances in Neural Information Processing Systems, volume 34, pages 29616–29629, 2021.
  * [52] Andrea Montanari, Feng Ruan, and Jun Yan.  Adapting to unknown noise distribution in matrix denoising.  arXiv preprint arXiv:1810.02954, 2018.
  * [53] Andrea Montanari and Ramji Venkataramanan.  Estimation of low-rank matrices via approximate message passing.  Ann. Statist., 49(1):321–345, 2021.
  * [54] Andrea Montanari and Yuchen Wu.  Fundamental limits of low-rank matrix estimation with diverging aspect ratios.  arXiv preprint arXiv:2211.00488, 2022.
  * [55] Andrea Montanari and Yuchen Wu.  Posterior sampling from the spiked models via diffusion processes.  arXiv preprint arXiv:2304.11449, 2023.
  * [56] Raj Rao Nadakuditi.  OptShrink: an algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage.  IEEE Trans. Inform. Theory, 60(5):3002–3018, 2014.
  * [57] Boaz Nadler.  Finite sample approximation results for principal component analysis: a matrix perturbation approach.  Ann. Statist., 36(6):2791–2817, 2008.
  * [58] Aleksandr Pak, Justin Ko, and Florent Krzakala.  Optimal algorithms for the inhomogeneous spiked wigner model.  In Advances in Neural Information Processing Systems, volume 36, pages 76409–76424, 2023.
  * [59] Dmitry Panchenko.  The Sherrington-Kirkpatrick model.  Springer Monographs in Mathematics. Springer, New York, 2013.
  * [60] Henrik Pedersen, Sebastian Kozerke, Steffen Ringgaard, Kay Nehrke, and Won Yong Kim.  k-t pca: Temporally constrained k-t blast reconstruction using principal component analysis.  Magnetic Resonance in Medicine, 62(3):706–716, 2009.
  * [61] S. Rangan.  Generalized approximate message passing for estimation with random linear mixing.  In IEEE International Symposium on Information Theory (ISIT), 2011\.
  * [62] Sundeep Rangan, Philip Schniter, and Alyson K Fletcher.  Vector approximate message passing.  IEEE Transactions on Information Theory, 65(10):6664–6684, 2019\.
  * [63] Galen Reeves.  Information-theoretic limits for the matrix tensor product.  IEEE Journal on Selected Areas in Information Theory, 1(3):777–798, 2020.
  * [64] R. Tyrrell Rockafellar.  Convex analysis.  Princeton Landmarks in Mathematics. Princeton University Press, Princeton, NJ, 1997.
  * [65] Charles M. Stein.  Estimation of the mean of a multivariate normal distribution.  The Annals of Statistics, 9(6):1135–1151, 1981.
  * [66] Pei-Chun Su and Hau-Tieng Wu.  Data-driven optimal shrinkage of singular values under high-dimensional noise with separable covariance structure.  arXiv preprint arXiv:2207.03466, 2022.
  * [67] O. Tamuz, T. Mazeh, and S. Zucker.  Correcting systematic effects in a large set of photometric light curves.  Monthly Notices of the Royal Astronomical Society, 356(4):1466–1470, 02 2005.
  * [68] Ramji Venkataramanan, Kevin Kögler, and Marco Mondelli.  Estimation in rotationally invariant generalized linear models via approximate message passing.  In International Conference on Machine Learning (ICML), 2022.
  * [69] Yihong Wu and Harrison H. Zhou.  Randomly initialized EM algorithm for two-component Gaussian mixture achieves near optimality in $O(\sqrt{n})$ iterations.  Math. Stat. Learn., 4(3-4):143–220, 2021.
  * [70] Kaylee Y Yang, Timothy LH Wee, and Zhou Fan.  Asymptotic mutual information in quadratic estimation problems over compact groups.  arXiv preprint arXiv:2404.10169, 2024.
  * [71] Xiaodong Yang, Buyu Lin, and Subhabrata Sen.  Fundamental limits of community detection from multi-view data: multi-layer, dynamic and partially labeled block models.  arXiv preprint arXiv:2401.08167, 2024.
  * [72] Anru R. Zhang, T. Tony Cai, and Yihong Wu.  Heteroskedastic PCA: algorithm, optimality, and applications.  Ann. Statist., 50(1):53–80, 2022.
  * [73] Cun-Hui Zhang and Stephanie S. Zhang.  Confidence intervals for low dimensional parameters in high dimensional linear models.  J. R. Stat. Soc. Ser. B. Stat. Methodol., 76(1):217–242, 2014.
  * [74] Lixin Zhang.  Spectral analysis of large dimentional random matrices.  PhD thesis, National University of Singapore, 2007.
  * [75] Yihan Zhang, Hong Chang Ji, Ramji Venkataramanan, and Marco Mondelli.  Spectral estimators for structured generalized linear models via approximate message passing.  arXiv preprint arXiv:2308.14507, 2023.
  * [76] Yihan Zhang, Marco Mondelli, and Ramji Venkataramanan.  Precise asymptotics for spectral methods in mixed generalized linear models.  arXiv preprint arXiv:2211.11368, 2022.
  * [77] Xinyi Zhong, Tianhao Wang, and Zhou Fan.  Approximate message passing for orthogonally invariant ensembles: Multivariate non-linearities and spectral initialization.  arXiv preprint arXiv:2110.02318, 2021.
  * [78] Yuchen Zhou and Yuxin Chen.  Deflated heteropca: Overcoming the curse of ill-conditioning in heteroskedastic pca.  arXiv preprint arXiv:2303.06198, 2023.

#### Notation.

All vectors are column vectors. The singular values of a matrix
$A\in{\mathbb{R}}^{n\times d}$ (where $n\geq d$ without loss of generality)
are denoted by $\sigma_{1}(A)\geq\cdots\geq\sigma_{d}(A)\geq 0$ and the
corresponding left/right singular vectors are denoted by
$u_{1}(A),\cdots,u_{d}(A)\in{\mathbb{S}}^{n-1}$ and
$v_{1}(A),\cdots,v_{d}(A)\in{\mathbb{S}}^{d-1}$. The (real) eigenvalues of a
symmetric matrix $B\in{\mathbb{R}}^{d\times d}$ are denoted by
$\lambda_{1}(B)\geq\cdots\geq\lambda_{d}(B)$ and the corresponding
eigenvectors are denoted by $v_{1}(B),\cdots,v_{d}(B)\in{\mathbb{S}}^{d-1}$
(which will not be confused with the right singular vectors, whenever they are
different, since we will never talk about both simultaneously for a square
asymmetric matrix). We generally put overlines on capital letters to indicate
a scalar random variable, e.g., $\overline{X}\in{\mathbb{R}}$, whose support
is denoted by $\operatorname{supp}(\overline{X})$. The limit/liminf/limsup in
probability are denoted by
$\operatorname*{p-lim},\operatorname*{p-liminf},\operatorname*{p-limsup}$. The
product distribution whose $i$-th ($i\in[k]$) marginal is given by $P_{i}$ is
denoted by $P_{1}\otimes\cdots\otimes P_{k}$, with the shorthand $P^{\otimes
k}$ when all $P_{i}$’s are equal to $P$. The gradient of
$f\colon{\mathbb{R}}^{n}\to{\mathbb{R}}$, or with abuse of notation, the
Jacobian matrix of $F\colon{\mathbb{R}}^{n}\to{\mathbb{R}}^{d}$ are denoted by
$\nabla f\in{\mathbb{R}}^{n},\nabla F\in{\mathbb{R}}^{d\times n}$. The partial
derivative of $f(x_{1},\cdots,x_{n})$ with respect to $x_{i}$ is denoted by
either $\frac{\partial}{\partial x_{i}}f(x_{1},\cdots,x_{n})$ or
$\partial_{i}f(x_{1},\cdots,x_{n})$. All $\log$ and $\exp$ are to the base
$e$. We use the standard notation of $\sup(S),\inf(S)$ for a subset
$S\subset{\mathbb{R}}$. We generally use $C>0$ to denote a sufficiently large
constant independent of $n,d$. Its dependence on other parameters will be
specified, though its value may change across passages. We use the standard
big O notation.

##  Appendix A Proof of Proposition 4.1

We eliminate $q_{u}$ and write a fixed point equation only involving $q_{v}$:

$\displaystyle q_{v}$
$\displaystyle=\mathbb{E}\left[\frac{\gamma\mathbb{E}\left[\frac{\alpha\gamma
q_{v}\overline{\Xi}^{-2}}{\alpha\gamma
q_{v}\overline{\Xi}^{-1}+1}\right]\overline{\Sigma}^{-2}}{\gamma\mathbb{E}\left[\frac{\alpha\gamma
q_{v}\overline{\Xi}^{-2}}{\alpha\gamma
q_{v}\overline{\Xi}^{-1}+1}\right]\overline{\Sigma}^{-1}+1}\right].$

Denote the RHS by $f(q_{v})$. Recall that we are only interested in non-
negative solutions $(q_{u},q_{v})$. So let us restrict attention on $f$ to the
domain ${\mathbb{R}}_{\geq 0}$. We have $f(0)=0$ and

$\displaystyle f^{\prime}(q_{u})$
$\displaystyle=\mathbb{E}\left[\frac{\gamma\overline{\Sigma}^{-2}}{\left(\gamma\mathbb{E}\left[\frac{\alpha\gamma
q_{v}\overline{\Xi}^{-2}}{\alpha\gamma
q_{v}\overline{\Xi}^{-1}+1}\right]\overline{\Sigma}^{-1}+1\right)^{2}}\right]\mathbb{E}\left[\frac{\alpha\gamma\overline{\Xi}^{-2}}{\left(\alpha\gamma
q_{v}\overline{\Xi}^{-1}+1\right)^{2}}\right]>0,$ $\displaystyle
f^{\prime}(0)$
$\displaystyle=\alpha\gamma^{2}\mathbb{E}\left[\overline{\Sigma}^{-2}\right]\mathbb{E}\left[\overline{\Xi}^{-2}\right],$
$\displaystyle f^{\prime\prime}(q_{v})$
$\displaystyle=-2\left(\mathbb{E}\left[\frac{\alpha\gamma\overline{\Xi}^{-2}}{\left(\alpha\gamma
q_{v}\overline{\Xi}^{-1}+1\right)^{2}}\right]^{2}\mathbb{E}\left[\frac{\gamma^{2}\overline{\Sigma}^{-3}}{\left(\gamma\mathbb{E}\left[\frac{\alpha\gamma
q_{v}\overline{\Xi}^{-2}}{\alpha\gamma
q_{v}\overline{\Xi}^{-1}+1}\right]\overline{\Sigma}^{-1}+1\right)^{3}}\right]\right.$
$\displaystyle\qquad\qquad\left.+\mathbb{E}\left[\frac{\gamma\overline{\Sigma}^{-2}}{\left(\gamma\mathbb{E}\left[\frac{\alpha\gamma
q_{v}\overline{\Xi}^{-2}}{\alpha\gamma
q_{v}\overline{\Xi}^{-1}+1}\right]\overline{\Sigma}^{-1}+1\right)^{2}}\right]\mathbb{E}\left[\frac{\alpha^{2}\gamma^{2}\overline{\Xi}^{-3}}{\left(\alpha\gamma
q_{v}\overline{\Xi}^{-1}+1\right)^{3}}\right]\right)<0,$
$\displaystyle\lim_{q_{v}\to\infty}f(q_{v})$
$\displaystyle=\mathbb{E}\left[\frac{\gamma\mathbb{E}\left[\overline{\Xi}^{-1}\right]\overline{\Sigma}^{-2}}{\gamma\mathbb{E}\left[\overline{\Xi}^{-1}\right]\overline{\Sigma}^{-1}+1}\right]\in(0,\infty).$

It then becomes evident that a non-trivial fixed point $q_{v}>0$ exists if and
only if $f^{\prime}(0)>1$ and in this case, the non-trivial fixed point is
unique.

Finally, by the first equation in 4.5, there is a non-trivial fixed point
$q_{u}$ if and only if there is a non-trivial fixed point $q_{v}$, which
completes the proof.

##  Appendix B Auxiliary Gaussian channel

We formally introduce here the auxiliary model mentioned in Section 4.
Consider a Gaussian channel with blocklength $n$, input $x^{*}$, output $Y$,
anisotropic Gaussian noise $\Sigma^{1/2}Z$ and SNR $\gamma$:

$\displaystyle Y$
$\displaystyle=\sqrt{\gamma}x^{*}+\Sigma^{1/2}Z\in{\mathbb{R}}^{n},$ (B.1)

where

$\displaystyle(x^{*},Z)$ $\displaystyle\sim P^{\otimes
n}\otimes{\mathcal{N}}(0_{n},I_{n}).$

By similar derivations as in Section 4, the posterior distribution of $x^{*}$
given $Y$ can be written as

$\displaystyle\,\mathrm{d}P(x\,|\,Y)$
$\displaystyle=\frac{1}{Z_{n}(\gamma)}\exp\left(H_{n}(x)\right)\,\mathrm{d}P^{\otimes
n}(x),$

where the Hamiltonian and the partition function are

$\displaystyle H_{n}(x)$
$\displaystyle\coloneqq\gamma{x^{*}}^{\top}\Sigma^{-1}x+\sqrt{\gamma}Z^{\top}\Sigma^{-1/2}x-\frac{\gamma}{2}x^{\top}\Sigma^{-1}x,$
$\displaystyle Z_{n}(\gamma)$
$\displaystyle\coloneqq\int_{{\mathbb{R}}^{n}}\exp\left(H_{n}(x)\right)\,\mathrm{d}P^{\otimes
n}(x).$

Define the free energy as

$\displaystyle F_{n}(\gamma)$
$\displaystyle\coloneqq\frac{1}{n}\mathbb{E}\left[\log Z_{n}(\gamma)\right].$

With $P={\mathcal{N}}(0,1)$, $Z_{n}(\gamma)$ becomes a Gaussian integral that
can be computed as below using Proposition G.1:

$\displaystyle Z_{n}(\gamma)$
$\displaystyle=\frac{1}{\sqrt{\det(\gamma\Sigma^{-1}+I_{n})}}\exp\left(\frac{1}{2}\left(\gamma\Sigma^{-1}x^{*}+\sqrt{\gamma}\Sigma^{-1/2}Z\right)^{\top}\left(\gamma\Sigma^{-1}+I_{n}\right)^{-1}\left(\gamma\Sigma^{-1}x^{*}+\sqrt{\gamma}\Sigma^{-1/2}Z\right)\right).$

Therefore, by Proposition G.2,

$\displaystyle\operatorname*{p-lim}_{n\to\infty}F_{n}(\gamma)$
$\displaystyle=-\frac{1}{2}\mathbb{E}\left[\log\left(\gamma\overline{\Sigma}^{-1}+1\right)\right]+\frac{1}{2}\gamma^{2}\mathbb{E}\left[\overline{\Sigma}^{-2}\left(\gamma\overline{\Sigma}^{-1}+1\right)^{-1}\right]+\frac{1}{2}\gamma\mathbb{E}\left[\overline{\Sigma}^{-1}\left(\gamma\overline{\Sigma}^{-1}+1\right)^{-1}\right]$
$\displaystyle=\frac{1}{2}\left(\gamma\mathbb{E}\left[\overline{\Sigma}^{-1}\right]-\mathbb{E}\left[\log\left(1+\gamma\overline{\Sigma}^{-1}\right)\right]\right).$
(B.2)

The above functional is nothing but $\psi_{\overline{\Sigma}}(\gamma)$
introduced in 4.15 which will play an important role in characterizing the
free energy of the original model 4.1.

##  Appendix C Proof of Theorem 4.4

Before diving into the proof, we make further notation adjustments for the
ease of applying the interpolation argument. Specifically, we will henceforth
assume $\gamma=1$ by incorporating the actual value of $\gamma$ into the prior
distributions $P,Q$,

$\displaystyle\int_{\mathbb{R}}x^{2}\,\mathrm{d}P(x)$
$\displaystyle=\gamma,\quad\int_{\mathbb{R}}x^{2}\,\mathrm{d}Q(x)=1.$

This is obviously equivalent to the previous setting. So we can drop the
dependence on $\gamma$ and write
$\operatorname{MMSE}_{n},{\mathcal{Z}}_{n},{\mathcal{F}}_{n}$ for
$\operatorname{MMSE}_{n}(\gamma),{\mathcal{Z}}_{n}(\gamma),{\mathcal{F}}_{n}(\gamma)$
defined in 4.2, 4.13 and 4.14.

We will also assume that $\Xi,\Sigma$ are diagonal. This is without loss of
generality since the Gaussianity of $P,Q,\widetilde{W}$ ensures that both the
prior distributions and the noise matrix are rotationally invariant.
Furthermore, we truncate $P,Q$ so that they are supported on $[-K,K]$ for a
constant $K>0$. The approximation error in the free energy due to truncation
can be made arbitrarily small if $K$ is sufficiently large, since the free
energy is pseudo-Lipschitz in the prior distribution with respect to the
Wasserstein-$2$ metric.

The proof follows an interpolation argument [9, 48, 49] with suitable
modifications to take care of the noise heteroscedasticity featured by the
covariances $\Xi,\Sigma$. To start with, define the interpolating models:

$\displaystyle Y_{t}$
$\displaystyle\coloneqq\sqrt{\frac{1-t}{n}}{u}^{*}{v^{*}}^{\top}+\Xi^{1/2}Z\Sigma^{1/2}\in{\mathbb{R}}^{n\times
d},$ $\displaystyle Y_{t}^{u}$ $\displaystyle\coloneqq\sqrt{\alpha
q_{1}(t)}{u}^{*}+\Xi^{1/2}Z^{u}\in{\mathbb{R}}^{n},$ $\displaystyle Y_{t}^{v}$
$\displaystyle\coloneqq\sqrt{q_{2}(t)}{v}^{*}+\Sigma^{1/2}Z^{v}\in{\mathbb{R}}^{d},$

where $q_{1}(t),q_{2}(t)\geq 0$ are to be determined and

$\displaystyle(u^{*},v^{*},Z,Z^{u},Z^{v})$ $\displaystyle\sim P^{\otimes
n}\otimes Q^{\otimes
d}\otimes{\mathcal{N}}(0_{nd},I_{nd})\otimes{\mathcal{N}}(0_{n},I_{n})\otimes{\mathcal{N}}(0_{d},I_{d}).$
(C.1)

By definition, $Y_{0}=Y$ is the model that we would like to understand, and
$Y_{t}^{u},Y_{t}^{v}$ are instances of Gaussian channels in B.1 whose free
energy we have already understood (see B.2). The idea is that $Y_{t}$ serves
as a path parametrized by $t\in[0,1]$ from the original model $Y$ to the
target models $(Y_{1}^{u},Y_{1}^{v})$. The crux of the interpolation argument
lies in showing that $Y_{t}$ and $(Y_{t}^{u},Y_{t}^{v})$ are equivalent (at
the level of free energy) along the path.

To study the interpolating models $(Y_{t},Y_{t}^{u},Y_{t}^{v})$, define the
Hamiltonian

$\displaystyle\begin{split}{\mathcal{H}}_{n,t}(\widetilde{u},\widetilde{v};q_{1},q_{2})&\coloneqq\sqrt{\frac{1-t}{n}}\widetilde{u}^{\top}Z\widetilde{v}+\frac{1-t}{n}\widetilde{u}^{\top}\widetilde{u}^{*}\widetilde{v}^{\top}\widetilde{v}^{*}-\frac{1-t}{2n}\left\|\widetilde{u}\right\|_{2}^{2}\left\|\widetilde{v}\right\|_{2}^{2}\\\
&\phantom{=}~{}+\alpha q_{1}\widetilde{u}^{\top}\widetilde{u}^{*}+\sqrt{\alpha
q_{1}}\widetilde{u}^{\top}Z^{u}-\frac{\alpha
q_{1}}{2}\left\|\widetilde{u}\right\|_{2}^{2}\\\
&\phantom{=}~{}+q_{2}\widetilde{v}^{\top}\widetilde{v}^{*}+\sqrt{q_{2}}\widetilde{v}^{\top}Z^{v}-\frac{q_{2}}{2}\left\|\widetilde{v}\right\|_{2}^{2}.\end{split}$
(C.2)

Then the posterior distribution of $({u}^{*},{v}^{*})$ given
$(Y_{t},Y_{t}^{u},Y_{t}^{v})$ is

$\displaystyle\,\mathrm{d}P({u},{v}\,|\,Y_{t},Y_{t}^{u},Y_{t}^{v})$
$\displaystyle=\frac{1}{{\mathcal{Z}}_{n,t}}\exp\left({\mathcal{H}}_{n,t}(\Xi^{-1/2}u,\Sigma^{-1/2}v;q_{1}(t),q_{2}(t))\right)\,\mathrm{d}P^{\otimes
n}(u)\,\mathrm{d}Q^{\otimes d}(v).$ (C.3)

Let the partition function be

$\displaystyle{\mathcal{Z}}_{n,t}$
$\displaystyle\coloneqq\int_{{\mathbb{R}}^{d}}\int_{{\mathbb{R}}^{n}}\exp\left({\mathcal{H}}_{n,t}(\Xi^{-1/2}u,\Sigma^{-1/2}v;q_{1}(t),q_{2}(t))\right)\,\mathrm{d}P^{\otimes
n}(u)\,\mathrm{d}Q^{\otimes d}(v)$
$\displaystyle=\int_{{\mathbb{R}}^{d}}\int_{{\mathbb{R}}^{n}}\exp\left({\mathcal{H}}_{n,t}(\widetilde{u},\widetilde{v};q_{1}(t),q_{2}(t))\right)\,\mathrm{d}\widetilde{P}(\widetilde{u})\,\mathrm{d}\widetilde{Q}(\widetilde{v}).$
(C.4)

Define the free energy as

$\displaystyle f_{n}(t)$
$\displaystyle\coloneqq\frac{1}{n}\mathbb{E}\left[\log{\mathcal{Z}}_{n,t}\right].$
(C.5)

The Gibbs bracket $\left\langle\cdot\right\rangle_{n,t}$ denotes the
expectation with respect to the posterior distribution in C.3:

$\displaystyle\left\langle g(\widetilde{u},\widetilde{v})\right\rangle_{n,t}$
$\displaystyle\coloneqq\frac{1}{{\mathcal{Z}}_{n,t}}\int_{{\mathbb{R}}^{d}}\int_{{\mathbb{R}}^{n}}g(\widetilde{u},\widetilde{v})\exp\left({\mathcal{H}}_{n,t}(\widetilde{u},\widetilde{v};q_{1}(t),q_{2}(t))\right)\,\mathrm{d}\widetilde{P}(\widetilde{u})\,\mathrm{d}\widetilde{Q}(\widetilde{v}),$
(C.6)

for any $g\colon{\mathbb{R}}^{n}\times{\mathbb{R}}^{d}\to{\mathbb{R}}$ such
that the expectation exists. That is,

$\displaystyle\left\langle g(\widetilde{u},\widetilde{v})\right\rangle_{n,t}$
$\displaystyle=\mathbb{E}\left[g(\widetilde{u}^{*},\widetilde{v}^{*})\,|\,Y_{t},Y_{t}^{u},Y_{t}^{v}\right],$

where we recall the notation

$\displaystyle\widetilde{u}^{*}\coloneqq\Xi^{-1/2}u^{*},\qquad\widetilde{v}^{*}\coloneqq\Sigma^{-1/2}v^{*}.$
(C.7)

We will also use the notation $\left\langle\cdot\right\rangle_{n}$ for the
Gibbs bracket with respect to the original posterior
$\,\mathrm{d}P(u,v\,|\,Y)$ in 4.11.

######  Lemma C.1.

Consider $f_{n}(t)$ defined in C.5 with $t\in\\{0,1\\}$. Assume that
$q_{1}(0),q_{2}(0)$ satisfy

$\displaystyle q_{1}(0)$ $\displaystyle\geq 0,\quad q_{2}(0)\geq
0,\quad\lim_{n\to\infty}q_{1}(0)=\lim_{n\to\infty}q_{2}(0)=0.$

Then we have

$\displaystyle f_{n}(0)$
$\displaystyle={\mathcal{F}}_{n}+{\mathcal{O}}(q_{1}(0)+q_{2}(0)),$ (C.8)
$\displaystyle\lim_{n\to\infty}f_{n}(1)$
$\displaystyle=\psi_{\overline{\Xi}}(\alpha
q_{1}(1)\gamma)+\alpha\psi_{\overline{\Sigma}}(q_{2}(1))+o_{K},$ (C.9)

where $\lim_{K\to\infty}o_{K}=0$.

###### Proof.

To show the first statement C.8, let us control $f_{n}(0)-{\mathcal{F}}_{n}$.
Denoting

$\displaystyle{\mathcal{H}}_{n,t}^{\prime}(\widetilde{u},\widetilde{v};q_{1},q_{2})$
$\displaystyle\coloneqq\alpha
q_{1}\widetilde{u}^{\top}\widetilde{u}^{*}+\sqrt{\alpha
q_{1}}\widetilde{u}^{\top}Z^{u}-\frac{\alpha
q_{1}}{2}\left\|\widetilde{u}\right\|_{2}^{2}+q_{2}\widetilde{v}^{\top}\widetilde{v}^{*}+\sqrt{q_{2}}\widetilde{v}^{\top}Z^{v}-\frac{q_{2}}{2}\left\|\widetilde{v}\right\|_{2}^{2}$

and recalling the Gibbs bracket notation $\left\langle\cdot\right\rangle_{n}$,
we have

$\displaystyle f_{n}(0)-{\mathcal{F}}_{n}$
$\displaystyle=\frac{1}{n}\mathbb{E}\left[\log\frac{{\mathcal{Z}}_{n,0}}{{\mathcal{Z}}_{n}}\right]=\frac{1}{n}\mathbb{E}\left[\log\left\langle\exp\left({\mathcal{H}}_{n,0}^{\prime}(\widetilde{u},\widetilde{v};q_{1}(0),q_{2}(0))\right)\right\rangle_{n}\right],$
(C.10)

where the outer expectation is over all randomness in
${u}^{*},{v}^{*},Z^{u},Z^{v}$. The second equality above follows since

$\displaystyle{\mathcal{H}}_{n,0}(\widetilde{u},\widetilde{v};q_{1}(0),q_{2}(0))$
$\displaystyle={\mathcal{H}}_{n}(\widetilde{u},\widetilde{v})+{\mathcal{H}}_{n,0}^{\prime}(\widetilde{u},\widetilde{v};q_{1}(0),q_{2}(0)).$

We will derive double-sided bounds on $f_{n}(0)-{\mathcal{F}}_{n}$.

To upper bound it, use Jensen’s inequality
$\mathbb{E}\left[\log(\cdot)\right]\leq\log\mathbb{E}\left[\cdot\right]$ on
the partial expectation over $Z^{u},Z^{v}$ in C.10:

$\displaystyle f_{n}(0)-{\mathcal{F}}_{n}$
$\displaystyle\leq\frac{1}{n}\mathop{\mathbb{E}}_{{u}^{*},{v}^{*}}\left[\log\mathop{\mathbb{E}}_{Z^{u},Z^{v}}\left[\left\langle\exp\left({\mathcal{H}}_{n,0}^{\prime}(\widetilde{u},\widetilde{v};q_{1}(0),q_{2}(0))\right)\right\rangle_{n}\right]\right]$
$\displaystyle=\frac{1}{n}\mathop{\mathbb{E}}_{{u}^{*},{v}^{*}}\left[\log\left\langle\mathop{\mathbb{E}}_{Z^{u},Z^{v}}\left[\exp\left({\mathcal{H}}_{n,0}^{\prime}(\widetilde{u},\widetilde{v};q_{1}(0),q_{2}(0))\right)\right]\right\rangle_{n}\right],$

where the equality is legit since $\left\langle\cdot\right\rangle_{n}$ does
not depend on $Z^{u},Z^{v}$. By the Gaussian integral formula (Proposition
G.1), the inner expectation equals

$\displaystyle\mathop{\mathbb{E}}_{Z^{u},Z^{v}}\left[\exp\left({\mathcal{H}}_{n,0}^{\prime}(\widetilde{u},\widetilde{v};q_{1}(0),q_{2}(0))\right)\right]$
$\displaystyle=\exp\left(\alpha
q_{1}(0)\widetilde{u}^{\top}\widetilde{u}^{*}+q_{2}(0)\widetilde{v}^{\top}\widetilde{v}^{*}\right).$

Replacing the Gibbs bracket with $\max$, we obtain an upper bound:

$\displaystyle f_{n}(0)-{\mathcal{F}}_{n}$
$\displaystyle\leq\frac{1}{n}\mathop{\mathbb{E}}_{{u}^{*},{v}^{*}}\left[\log\max_{(\widetilde{u},\widetilde{v})\in\Xi^{-1/2}[-K,K]^{n}\times\Sigma^{-1/2}[-K,K]^{d}}\exp\left(\alpha
q_{1}(0)\widetilde{u}^{\top}\widetilde{u}^{*}+q_{2}(0)\widetilde{v}^{\top}\widetilde{v}^{*}\right)\right]$
$\displaystyle\leq\alpha
q_{1}(0)\left\|\Xi\right\|_{2}^{-1}K^{2}+q_{2}(0)\left\|\Sigma\right\|_{2}^{-1}K^{2}\frac{d}{n}$
$\displaystyle\leq\frac{2\alpha
q_{1}(0)K^{2}}{\inf\operatorname{supp}(\overline{\Xi})}+\frac{2\alpha
q_{2}(0)K^{2}}{\inf\operatorname{supp}(\overline{\Sigma})},$ (C.11)

where the last inequality holds for all sufficiently large $n$ by 3.3 and
$n\asymp d$.

To lower bound $f_{n}-{\mathcal{F}}_{n}$, we use Jensen’s inequality again but
this time on the Gibbs bracket in C.10:

$\displaystyle f_{n}(0)-{\mathcal{F}}_{n}$
$\displaystyle\geq\frac{1}{n}\mathbb{E}\left[\left\langle{\mathcal{H}}_{n,0}^{\prime}(\widetilde{u},\widetilde{v};q_{1}(0),q_{2}(0))\right\rangle_{n}\right]=\frac{1}{n}\mathop{\mathbb{E}}_{{u}^{*},{v}^{*}}\left[\left\langle\mathop{\mathbb{E}}_{Z^{u},Z^{v}}\left[{\mathcal{H}}_{n,0}^{\prime}(\widetilde{u},\widetilde{v};q_{1}(0),q_{2}(0))\right]\right\rangle_{n}\right].$
(C.12)

Since
$(Z^{u},Z^{v})\sim{\mathcal{N}}(0_{n},I_{n})\otimes{\mathcal{N}}(0_{d},I_{d})$,
the inner expectation equals

$\displaystyle\mathop{\mathbb{E}}_{Z^{u},Z^{v}}\left[{\mathcal{H}}_{n,0}^{\prime}(\widetilde{u},\widetilde{v};q_{1}(0),q_{2}(0))\right]$
$\displaystyle=\alpha
q_{1}(0)\widetilde{u}^{\top}\widetilde{u}^{*}-\frac{\alpha
q_{1}(0)}{2}\left\|\widetilde{u}\right\|_{2}^{2}+q_{2}(0)\widetilde{v}^{\top}\widetilde{v}^{*}-\frac{q_{2}(0)}{2}\left\|\widetilde{v}\right\|_{2}^{2}.$

So

$\max_{(\widetilde{u},\widetilde{v})\in\Xi^{-1/2}[-K,K]^{n}\times\Sigma^{-1/2}[-K,K]^{d}}\left|\mathop{\mathbb{E}}_{Z^{u},Z^{v}}\left[{\mathcal{H}}_{n,0}^{\prime}(\widetilde{u},\widetilde{v};q_{1}(0),q_{2}(0))\right]\right|\\\
\leq\frac{3}{2}\alpha
q_{1}(0)\left\|\Xi\right\|_{2}^{-1}nK^{2}+\frac{3}{2}q_{2}(0)\left\|\Sigma\right\|_{2}^{-1}dK^{2}.$
(C.13)

Using this, we obtain a lower bound on $f_{n}-{\mathcal{F}}_{n}$ by replacing
the Gibbs bracket on the RHS of C.12 with $-\max\left|\cdot\right|$:

$\displaystyle f_{n}(0)-{\mathcal{F}}_{n}$
$\displaystyle\geq-\frac{1}{n}\mathop{\mathbb{E}}_{{u}^{*},{v}^{*}}\left[\max_{(\widetilde{u},\widetilde{v})\in\Xi^{-1/2}[-K,K]^{n}\times\Sigma^{-1/2}[-K,K]^{d}}\left|\mathop{\mathbb{E}}_{Z^{u},Z^{v}}\left[{\mathcal{H}}_{n,0}^{\prime}(\widetilde{u},\widetilde{v};q_{1}(0),q_{2}(0))\right]\right|\right]$
$\displaystyle\geq-\frac{3}{2}\alpha
q_{1}(0)\left\|\Xi\right\|_{2}^{-1}K^{2}-\frac{3}{2}q_{2}(0)\left\|\Sigma\right\|_{2}^{-1}K^{2}\frac{d}{n}$
$\displaystyle\geq-\frac{3\alpha
q_{1}(0)K^{2}}{\inf\operatorname{supp}(\overline{\Xi})}-\frac{3\alpha
q_{2}(0)K^{2}}{\inf\operatorname{supp}(\overline{\Sigma})},$ (C.14)

where the last inequality holds for all sufficiently large $n$ by 3.3 and
$n\asymp d$. Combining C.14 and C.11 gives the first result C.8.

We then prove the second statement C.9. Since $f_{n}$ is pseudo-Lipschitz as a
function of the priors, up to a term $o_{K}$ that vanishes as $K\to\infty$
uniformly over $n$, it suffices to ignore the truncation at $K$ and assume
$P={\mathcal{N}}(0,\gamma),Q={\mathcal{N}}(0,1)$. From the definition C.4 of
${\mathcal{Z}}_{n,t}$, we have

$\displaystyle{\mathcal{Z}}_{n,1}=\int_{{\mathbb{R}}^{d}}\int_{{\mathbb{R}}^{n}}\exp\left({\mathcal{H}}_{n,1}^{\prime}(\Xi^{-1/2}u,\Sigma^{-1/2}v;q_{1}(1),q_{2}(1))\right)\,\mathrm{d}P^{\otimes
n}(u)\,\mathrm{d}Q^{\otimes d}(v)$
$\displaystyle=\sqrt{\frac{\det(\Xi)}{\gamma^{n}\det(\Xi/\gamma+\alpha
q_{1}(1)I_{n})}}$ $\displaystyle\times\exp\left(\frac{1}{2}\left(\alpha
q_{1}(1)\widetilde{u}^{*}+\sqrt{\alpha
q_{1}(1)}Z^{u}\right)^{\top}\left(\Xi/\gamma+\alpha
q_{1}(1)I_{n}\right)^{-1}\left(\alpha q_{1}(1)\widetilde{u}^{*}+\sqrt{\alpha
q_{1}(1)}Z^{u}\right)\right)$
$\displaystyle\times\sqrt{\frac{\det(\Sigma)}{\det(\Sigma+q_{2}(1)I_{d})}}\exp\left(\frac{1}{2}\left(q_{2}(1)\widetilde{v}^{*}+\sqrt{q_{2}(1)}Z^{v}\right)^{\top}\left(\Sigma+q_{2}(1)I_{d}\right)^{-1}\left(q_{2}(1)\widetilde{v}^{*}+\sqrt{q_{2}(1)}Z^{v}\right)\right),$

where in the second equality, we use the Gaussian integral formula
(Proposition G.1). Therefore

$\displaystyle\lim_{n\to\infty}f_{n}(1)$
$\displaystyle=\lim_{n\to\infty}\frac{1}{n}\mathbb{E}\left[\log{\mathcal{Z}}_{n,1}\right]$
$\displaystyle=\frac{1}{2}\alpha
q_{1}(1)\gamma\mathbb{E}\left[\overline{\Xi}^{-1}\right]-\frac{1}{2}\mathbb{E}\left[\log\left(1+\alpha
q_{1}(1)\gamma\overline{\Xi}^{-1}\right)\right]$
$\displaystyle\quad+\frac{\alpha}{2}q_{2}(1)\mathbb{E}\left[\overline{\Sigma}^{-1}\right]-\frac{\alpha}{2}\mathbb{E}\left[\log\left(1+q_{2}(1)\overline{\Sigma}^{-1}\right)\right],$

verifying the identity C.9. In the second equality, we have used Proposition
G.2. ∎

######  Lemma C.2 (Free energy variation).

For all $t\in(0,1)$,

$\displaystyle f_{n}^{\prime}(t)$
$\displaystyle=\frac{\alpha}{2}q_{1}^{\prime}(t)q_{2}^{\prime}(t)-\frac{1}{2}\mathbb{E}\left[\left\langle\left(\frac{\widetilde{u}^{\top}\widetilde{u}^{*}}{n}-q_{2}^{\prime}(t)\right)\left(\frac{\widetilde{v}^{\top}\widetilde{v}^{*}}{n}-\alpha
q_{1}^{\prime}(t)\right)\right\rangle_{n,t}\right].$ (C.15)

###### Proof.

From the definitions C.5 and C.4, we compute

$\displaystyle f_{n}^{\prime}(t)$
$\displaystyle=\frac{1}{n}\mathbb{E}\left[\frac{1}{{\mathcal{Z}}_{n,t}}\frac{\partial}{\partial
t}{\mathcal{Z}}_{n,t}\right]$
$\displaystyle=\frac{1}{n}\mathbb{E}\left[\frac{1}{{\mathcal{Z}}_{n,t}}\int_{{\mathbb{R}}^{d}}\int_{{\mathbb{R}}^{n}}\left(\frac{\partial}{\partial
t}{\mathcal{H}}_{n,t}(\widetilde{u},\widetilde{v};q_{1}(t),q_{2}(t))\right)\exp\left({\mathcal{H}}_{n,t}(\widetilde{u},\widetilde{v};q_{1}(t),q_{2}(t))\right)\,\mathrm{d}\widetilde{P}(\widetilde{u})\,\mathrm{d}\widetilde{Q}(\widetilde{v})\right]$
$\displaystyle=\frac{1}{n}\mathbb{E}\left[\left\langle\frac{\partial}{\partial
t}{\mathcal{H}}_{n,t}(\widetilde{u},\widetilde{v};q_{1}(t),q_{2}(t))\right\rangle_{n,t}\right].$

By the definition C.2, the time derivative of the Hamiltonian is

$\displaystyle\begin{split}\frac{\partial}{\partial
t}{\mathcal{H}}_{n,t}(\widetilde{u},\widetilde{v};q_{1}(t),q_{2}(t))&=-\frac{1}{2\sqrt{(1-t)n}}\widetilde{u}^{\top}Z\widetilde{v}-\frac{1}{n}\widetilde{u}^{\top}\widetilde{u}^{*}\widetilde{v}^{\top}\widetilde{v}^{*}+\frac{1}{2n}\left\|\widetilde{u}\right\|_{2}^{2}\left\|\widetilde{v}\right\|_{2}^{2}\\\
&\quad+\alpha
q_{1}^{\prime}(t)\widetilde{u}^{\top}\widetilde{u}^{*}+\frac{\sqrt{\alpha}}{2\sqrt{q_{1}(t)}}q_{1}^{\prime}(t)\widetilde{u}^{\top}Z^{u}-\frac{\alpha}{2}q_{1}^{\prime}(t)\left\|\widetilde{u}\right\|_{2}^{2}\\\
&\quad+q_{2}^{\prime}(t)\widetilde{v}^{\top}\widetilde{v}^{*}+\frac{1}{2\sqrt{q_{2}(t)}}q_{2}^{\prime}(t)\widetilde{v}^{\top}Z^{v}-\frac{1}{2}q_{2}^{\prime}(t)\left\|\widetilde{v}\right\|_{2}^{2}.\end{split}$
(C.16)

The expectation of
$\left\langle\partial_{t}{\mathcal{H}}_{n,t}\right\rangle_{n,t}$ can be
computed using the Stein’s lemma (Proposition G.7). Indeed, let us consider
the term

$\displaystyle\begin{split}\mathbb{E}\left[\left\langle\widetilde{u}^{\top}Z\widetilde{v}\right\rangle_{n,t}\right]&=\sum_{i=1}^{n}\sum_{j=1}^{d}\mathbb{E}\left[\left\langle\widetilde{u}_{i}\widetilde{v}_{j}\right\rangle_{n,t}Z_{i,j}\right]=\sum_{i=1}^{n}\sum_{j=1}^{d}\mathbb{E}\left[\frac{\partial\left\langle\widetilde{u}_{i}\widetilde{v}_{j}\right\rangle_{n,t}}{\partial
Z_{i,j}}\right]\\\
&=\sum_{i=1}^{n}\sum_{j=1}^{d}\sqrt{\frac{1-t}{n}}\mathbb{E}\left[\left\langle\widetilde{u}_{i}^{2}\widetilde{v}_{j}^{2}\right\rangle_{n,t}\right]-\sum_{i=1}^{n}\sum_{j=1}^{d}\sqrt{\frac{1-t}{n}}\mathbb{E}\left[\left\langle\widetilde{u}_{i}\widetilde{v}_{j}\right\rangle_{n,t}^{2}\right]\\\
&=\sqrt{\frac{1-t}{n}}\mathbb{E}\left[\left\langle\left\|\widetilde{u}\right\|_{2}^{2}\left\|\widetilde{v}\right\|_{2}^{2}\right\rangle_{n,t}\right]-\sqrt{\frac{1-t}{n}}\mathbb{E}\left[\left\langle\widetilde{u}^{\top}\widetilde{u}^{*}\widetilde{v}^{\top}\widetilde{v}^{*}\right\rangle_{n,t}\right],\end{split}$
(C.17)

where the last step is by the Nishimori identity (Proposition G.4). So the
first line of C.16 upon taken the Gibbs bracket and the expectation becomes

$\displaystyle-\frac{1}{2n}\mathbb{E}\left[\left\langle\widetilde{u}^{\top}\widetilde{u}^{*}\widetilde{v}^{\top}\widetilde{v}^{*}\right\rangle_{n,t}\right].$

Similar cancellations happen for the second and third lines of C.16. Putting
them together, we obtain

$\displaystyle f^{\prime}_{n}(t)$
$\displaystyle=-\frac{1}{2n^{2}}\mathbb{E}\left[\left\langle\widetilde{u}^{\top}\widetilde{u}^{*}\widetilde{v}^{\top}\widetilde{v}^{*}\right\rangle_{n,t}\right]+\frac{\alpha}{2n}q_{1}^{\prime}(t)\mathbb{E}\left[\left\langle\widetilde{u}^{\top}\widetilde{u}^{*}\right\rangle_{n,t}\right]+\frac{1}{2n}q_{2}^{\prime}(t)\mathbb{E}\left[\left\langle\widetilde{v}^{\top}\widetilde{v}^{*}\right\rangle_{n,t}\right],$

which is the same as C.15 with the parentheses opened up. ∎

In what follows, our strategy is:

  1. 1.

Show that
$\left\langle\widetilde{u}^{\top}\widetilde{u}^{*}\right\rangle_{n,t}$ is
concentrated around its mean
$\mathbb{E}\left[\left\langle\widetilde{u}^{\top}\widetilde{u}^{*}\right\rangle_{n,t}\right]$;

  2. 2.

Choose $q_{2}(t)$ to be the solution to

$\displaystyle q_{2}^{\prime}(t)$
$\displaystyle=\frac{1}{n}\mathbb{E}\left[\left\langle\widetilde{u}^{\top}\widetilde{u}^{*}\right\rangle_{n,t}\right].$

Once Items 1 and 2 are done, we then have

$\displaystyle{\mathcal{F}}_{n}$ $\displaystyle\approx
f_{n}(0)=f_{n}(1)-\int_{0}^{1}f_{n}^{\prime}(t)\,\mathrm{d}t$
$\displaystyle\approx\psi_{\overline{\Xi}}(\alpha
q_{1}(1)\gamma)+\alpha\psi_{\overline{\Sigma}}(q_{2}(1))-\int_{0}^{1}\frac{\alpha}{2}q_{1}^{\prime}(t)q_{2}^{\prime}(t)\,\mathrm{d}t$
$\displaystyle\quad+\int_{0}^{1}\frac{1}{2}\mathbb{E}\left[\left\langle\left(\frac{\widetilde{u}^{\top}\widetilde{u}^{*}}{n}-q_{2}^{\prime}(t)\right)\left(\frac{\widetilde{v}^{\top}\widetilde{v}^{*}}{n}-\alpha
q_{1}^{\prime}(t)\right)\right\rangle_{n,t}\right]\,\mathrm{d}t$
$\displaystyle\approx\psi_{\overline{\Xi}}(\alpha
q_{1}(1)\gamma)+\alpha\psi_{\overline{\Sigma}}(q_{2}(1))-\int_{0}^{1}\frac{\alpha}{2}q_{1}^{\prime}(t)q_{2}^{\prime}(t)\,\mathrm{d}t,$

where the first line above uses C.8 in Lemma C.1; the second line uses C.9 in
Lemma C.1 and Lemma C.2; the third line uses Items 1 and 2. This will almost
lead to the desired characterization of the free energy ${\mathcal{F}}_{n}$ in
Theorem 4.4:

$\displaystyle\sup_{q_{v}\geq 0}\inf_{q_{u}\geq
0}\psi_{\overline{\Xi}}(\alpha\gamma
q_{v})+\alpha\psi_{\overline{\Sigma}}(q_{u})-\frac{\alpha}{2}q_{u}q_{v}.$

Consider the function

$\displaystyle\phi_{t}(q_{1},q_{2})$
$\displaystyle=\frac{1}{n}\log{\mathcal{Z}}_{n,t}.$

Note that $\phi_{t}(q_{1},q_{2})$ also depends on
$n,u^{*},v^{*},Z,Z^{u},Z^{v}$, and
$\mathbb{E}\left[\phi_{t}(q_{1},q_{2})\right]=f_{n}(t)$ where the expectation
is over $(u^{*},v^{*},Z,Z^{u},Z^{v})$ distributed according to C.1.

######  Lemma C.3 (Free energy concentration).

Fix a constant $M>0$. There exists a constant $C>0$ depending only on
$K,M,\alpha,\overline{\Xi},\overline{\Sigma}$ such that for any $t\in[0,1]$,
$0\leq q_{1}(t),q_{2}(t)\leq M$ and sufficiently large $n$,

$\displaystyle\mathbb{E}\left[\left|\phi_{t}(q_{1},q_{2})-\mathbb{E}\left[\phi_{t}(q_{1},q_{2})\right]\right|\right]$
$\displaystyle\leq\frac{C}{\sqrt{n}}.$

###### Proof.

Fix $u^{*},v^{*}$. Consider $\phi_{t}(q_{1},q_{2})$ as a function of
$(Z,Z^{u},Z^{v})$. Then

$\displaystyle\left\|\nabla_{(Z,Z^{u},Z^{v})}\phi_{t}(q_{1},q_{2})\right\|_{2}^{2}$
$\displaystyle=\left\|\nabla_{Z}\phi_{t}(q_{1},q_{2})\right\|_{2}^{2}+\left\|\nabla_{Z^{u}}\phi_{t}(q_{1},q_{2})\right\|_{2}^{2}+\left\|\nabla_{Z^{v}}\phi_{t}(q_{1},q_{2})\right\|_{2}^{2}$
$\displaystyle=\sum_{i=1}^{n}\sum_{j=1}^{d}\left(\frac{1}{n}\mathbb{E}\left[\left\langle\sqrt{\frac{1-t}{n}}\widetilde{u}_{i}\widetilde{v}_{j}\right\rangle_{n,t}\right]\right)^{2}+\sum_{i=1}^{n}\left(\frac{1}{n}\mathbb{E}\left[\left\langle\sqrt{\alpha
q_{1}}\widetilde{u}_{i}\right\rangle_{n,t}\right]\right)^{2}+\sum_{j=1}^{d}\left(\frac{1}{n}\mathbb{E}\left[\left\langle\sqrt{q_{2}}\widetilde{v}_{j}\right\rangle_{n,t}\right]\right)^{2}$
$\displaystyle\leq\frac{1}{n^{3}}\mathbb{E}\left[\left\langle\left\|\widetilde{u}\right\|_{2}^{2}\left\|\widetilde{v}\right\|_{2}^{2}\right\rangle_{n,t}\right]+\frac{\alpha
q_{1}}{n^{2}}\mathbb{E}\left[\left\langle\left\|\widetilde{u}\right\|_{2}^{2}\right\rangle_{n,t}\right]+\frac{q_{2}}{n^{2}}\mathbb{E}\left[\left\langle\left\|\widetilde{v}\right\|_{2}^{2}\right\rangle_{n,t}\right]$
$\displaystyle\leq\frac{\alpha}{n}\left\|\Xi\right\|_{2}^{-1}\left\|\Sigma\right\|_{2}^{-1}K^{4}+\frac{\alpha}{n}M\left\|\Xi\right\|_{2}^{-1}K^{2}+\frac{\alpha}{n}M\left\|\Sigma\right\|_{2}^{-1}K^{2}\leq\frac{C}{n},$

where $C>0$ is a constant depending only on
$\alpha,M,\overline{\Xi},\overline{\Sigma},K$. The penultimate line is by
Cauchy–Schwarz and the last line holds for all sufficiently large $n$ by 3.3
and $n\asymp d$. Then by the Gaussian Poincaré inequality (Proposition G.8),

$\displaystyle\mathop{\mathbb{E}}_{Z,Z^{u},Z^{v}}\left[\left|\phi_{t}(q_{1},q_{2})-\mathop{\mathbb{E}}_{Z,Z^{u},Z^{v}}\left[\phi_{t}(q_{1},q_{2})\right]\right|\right]$
$\displaystyle\leq\sqrt{\mathop{\mathrm{Var}}_{Z,Z^{u},Z^{v}}\left[\phi_{t}(q_{1},q_{2})\right]}\leq\frac{C}{\sqrt{n}}.$
(C.18)

The above result holds for any fixed $u^{*},v^{*}$. We then verify that
$\mathop{\mathbb{E}}_{Z,Z^{u},Z^{v}}\left[\phi_{t}(q_{1},q_{2})\right]$ has
bounded difference as a function of $u^{*},v^{*}$. We do so by bounding the
derivatives of
$\mathop{\mathbb{E}}_{Z,Z^{u},Z^{v}}\left[\phi_{t}(q_{1},q_{2})\right]$ with
respect to $u_{i}^{*},v_{j}^{*}$ for any $i\in[n],j\in[d]$. We have

$\displaystyle\frac{\partial}{\partial
u_{i}^{*}}\mathop{\mathbb{E}}_{Z,Z^{u},Z^{v}}\left[\phi_{t}(q_{1},q_{2})\right]$
$\displaystyle=\frac{1}{n}\mathop{\mathbb{E}}_{Z,Z^{u},Z^{v}}\left[\left\langle\frac{\partial}{\partial
u_{i}^{*}}{\mathcal{H}}_{n,t}(\widetilde{u},\widetilde{v};q_{1},q_{2})\right\rangle_{n,t}\right],$
(C.19)

and a similar expression holds for the derivative with respect to $v_{j}^{*}$.
Recall the definition of ${\mathcal{H}}_{n,t}$ from C.2. We have

$\displaystyle\frac{\partial}{\partial
u^{*}_{i}}{\mathcal{H}}_{n,t}(\Xi^{-1/2}u,\Sigma^{-1/2}v;q_{1},q_{2})$
$\displaystyle=\frac{1-t}{n}v^{\top}\Sigma^{-1}v^{*}(\Xi^{-1})_{i,i}u_{i}+\alpha
q_{1}(\Xi^{-1})_{i,i}u_{i},$

where we have used the fact that $\Xi$ is diagonal. Therefore,

$\displaystyle\left|\frac{\partial}{\partial
u^{*}_{i}}{\mathcal{H}}_{n,t}(\Xi^{-1/2}u,\Sigma^{-1/2}v;q_{1},q_{2})\right|$
$\displaystyle\leq\frac{1}{n}\cdot
dK^{2}\left\|\Sigma\right\|_{2}^{-1}\cdot\left\|\Xi\right\|_{2}^{-1}K+\alpha
q_{1}\cdot\left\|\Xi\right\|_{2}^{-1}K\leq C,$

where $C>0$ is a constant depending only on
$\alpha,M,K,\overline{\Xi},\overline{\Sigma}$. The last inequality holds for
all sufficiently large $n$. A similar bound holds for
$\frac{\partial}{\partial
v^{*}_{j}}{\mathcal{H}}_{n,t}(\Xi^{-1/2}u,\Sigma^{-1/2}v;q_{1},q_{2})$. This,
by C.19, implies that
$\mathop{\mathbb{E}}_{Z,Z^{u},Z^{v}}\left[\phi_{t}(q_{1},q_{2})\right]$ as a
function of $(u^{*},v^{*})$ satisfies the bounded difference property with
$c_{i}=C/n$ (see Proposition G.9). So by Proposition G.9,

$\displaystyle\mathop{\mathbb{E}}_{u^{*},v^{*}}\left[\left|\mathop{\mathbb{E}}_{Z,Z^{u},Z^{v}}\left[\phi_{t}(q_{1},q_{2})\right]-\mathop{\mathbb{E}}_{u^{*},v^{*}}\left[\mathop{\mathbb{E}}_{Z,Z^{u},Z^{v}}\left[\phi_{t}(q_{1},q_{2})\right]\right]\right|\right]$
$\displaystyle\leq\sqrt{\mathop{\mathrm{Var}}_{u^{*},v^{*}}\left[\mathop{\mathbb{E}}_{Z,Z^{u},Z^{v}}\left[\phi_{t}(q_{1},q_{2})\right]\right]}\leq\frac{C}{\sqrt{n}}.$
(C.20)

Finally using C.18 and C.20 and the triangle inequality,

$\displaystyle\mathbb{E}\left[\left|\phi_{t}(q_{1},q_{2})-\mathbb{E}\left[\phi_{t}(q_{1},q_{2})\right]\right|\right]$
$\displaystyle\leq\mathop{\mathbb{E}}_{u^{*},v^{*}}\left[\mathop{\mathbb{E}}_{Z,Z^{u},Z^{v}}\left[\left|\phi_{t}(q_{1},q_{2})-\mathop{\mathbb{E}}_{Z,Z^{u},Z^{v}}\left[\phi_{t}(q_{1},q_{2})\right]\right|\right]+\left|\mathop{\mathbb{E}}_{Z,Z^{u},Z^{v}}\left[\phi_{t}(q_{1},q_{2})\right]-\mathop{\mathbb{E}}_{u^{*},v^{*}}\left[\mathop{\mathbb{E}}_{Z,Z^{u},Z^{v}}\left[\phi_{t}(q_{1},q_{2})\right]\right]\right|\right]$
$\displaystyle\leq\frac{C}{\sqrt{n}},$

concluding the proof. ∎

Suppose $a,b\geq 0$ are constants. With $q_{1}=s_{n}a^{2},q_{2}=s_{n}b^{2}$,
we can write ${\mathcal{H}}_{n,t}$ in C.2 as

$\displaystyle{\mathcal{H}}_{n,t}(\widetilde{u},\widetilde{v};s_{n}a^{2},s_{n}b^{2})$
$\displaystyle={\mathcal{H}}_{n,t}(\widetilde{u},\widetilde{v})+{\mathcal{H}}_{n,a}^{u}(\widetilde{u})+{\mathcal{H}}_{n,b}^{v}(\widetilde{v}),$

where

$\displaystyle{\mathcal{H}}_{n,t}(\widetilde{u},\widetilde{v})$
$\displaystyle\coloneqq\sqrt{\frac{1-t}{n}}\widetilde{u}^{\top}Z\widetilde{v}+\frac{1-t}{n}\widetilde{u}^{\top}\widetilde{u}^{*}\widetilde{v}^{\top}\widetilde{v}^{*}-\frac{1-t}{2n}\left\|\widetilde{u}\right\|_{2}^{2}\left\|\widetilde{vv}\right\|_{2}^{2},$
$\displaystyle{\mathcal{H}}_{n,a}^{u}(\widetilde{u})$
$\displaystyle\coloneqq\alpha
s_{n}a^{2}\widetilde{u}^{\top}\widetilde{u}^{*}+a\sqrt{\alpha
s_{n}}\widetilde{u}^{\top}Z^{u}-\frac{\alpha
s_{n}a^{2}}{2}\left\|\widetilde{u}\right\|_{2}^{2},$
$\displaystyle{\mathcal{H}}_{n,b}^{v}(\widetilde{v})$ $\displaystyle\coloneqq
s_{n}b^{2}\widetilde{v}^{\top}\widetilde{v}^{*}+b\sqrt{s_{n}}\widetilde{v}^{\top}Z^{v}-\frac{s_{n}b^{2}}{2}\left\|\widetilde{v}\right\|_{2}^{2}.$

Fix $A\geq 2$. Define

$\displaystyle\phi^{u}(a)$
$\displaystyle=\frac{1}{ns_{n}}\log\int_{{\mathbb{R}}^{n}}\exp\left({\mathcal{H}}_{n,a}^{u}(\widetilde{u})\right)\,\mathrm{d}\widetilde{P}(\widetilde{u}),$
$\displaystyle\phi^{v}(b)$
$\displaystyle=\frac{1}{ns_{n}}\log\int_{{\mathbb{R}}^{d}}\exp\left({\mathcal{H}}_{n,b}^{v}(\widetilde{v})\right)\,\mathrm{d}\widetilde{Q}(\widetilde{v}),$
$\displaystyle\xi_{n}^{u}(s_{n})$ $\displaystyle\coloneqq\sup_{1/2\leq a\leq
A+1/2}\mathbb{E}\left[\left|\phi^{u}(a)-\mathbb{E}\left[\phi^{u}(a)\right]\right|\right],$
$\displaystyle\xi_{n}^{v}(s_{n})$ $\displaystyle\coloneqq\sup_{1/2\leq b\leq
A+1/2}\mathbb{E}\left[\left|\phi^{v}(b)-\mathbb{E}\left[\phi^{v}(b)\right]\right|\right],$
(C.21) $\displaystyle\Upsilon^{u}(\widetilde{u})$
$\displaystyle\coloneqq\frac{1}{ns_{n}}\frac{\partial}{\partial
a}{\mathcal{H}}_{n,a}^{u}(\widetilde{u}),$
$\displaystyle\Upsilon^{v}(\widetilde{v})$
$\displaystyle\coloneqq\frac{1}{ns_{n}}\frac{\partial}{\partial
b}{\mathcal{H}}_{n,b}^{b}(\widetilde{v}).$

Denote by $\left\langle g(u^{*},v^{*})\right\rangle_{n,a,b}$ the conditional
expectation of $g(u^{*},v^{*})$ given

$\displaystyle Y_{t},\quad Y_{t}^{u}(a)\coloneqq a\sqrt{\alpha
s_{n}}u^{*}+\Xi^{1/2}Z^{u},\quad Y_{t}^{v}(b)\coloneqq
b\sqrt{s_{n}}v^{*}+\Sigma^{1/2}Z^{v},$ (C.22)

where the expectation is with respect to the distribution in C.1.

######  Corollary C.4.

Let $s_{n}=n^{-1/32}$ and $A\leq\sqrt{M/s_{n}}-1/2$ for a constant $M>0$
independent of $n$. Then there exists a constant $C>0$ depending only on
$K,\alpha,\overline{\Xi}$ such that for all sufficiently large $n$,

$\displaystyle\xi_{n}^{u}(s_{n})$ $\displaystyle\leq\frac{C}{s_{n}\sqrt{n}}.$

###### Proof.

Note that if $b=0,t=1$, it holds that
${\mathcal{H}}_{n,1}(\widetilde{u},\widetilde{v};s_{n}a^{2},0)={\mathcal{H}}_{n,a}^{u}(\widetilde{u})$
and $\phi_{1}(s_{n}a^{2},0)=s_{n}\phi^{u}(a)$. The conclusion then follows
immediately from Lemma C.3 since by the assumption on $A$,
$q_{1}=s_{n}a^{2}\in[0,M]$ for any $1/2\leq a\leq A+1/2$. ∎

######  Lemma C.5.

Let $s_{n}=n^{-1/32}$. For all $A\geq 2$,

$\displaystyle\frac{1}{A-1}\int_{1}^{A}\mathbb{E}\left[\left\langle\left|\Upsilon^{u}(\widetilde{u})-\mathbb{E}\left[\left\langle\Upsilon^{u}(\widetilde{u})\right\rangle_{n,a,b}\right]\right|\right\rangle_{n,a,b}\right]\,\mathrm{d}a$
$\displaystyle\leq
C\left(\frac{1}{\sqrt{ns_{n}}}+\sqrt{\xi_{n}^{u}(s_{n})}\right),$

where $C>0$ only depends on $\alpha,\overline{\Xi},K$.

###### Proof.

By the triangle inequality,

$\displaystyle\begin{split}\mathbb{E}\left[\left\langle\left|\Upsilon^{u}(\widetilde{u})-\mathbb{E}\left[\left\langle\Upsilon^{u}(\widetilde{u})\right\rangle_{n,a,b}\right]\right|\right\rangle_{n,a,b}\right]&\leq\mathbb{E}\left[\left\langle\left|\Upsilon^{u}(\widetilde{u})-\left\langle\Upsilon^{u}(\widetilde{u})\right\rangle_{n,a,b}\right|\right\rangle_{n,a,b}\right]\\\
&\quad+\mathbb{E}\left[\left\langle\left|\left\langle\Upsilon^{u}(\widetilde{u})\right\rangle_{n,a,b}-\mathbb{E}\left[\left\langle\Upsilon^{u}(\widetilde{u})\right\rangle_{n,a,b}\right]\right|\right\rangle_{n,a,b}\right].\end{split}$
(C.23)

We will bound the two terms on the RHS separately. We first bound

$\displaystyle\frac{1}{A-1}\int_{1}^{A}\mathbb{E}\left[\left\langle\left|\Upsilon^{u}(\widetilde{u})-\left\langle\Upsilon^{u}(\widetilde{u})\right\rangle_{n,a,b}\right|\right\rangle_{n,a,b}\right]\,\mathrm{d}a$
$\displaystyle\leq\left(\frac{1}{A-1}\int_{1}^{A}\mathbb{E}\left[\left\langle\left(\Upsilon^{u}(\widetilde{u})-\left\langle\Upsilon^{u}(\widetilde{u})\right\rangle_{n,a,b}\right)^{2}\right\rangle_{n,a,b}\right]\,\mathrm{d}a\right)^{1/2}.$
(C.24)

The first two derivatives of $\phi^{u}$ are

$\displaystyle\begin{split}(\phi^{u})^{\prime}(a)&=\frac{1}{ns_{n}}\frac{\int_{{\mathbb{R}}^{n}}\exp\left({\mathcal{H}}_{n,a}^{u}(\widetilde{u})\right)\frac{\partial}{\partial
a}{\mathcal{H}}_{n,a}^{u}(\widetilde{u})\,\mathrm{d}\widetilde{P}(\widetilde{u})}{\int_{{\mathbb{R}}^{n}}\exp\left({\mathcal{H}}_{n,a}^{u}(\widetilde{u})\right)\,\mathrm{d}\widetilde{P}(\widetilde{u})}=\left\langle\Upsilon^{u}(\widetilde{u})\right\rangle_{n,a,b},\end{split}$
(C.25a)
$\displaystyle\begin{split}(\phi^{u})^{\prime\prime}(a)&=\frac{1}{ns_{n}}\left[\left\langle\left(\frac{\partial}{\partial
a}{\mathcal{H}}_{n,a}^{u}(\widetilde{u})\right)^{2}\right\rangle_{n,a,b}-\left\langle\frac{\partial}{\partial
a}{\mathcal{H}}_{n,a}^{u}(\widetilde{u})\right\rangle_{n,a,b}^{2}+\left\langle\frac{\partial^{2}}{\partial
a^{2}}{\mathcal{H}}_{n,a}^{u}(\widetilde{u})\right\rangle_{n,a,b}\right]\\\
&=ns_{n}\left[\left\langle\Upsilon^{u}(\widetilde{u})^{2}\right\rangle_{n,a,b}-\left\langle\Upsilon^{u}(\widetilde{u})\right\rangle_{n,a,b}^{2}\right]+\frac{\alpha}{n}\left[2\left\langle\widetilde{u}^{\top}\widetilde{u}^{*}\right\rangle_{n,a,b}-\left\langle\left\|\widetilde{u}\right\|_{2}^{2}\right\rangle_{n,a,b}\right].\end{split}$
(C.25b)

Since there exists $C>0$ depending only on $\alpha,\overline{\Xi},K$ such that
for all sufficiently large $n$,

$\displaystyle\left|\frac{\alpha}{n}\left[2\left\langle\widetilde{u}^{\top}\widetilde{u}^{*}\right\rangle_{n,a,b}-\left\langle\left\|\widetilde{u}\right\|_{2}^{2}\right\rangle_{n,a,b}\right]\right|$
$\displaystyle\leq\frac{4\alpha
K^{2}}{\inf\operatorname{supp}(\overline{\Xi})}\eqqcolon C,$ (C.26)

the second result C.25b above implies that for all sufficiently large $n$,

$\displaystyle\left\langle\left(\Upsilon^{u}(\widetilde{u})-\left\langle\Upsilon^{u}(\widetilde{u})\right\rangle_{n,a,b}\right)^{2}\right\rangle_{n,a,b}$
$\displaystyle\leq\frac{1}{ns_{n}}\left((\phi^{u})^{\prime\prime}(a)+C\right).$

Consequently,

$\displaystyle\int_{1}^{A}\mathbb{E}\left[\left\langle\left(\Upsilon^{u}(\widetilde{u})-\left\langle\Upsilon^{u}(\widetilde{u})\right\rangle_{n,a,b}\right)^{2}\right\rangle_{n,a,b}\right]\,\mathrm{d}a$
$\displaystyle\leq\frac{1}{ns_{n}}\left(\mathbb{E}\left[(\phi^{u})^{\prime}(A)\right]-\mathbb{E}\left[(\phi^{u})^{\prime}(1)\right]+C(A-1)\right).$
(C.27)

To proceed, we compute for any $a$,

$\displaystyle\mathbb{E}\left[(\phi^{u})^{\prime}(a)\right]$
$\displaystyle=\mathbb{E}\left[\left\langle\Upsilon^{u}(\widetilde{u})\right\rangle_{n,a,b}\right],$
(C.28)

which is by C.25a. By definition,

$\displaystyle\Upsilon^{u}(\widetilde{u})$
$\displaystyle=\frac{1}{n}\left(2\alpha
a\widetilde{u}^{\top}\widetilde{u}^{*}+\sqrt{\alpha/s_{n}}\widetilde{u}^{\top}Z^{u}-\alpha
a\left\|\widetilde{u}\right\|_{2}^{2}\right),$

whose expectation is therefore given by

$\displaystyle\mathbb{E}\left[\left\langle\Upsilon^{u}(\widetilde{u})\right\rangle_{n,a,b}\right]$
$\displaystyle=\frac{2\alpha
a}{n}\mathbb{E}\left[\left\langle\widetilde{u}^{\top}\widetilde{u}^{*}\right\rangle_{n,a,b}\right]+\frac{\sqrt{\alpha}}{n\sqrt{s_{n}}}\mathbb{E}\left[\left\langle\widetilde{u}^{\top}Z^{u}\right\rangle_{n,a,b}\right]-\frac{\alpha
a}{n}\mathbb{E}\left[\left\langle\left\|\widetilde{u}\right\|_{2}^{2}\right\rangle_{n,a,b}\right].$
(C.29)

Using Stein’s lemma (Proposition G.7), the middle term is equal to

$\displaystyle\frac{\sqrt{\alpha}}{n\sqrt{s_{n}}}\mathbb{E}\left[\left\langle\widetilde{u}^{\top}Z^{u}\right\rangle_{n,a,b}\right]$
$\displaystyle=\frac{\sqrt{\alpha}}{n\sqrt{s_{n}}}\sum_{i=1}^{n}\mathbb{E}\left[Z^{u}_{i}\left\langle\widetilde{u}_{i}\right\rangle_{n,a,b}\right]=\frac{\sqrt{\alpha}}{n\sqrt{s_{n}}}\sum_{i=1}^{n}\mathbb{E}\left[\frac{\partial}{\partial
Z^{u}_{i}}\left\langle\widetilde{u}_{i}\right\rangle_{n,a,b}\right]$
$\displaystyle=\frac{a\alpha}{n}\mathbb{E}\left[\left\langle\left\|\widetilde{u}\right\|_{2}^{2}\right\rangle_{n,a,b}\right]-\frac{a\alpha}{n}\mathbb{E}\left[\left\langle\widetilde{u}^{\top}\widetilde{u}^{*}\right\rangle_{n,a,b}\right].$

Therefore,

$\displaystyle\mathbb{E}\left[\left\langle\Upsilon^{u}(\widetilde{u})\right\rangle_{n,a,b}\right]$
$\displaystyle=\frac{\alpha
a}{n}\mathbb{E}\left[\left\langle\widetilde{u}^{\top}\widetilde{u}^{*}\right\rangle_{n,a,b}\right],$
(C.30)

and

$\displaystyle\left|\mathbb{E}\left[\left\langle\Upsilon^{u}(\widetilde{u})\right\rangle_{n,a,b}\right]\right|$
$\displaystyle\leq aC,$ (C.31)

where $C$ depends only on $\alpha,\overline{\Xi},K$.

Using the last inequality, we can further upper bound the RHS of C.27 by
$\frac{CA}{ns_{n}}$ for some $C$ depending only on $\alpha,\overline{\Xi},K$.
Putting this back to C.24, we get

$\displaystyle\frac{1}{A-1}\int_{1}^{A}\mathbb{E}\left[\left\langle\left|\Upsilon^{u}(\widetilde{u})-\left\langle\Upsilon^{u}(\widetilde{u})\right\rangle_{n,a,b}\right|\right\rangle_{n,a,b}\right]\,\mathrm{d}a$
$\displaystyle\leq\sqrt{\frac{CA}{ns_{n}(A-1)}}\leq\sqrt{\frac{2C}{ns_{n}}},$
(C.32)

since $A\geq 2$.

Now it remains to bound the second term on the RHS of C.23 which can be
written as

$\displaystyle\mathbb{E}\left[\left\langle\left|\left\langle\Upsilon^{u}(\widetilde{u})\right\rangle_{n,a,b}-\mathbb{E}\left[\left\langle\Upsilon^{u}(\widetilde{u})\right\rangle_{n,a,b}\right]\right|\right\rangle_{n,a,b}\right]$
$\displaystyle=\mathbb{E}\left[\left|(\phi^{u})^{\prime}(a)-\mathbb{E}\left[(\phi^{u})^{\prime}(a)\right]\right|\right]$
(C.33)

using C.25a. To further bound the RHS, consider the following two functions

$\displaystyle a$ $\displaystyle\mapsto\phi^{u}(a)+\frac{2\alpha
K^{2}a^{2}}{\inf\operatorname{supp}(\overline{\Xi})},\quad
a\mapsto\mathbb{E}\left[\phi^{u}(a)\right]+\frac{2\alpha
K^{2}a^{2}}{\inf\operatorname{supp}(\overline{\Xi})}.$

They are both differentiable and convex for all sufficiently large $n$ since
their second derivatives are non-negative by C.25b and C.26. Applying
Proposition G.10 with the above two functions, taking the expectation and
using the triangle inequality, we have that for any $1\leq a\leq
A,0<a^{\prime}\leq 1/2$,

$\displaystyle\mathbb{E}\left[\left|(\phi^{u})^{\prime}(a)-\mathbb{E}\left[(\phi^{u})^{\prime}(a)\right]\right|\right]$
$\displaystyle\leq\mathbb{E}\left[(\phi^{u})^{\prime}(a+a^{\prime})\right]-\mathbb{E}\left[(\phi^{u})^{\prime}(a-a^{\prime})\right]+3\xi_{n}^{u}(s_{n})/a^{\prime}+\frac{8\alpha
K^{2}a^{\prime}}{\inf\operatorname{supp}(\overline{\Xi})}.$ (C.34)

Then

$\displaystyle\int_{1}^{A}\mathbb{E}\left[(\phi^{u})^{\prime}(a+a^{\prime})\right]-\mathbb{E}\left[(\phi^{u})^{\prime}(a-a^{\prime})\right]\,\mathrm{d}a$
$\displaystyle=\left(\mathbb{E}\left[\phi^{u}(A+a^{\prime})\right]-\mathbb{E}\left[\phi^{u}(1+a^{\prime})\right]\right)-\left(\mathbb{E}\left[\phi^{u}(A-a^{\prime})\right]-\mathbb{E}\left[\phi^{u}(1-a^{\prime})\right]\right)$
$\displaystyle=\left(\mathbb{E}\left[\phi^{u}(A+a^{\prime})\right]-\mathbb{E}\left[\phi^{u}(A-a^{\prime})\right]\right)-\left(\mathbb{E}\left[\phi^{u}(1+a^{\prime})\right]-\mathbb{E}\left[\phi^{u}(1-a^{\prime})\right]\right)$
$\displaystyle=\int_{-a^{\prime}}^{a^{\prime}}\mathbb{E}\left[(\phi^{u})^{\prime}(A+a)\right]\,\mathrm{d}a-\int_{-a^{\prime}}^{a^{\prime}}\mathbb{E}\left[(\phi^{u})^{\prime}(1+a)\right]\,\mathrm{d}a$
$\displaystyle\leq 4a^{\prime}(A+a^{\prime})C\leq 8a^{\prime}AC,$

where the last step is by C.28 and C.31, and $C$ depends only on
$\alpha,\overline{\Xi},K$. Using this in C.34, we obtain

$\displaystyle\int_{1}^{A}\mathbb{E}\left[\left|(\phi^{u})^{\prime}(a)-\mathbb{E}\left[(\phi^{u})^{\prime}(a)\right]\right|\right]\,\mathrm{d}a$
$\displaystyle\leq CA\left(a^{\prime}+\xi_{n}^{u}(s_{n})/a^{\prime}\right),$

and the RHS is minimized by $a^{\prime}=\sqrt{\xi_{n}^{u}(s_{n})}$ which lies
in the interval $(0,1/2]$ for all sufficiently large $n$ due to Corollary C.4.
Using this result in C.33 and integrating over $a$, we have

$\displaystyle\frac{1}{A-1}\int_{1}^{A}\mathbb{E}\left[\left\langle\left|\left\langle\Upsilon^{u}(\widetilde{u})\right\rangle_{n,a,b}-\mathbb{E}\left[\left\langle\Upsilon^{u}(\widetilde{u})\right\rangle_{n,a,b}\right]\right|\right\rangle_{n,a,b}\right]\,\mathrm{d}a$
$\displaystyle\leq\frac{CA}{A-1}\cdot
2\sqrt{\xi_{n}^{u}(s_{n})}\leq{4C}\sqrt{\xi_{n}^{u}(s_{n})},$ (C.35)

since $A\geq 2$.

Finally, combining C.32, C.35 and C.23 proves the lemma. ∎

######  Lemma C.6.

There exists $C>0$ depending only on $\alpha,K,\overline{\Xi}$ such that for
any $A\geq 2$,

$\displaystyle\frac{1}{A-1}\int_{1}^{A}\frac{1}{n^{2}}\mathbb{E}\left[\left\langle\left((\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(2)}-\mathbb{E}\left[\left\langle(\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(2)}\right\rangle_{n,a,b}\right]\right)^{2}\right\rangle_{n,a,b}\right]\,\mathrm{d}a$
$\displaystyle\leq
C\left(\frac{1}{\sqrt{ns_{n}}}+\sqrt{\xi_{n}^{u}(s_{n})}\right),$

where
$\widetilde{u}^{(1)}=\Xi^{-1/2}u^{(1)},\widetilde{u}^{(2)}=\Xi^{-1/2}u^{(2)}$
with $u^{(1)},u^{(2)}$ being two i.i.d. copies from the conditional law of
$u^{*},v^{*}$ given C.22.

###### Proof.

Since
$\operatorname{supp}(P),\operatorname{supp}(Q),\left\|\Xi\right\|_{2}^{-1}$
are all bounded, by the triangle inequality,

$\displaystyle\frac{1}{n}\left|\mathbb{E}\left[\left\langle\Upsilon^{u}(\widetilde{u}^{(1)})(\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(2)}\right\rangle_{n,a,b}\right]-\mathbb{E}\left[\left\langle\Upsilon^{u}(\widetilde{u})\right\rangle_{n,a,b}\right]\mathbb{E}\left[\left\langle(\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(2)}\right\rangle_{n,a,b}\right]\right|$
(C.36)
$\displaystyle\leq\left\|\Xi\right\|_{2}^{-1}K^{2}\mathbb{E}\left[\left\langle\left|\Upsilon^{u}(\widetilde{u})-\mathbb{E}\left[\left\langle\Upsilon^{u}(\widetilde{u})\right\rangle_{n,a,b}\right]\right|\right\rangle_{n,a,b}\right].$
(C.37)

On the other hand, let us compute C.36 on the LHS of the above inequality.
Recall from C.30 in the proof of Lemma C.5 that

$\displaystyle\mathbb{E}\left[\left\langle\Upsilon^{u}(\widetilde{u}^{(1)})\right\rangle_{n,a,b}\right]=\frac{\alpha
a}{n}\mathbb{E}\left[\left\langle(\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(2)}\right\rangle_{n,a,b}\right].$
(C.38)

Similar to C.29, we have

$\mathbb{E}\left[\left\langle\Upsilon^{u}(\widetilde{u}^{(1)})(\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(2)}\right\rangle_{n,a,b}\right]=\frac{2\alpha
a}{n}\mathbb{E}\left[\left\langle(\widetilde{u}^{(1)})^{\top}\widetilde{u}^{*}(\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(2)}\right\rangle_{n,a,b}\right]\\\
+\frac{\sqrt{\alpha}}{n\sqrt{s_{n}}}\mathbb{E}\left[\left\langle(\widetilde{u}^{(1)})^{\top}Z^{u}(\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(2)}\right\rangle_{n,a,b}\right]-\frac{\alpha
a}{n}\mathbb{E}\left[\left\langle\left\|\widetilde{u}^{(1)}\right\|_{2}^{2}(\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(2)}\right\rangle_{n,a,b}\right].$
(C.39)

Using Stein’s lemma and following similar derivations leading to C.30, the
second term on the RHS equals

$\frac{\sqrt{\alpha}}{n\sqrt{s_{n}}}\mathbb{E}\left[\left\langle(\widetilde{u}^{(1)})^{\top}Z^{u}(\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(2)}\right\rangle_{n,a,b}\right]=\frac{a\alpha}{n}\mathbb{E}\left[\left\langle\left\|\widetilde{u}^{(1)}\right\|_{2}^{2}(\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(2)}\right\rangle_{n,a,b}\right]\\\
-\frac{2a\alpha}{n}\mathbb{E}\left[\left\langle(\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(3)}(\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(2)}\right\rangle_{n,a,b}\right]+\frac{a\alpha}{n}\mathbb{E}\left[\left\langle\left((\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(2)}\right)^{2}\right\rangle_{n,a,b}\right].$
(C.40)

Therefore,

$\displaystyle\mathbb{E}\left[\left\langle\Upsilon^{u}(\widetilde{u}^{(1)})(\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(2)}\right\rangle_{n,a,b}\right]$
$\displaystyle=\frac{\alpha
a}{n}\mathbb{E}\left[\left\langle\left((\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(2)}\right)^{2}\right\rangle_{n,a,b}\right].$
(C.41)

Putting C.38 and C.41 together and using Nishimori identity (Proposition G.4),
we have that C.36 equals

$\frac{\alpha
a}{n^{2}}\left|\mathbb{E}\left[\left\langle\left((\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(2)}\right)^{2}\right\rangle_{n,a,b}\right]-\mathbb{E}\left[\left\langle(\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(2)}\right\rangle_{n,a,b}\right]^{2}\right|\\\
=\frac{\alpha
a}{n^{2}}\left|\mathbb{E}\left[\left\langle\left((\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(2)}-\mathbb{E}\left[\left\langle(\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(2)}\right\rangle_{n,a,b}\right]\right)^{2}\right\rangle_{n,a,b}\right]\right|.$
(C.42)

So by the inequality C.37, for any $a\geq 1$,

$\displaystyle\frac{1}{n^{2}}\mathbb{E}\left[\left\langle\left((\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(2)}-\mathbb{E}\left[\left\langle(\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(2)}\right\rangle_{n,a,b}\right]\right)^{2}\right\rangle_{n,a,b}\right]$
$\displaystyle\leq\frac{\left\|\Xi\right\|_{2}^{-1}K^{2}}{\alpha}\mathbb{E}\left[\left\langle\left|\Upsilon^{u}(\widetilde{u})-\mathbb{E}\left[\left\langle\Upsilon^{u}(\widetilde{u})\right\rangle_{n,a,b}\right]\right|\right\rangle_{n,a,b}\right].$

Integrating over $a\in[1,A]$ and invoking Lemma C.5 concludes the proof. ∎

######  Lemma C.7 (Overlap concentration).

Let $R_{1},R_{2}\colon[0,1]\times{\mathbb{R}}_{>0}^{2}\to{\mathbb{R}}_{\geq
0}$ be two continuous bounded functions such that their partial derivatives
with respect to the second and third arguments are continuous and non-
negative. Let $s_{n}=n^{-1/32}$. For
$(\varepsilon_{1},\varepsilon_{2})\in[1,2]^{2}$, slightly abusing notation,
let
$q_{1}(\cdot,\varepsilon_{1},\varepsilon_{2}),q_{2}(\cdot,\varepsilon_{1},\varepsilon_{2})$
be the unique solution to

$\displaystyle\begin{cases}q_{1}(0)=s_{n}\varepsilon_{1}\\\
q_{2}(0)=s_{n}\varepsilon_{2}\end{cases},\qquad\begin{cases}q_{1}^{\prime}(t)=R_{1}(t,q_{1}(t),q_{2}(t))\\\
q_{2}^{\prime}(t)=R_{2}(t,q_{1}(t),q_{2}(t))\end{cases}.$ (C.43)

Then there exists a constant $C>0$ depending only on $K$, $\alpha$,
$\left\|R_{1}\right\|_{\infty}$, $\left\|R_{2}\right\|_{\infty}$,
$\overline{\Xi}$ such that for every $t\in[0,1]$,

$\displaystyle\int_{1}^{2}\int_{1}^{2}\frac{1}{n^{2}}\mathbb{E}\left[\left\langle\left(\widetilde{u}^{\top}\widetilde{u}^{*}-\mathbb{E}\left[\left\langle\widetilde{u}^{\top}\widetilde{u}^{*}\right\rangle_{n,t}\right]\right)^{2}\right\rangle_{n,t}\right]\,\mathrm{d}\varepsilon_{1}\,\mathrm{d}\varepsilon_{2}\leq
Cn^{-1/8}.$

###### Proof.

The existence and uniqueness of the solution to the Cauchy problem C.43 is a
direct consequence of the Cauchy–Lipschitz theorem [33, Theorem 3.1, Chapter
V]. For any $t\in[0,1]$, the function
$Q_{t}(\varepsilon_{1},\varepsilon_{2})=(q_{1}(t,\varepsilon_{1},\varepsilon_{2}),q_{2}(t,\varepsilon_{1},\varepsilon_{2}))$
is a $C^{1}$-diffeomorphism. Its Jacobian determinant is given by the
Liouville formula [33, Corollary 3.1, Chapter V] and can be lower bounded as

$\displaystyle J(\varepsilon_{1},\varepsilon_{2})$
$\displaystyle\coloneqq\det(\nabla Q_{t}(\varepsilon_{1},\varepsilon_{2}))$
$\displaystyle=s_{n}^{2}\exp\left(\int_{0}^{t}\partial_{2}R_{1}(s,Q_{s}(\varepsilon_{1},\varepsilon_{2}))\,\mathrm{d}s+\int_{0}^{t}\partial_{3}R_{2}(s,Q_{s}(\varepsilon_{1},\varepsilon_{2}))\,\mathrm{d}s\right)\geq
s_{n}^{2},$ (C.44)

since the partial derivatives are non-negative by assumptions.

We then view the RHS below as a function of $q_{1},q_{2}$ and denote it by

$\displaystyle V(q_{1},q_{2})$
$\displaystyle=\mathbb{E}\left[\left\langle\left(\widetilde{u}^{\top}\widetilde{u}^{*}-\mathbb{E}\left[\left\langle\widetilde{u}^{\top}\widetilde{u}^{*}\right\rangle_{n,t}\right]\right)^{2}\right\rangle_{n,t}\right].$
(C.45)

Denote $\Omega_{t}=Q_{t}([1,2]^{2})/s_{n}$ and
$M\coloneqq\max\left\\{\left\|R_{1}\right\|_{\infty},\left\|R_{2}\right\|_{\infty}\right\\}+2$.
Since $R_{1},R_{2}\geq 0$ by assumptions, $q_{1},q_{2}$ are non-decreasing in
$t$ by C.43. So for $i\in\\{1,2\\}$, and any $t\in[0,1]$ and
$(\varepsilon_{1},\varepsilon_{2})\in[1,2]^{2}$,

$\displaystyle q_{i}(t,\varepsilon_{1},\varepsilon_{2})/s_{n}$
$\displaystyle\geq
q_{i}(0,\varepsilon_{1},\varepsilon_{2})/s_{n}=\varepsilon_{i}\geq 1,$
$\displaystyle q_{i}(t,\varepsilon_{1},\varepsilon_{2})/s_{n}$
$\displaystyle\leq
q_{i}(1,\varepsilon_{1},\varepsilon_{2})/s_{n}=s_{n}^{-1}\int_{0}^{1}q_{i}^{\prime}(t,\varepsilon_{1},\varepsilon_{2})\,\mathrm{d}t+s_{n}^{-1}q_{i}(0,\varepsilon_{1},\varepsilon_{2})\leq
s_{n}^{-1}(\left\|R_{i}\right\|_{\infty}+2).$

We obtain the relation $\Omega_{t}\subset[1,M/s_{n}]^{2}$ for any $t\in[0,1]$.

Next, using the change of variable
$(r_{1},r_{2})=Q_{t}(\varepsilon_{1},\varepsilon_{2})/s_{n}$, we have

$\displaystyle\int_{1}^{2}\int_{1}^{2}\frac{1}{n^{2}}\mathbb{E}\left[\left\langle\left(\widetilde{u}^{\top}\widetilde{u}^{*}-\mathbb{E}\left[\left\langle\widetilde{u}^{\top}\widetilde{u}^{*}\right\rangle_{n,t}\right]\right)^{2}\right\rangle_{n,t}\right]\,\mathrm{d}\varepsilon_{1}\,\mathrm{d}\varepsilon_{2}$
$\displaystyle=\frac{1}{n^{2}}\int_{1}^{2}\int_{1}^{2}V(q_{1}(t,\varepsilon_{1},\varepsilon_{2}),q_{2}(t,\varepsilon_{1},\varepsilon_{2}))\,\mathrm{d}\varepsilon_{1}\,\mathrm{d}\varepsilon_{2}$
$\displaystyle=\frac{1}{n^{2}}\int_{\Omega_{t}}\frac{V(s_{n}r_{1},s_{n}r_{2})s_{n}^{2}}{J(Q_{t}^{-1}(s_{n}r_{1},s_{n}r_{2}))}\,\mathrm{d}(r_{1},r_{2})$
$\displaystyle\leq\frac{1}{n^{2}}\int_{1}^{M/s_{n}}\int_{1}^{M/s_{n}}V(s_{n}r_{1},s_{n}r_{2})\,\mathrm{d}r_{1}\,\mathrm{d}r_{2},$
(C.46)

where the last step is by C.44. Further applying the change of variable
$r_{1}=a^{2}$, we have that for all $r_{2}\geq 1$,

$\displaystyle\frac{1}{n^{2}}\int_{1}^{M/s_{n}}V(s_{n}r_{1},s_{n}r_{2})\,\mathrm{d}r_{1}$
$\displaystyle=\frac{1}{n^{2}}\int_{1}^{\sqrt{M/s_{n}}}V(s_{n}a^{2},s_{n}r_{2})2a\,\mathrm{d}a\leq\frac{2\sqrt{M/s_{n}}}{n^{2}}\int_{1}^{\sqrt{M/s_{n}}}V(s_{n}a^{2},s_{n}r_{2})\,\mathrm{d}a.$
(C.47)

Recalling $V$ from C.45, we recognize that

$\displaystyle V(s_{n}a^{2},s_{n}r_{2})$
$\displaystyle=\mathbb{E}\left[\left\langle\left(\widetilde{u}^{\top}\widetilde{u}^{*}-\mathbb{E}\left[\left\langle\widetilde{u}^{\top}\widetilde{u}^{*}\right\rangle_{n,a,\sqrt{r_{2}}}\right]\right)^{2}\right\rangle_{n,a,\sqrt{r_{2}}}\right]$
$\displaystyle=\mathbb{E}\left[\left\langle\left((\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(2)}-\mathbb{E}\left[\left\langle(\widetilde{u}^{(1)})^{\top}\widetilde{u}^{(2)}\right\rangle_{n,a,\sqrt{r_{2}}}\right]\right)^{2}\right\rangle_{n,a,\sqrt{r_{2}}}\right],$

where $\left\langle\cdot\right\rangle_{n,a,b}$ is defined in C.22 and the
second equality is by Nishimori identity (Proposition G.4). Since
$\sqrt{M/s_{n}}\geq 2$ for all sufficiently large $n$, applying Lemma C.6, we
get

$\displaystyle\frac{1}{\sqrt{M/s_{n}}-1}\int_{1}^{\sqrt{M/s_{n}}}\frac{1}{n^{2}}V(s_{n}a^{2},s_{n}r_{2})\,\mathrm{d}a$
$\displaystyle\leq
C\left(\frac{1}{\sqrt{ns_{n}}}+\sqrt{\xi_{n}^{u}(s_{n})}\right),$

where $C>0$ depends only on $\alpha,K,\overline{\Xi}$ and $\xi_{n}^{u}(s_{n})$
is given in C.21 with $A=\sqrt{M/s_{n}}$. Using this back in C.47 and then in
C.46, we obtain

$\displaystyle\int_{1}^{2}\int_{1}^{2}\frac{1}{n^{2}}\mathbb{E}\left[\left\langle\left(\widetilde{u}^{\top}\widetilde{u}^{*}-\mathbb{E}\left[\left\langle\widetilde{u}^{\top}\widetilde{u}^{*}\right\rangle_{n,t}\right]\right)^{2}\right\rangle_{n,t}\right]\,\mathrm{d}\varepsilon_{1}\,\mathrm{d}\varepsilon_{2}$
$\displaystyle\leq
2C(M/s_{n})^{2}(1/\sqrt{ns_{n}}+\sqrt{\xi_{n}^{u}(s_{n})}).$

Corollary C.4 guarantees that $\xi_{n}^{u}(s_{n})\leq\frac{C}{s_{n}\sqrt{n}}$
for some $C>0$ depending only on $\alpha,K,\overline{\Xi}$. By the choice of
$s_{n}$, we can finally upper bound (up to a positive constant depending only
on $\alpha,K,\overline{\Xi},M$) the RHS above by

$\displaystyle\frac{1}{s_{n}^{2}}\left(\frac{1}{\sqrt{ns_{n}}}+\frac{1}{\sqrt{s_{n}\sqrt{n}}}\right)\leq\frac{1}{s_{n}^{2.5}}\cdot\frac{2}{n^{1/4}}=2\cdot
n^{-11/64}\leq 2\cdot n^{-1/8},$

for all sufficiently large $n$, which completes the proof. ∎

Recalling $\left\langle\cdot\right\rangle_{n,t}$ defined in C.6, let us
identify
$\mathbb{E}\left[\left\langle\widetilde{u}^{\top}\widetilde{u}^{*}\right\rangle_{n,t}\right]$
as a function of $(t,q_{1},q_{2})$:

$\displaystyle\frac{1}{n}\mathbb{E}\left[\left\langle\widetilde{u}^{\top}\widetilde{u}^{*}\right\rangle_{n,t}\right]$
$\displaystyle=\Delta(t,q_{1},q_{2}).$ (C.48)

Note that $\Delta$ is continuous, non-negative (by Nishimori identity) on
$[0,1]\times{\mathbb{R}}_{\geq 0}^{2}$ and bounded by $K^{2}$. Its partial
derivatives with respect to the second and third arguments are continuous and
non-negative, since the correlation between $\widetilde{u}^{*}$ and
$\left\langle\widetilde{u}\right\rangle_{n,t}$ is a non-decreasing function of
the SNRs $q_{1},q_{2}$.

######  Lemma C.8 (Fundamental sum rule).

In the setting of Lemma C.7, for
$(\varepsilon_{1},\varepsilon_{2})\in[1,2]^{2}$, let
$q_{1}(t,\varepsilon_{1},\varepsilon_{2}),q_{2}(t,\varepsilon_{1},\varepsilon_{2})$
be the solution to C.43 with $R_{2}=\Delta$ defined in C.48. Then we have

$\displaystyle{\mathcal{F}}_{n}$
$\displaystyle=\int_{1}^{2}\int_{1}^{2}\int_{0}^{1}\psi_{\overline{\Xi}}(\alpha\gamma
q_{1}(1,\varepsilon_{1},\varepsilon_{2}))+\alpha\psi_{\overline{\Sigma}}(q_{2}(1,\varepsilon_{1},\varepsilon_{2}))-\frac{\alpha}{2}q_{1}^{\prime}(t,\varepsilon_{1},\varepsilon_{2})q_{2}^{\prime}(t,\varepsilon_{1},\varepsilon_{2})\,\mathrm{d}t\,\mathrm{d}\varepsilon_{1}\,\mathrm{d}\varepsilon_{2}+o(1).$

###### Proof.

Fix $(\varepsilon_{1},\varepsilon_{2})\in[1,2]^{2}$. By the choice
$R_{2}=\Delta$,

$\displaystyle q_{2}^{\prime}(t,\varepsilon_{1},\varepsilon_{2})$
$\displaystyle=\Delta(t,q_{1}(t,\varepsilon_{1},\varepsilon_{2}),q_{2}(t,\varepsilon_{1},\varepsilon_{2}))=\frac{1}{n}\mathbb{E}\left[\left\langle\widetilde{u}^{\top}\widetilde{u}^{*}\right\rangle_{n,t}\right].$

Plugging this into C.15, integrating over
$(\varepsilon_{1},\varepsilon_{2})\in[1,2]^{2}$ and applying Lemma C.7, we
have

$\displaystyle f_{n}^{\prime}(t)$
$\displaystyle=\frac{\alpha}{2}q_{1}^{\prime}(t)q_{2}^{\prime}(t)+o(1),$

where $o(1)\to 0$ as $n\to\infty$, uniformly over $t$. By Lemma C.1, we
conclude

$\displaystyle{\mathcal{F}}_{n}$
$\displaystyle=\int_{1}^{2}\int_{1}^{2}f_{n}(0)\,\mathrm{d}\varepsilon_{1}\,\mathrm{d}\varepsilon_{2}+o(1)=\int_{1}^{2}\int_{1}^{2}\left(f_{n}(1)-\int_{0}^{1}f_{n}^{\prime}(t)\,\mathrm{d}t\right)\,\mathrm{d}\varepsilon_{1}\,\mathrm{d}\varepsilon_{2}+o(1)$
$\displaystyle=\int_{1}^{2}\int_{1}^{2}\int_{0}^{1}\psi_{\overline{\Xi}}(\alpha\gamma
q_{1}(1,\varepsilon_{1},\varepsilon_{2}))+\alpha\psi_{\overline{\Sigma}}(q_{2}(1,\varepsilon_{1},\varepsilon_{2}))-\frac{\alpha}{2}q_{1}^{\prime}(t,\varepsilon_{1},\varepsilon_{2})q_{2}^{\prime}(t,\varepsilon_{1},\varepsilon_{2})\,\mathrm{d}t\,\mathrm{d}\varepsilon_{1}\,\mathrm{d}\varepsilon_{2}+o(1),$

as desired. Here the first and last equalities are by C.8 and C.9,
respectively. ∎

Finally, we prove a pair of matching upper and lower bounds, completing the
proof of Theorem 4.4.

######  Lemma C.9 (Lower bound).

Let $s_{n}=n^{-1/32}$ and $R_{2}=\Delta$. Then

$\displaystyle\liminf_{n\to\infty}{\mathcal{F}}_{n}$
$\displaystyle\geq\sup_{q_{v}\geq 0}\inf_{q_{u}\geq
0}{\mathcal{F}}(q_{u},q_{v}).$

###### Proof.

Fix an arbitrary $q_{v}\geq 0$. Let $R_{1}=q_{v}$. Then
$q_{1}(t,\varepsilon_{1},\varepsilon_{2})=s_{n}\varepsilon_{1}+tq_{v}$ and
Lemma C.8 gives

$\displaystyle{\mathcal{F}}_{n}$
$\displaystyle=\int_{1}^{2}\int_{1}^{2}\int_{0}^{1}\psi_{\overline{\Xi}}(\alpha\gamma(s_{n}\varepsilon_{1}+q_{v}))+\alpha\psi_{\overline{\Sigma}}(q_{2}(1,\varepsilon_{1},\varepsilon_{2}))-\frac{\alpha}{2}q_{v}q_{2}^{\prime}(t,\varepsilon_{1},\varepsilon_{2})\,\mathrm{d}t\,\mathrm{d}\varepsilon_{1}\,\mathrm{d}\varepsilon_{2}+o(1)$
$\displaystyle=\int_{1}^{2}\int_{1}^{2}\psi_{\overline{\Xi}}(\alpha\gamma
q_{v})+\alpha\psi_{\overline{\Sigma}}(q_{2}(1,\varepsilon_{1},\varepsilon_{2}))-\frac{\alpha}{2}q_{v}q_{2}(1,\varepsilon_{1},\varepsilon_{2})\,\mathrm{d}\varepsilon_{1}\,\mathrm{d}\varepsilon_{2}+o(1)$
$\displaystyle\geq\inf_{q_{2}\geq 0}\psi_{\overline{\Xi}}(\alpha\gamma
q_{v})+\alpha\psi_{\overline{\Sigma}}(q_{2})-\frac{\alpha}{2}q_{v}q_{2}+o(1)=\inf_{q_{u}\geq
0}{\mathcal{F}}(q_{u},q_{v})+o(1),$

where the second line holds since $\psi_{\overline{\Xi}}$ is Lipschitz and
$q_{2}(0,\varepsilon_{1},\varepsilon_{2})=s_{n}\varepsilon_{2}=o(1)$. This
completes the proof since the above lower bound holds for all $q_{v}\geq 0$. ∎

######  Lemma C.10 (Upper bound).

Let $s_{n}=n^{-1/32}$ and $R_{2}=\Delta$. Then

$\displaystyle\limsup_{n\to\infty}{\mathcal{F}}_{n}$
$\displaystyle\leq\sup_{q_{v}\geq 0}\inf_{q_{u}\geq
0}{\mathcal{F}}(q_{u},q_{v}).$

###### Proof.

We apply Lemma C.8 with

$\displaystyle R_{1}(t,q_{1},q_{2})$
$\displaystyle=2\alpha\psi_{\overline{\Sigma}}^{\prime}\left(\Delta(t,q_{1},q_{2})\right).$
(C.49)

Since $\psi_{\overline{\Xi}}$ is Lipschitz and convex,

$\displaystyle\psi_{\overline{\Xi}}(\alpha\gamma
q_{1}(1,\varepsilon_{1},\varepsilon_{2}))$
$\displaystyle=\psi_{\overline{\Xi}}\left(\alpha\gamma(q_{1}(1,\varepsilon_{1},\varepsilon_{2})-q_{1}(0,\varepsilon_{1},\varepsilon_{2}))\right)+o(1)=\psi_{\overline{\Xi}}\left(\alpha\gamma\int_{0}^{1}q_{1}^{\prime}(t,\varepsilon_{1},\varepsilon_{2})\,\mathrm{d}t\right)+o(1)$
$\displaystyle\leq\int_{0}^{1}\psi_{\overline{\Xi}}(\alpha\gamma
q_{1}^{\prime}(t,\varepsilon_{1},\varepsilon_{2}))\,\mathrm{d}t+o(1),$

and similarly

$\displaystyle\psi_{\overline{\Sigma}}(q_{2}(1,\varepsilon_{1},\varepsilon_{2}))$
$\displaystyle\leq\int_{0}^{1}\psi_{\overline{\Sigma}}(q_{2}^{\prime}(t,\varepsilon_{1},\varepsilon_{2}))\,\mathrm{d}t+o(1).$

Now Lemma C.8 implies

$\displaystyle{\mathcal{F}}_{n}$
$\displaystyle\leq\int_{1}^{2}\int_{1}^{2}\int_{0}^{1}\psi_{\overline{\Xi}}(\alpha\gamma
q_{1}^{\prime}(t,\varepsilon_{1},\varepsilon_{2}))+\alpha\psi_{\overline{\Sigma}}(q_{2}^{\prime}(t,\varepsilon_{1},\varepsilon_{2}))-\frac{\alpha}{2}q_{1}^{\prime}(t,\varepsilon_{1},\varepsilon_{2})q_{2}^{\prime}(t,\varepsilon_{1},\varepsilon_{2})\,\mathrm{d}t\,\mathrm{d}\varepsilon_{1}\,\mathrm{d}\varepsilon_{2}+o(1)$
$\displaystyle=\int_{1}^{2}\int_{1}^{2}\int_{0}^{1}{\mathcal{G}}(q_{2}^{\prime}(t,\varepsilon_{1},\varepsilon_{2}),q_{1}^{\prime}(t,\varepsilon_{1},\varepsilon_{2}))\,\mathrm{d}t\,\mathrm{d}\varepsilon_{1}\,\mathrm{d}\varepsilon_{2}+o(1),$

where

$\displaystyle{\mathcal{G}}(q_{u},q_{v})$
$\displaystyle\coloneqq\psi_{\overline{\Xi}}(\alpha\gamma
q_{v})+\alpha\psi_{\overline{\Sigma}}(q_{u})-\frac{\alpha}{2}q_{u}q_{v}.$

With the choice of $R_{1}$ in C.49 and $R_{2}=\Delta$, the ODE in Lemma C.7
gives

$\displaystyle q_{1}^{\prime}(t,\varepsilon_{1},\varepsilon_{2})$
$\displaystyle=2\alpha\psi_{\overline{\Sigma}}^{\prime}(q_{2}^{\prime}(t,\varepsilon_{1},\varepsilon_{2})),$

which corresponds to the criticality condition of ${\mathcal{G}}$ with respect
to $q_{u}$:

$\displaystyle\partial_{1}{\mathcal{G}}(q_{2}^{\prime}(t,\varepsilon_{1},\varepsilon_{2}),q_{1}^{\prime}(t,\varepsilon_{1},\varepsilon_{2}))$
$\displaystyle=0.$

Since $\psi_{\overline{\Sigma}}$ is convex and $-\frac{\alpha}{2}q_{u}q_{v}$
is linear in $q_{2}$, we have that ${\mathcal{G}}$ is convex in $q_{u}$.
Therefore

$\displaystyle{\mathcal{G}}(q_{2}^{\prime}(t,\varepsilon_{1},\varepsilon_{2}),q_{1}^{\prime}(t,\varepsilon_{1},\varepsilon_{2}))$
$\displaystyle=\inf_{q_{u}\geq
0}{\mathcal{G}}(q_{u},q_{1}^{\prime}(t,\varepsilon_{1},\varepsilon_{2}))=\inf_{q_{u}\geq
0}{\mathcal{F}}(q_{u},q_{1}^{\prime}(t,\varepsilon_{1},\varepsilon_{2}))\leq\sup_{q_{v}\geq
0}\inf_{q_{u}\geq 0}{\mathcal{F}}(q_{u},q_{v}),$

which completes the proof. ∎

##  Appendix D Proofs of Theorem 4.2 and Corollary 4.3

###  D.1 Proof of 4.7

We compute the derivative of ${\mathcal{F}}_{n}(\gamma)$:

$\displaystyle{\mathcal{F}}_{n}^{\prime}(\gamma)$
$\displaystyle=\frac{1}{n}\mathbb{E}\left[\frac{{\mathcal{Z}}_{n}^{\prime}(\gamma)}{{\mathcal{Z}}_{n}(\gamma)}\right]=\frac{1}{n}\mathbb{E}\left[\left\langle\frac{1}{2\sqrt{\gamma
n}}\widetilde{u}^{\top}Z\widetilde{v}+\frac{1}{n}\widetilde{u}^{\top}\widetilde{u}^{*}\widetilde{v}^{\top}\widetilde{v}^{*}-\frac{1}{2n}\widetilde{u}^{\top}\widetilde{u}\widetilde{v}^{\top}\widetilde{v}\right\rangle_{n}\right]$
$\displaystyle=\frac{1}{2n^{2}}\mathbb{E}\left[\left\langle\widetilde{u}^{\top}\widetilde{u}^{*}\widetilde{v}^{\top}\widetilde{v}^{*}\right\rangle_{n}\right],$
(D.1)

where the last step follows similar calculations in C.17. Since
${\mathcal{F}}_{n}(\gamma)\to\sup_{q_{v}}\inf_{q_{u}}{\mathcal{F}}(q_{u},q_{v})$
as $n\to\infty$, we have
${\mathcal{F}}_{n}^{\prime}(\gamma)\to\frac{\partial}{\partial\gamma}\sup_{q_{v}}\inf_{q_{u}}{\mathcal{F}}(q_{u},q_{v})$.
To compute the RHS, note that

$\displaystyle\sup_{q_{v}\geq 0}\inf_{q_{u}\geq 0}{\mathcal{F}}(q_{u},q_{v})$
$\displaystyle=\sup_{q_{v}\geq 0}\left\\{\psi_{\overline{\Xi}}(\alpha\gamma
q_{v})+\inf_{q_{u}\geq 0}\left\\{\alpha\psi_{\overline{\Sigma}}(\gamma
q_{u})-\frac{\alpha}{2}\gamma q_{u}q_{v}\right\\}\right\\},$

and the value of the infimum does not depend on $\gamma$. Therefore, we have

$\displaystyle\frac{\partial}{\partial\gamma}\sup_{q_{v}\geq 0}\inf_{q_{u}\geq
0}{\mathcal{F}}(q_{u},q_{v})$
$\displaystyle=\psi_{\overline{\Xi}}^{\prime}(\alpha\gamma q_{v}^{*})\alpha
q_{v}^{*}=\frac{\alpha q_{u}^{*}q_{v}^{*}}{2},$ (D.2)

where first equality is by the envelope theorem from [47, Corollary 4] and the
last equality follows since the extremizers $q_{u}^{*},q_{v}^{*}$ solve a pair
of equations in 4.16.

On the other hand, we relate ${\mathcal{F}}_{n}^{\prime}(\gamma)$ to the MMSE
4.2 as follows:

$\displaystyle\operatorname{MMSE}_{n}(\gamma)$
$\displaystyle=\frac{1}{nd}\mathbb{E}\left[\left\|{\widetilde{u}^{*}}(\widetilde{v}^{*})^{\top}-\left\langle\widetilde{u}\widetilde{v}^{\top}\right\rangle_{n}\right\|_{\mathrm{F}}^{2}\right]=\frac{1}{nd}\mathbb{E}\left[\left\|\widetilde{u}^{*}\right\|_{2}^{2}\left\|\widetilde{v}^{*}\right\|_{2}^{2}+\left\|\left\langle\widetilde{u}\widetilde{v}^{\top}\right\rangle_{n}\right\|_{\mathrm{F}}^{2}-2(\widetilde{u}^{*})^{\top}\left\langle\widetilde{u}\widetilde{v}^{\top}\right\rangle_{n}\widetilde{v}^{*}\right]$
$\displaystyle=\frac{\operatorname{Tr}(\Xi^{-1})}{n}\frac{\operatorname{Tr}(\Sigma^{-1})}{d}-\frac{1}{nd}\mathbb{E}\left[\left\langle\widetilde{u}^{\top}\widetilde{u}^{*}\widetilde{v}^{\top}\widetilde{v}^{*}\right\rangle_{n}\right],$
(D.3)

where the last step is by Nishimori identity (Proposition G.4). Combining the
above with D.1 and D.2, we conclude

$\displaystyle\operatorname{MMSE}_{n}(\gamma)$
$\displaystyle\to\mathbb{E}\left[\overline{\Xi}^{-1}\right]\mathbb{E}\left[\overline{\Sigma}^{-1}\right]-q_{u}^{*}q_{v}^{*},$

as claimed.

###  D.2 Proof of 4.8

Recall $Y$ from 4.1 and define for some $\gamma^{\prime}\geq 0$,

$\displaystyle
Y^{\prime}\coloneqq\sqrt{\frac{\gamma^{\prime}}{n}}u^{*}{u^{*}}^{\top}+\Xi^{1/2}Z^{\prime}\Xi^{1/2},$
(D.4)

where $Z^{\prime}\in{\mathbb{R}}^{n\times n}$ is a symmetric random matrix
independent of $u^{*},v^{*}$ with
$Z^{\prime}_{i,i}\overset{\mathrm{i.i.d.}}{\sim}{\mathcal{N}}(0,2)$ and
$Z^{\prime}_{i,j}\overset{\mathrm{i.i.d.}}{\sim}{\mathcal{N}}(0,1)$ for all
$1\leq i<j\leq n$. By similar derivations as before, the free energy
associated with $(Y,Y^{\prime})$ is given by

$\displaystyle{\mathsf{F}}_{n}(\gamma,\gamma^{\prime})$
$\displaystyle=\frac{1}{n}{\mathbb{E}}\Bigg{[}\log\int_{{\mathbb{R}}^{d}}\int_{{\mathbb{R}}^{n}}\exp\bigg{(}{\mathcal{H}}_{n}(\Xi^{-1/2}u,\Sigma^{-1/2}v)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+\frac{1}{2}\sqrt{\frac{\gamma^{\prime}}{n}}u^{\top}\Xi^{-1}Y^{\prime}\Xi^{-1}u-\frac{\gamma^{\prime}}{4n}(u^{\top}\Xi^{-1}u)^{2}\bigg{)}\,\mathrm{d}P^{\otimes
n}(u)\,\mathrm{d}Q^{\otimes n}(v)\Bigg{]},$

where ${\mathcal{H}}_{n}$ is given in 4.12. Denote by
$\left\llangle\cdot\right\rrangle_{n}$ the Gibbs bracket with respect to the
corresponding Hamiltonian. Let

$\displaystyle{\mathsf{F}}(\gamma,\gamma^{\prime})$
$\displaystyle\coloneqq\sup_{q_{u},q_{v}\geq
0}\frac{\gamma^{\prime}}{4}q_{u}^{2}+\frac{\alpha\gamma}{2}q_{u}q_{v}-\psi_{\overline{\Xi}}^{*}\left(\frac{q_{u}}{2}\right)-\alpha\psi_{\overline{\Sigma}}^{*}\left(\frac{q_{v}}{2}\right),$
(D.5)

where $f^{*}$ denotes the monotone conjugate of a convex non-decreasing