We did not manage to prove the universality conjecture, i.e. the independence
of the MMSE on the law of $X$, and its coincidence with the orthogonally
invariant result for a Gaussian noise, in its full generality; in this section
and in the Appendix E we shall content ourselves with a sketch of a possible
roadmap towards a proof of the weak form of the conjecture for arbitrary $D$
(assuming of course it is true), in the Gaussian noise case; the global
structure of the reasoning mirrors the proposition 3.3 in MoWe22 , while the
details are rather different.

At a first level of description the main idea is to define a set of basic
functions $\\{b_{G}\\}_{G\in\mathcal{A}^{(D)}}$, indexed by the marked graphs
with at most $D$ edges introduced in Sec. V.3, whose linear combinations span
all possible polynomial estimators that are equivariant under the permutation
group. They are divided into the “good” and the “bad” ones, with the index set
decomposed as the disjoint union $\mathcal{A}^{(D)}=\mathcal{A}^{(D,{\rm
g})}\cup\mathcal{A}^{(D,{\rm b})}$. This distinction is made in such a way
that the linear combinations of $\\{b_{G}\\}_{G\in\mathcal{A}^{(D,{\rm g})}}$
correspond to the estimators that are equivariant under the orthogonal group.
The computation of ${\rm MMSE}^{(n,D)}$ goes as usual through the minimization
of the quadratic plus linear form with the matrix $\mathcal{M}^{(n,D)}$ and
vector $\mathcal{R}^{(n,D)}$, indexed by $\mathcal{A}^{(D)}$, with elements
$\mathcal{M}^{(n,D)}_{G,G^{\prime}}=\mathbb{E}[{\rm
tr}(b_{G}(Y)b_{G^{\prime}}(Y))]$ and $\mathcal{R}^{(n,D)}_{G}=\mathbb{E}[{\rm
tr}(S\,b_{G}(Y))]$. These can be decomposed in blocks of good and bad
components,

$\mathcal{M}^{(n,D)}=\begin{pmatrix}\mathcal{M}^{(n,D,{\rm
gg})}&\mathcal{M}^{(n,D,{\rm gb})}\\\ (\mathcal{M}^{(n,D,{\rm
gb})})^{T}&\mathcal{M}^{(n,D,{\rm bb})}\end{pmatrix}\
,\qquad\mathcal{R}^{(n,D)}=\begin{pmatrix}\mathcal{R}^{(n,D,{\rm g})}\\\
\mathcal{R}^{(n,D,{\rm b})}\end{pmatrix}\ .$ (129)

Let us denote $\mathcal{M}^{(\infty,D)}$ and $\mathcal{R}^{(\infty,D)}$ their
large $n$ limit. The proof would be complete if one managed to show that

$\mathcal{M}^{(\infty,D)}\ \ \text{is
invertible,}\qquad\mathcal{M}^{(\infty,D,{\rm gb})}=0\ ,\qquad\text{and}\
\mathcal{R}^{(\infty,D,{\rm b})}=0\ .$ (130)

As a matter of fact one could then exchange the large $n$ limit and the
minimization (recall the discussion of Sec. III.2) to write

$\lim_{n\to\infty}{\rm
MMSE}^{(n,D)}=\mu_{S,2}-(\mathcal{R}^{(\infty,D)})^{T}(\mathcal{M}^{(\infty,D)})^{-1}\mathcal{R}^{(\infty,D)}=\mu_{S,2}-(\mathcal{R}^{(\infty,D,{\rm
g})})^{T}(\mathcal{M}^{(\infty,D,{\rm gg})})^{-1}\mathcal{R}^{(\infty,D,{\rm
g})}\ .$ (131)

The last expression only involves the good components, i.e. the orthogonally
invariant ones, and as such reduces precisely to the computation performed in
Sec. IV. The latter involved the limit law $\mu_{Y}$, which is universal for
the addition of an arbitrary Wishart matrix and a GOE distributed one. The
fluctuation moments of the centered traces introduced in Eq. (2) have instead
a non-universal limit covariance, but the reasoning of Sec. IV.3 only relied
on the existence of a Gaussian limit, not on the explicit form of the
covariance. One could thus conclude at this point in favor of the universality
of $\lim_{n\to\infty}{\rm MMSE}^{(n,D)}$. Note that for $D=2$ this is
precisely what we obtained above: the limit matrix in Eq. (119) has the
correct invertible block diagonal structure, the good components being the
first three ones, and the bad row vanishes in $\mathcal{R}^{(\infty,2)}$. We
defer to Appendix E some additional details on a possible implementation of
this strategy, and the difficulties it involves.

##  VI Numerical results on finite-size matrices

We present in this Section the results of numerical simulations on the
denoising of matrices of finite size $n$. We used as a signal distribution the
Wishart ensemble defined in Eq. (107), with matrix elements $X_{i,\mu}$ drawn
either from the Gaussian law with zero mean and variance 1, or from the
Rademacher-Bernouilli law where $X_{i,\mu}=0$ with probability $\rho$,
$1/\sqrt{\rho}$ with probability $\rho/2$ and $-1/\sqrt{\rho}$ with the same
probability, in such a way that the first two moments of $X_{i,\mu}$ coincide
in these two ensembles. The noise was extracted from the GOE, with a width
parametrized by $\Delta$. The BABP denoiser $\widehat{S}(Y)$ was implemented
by a numerical diagonalization of the matrix $Y$, the function
$\mathcal{D}_{\rm BABP}$ being applied on each of its eigenvalues $\lambda$;
the computation of this scalar function was performed with Cardan’s formula to
solve the cubic equation (98) with $z=\lambda$. As explained in Sec. IV.4 the
relevant root can be unambiguously selected on the support of $\mu_{Y}$ by
choosing the one with a positive imaginary part; when $\lambda$ is outside the
support (which can happen at finite $n$) we took the unique branch that
ensures both the continuity of $h_{Y}(\lambda)$ on $\mathbb{R}$ and the
asymptotic behavior $h_{Y}(\lambda)\sim 1/\lambda$ as $\lambda\to\pm\infty$.

In the left panel of Fig. 7 we plot the MSE achieved by the BABP denoiser as a
function of the noise intensity $\Delta$ for a given aspect ratio $\alpha$,
and compare it to the analytical prediction ${\rm MSE}_{\rm BABP}$ from Eq.
(76). The agreement, which can be expected to hold exactly only in the large
$n$ limit, is very good already for matrices of size $n=100$, except at large
noise $\Delta$ where some deviations are visible. These finite-size effects
are further studied in the right panel of Fig. 7, where the MSE is now plotted
as a function of $1/n$ for one choice of $(\alpha,\Delta)$. The data
extrapolates in a very accurate way to the analytical prediction in the
$n\to\infty$ limit, with dominant corrections of order $1/n$. Recall that the
analytical formula (76) was obtained by taking the $n\to\infty$ limit after
the $D\to\infty$ limit, hence this numerical result supports the hypothesis
that the same result is obtained when the limits are taken in the reverse
order. Furthermore, the two ensembles of Wishart matrices (Gaussian and
Rademacher-Bernoulli) exhibit different finite-size corrections, but the
extrapolation of the MSE in the large $n$ limit coincide with ${\rm MSE}_{\rm
BABP}$ for both of them. Of course this is not a compelling evidence for the
universality conjecture stated above, since another estimator could in
principe improve the asymptotic MSE for the Rademacher-Bernoulli ensemble. Yet
this shows that the properties of the eigenvectors of the matrices $S$ in both
ensembles are sufficiently similar for the computations recalled in (72) to
give the same asymptotic result.

Figure 7: The MSE reached by the BABP denoiser on Wishart matrices with aspect
ratio $\alpha=1$ and $X_{i,\mu}$ drawn as a standard Gaussian (points denoted
G) or as a Rademacher-Bernouilli random variable with $\rho=0.2$ (RB). Left
panel: the MSE is plotted as a function of $\Delta$, the solid line is the
analytical prediction from Eq. (76), the symbols correspond to matrices of
size $n=100$, each point is an average over 200 repetitions, error bars are of
the order of the symbol size. Right panel: study of the finite-size effects
for $\Delta=8$. The symbols have been obtained from simulations with matrices
of size $n$ between 100 and 300, each point corresponding to an average over
2000 repetitions. The dashed lines are fits of the data of the form $a+b/n$,
with the asymptotic value $a$ fixed to the analytical prediction ${\rm
MSE}_{\rm BABP}$.

In figure 8 we report some results obtained with the second order permutation
invariant denoiser studied in Sec. V.4, that estimates the off-diagonal
elements of $S$ with the linear combinations of the four functions of Eq.
(116). For each finite size $n$ the optimal coefficients have been computed by
solving numerically the linear system of dimension 4
$\mathcal{M}^{(n,2)}c=\mathcal{R}^{(n,2)}$, with the exact expressions at
finite $n$ of $\mathcal{M}^{(n,2)}$ and $\mathcal{R}^{(n,2)}$ given in
Appendix D. In the left panel we check, for one choice of $(\alpha,\Delta)$,
that the expression of ${\rm MMSE}_{\rm od}^{(n,2)}$ computed from (117)
coincides, as it should, with the MSE obtained by the optimal second order
denoiser. This plot also illustrates the universality of the large $n$ limit
that converges to ${\rm MMSE}^{(2)}$, while the finite $n$ corrections depend
on the details of the law of $X$. The right panel of Fig. 8 presents a
comparison between this denoiser and the BABP one (to allow for a fair
comparison we computed the off-diagonal part of the MSE for both denoisers).
What is to be noticed on this plot is the fact that for small sizes the
optimal denoiser of degree 2 reaches a smaller MSE than the BABP one. As a
matter of fact the latter is expected to be optimal only in the infinite size
limit, finite degree estimators with $n$-dependent coefficients (and possibly
terms breaking the orthogonal invariance as $b_{4}$ for RB matrices) may well
be more accurate at finite $n$. This happens here only for rather small
matrices ($n\lesssim 40$ in the RB case, $n\lesssim 25$ in the Gaussian one),
and with an ad-hoc choice of the parameters $(\alpha,\Delta)$ such that ${\rm
MMSE}^{(2)}$ is only slightly larger than ${\rm MSE}_{\rm BABP}$. One can
however expect that including polynomial terms of degree higher than 2 in the
exact finite $n$ treatment would widen the range of $(n,\alpha,\Delta)$ for
which this phenomenon happens.

Figure 8: A numerical study of the second order denoiser of Sec. V.4, for
$\alpha=1$ and $\Delta=0.8$. Left panel: the two solid lines are the
predictions MMSE${}_{\rm od}^{(n,2)}$ of Eq. (117), with the non-universal
parameters corresponding to $X_{i,\mu}$ of laws either Gaussian (G) or
Rademacher-Bernoulli (RN) with $\rho=0.2$. The symbols represent the results
of numerical simulations with matrices of size between $n=50$ and $n=225$,
averaged over $10^{4}$ repetitions for the smallest sizes, $10^{3}$ for the
largest ones (error bars are of the order of the symbol size). Right panel:
the solid lines are the same as in the left panel, the symbols are obtained
with the BABP estimator, each point is an average over $10^{4}$ samples.

##  VII Conclusions

Let us conclude this manuscript with a few remarks and perspectives for
possible future work. We focused our attention on the additive noise model
$Y=S+Z$, but the multiplicative noise one $Y=\sqrt{S}Z\sqrt{S}$ could be
handled along the same lines. In particular the coefficients of the optimal
scalar polynomial estimators of degree at most $D$ will again be given as the
solution of the matrix equation
$\widehat{\mathcal{M}}^{(D)}c=\widehat{\mathcal{R}}^{(D)}$, where
$\widehat{\mathcal{M}}$ will be the Hankel matrix of
$\mu_{Y}=\mu_{S}\boxtimes\mu_{Z}$, the multiplicative free convolution
NiSp_book ; PoBo_book  replacing the additive one, and where the coefficients
of the right hand side will read $\widehat{\mathcal{R}}_{p}=\tau(S(SZ)^{p})$.
Free-probability manipulations analogous to those performed in Sec. IV.2.3 and
IV.2.4 should allow to recover, in the large $D$ limit, the expression of the
denoiser given in BuAlBoPo16 .

Another direction that could be investigated concerns the case of non-
symmetric or rectangular matrices, for which TrErKrMaZd22 ; PoMa23b ; LaMeGa23
derived an estimator that acts on the singular values of the observed matrix,
keeping its singular vectors untouched. The low-degree polynomial approach
presented here could most probably be adapted to this case; this rectangular
setting is also relevant for the estimation of cross-correlations in datasets
BeBoPo23 .

One could also investigate further the case of an orthogonally invariant
signal perturbed by a non-Gaussian Wigner noise: the matrix elements
$(\widehat{S}(Y))_{i,j}$ of the asymptotically optimal polynomial estimators
of degree at most $D$ will certainly contain a linear combination of
$(Y^{p})_{i,j}$ and $Y_{i,j}^{r}$, with $p,r\leq D$ and $r$ odd, as was
demonstrated for $D=3$ in Sec. V.5. Whether these are the only relevant terms
for generic $D$ is not clear to us at the moment.

At variance with the derivation of BuAlBoPo16  we did not directly deal with
the overlaps between the eigenvectors of $S$ and those of $Y=S+Z$; however the
free-probability computation of $\widehat{\mathcal{R}}_{p}=\tau(S(S+Z)^{p})$
performed in Sec. IV.2.3 can actually be interpreted in terms of these
overlaps, as explained in Section 6.2 of Bi03  (see also Section 19.4 of
PoBo_book  for a related re-interpretation of this result).

A drawback of our approach is the order in which the limits $n\to\infty$ and
$D\to\infty$ are taken, one could indeed worry that estimators that are not
polynomials, or polynomials of growing order $D(n)$, reach a better accuracy.
For orthogonally invariant priors there are however some reasons (besides the
coincidence with the result of BuAlBoPo16 ) to believe that this is not an
issue: recall indeed that, at any finite $n$, the optimal estimator is the
posterior mean, and is thus equivariant. The generic form of such functions
was given in Eq. (31); assuming some regularity for its dependency on the
first argument of $\widehat{f}$ the latter can be approximated within an
arbitrary precision by a polynomial of fixed order. Moreover the traces of
growing powers of $Y$ are more and more concentrated on the right edge of the
support of $\mu_{Y}$, and hence should not spoil the optimality result
obtained by taking $D\to\infty$ after $n\to\infty$.

The status of this interversion of limit is less clear when the priors are
only approximately orthogonally invariant. As a matter of fact the
universality we conjecture for the denoising of Wishart matrices perturbed by
GOE noise seems to contradict some of the findings of CaMe23a ; CaMe23b . This
paper presents indeed numerical results on the denoising of Wishart matrices
with a Rademacher-Bernouilli law for the entries of $X$ by an estimator that
reaches a smaller MSE than the BABP one, in a certain range of the parameters
$(\alpha,\Delta)$. Two possible explanations of this apparent contradiction
are that either the universality conjecture is plainly wrong, or the numerical
results from CaMe23a ; CaMe23b  suffer from finite-size effects and are not
representative of the asymptotic behavior as $n\to\infty$. There is however a
third, more far-fetched but richer scenario: it would be that the weak form of
the conjecture is true but the strong one is wrong, i.e. that for all finite
$D$ the asymptotic MMSE among polynomial estimators of degree $D$ is
universal, but this is not the case for unbounded degree estimators. Since
low-degree polynomials are expected to emulate all polynomial-time algorithms,
and since the procedure of CaMe23a ; CaMe23b  has an exponential computational
cost, this third explanation would suggest that matrix denoising with some
priors that break the orthogonal invariance can have a hard phase in which the
information-theoretic optimal accuracy cannot be reached in polynomial time.

##  Appendix A Some free probability formulas

We collect in this Appendix some definitions on transforms of probability
measures, and a computation for the moments of the free convolution of a
Marcenko-Pastur distribution with a semi-circular law that was stated without
proof in Eq. (100).

We consider a probability measure $\mu$ on the reals, with bounded support.
Its moments are denoted $\mu_{p}=\int\mu({\rm d}\lambda)\lambda^{p}$. One
defines its Cauchy transform $g(z)$, for $z\in\mathbb{C}\setminus\mathbb{R}$,
by $g(z)=\int\mu({\rm d}\lambda)\frac{1}{z-\lambda}$. The $R$-transform of
$\mu$ is the function $R(g)=z(g)-\frac{1}{g}$, where $z(g)$ is the functional
inverse of $g(z)$. This transformation is invertible, i.e. the knowledge of
$R$ caracterizes the probability distribution $\mu$. The free cumulants
$\kappa_{p}$ of $\mu$ are the coefficients of the series expansion of $R(g)$
around $g=0$, namely $R(g)=\sum_{p=0}^{\infty}\kappa_{p+1}g^{p}$. When
manipulating different probability measures $\mu_{S}$, $\mu_{Z}$ and $\mu_{Y}$
we add the corresponding index to the various transforms, and denote for
instance $g_{Y}$, $R_{Y}$, $\mu_{Y,p}$ and $\kappa_{Y,p}$ the various
quantities associated to $\mu_{Y}$.

The free convolution is an operation taking two probability measures and
forming a third one, in such a way that the $R$-transform (and hence the free
cumulants) is additive. Explicitly, $\mu_{Y}=\mu_{S}\boxplus\mu_{Z}$ is
defined by $R_{Y}(g)=R_{S}(g)+R_{Z}(g)$, or
$\kappa_{Y,p}=\kappa_{S,p}+\kappa_{Z,p}$ for all $p$.

The free cumulant of order $p$ of a measure $\mu_{Y}$, $\kappa_{Y,p}$, can be
expressed in terms of the moments $\\{\mu_{Y,p^{\prime}}\\}_{p^{\prime}\leq
p}$, and reciprocally $\mu_{Y,p}$ is a function of
$\\{\kappa_{Y,p^{\prime}}\\}_{p^{\prime}\leq p}$. Explicit forms of these
relationships can be obtained for small values of $p$ by truncating the formal
expansions of the power series $R_{Y}(g)$ and $g_{Y}(z)$ (around $g=0$ and
$z=\infty$), imposing the functional inverse relationship between $g_{Y}$ and
$R_{Y}$ order by order. A more systematic approach is based on the Lagrange
inversion theorem, that allows to compute the coefficients of the functional
inverse of a formal power series in terms of the coefficients of the series
itself. Its application in the present context yields a simple expression of
the moments in terms of the $R$-transform (see Prop. 16.20 in NiSp_book ):

$\mu_{Y,p}=\frac{1}{p+1}[g^{p}](1+gR_{Y}(g))^{p+1}\ ,$ (132)

where $[g^{p}]f(g)$ denotes the coefficient of the $p$-th term in the formal
power series $f(g)$.

We will now apply this formula to the case considered in Sec. IV.4, where
$\mu_{Y}$ is the free convolution of a Marcenko-Pastur distribution with a
semi-circular law (its $R$-transform was given in (96)). This yields

$\displaystyle\mu_{Y,p}$ $\displaystyle=\frac{1}{p+1}[g^{p}]\left(1+\Delta
g^{2}+g^{2}\frac{1}{1-\sqrt{\alpha}g}\right)^{p+1}$ (133)
$\displaystyle=\frac{1}{p+1}\sum_{\begin{subarray}{c}q,r\geq 0\\\ q+r\leq
p+1\end{subarray}}\binom{p+1}{q,r}\Delta^{q}[g^{p}]g^{2q+2r}(1-\sqrt{\alpha}g)^{-r}$
(134) $\displaystyle=\frac{1}{p+1}\sum_{\begin{subarray}{c}q,r\geq 0\\\
q+r\leq p+1\\\ 2(q+r)\leq
p\end{subarray}}\binom{p+1}{q,r}\Delta^{q}[g^{p-2q-2r}](1-\sqrt{\alpha}g)^{-r}\
,$ (135)

with a multinomial expansion to go from the first to the second line. For
$r=0$, $[g^{k}](1-\sqrt{\alpha}g)^{-r}=\delta_{k,0}$, while for $r\geq 1$,

$\displaystyle[g^{k}](1-\sqrt{\alpha}g)^{-r}$
$\displaystyle=(-1)^{k}\alpha^{\frac{k}{2}}\frac{1}{k!}(-r)(-r-1)\dots(-r-k+1)$
(136) $\displaystyle=\alpha^{\frac{k}{2}}\frac{1}{k!}r(r+1)\dots(r+k-1)$ (137)
$\displaystyle=\alpha^{\frac{k}{2}}\binom{r+k-1}{r-1}\ .$ (138)

Separating the contributions of $r=0$ and $r\geq 1$ in (135), and noting than
in the latter case the condition $2(q+r)\leq p$ implies $q+r\leq p+1$, one
obtains (100).

##  Appendix B Equivariant functions

###  B.1 Orthogonal group, direct reasoning

In this appendix we shall characterize the functions $f:M_{n}^{\rm
sym}(\mathbb{R})\to M_{n}^{\rm sym}(\mathbb{R})$, where $M_{n}^{\rm
sym}(\mathbb{R})$ is the set of $n\times n$ real symmetric matrices, that are
equivariant under conjugation by orthogonal matrices, i.e. such that
$f(OYO^{T})=Of(Y)O^{T}$ for all $O\in\mathcal{O}_{n}$ and $Y\in M_{n}^{\rm
sym}(\mathbb{R})$, where $\mathcal{O}_{n}=\\{O\in
M_{n}(\mathbb{R}):OO^{T}=O^{T}O=\mathds{1}_{n}\\}$ is the orthogonal group.

In a first step we shall use the fact that any real symmetric matrix $Y$ is
diagonalizable as $Y=P\,{\rm diag}(\lambda_{1},\dots,\lambda_{n})P^{T}$, with
$P\in\mathcal{O}_{n}$; for an equivariant function this yields $f(Y)=P\,f({\rm
diag}(\lambda_{1},\dots,\lambda_{n}))\,P^{T}$, it is thus enough to
characterize $f$ when its argument is a diagonal matrix. It turns out that the
image of a diagonal matrix by an equivariant $f$ is also diagonal (we follow
here the proof of this fact given in Appendix A.1 of PoMa23 ): let us define,
for $k\in[n]$, the matrix $I(k)={\rm diag}(1,\dots,1,-1,1,\dots 1)$, with the
negative element in the $k$-th row. Since $I(k)=I(k)^{T}$ and
$I(k)^{2}=\mathds{1}_{n}$, this matrix is orthogonal, and for a diagonal
matrix $D$ one has $I(k)DI(k)=D$. As a consequence for an equivariant function
$f$ one obtains $f(D)=f(I(k)DI(k))=I(k)f(D)I(k)$. Taking the matrix elements
of this equation gives
$f(D)_{i,j}=(-1)^{\mathds{1}(i=k)+\mathds{1}(j=k)}f(D)_{i,j}$. In particular
for $i\neq j=k$, $f(D)_{i,k}=-f(D)_{i,k}=0$. Since this is valid for all $k$,
one concludes that an equivariant $f$ sends diagonal matrices to diagonal
matrices.

We can thus write $f({\rm diag}(\lambda_{1},\dots,\lambda_{n}))={\rm
diag}(f_{1}(\lambda_{1},\dots,\lambda_{n}),\dots,f_{n}(\lambda_{1},\dots,\lambda_{n}))$,
and it remains to characterize the $n$ functions
$f_{i}:\mathbb{R}^{n}\to\mathbb{R}$. For $\pi\in\mathcal{S}_{n}$ a permutation
of $[n]$, consider the matrix $\Pi$ with matrix elements
$\Pi_{i,j}=\mathds{1}(j=\pi(i))$, which is easily seen to belong to
$\mathcal{O}_{n}$. One has $\Pi\,{\rm
diag}(\lambda_{1},\dots,\lambda_{n})\,\Pi^{T}={\rm
diag}(\lambda_{\pi(1)},\dots,\lambda_{\pi(n)})$, hence the functions $f_{i}$
obey the equations

$f_{i}(\lambda_{\pi(1)},\dots,\lambda_{\pi(n)})=f_{\pi(i)}(\lambda_{1},\dots,\lambda_{n})\quad\forall\pi\in\mathcal{S}_{n}\
.$ (139)

Consider first the $(n-1)!$ permutations that fixes the point $1$, i.e. such
that $\pi(1)=1$. Equation (139) with $i=1$ becomes

$f_{1}(\lambda_{1},\lambda_{2}\dots,\lambda_{n})=f_{1}(\lambda_{1},\lambda_{\pi(2)}\dots,\lambda_{\pi(n)})\
,$ (140)

which shows that the dependency of $f_{1}$ on its last $n-1$ arguments is
symmetric under the permutations. The theory of symmetric functions allows
then to conclude on the existence of a function $\widetilde{f}$ such that

$f_{1}(\lambda_{1},\lambda_{2},\dots,\lambda_{n})=\widetilde{f}\left(\lambda_{1};\sum_{j=2}^{n}\lambda_{j},\sum_{j=2}^{n}\lambda_{j}^{2},\dots,\sum_{j=2}^{n}\lambda_{j}^{n-1}\right)\
.$ (141)

Indeed, such a symmetric dependency on $(\lambda_{2},\dots,\lambda_{n})$ can
be interpreted as a dependency on the multiset
$\\{\lambda_{2},\dots,\lambda_{n}\\}$, with multiplicities, or equivalently on
the roots of the univariate polynomial
$P(t)=(t-\lambda_{2})(t-\lambda_{3})\dots(t-\lambda_{n})$. There are several
equivalent ways to encode the coefficients of this polynomial, see for
instance Stanley_vol2  for more details, one being in terms of the powersums,
which justifies the expression (141). Defining a new function $\widehat{f}$ as
$\widehat{f}(\lambda;x_{1},\dots,x_{n-1})=\widetilde{f}(\lambda;x_{1}-\lambda,x_{2}-\lambda^{2},\dots,x_{n-1}-\lambda^{n-1})$,
we arrive at

$f_{1}(\lambda_{1},\lambda_{2},\dots,\lambda_{n})=\widehat{f}\left(\lambda_{1};\sum_{j=1}^{n}\lambda_{j},\sum_{j=1}^{n}\lambda_{j}^{2},\dots,\sum_{j=1}^{n}\lambda_{j}^{n-1}\right)\
.$ (142)

Consider now an index $k\neq 1$, and a permutation $\pi$ such that $\pi(1)=k$;
from (139) with $i=1$ we get

$\displaystyle f_{k}(\lambda_{1},\dots,\lambda_{n})$
$\displaystyle=f_{1}(\lambda_{k},\lambda_{\pi(2)},\dots,\lambda_{\pi(n)})$
(143) $\displaystyle=\widetilde{f}\left(\lambda_{k};\sum_{j\neq
k}\lambda_{j},\sum_{j\neq k}\lambda_{j}^{2},\dots,\sum_{j\neq
k}\lambda_{j}^{n-1}\right)$ (144)
$\displaystyle=\widehat{f}\left(\lambda_{k};\sum_{j=1}^{n}\lambda_{j},\sum_{j=1}^{n}\lambda_{j}^{2},\dots,\sum_{j=1}^{n}\lambda_{j}^{n-1}\right)\
.$ (145)

We can conclude at this point that an arbitrary equivariant function $f$ acts
on diagonal matrices as $f(D)=\widehat{f}(D;{\rm Tr}(D),\dots,{\rm
Tr}(D^{n-1}))$, with the convention that the function $\widehat{f}$ acts
componentwise through its first argument. Coming back to arbitrary matrices we
write $Y=PDP^{T}$ with $P\in\mathcal{O}_{n}$ and $D$ diagonal, use the
equivariance of $f$ to deduce $f(Y)=Pf(D)P^{T}$, and the invariance of the
trace under orthogonal conjugation to write ${\rm Tr}(Y^{j})={\rm Tr}(D^{j})$.
This yields $f(Y)=\widehat{f}(Y;{\rm Tr}(Y),\dots,{\rm Tr}(Y^{n-1}))$, under
the convention that a scalar function acts on $Y$ by being applied
componentwise on its eigenvalues without altering its eigenvectors, and
concludes the justification of the statement (31) of the main text, namely
that all equivariant functions $f$ are of the form $f(Y)=\widehat{f}(Y;{\rm
Tr}(Y),\dots,{\rm Tr}(Y^{n-1}))$, for an arbitrary function
$\widehat{f}:\mathbb{R}^{n}\to\mathbb{R}$.

As a side remark let us mention that a similar reasoning can be made to
characterize the invariant functions $f:M_{n}^{\rm
sym}(\mathbb{R})\to\mathbb{R}$, such that $f(OYO^{T})=f(Y)$ for all
$O\in\mathcal{O}_{n}$ and $Y\in M_{n}^{\rm sym}(\mathbb{R})$: these are of the
form $f(Y)=F({\rm Tr}(Y),{\rm Tr}(Y^{2}),\dots,{\rm Tr}(Y^{n}))$, for an
arbitrary function $F:\mathbb{R}^{n}\to\mathbb{R}$.

###  B.2 Orthogonal group, via Weingarten calculus

We present here an alternative justification of a weaker version of the
statement made above on the structure of the equivariant functions, based on
the symmetrization defined in Eq. (23), which in the present context amounts
to take an arbitrary function $f:M_{n}^{\rm sym}(\mathbb{R})\to M_{n}^{\rm
sym}(\mathbb{R})$, and symmetrize it with the action considered here, defining

$f^{\rm eq}(Y)=\mathbb{E}_{O}[O^{T}f(OYO^{T})O]\ ,$ (146)

where the expectation is over the Haar measure of the orthogonal group, that
is invariant under multiplication by a fixed element, as required in Sec.
III.3. One can check indeed that $f^{\rm eq}$ is equivariant for arbitrary $f$
(provided the expectation above exists), and since $f^{\rm eq}=f$ if $f$ is
equivariant, the map $f\to f^{\rm eq}$ is a linear projection in the space of
functions from $M_{n}^{\rm sym}(\mathbb{R})$ to itself, hence all equivariant
functions can be obtained by such a symmetrization.

It might be difficult to compute this expectation for generic $f$, or even to
ensure that the expectation is well-defined; we shall bypass this difficulty
by treating a special case, assuming that each matrix element of $f(Y)$ is a
polynomial of total degree $m$ in the matrix elements of $Y$. The computation
of $f^{\rm eq}$ amounts then to compute the expectation of a product of $2m+2$
matrix elements of an Haar distributed orthogonal matrix $O$, a problem which
has been solved in CoSn06 , with the following result:

$\mathbb{E}_{O}[O_{i_{1},j_{1}}\dots
O_{i_{2m+2},j_{2m+2}}]=\sum_{\mathfrak{p},\mathfrak{q}}{\rm
Wg}(n,\mathfrak{p},\mathfrak{q})\prod_{\\{a,b\\}\in\mathfrak{p}}\delta_{i_{a},i_{b}}\prod_{\\{a,b\\}\in\mathfrak{q}}\delta_{j_{a},j_{b}}\
,$ (147)

where the sum is over two independent pairings $\mathfrak{p},\mathfrak{q}$ of
$[2m+2]$, a pairing being a partition of $[2m+2]$ in which all blocks denoted
$\\{a,b\\}$ are of size 2, the Kronecker deltas imposing the equality of the
row (resp. column) indices $i$ (resp. $j$) in the blocks of the pairing
$\mathfrak{p}$ (resp. $\mathfrak{q})$. We shall not need here an explicit
expression of the so-called Weingarten function ${\rm
Wg}(n,\mathfrak{p},\mathfrak{q})$, the rest of the computation only relies on
the fact that the non-zero terms in the expectation are necessarily of this
form.

Let us now write a polynomial $f$ of degree $m$ as

$f(Y)_{i_{1},i_{2}}=\sum_{i_{3},\dots,i_{2m+2}}F_{i_{1},i_{2},i_{3},\dots,i_{2m+2}}Y_{i_{3},i_{4}}\dots
Y_{i_{2m+1},i_{2m+2}}\ ,$ (148)

where here and in the following we keep implicit the range $[n]$ of the
summations over indices $i_{a},j_{a}$, and where $F$ is an arbitrary array,
invariant under the exchange of its first two arguments for the image of $f$
to be a symmetric matrix. Inserting this definition in (146) and computing the
expectation thanks to (147) yields:

$\displaystyle f^{\rm eq}(Y)_{j_{1},j_{2}}$
$\displaystyle=\mathbb{E}_{O}[(O^{T}f(OYO^{T})O)_{j_{1},j_{2}}]$ (149)
$\displaystyle=\sum_{i_{1},i_{2}}\mathbb{E}_{O}[O_{i_{1},j_{1}}O_{i_{2},j_{2}}f(OYO^{T})_{i_{1},i_{2}}]$
(150)
$\displaystyle=\sum_{i_{1},\dots,i_{2m+2}}F_{i_{1},\dots,i_{2m+2}}\mathbb{E}_{O}[O_{i_{1},j_{1}}O_{i_{2},j_{2}}(OYO^{T})_{i_{3},i_{4}}\dots(OYO^{T})_{i_{2m+1},i_{2m+2}}]$
(151) $\displaystyle=\sum_{\begin{subarray}{c}i_{1},\dots,i_{2m+2}\\\
j_{3},\dots,j_{2m+2}\end{subarray}}F_{i_{1},\dots,i_{2m+2}}\mathbb{E}_{O}[O_{i_{1},j_{1}}\dots
O_{i_{2m+2},j_{2m+2}}]Y_{j_{3},j_{4}}\dots Y_{j_{2m+1},j_{2m+2}}$ (152)
$\displaystyle=\sum_{\mathfrak{p},\mathfrak{q}}{\rm
Wg}(n,\mathfrak{p},\mathfrak{q})\left(\sum_{i_{1},\dots,i_{2m+2}}F_{i_{1},\dots,i_{2m+2}}\prod_{\\{a,b\\}\in\mathfrak{p}}\delta_{i_{a},i_{b}}\right)$
(153) $\displaystyle\hskip
113.81102pt\left(\sum_{j_{3},\dots,j_{2p+2}}\left(\prod_{\\{a,b\\}\in\mathfrak{q}}\delta_{j_{a},j_{b}}\right)Y_{j_{3},j_{4}}\dots
Y_{j_{2m+1},j_{2m+2}}\right)\ .$ (154)

A moment of thought reveals that the last parenthesis is equal to
$(Y^{p})_{j_{1},j_{2}}({\rm Tr}(Y))^{q_{1}}\dots({\rm Tr}(Y^{m}))^{q_{m}}$,
where the non-negative integers $p,q_{1},\dots,q_{m}$ depend on the pairing
$\mathfrak{q}$, and add up to $m$, hence the dependency of $f^{\rm eq}$ on $Y$
is a linear combination of such terms. This is most easily seen through a
graphical reasoning illustrated on one example in Fig. 9, by representing the
pairing $\mathfrak{q}$ along with the canonical pairing
$\\{3,4\\},\dots,\\{2p+1,2p+2\\}$ complemented by two “external legs”
corresponding to the indices $1$ and $2$. The combination of these two
combinatorial objects yields a single path between the two external lengths of
length $p$, where length is counted in terms of the number of blocks crossed
in the canonical pairing, and an union of cycles, $q_{l}$ corresponding to the
number of cycles of length $l$. In this reasoning one uses crucially the
symmetry of $Y$, in its absence one would get arbitrary words in $Y$ and
$Y^{T}$ instead of simple matrix powers of $Y$. To fully recover the result of
Appendix B.1 one could use the theory of symmetric polynomials to express, for
$l>n$, ${\rm Tr}(Y^{l})$ as a polynomial in the ${\rm Tr}(Y),\dots,{\rm
Tr}(Y^{n})$, we skip these details which shall not be useful for us.

$j_{1}$$j_{2}$$Y$$Y$$Y$ Figure 9: An illustration of the contribution of Eq.
(154). In this figure $m=3$, the pairing
$\mathfrak{q}=\\{\\{1,6\\},\\{2,3\\},\\{4,5\\},\\{7,8\\}\\}$ is represented
above the horizontal axis, while below are the representation of the matrices
$Y$ as the blocks of the canonical pairing
$\\{\\{3,4\\},\\{5,6\\},\\{7,8\\}\\}$ and two “external legs”. This yields
here $(Y^{2})_{j_{1},j_{2}}{\rm Tr}(Y)$.

Note that a similar computation allows to characterize the polynomials
$f:M_{n}^{\rm sym}(\mathbb{R})\to\mathbb{R}$ that are invariant in the sense
that $f(OYO^{T})=f(Y)$ for all $O\in\mathcal{O}_{n}$ and $Y\in M_{n}^{\rm
sym}(\mathbb{R})$: starting from an arbitrary polynomial $f$ one computes its
symmetrized version $f^{\rm inv}(Y)=\mathbb{E}_{O}[f(OYO^{T})]$, that is
indeed invariant. Following the same steps as above shows that $f^{\rm inv}$
is a polynomial in the traces of the powers of $Y$ (one obtains the same
graphical representation without the external legs).

Let us finally mention that BlVi22  presents another method, attributed to
Malgrange, to obtain the equivariant functions $f:V\to W$ for generic group
actions on $V$ and $W$, via the determination of the invariant functions on
$V\times W^{*}$, with $W^{*}$ the dual of $W$. It should be possible to use it
to rederive the results presented here for the conjugation action of the
orthogonal group.

###  B.3 Symmetric group

We justify in this Appendix the claims made in Sec. V.3 on the form of the
polynomials equivariant under the symmetric group of permutations, adapting
the proof of Appendix B.2 (a similar reasoning can be found in PuLiKiMaLi23 ).
Consider indeed an arbitrary polynomial function as in (148), and its
symmetrization

$f^{\rm
eq}(Y)=\mathbb{E}_{\sigma}[O(\sigma)^{T}f(O(\sigma)YO(\sigma)^{T})O(\sigma)]\
,$ (155)

where the average is under the uniform measure on the $n!$ permutations of
$\mathcal{S}_{n}$, and $O(\sigma)$ is the associated matrix with elements
$O(\sigma)_{i,j}=\delta_{i,\sigma(j)}$. One can perform exactly the same steps
as in Eqs. (149-152), with the replacement of the Haar measure on
$\mathcal{O}_{n}$ by the one on $\mathcal{S}_{n}$. The computation thus boils
down to find the equivalent of (147), i.e. the Weingarten function of the
symmetric group instead of the orthogonal one, an exercise performed in
Section 3 of CoMaNo22  and that we reproduce here. This function turns out to
be much simpler, since

$\mathbb{E}_{\sigma}[O(\sigma)_{i_{1},j_{1}}\dots
O(\sigma)_{i_{2m+2},j_{2m+2}}]=\mathbb{P}_{\sigma}[i_{a}=\sigma(j_{a})\
\forall a\in[2m+2]]=\mathbb{P}_{\sigma}[i=\sigma\circ j]\ ,$ (156)

where in the second step one interprets the multiplet of indices
$(i_{1},\dots,i_{2m+2})$ as a function $i:[2m+2]\to[n]$, and similarly for
$j$. This expectation is thus the ratio between the number of permutations
$\sigma$ such that $i=\sigma\circ j$ and the total number $n!$ of
permutations. To complete the computation let us define the kernel $\ker(i)$
(resp. $\ker(j)$) as the partition of $[2m+2]$ induced by the equivalence
relation $a\sim b$ iff $i_{a}=i_{b}$ (resp. $j_{a}=j_{b}$), in other words the
blocks of $\ker(i)$ encode the patterns of equalities in the indices
$(i_{1},\dots,i_{2m+2})$. A moment of thought reveals that one can find a
permutation $\sigma$ with $i=\sigma\circ j$ if and only if $\ker(i)=\ker(j)$,
and that in this case there are $(n-|\ker(i)|)!$ satisfying permutations, with
$|\ker(i)|$ the number of blocks in the partition: the value of $\sigma$ on
the image of each block of $j$ is fixed to match the corresponding values in
$i$, and $\sigma$ is a bijection from the indices absent from $j$ to those
absent from $i$. In formula this reads

$\mathbb{E}_{\sigma}[O(\sigma)_{i_{1},j_{1}}\dots
O(\sigma)_{i_{2m+2},j_{2m+2}}]=\mathds{1}(\ker(i)=\ker(j))\frac{(n-|\ker(i)|)!}{n!}\
,$ (157)

which inserted in (152) yields

$f^{\rm
eq}(Y)_{j_{1},j_{2}}=\sum_{i_{1},\dots,i_{2m+2}}F_{i_{1},\dots,i_{2m+2}}\frac{(n-|\ker(i)|)!}{n!}\sum_{\begin{subarray}{c}j_{3},\dots,j_{2p+2}\\\
\ker(j)=\ker(i)\end{subarray}}Y_{j_{3},j_{4}}\dots Y_{j_{2m+1},j_{2m+2}}\ .$
(158)

One can conclude from here that such a function is indeed a linear combination
of the $\widehat{b}_{G}$ defined in (110): for each choice of $i$ consider the
graph $G$ with $|\ker(i)|$ vertices associated to the blocks of $\ker(i)$, $m$
edges between the vertices $B(i_{2p+1})$ and $B(i_{2p+2})$, for $p=1,\dots,m$,
where $B(i_{a})$ denotes the block containing $i_{a}$, and mark the vertices
$B(i_{1})$ and $B(i_{2})$ (which can be equal). The summation over $j$ in
(158) is then equal to $(\widehat{b}^{\prime}_{G})_{j_{1},j_{2}}$. Moreover
the invariance of the array $F$ under the exchange of its first two indices
imply the symmetrization of (111).

###  B.4 Inversion symmetry

We turn now to the consequences of the inversion symmetry that appears when
the laws of $X_{1,1}$ and $B_{1,2}$ are invariant under a change of sign.
Consider the group $\\{-1,1\\}^{n}$, with componentwise multiplication, and
let it act on $M_{n}^{\rm sym}(\mathbb{R})$ as $(\varepsilon\cdot
Y)_{i,j}=\varepsilon_{i}\varepsilon_{j}Y_{i,j}$ for
$\varepsilon=(\varepsilon_{1},\dots,\varepsilon_{n})\in\\{-1,1\\}^{n}$, which
can be seen as the conjugation with the orthogonal matrix ${\rm
diag}(\varepsilon_{1},\dots,\varepsilon_{n})$. The assumptions on the laws of
$X_{1,1}$ and $B_{1,2}$ imply that this action is a symmetry for the inference
problem. We can therefore symmetrize the permutation invariant polynomials in
order to further reduce the set of relevant basic functions. Since
$\varepsilon^{-1}=\varepsilon$, the symmetrization of
$\widehat{b}^{\prime}_{G}$ of (110) yields

$\mathbb{E}_{\varepsilon}\left[(\varepsilon\cdot\widehat{b}^{\prime}_{G}(\varepsilon\cdot
Y))_{i,j}\right]=\sum_{\begin{subarray}{c}\phi\in\mathcal{I}(V,[n])\\\
\phi(v)=i,\phi(w)=j\end{subarray}}\left(\prod_{e=\\{a,b\\}\in
E}Y_{\phi(a),\phi(b)}\right)\mathbb{E}_{\varepsilon}\left[\varepsilon_{i}\varepsilon_{j}\prod_{e=\\{a,b\\}\in
E}\varepsilon_{\phi(a)}\varepsilon_{\phi(b)}\right]\ ,$ (159)

where $\mathbb{E}_{\varepsilon}$ denotes the uniform average over the $2^{n}$
elements of $\\{-1,1\\}^{n}$. The expectation in the right hand side vanishes
unless for all $k\in[n]$, the power of $\varepsilon_{k}$ is even. If $k$ is
not in the image of $\phi$ this power is zero; if it is, the power arising
from the product over the edges is the degree of the vertex $a=\phi^{-1}(k)$
(here the injectivity of $\phi$ is crucial), while each of the marked vertices
(possibly equal) add 1 to the corresponding powers. Hence the average over
$\varepsilon$ is independent of $\phi$, and is equal to 1 if the degree
constraints stated at the end of Sec. V.3 are fulfilled, 0 otherwise.

##  Appendix C Gaussian integration by parts

This Appendix presents an alternative derivation of the expression (50) of
$\widehat{\mathcal{R}}_{p}$ in the Gaussian noise case, avoiding the free
probability combinatorial interpretation of the main text.

Before starting the derivation let us recall two well-known facts:

  * •

if $G=(G_{1},\dots,G_{N})$ is a vector of independent centered Gaussian
variables, an integration by parts reveals that

$\mathbb{E}[G_{i}f(G)]=\mathbb{E}[G_{i}^{2}]\mathbb{E}\left[\frac{\partial
f}{\partial G_{i}}\right]\ ,$ (160)

whenever the function $f$ is smooth enough and behaves at infinity in such a
way that both members of this equality exist.

  * •

if $A(t)$ is a matrix dependent on a parameter $t$, the derivatives of its
matrix powers read

$\frac{{\rm d}}{{\rm d}t}A(t)^{p}=\sum_{k=1}^{p}A(t)^{k-1}\frac{{\rm d}A}{{\rm
d}t}A(t)^{p-k}\ .$ (161)

We have to compute the value of $\frac{1}{n}\mathbb{E}[{\rm Tr}(ZY^{p})]$,
where $Y=S+Z$, with $S$ and $Z$ two independent symmetric matrices of size
$n$, and specifically in the Gaussian noise case the matrix elements
$\\{Z_{i,j}\\}_{i\leq j}$ are independent centered Gaussian random variables,
with $\mathbb{E}[(Z_{i,i})^{2}]=\frac{2\Delta}{n}$, and
$\mathbb{E}[(Z_{i,j})^{2}]=\frac{\Delta}{n}$ for $i<j$. The matrix elements
$Z_{i,j}$ for $i>j$ are given as $Z_{i,j}=Z_{j,i}$, to ensure the symmetry. We
thus write the above expression as

$\displaystyle\frac{1}{n}\mathbb{E}[{\rm Tr}(ZY^{p})]=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}[Z_{i,i}(Y^{p})_{i,i}]+\frac{2}{n}\sum_{1\leq
i<j\leq n}\mathbb{E}[Z_{i,j}(Y^{p})_{i,j}]$ (162) $\displaystyle=$
$\displaystyle\frac{2\Delta}{n^{2}}\sum_{i=1}^{n}\mathbb{E}\left[\frac{\partial}{\partial
Z_{i,i}}(Y^{p})_{i,i}\right]+\frac{2\Delta}{n^{2}}\sum_{1\leq i<j\leq
n}\mathbb{E}\left[\frac{\partial}{\partial Z_{i,j}}(Y^{p})_{i,j}\right]$ (163)
$\displaystyle=$
$\displaystyle\frac{2\Delta}{n^{2}}\sum_{k=1}^{p}\left(\sum_{i=1}^{n}\mathbb{E}\left[\left(Y^{k-1}\frac{\partial
Y}{\partial Z_{i,i}}Y^{p-k}\right)_{i,i}\right]+\sum_{1\leq i<j\leq
n}\mathbb{E}\left[\left(Y^{k-1}\frac{\partial Y}{\partial
Z_{i,j}}Y^{p-k}\right)_{i,j}\right]\right)\ ,$ (164)

by applying the two facts recalled above. Since $\left(\frac{\partial
Y}{\partial Z_{i,i}}\right)_{l,m}=\delta_{l,i}\delta_{m,i}$ and for $i<j$
$\left(\frac{\partial Y}{\partial
Z_{i,j}}\right)_{l,m}=\delta_{l,i}\delta_{m,j}+\delta_{l,j}\delta_{m,i}$,
because of the symmetry requirement, one has

$\displaystyle\frac{1}{n}\mathbb{E}[{\rm Tr}(ZY^{p})]=$
$\displaystyle\frac{2\Delta}{n^{2}}\sum_{k=1}^{p}\left(\sum_{i=1}^{n}\mathbb{E}\left[(Y^{k-1})_{i,i}(Y^{p-k})_{i,i}\right]+\sum_{1\leq
i<j\leq n}\mathbb{E}\left[(Y^{k-1})_{i,i}(Y^{p-k})_{j,j}\right]+\sum_{1\leq
i<j\leq n}\mathbb{E}\left[(Y^{k-1})_{i,j}(Y^{p-k})_{j,i}\right]\right)\ ,$
$\displaystyle=$
$\displaystyle\frac{\Delta}{n^{2}}\sum_{k=1}^{p}\left(\sum_{i,j=1}^{n}\mathbb{E}\left[(Y^{k-1})_{i,i}(Y^{p-k})_{j,j}\right]+\sum_{i,j=1}^{n}\mathbb{E}\left[(Y^{k-1})_{i,j}(Y^{p-k})_{j,i}\right]\right)$
(165) $\displaystyle=$ $\displaystyle\Delta\sum_{k=1}^{p}\mathbb{E}[{\rm
tr}(Y^{k-1}){\rm tr}(Y^{p-k})]+\frac{1}{n}\Delta p\mathbb{E}[{\rm
tr}(Y^{p-1})]\ ,$ (166)

with ${\rm tr}(\bullet)={\rm Tr}(\bullet)/n$ the normalized trace. In the
large $n$ limit the assumptions of convergence of the average moments and
concentration expressed in Eqs. (1,2) imply that the second term above is of
order $1/n$ and that the first one factorizes, yielding

$\lim_{n\to\infty}\frac{1}{n}\mathbb{E}[{\rm
Tr}(ZY^{p})]=\Delta\sum_{k=1}^{p}\mu_{Y,k-1}\mu_{Y,p-k}\ .$ (167)

This agrees indeed with (50), after the change of variable $j=k-1$.

##  Appendix D Degree 2 denoisers for arbitrary Wishart matrices perturbed by
Wigner noise

We present in this Appendix the details on the computation whose results were
stated in Sec. V.4.

###  D.1 Off-diagonal terms

We start with the determination of the matrix $\mathcal{M}^{(n,2)}$ defined in
Eq. (118), with the functions $b_{1},\dots,b_{4}$ of Eq. (116). A first step
is to reduce its computation to those of joint moments of the matrix elements
of $Y$, exploiting its invariance under permutation to group together
equivalent terms:

$\displaystyle\mathcal{M}_{1,1}^{(n,2)}=$
$\displaystyle(n-1)\mathbb{E}[Y_{1,2}^{2}]\ ,$
$\displaystyle\mathcal{M}_{1,2}^{(n,2)}=$ $\displaystyle
2(n-1)\mathbb{E}[Y_{1,1}Y_{1,2}^{2}]+(n-1)(n-2)\mathbb{E}[Y_{1,2}Y_{2,3}Y_{3,1}]\
,$ $\displaystyle\mathcal{M}_{1,3}^{(n,2)}=$ $\displaystyle
2(n-1)\mathbb{E}[Y_{1,1}Y_{1,2}^{2}]\ ,$
$\displaystyle\mathcal{M}_{1,4}^{(n,2)}=$
$\displaystyle\sqrt{n}(n-1)\mathbb{E}[Y_{1,1}Y_{1,2}^{2}]\ ,$
$\displaystyle\mathcal{M}_{2,2}^{(n,2)}=$ $\displaystyle
2(n-1)\mathbb{E}[Y_{1,1}^{2}Y_{1,2}^{2}]+2(n-1)\mathbb{E}[Y_{1,1}Y_{1,2}^{2}Y_{2,2}]+4(n-1)(n-2)\mathbb{E}[Y_{1,1}Y_{1,2}Y_{2,3}Y_{3,1}]$
$\displaystyle+(n-1)(n-2)\mathbb{E}[Y_{1,2}^{2}Y_{2,3}^{2}]+(n-1)(n-2)(n-3)\mathbb{E}[Y_{1,2}Y_{2,3}Y_{3,4}Y_{4,1}]\
,$ $\displaystyle\mathcal{M}_{2,3}^{(n,2)}=$ $\displaystyle
2(n-1)\mathbb{E}[Y_{1,1}^{2}Y_{1,2}^{2}]+2(n-1)\mathbb{E}[Y_{1,1}Y_{1,2}^{2}Y_{2,2}]+3(n-1)(n-2)\mathbb{E}[Y_{1,1}Y_{1,2}Y_{2,3}Y_{3,1}]\
,$ $\displaystyle\mathcal{M}_{2,4}^{(n,2)}=$
$\displaystyle\sqrt{n}(n-1)\mathbb{E}[Y_{1,1}^{2}Y_{1,2}^{2}]+\sqrt{n}(n-1)\mathbb{E}[Y_{1,1}Y_{1,2}^{2}Y_{2,2}]+\sqrt{n}(n-1)(n-2)\mathbb{E}[Y_{1,1}Y_{1,2}Y_{2,3}Y_{3,1}]\
,$ $\displaystyle\mathcal{M}_{3,3}^{(n,2)}=$ $\displaystyle
2(n-1)\mathbb{E}[Y_{1,1}^{2}Y_{1,2}^{2}]+2(n-1)\mathbb{E}[Y_{1,1}Y_{1,2}^{2}Y_{2,2}]+(n-1)(n-2)\mathbb{E}[Y_{1,1}^{2}]\mathbb{E}[Y_{1,2}^{2}]\
,$ $\displaystyle\mathcal{M}_{3,4}^{(n,2)}=$
$\displaystyle\sqrt{n}(n-1)\mathbb{E}[Y_{1,1}^{2}Y_{1,2}^{2}]+\sqrt{n}(n-1)\mathbb{E}[Y_{1,1}Y_{1,2}^{2}Y_{2,2}]\
,$ $\displaystyle\mathcal{M}_{4,4}^{(n,2)}=$
$\displaystyle\frac{1}{2}n(n-1)\mathbb{E}[Y_{1,1}^{2}Y_{1,2}^{2}]+\frac{1}{2}n(n-1)\mathbb{E}[Y_{1,1}Y_{1,2}^{2}Y_{2,2}]\
.$ (168)

Let us explain the main ideas underlying the reasonings followed to obtain
these expressions, taking for instance the term $\mathcal{M}_{1,2}^{(n,2)}$;
starting from its definition,

$\mathcal{M}_{1,2}^{(n,2)}=(n-1)\mathbb{E}[Y_{1,2}(Y^{2})_{1,2}]=(n-1)\sum_{k=1}^{n}\mathbb{E}[Y_{1,2}Y_{1,k}Y_{k,2}]\
.$ (169)

Using the permutation invariance the two terms with $k\in\\{1,2\\}$ are both
equal to $\mathbb{E}[Y_{1,1}Y_{1,2}^{2}]$, while the $n-2$ summands with
$k\geq 3$ all yield $\mathbb{E}[Y_{1,2}Y_{2,3}Y_{3,1}]$. All the other lines
in (168) are obtained with similar combinatorial reasonings, enumerating the
number of contributions with equivalent patterns of indices under the
permutation symmetry, and exploiting the following independence property: if
$i_{1},\dots,i_{2p}$ are indices in $I\subset[n]$, and $j_{1},\dots,j_{2q}$
are in another subset $J$ disjoint from $I$, then
$\mathbb{E}[Y_{i_{1},i_{2}}\dots Y_{i_{2p-1},i_{2p}}Y_{j_{1},j_{2}}\dots
Y_{j_{2p-1},j_{2p}}]=\mathbb{E}[Y_{i_{1},i_{2}}\dots
Y_{i_{2p-1},i_{2p}}]\mathbb{E}[Y_{j_{1},j_{2}}\dots Y_{j_{2p-1},j_{2p}}]$.
Indeed the noise random variables $B_{k,l}$ are independent from one edge to
the other, while the signal matrix elements $S_{k,l}$ are correlated but only
through their dependency on $X_{k,\mu}$ and $X_{l,\mu}$. As a consequence the
matrix elements $S_{i_{a},i_{b}}$ only depend on $\\{X_{i,\mu}\\}_{i\in I}$
and the $S_{j_{a},j_{b}}$ on the $\\{X_{j,\mu}\\}_{j\in J}$, hence the
independence at the origin of this factorization.

Similar considerations yield the following expressions for the elements of the
vector $\mathcal{R}^{(n,2)}$ also defined in Eq. (118):

$\displaystyle\mathcal{R}_{1}^{(n,2)}=$
$\displaystyle(n-1)\mathbb{E}[S_{1,2}Y_{1,2}]\ ,$
$\displaystyle\mathcal{R}_{2}^{(n,2)}=$ $\displaystyle
2(n-1)\mathbb{E}[S_{1,2}Y_{1,1}Y_{1,2}]+(n-1)(n-2)\mathbb{E}[S_{1,2}Y_{2,3}Y_{3,1}]\
,$ $\displaystyle\mathcal{R}_{3}^{(n,2)}=$ $\displaystyle
2(n-1)\mathbb{E}[S_{1,2}Y_{1,1}Y_{1,2}]\ ,$
$\displaystyle\mathcal{R}_{4}^{(n,2)}=$
$\displaystyle\sqrt{n}(n-1)\mathbb{E}[S_{1,2}Y_{1,1}Y_{1,2}]\ .$ (170)

To complete the computation it remains to determine more explicitly the joint
moments of the matrix elements of $Y$. To do so one expands each $Y_{i,j}$ as
$S_{i,j}+Z_{i,j}$; since the $Z$ are independent from the $S$, and independent
from one edge to the other, one can easily perform the average on the former.
To average over the $S$ it is necessary to express them in terms of the
$X_{i,\mu}$, according to the definition (107), and then perform the averages
exploiting the i.i.d. character of $X_{i,\mu}$, with again some combinatorial
considerations to count the number of equivalent terms. The first order terms
trivially vanish, $\mathbb{E}[Y_{1,1}]=\mathbb{E}[Y_{1,2}]=0$. The moments of
order 2 that appear in (168) are

$\mathbb{E}[Y_{1,1}^{2}]=\frac{1}{n}(\mathbb{E}[X_{1,1}^{4}]-1+\Delta\mathbb{E}[B_{1,1}^{2}])\
,\qquad\mathbb{E}[Y_{1,2}^{2}]=\frac{1}{n}(1+\Delta)\ .$ (171)

Similarly one finds for the required third order terms:

$\mathbb{E}[Y_{1,1}Y_{1,2}^{2}]=\frac{1}{n^{2}}\sqrt{\alpha}(\mathbb{E}[X_{1,1}^{4}]-1)\
,\qquad\mathbb{E}[Y_{1,2}Y_{2,3}Y_{3,1}]=\frac{1}{n^{2}}\sqrt{\alpha}\ ,$
(172)

and finally for the fourth order terms:

$\displaystyle\mathbb{E}[Y_{1,1}^{2}Y_{1,2}^{2}]$
$\displaystyle=\frac{1}{n^{2}}(1+\Delta)(\mathbb{E}[X_{1,1}^{4}]-1+\Delta\mathbb{E}[B_{1,1}^{2}])+\frac{1}{n^{3}}\alpha(\mathbb{E}[X_{1,1}^{6}]-3\mathbb{E}[X_{1,1}^{4}]+2)\
,$ (173) $\displaystyle\mathbb{E}[Y_{1,1}Y_{1,2}^{2}Y_{2,2}]$
$\displaystyle=\frac{1}{n^{3}}\alpha(\mathbb{E}[X_{1,1}^{4}]-1)^{2}\ ,$ (174)
$\displaystyle\mathbb{E}[Y_{1,1}Y_{1,2}Y_{2,3}Y_{3,1}]$
$\displaystyle=\frac{1}{n^{3}}\alpha(\mathbb{E}[X_{1,1}^{4}]-1)\ ,$ (175)
$\displaystyle\mathbb{E}[Y_{1,2}^{2}Y_{2,3}^{2}]$
$\displaystyle=\frac{1}{n^{2}}(1+\Delta)^{2}+\frac{1}{n^{3}}\alpha(\mathbb{E}[X_{1,1}^{4}]-1)\
,$ (176) $\displaystyle\mathbb{E}[Y_{1,2}Y_{2,3}Y_{3,4}Y_{4,1}]$
$\displaystyle=\frac{1}{n^{3}}\alpha\ .$ (177)

The expressions for the joint moments between $S$ and $Y$ that appear in the
expression (170) of $\mathcal{R}^{(n,2)}$ are established through a similar
expansion and read

$\mathbb{E}[S_{1,2}Y_{1,2}]=\frac{1}{n}\
,\qquad\mathbb{E}[S_{1,2}Y_{1,1}Y_{1,2}]=\frac{1}{n^{2}}\sqrt{\alpha}(\mathbb{E}[X_{1,1}^{4}]-1)\
,\qquad\mathbb{E}[S_{1,2}Y_{2,3}Y_{3,1}]=\frac{1}{n^{2}}\sqrt{\alpha}\ .$
(178)

Collecting all these relations gives explicit expressions for the matrix
$\mathcal{M}^{(n,2)}$ and vector $\mathcal{R}^{(n,2)}$ in terms of the
parameters $\alpha$, $\Delta$, $n$, and the non-universal quantities
$\mathbb{E}[X_{1,1}^{4}]$, $\mathbb{E}[X_{1,1}^{6}]$ and
$\mathbb{E}[B_{1,1}^{2}]$. Taking their large $n$ limit then yields the
expression (119) of the main text; with the help of a symbolic computing
software we checked that the solution of
$\mathcal{M}^{(n,2)}c=\mathcal{R}^{(n,2)}$ has $c_{4}=0$ when plugging the
values $\mathbb{E}[X_{1,1}^{4}]=3$, $\mathbb{E}[X_{1,1}^{6}]=15$ and
$\mathbb{E}[B_{1,1}^{2}]=2$, restoring as it should the orthogonal invariance
of the optimal estimator when the priors on the signal and noise are invariant
under the orthogonal group. The off-diagonal contribution to the MMSE reads
${\rm MMSE}_{\rm
od}^{(n,2)}=\left(1-\frac{1}{n}\right)-c^{T}\mathcal{R}^{(n,2)}$, with $c$ the
solution of $\mathcal{M}^{(n,2)}c=\mathcal{R}^{(n,2)}$. With a symbolic
expansion in the large $n$ limit we obtained the $1/n$ term stated in Eq.
(123), as well as the $1/n$ corrections on the coefficients $c_{1}$ and
$c_{2}$, and the dominant contribution for the coefficients of the optimal
estimator that vanished in the large $n$ limit:

$\displaystyle c_{3}^{(n)}$
$\displaystyle=-\frac{1}{n}\frac{\alpha^{\frac{3}{2}}\Delta(\mathbb{E}[X_{1,1}^{4}]-1)}{(1+\Delta)((1+\Delta)^{3}+\alpha\Delta)(\mathbb{E}[X_{1,1}^{4}]-1+\Delta\mathbb{E}[B_{1,1}^{2}])}+O\left(\frac{1}{n^{2}}\right)\
,$ (179) $\displaystyle c_{4}^{(n)}$
$\displaystyle=\frac{1}{\sqrt{n}}\frac{2\sqrt{\alpha}\Delta^{2}(\mathbb{E}[X_{1,1}^{4}]-1-\mathbb{E}[B_{1,1}^{2}])}{((1+\Delta)^{3}+\alpha\Delta)(\mathbb{E}[X_{1,1}^{4}]-1+\Delta\mathbb{E}[B_{1,1}^{2}])}+O\left(\frac{1}{n^{3/2}}\right)\
.$ (180)

###  D.2 Diagonal terms

We shall now explain the main steps that brought us to the expansion (124) for
the contributions of the diagonal terms to the MMSE. According to the rules
explained in Sec. V.3 there are 9 graphs, represented on Fig. 10 that
contribute: these are indeed the unlabelled graphs with one marked vertex
(since we consider the diagonal terms of the estimator), at most two edges
(since we only consider polynomials of degree at most two), no isolated
unmarked vertex, and all vertices with even degrees (assuming the inversion
symmetry).

= ,= ,= ,= ,= ,= ,= ,= ,= . Figure 10: The nine graphs in $\mathcal{A}_{\rm
d}^{(2)}$ that contribute to the diagonal terms of the estimators at order 2,
in presence of the inversion symmetry. White circles represent the marked
vertex, black ones the unmarked vertices.

These estimators $b_{\beta}(Y)$ have by definition $(b_{\beta}(Y))_{i,j}=0$
for $i\neq j$, and their respective values for $(b_{\beta}(Y))_{i,i}$ read:

$1\ ,\quad Y_{i,i}\ ,\quad\sum_{j\neq i}Y_{i,j}^{2}\ ,\quad Y_{i,i}^{2}\
,\quad\sum_{j\neq i}Y_{j,j}\ ,\quad\sum_{j\neq i}Y_{j,j}^{2}\ ,\quad
Y_{i,i}\sum_{j\neq i}Y_{j,j}\ ,\quad\sum_{\begin{subarray}{c}j,k\neq i\\\
j\neq k\end{subarray}}Y_{j,k}^{2}\ ,\quad\sum_{\begin{subarray}{c}j,k\neq i\\\
j\neq k\end{subarray}}Y_{j,j}Y_{k,k}\ .$ (181)

One is thus led to consider a matrix $\mathcal{M}^{(n,2,{\rm d})}$ and a
vector $\mathcal{R}^{(n,2,{\rm d})}$, of dimension 9, defined as
$\mathcal{M}_{\beta,\beta^{\prime}}^{(n,2,{\rm
d})}=\mathbb{E}[(b_{\beta}(Y))_{1,1}(b_{\beta^{\prime}}(Y))_{1,1}]$ and
$\mathcal{R}_{\beta}^{(n,2,{\rm d})}=\mathbb{E}[S_{1,1}(b_{\beta}(Y))_{1,1}]$.
Their matrix elements can be expressed in terms of joint moments of $Y$ with
the same kind of computations as those detailed in the off-diagonal case.
There are some additional moments that appear:

$\displaystyle\mathbb{E}[Y_{1,1}^{3}]=$
$\displaystyle\frac{1}{n^{2}}\sqrt{\alpha}(\mathbb{E}[X_{1,1}^{6}]-3\,\mathbb{E}[X_{1,1}^{4}]+2)\
,$ (182) $\displaystyle\mathbb{E}[Y_{1,1}^{4}]=$
$\displaystyle\frac{1}{n^{2}}(3(\mathbb{E}[X_{1,1}^{4}]-1)^{2}+6\Delta(\mathbb{E}[X_{1,1}^{4}]-1)\mathbb{E}[B_{1,1}^{2}]+\Delta^{2}\mathbb{E}[B_{1,1}^{4}])$
(183)
$\displaystyle+\frac{1}{n^{3}}\alpha(\mathbb{E}[X_{1,1}^{8}]-3\,\mathbb{E}[X_{1,1}^{4}]^{2}-4\,\mathbb{E}[X_{1,1}^{6}]+12\,\mathbb{E}[X_{1,1}^{4}]-6)\
,$ (184) $\displaystyle\mathbb{E}[Y_{1,2}^{4}]=$
$\displaystyle\frac{1}{n^{2}}(3+6\Delta+\Delta^{2}\mathbb{E}[B_{1,1}^{4}])+\frac{1}{n^{3}}\alpha(\mathbb{E}[X_{1,1}^{4}]^{2}-3)\
,$ (185) $\displaystyle\mathbb{E}[S_{1,1}Y_{1,1}]=$
$\displaystyle\frac{1}{n}(\mathbb{E}[X_{1,1}^{4}]-1)\ ,$ (186)
$\displaystyle\mathbb{E}[S_{1,1}Y_{1,1}^{2}]=$
$\displaystyle\frac{1}{n^{2}}\sqrt{\alpha}(\mathbb{E}[X_{1,1}^{6}]-3\,\mathbb{E}[X_{1,1}^{4}]+2)\
,$ (187) $\displaystyle\mathbb{E}[S_{1,1}Y_{1,2}^{2}]=$
$\displaystyle\frac{1}{n^{2}}\sqrt{\alpha}(\mathbb{E}[X_{1,1}^{4}]-1)\ .$
(188)

Following the steps encountered several times before one computes the optimal
diagonal MSE as

${\rm MMSE}_{\rm
d}^{(n,2)}=\mathbb{E}[S_{1,1}^{2}]-c^{T}\mathcal{R}^{(n,2,{\rm
d})}=\frac{1}{n}(\mathbb{E}[X_{1,1}^{4}]-1)-c^{T}\mathcal{R}^{(n,2,{\rm d})}\
,$ (189)

where $c$ is the 9-dimensional vector solution of $\mathcal{M}^{(n,2,{\rm
d})}c=\mathcal{R}^{(n,2,{\rm d})}$. With the help of a symbolic computing
software we thus obtained the expansion stated in Eq. (124). We spare the
reader the details of this nine-dimensional optimization problem since
actually the leading $1/n$ order can be determined in a simpler way, working
with the reduced set of diagonal basic functions $(b_{1}(Y))_{i,i}=1$,
$(b_{2}(Y))_{i,i}=Y_{i,i}$, $(b_{3}(Y))_{i,i}=(Y^{2})_{i,i}$ (the discarded
six-dimensional variational space only contribute to the order $1/n^{2}$ in
the expansion of ${\rm MMSE}_{\rm d}^{(n,2)}$). The elements of the
corresponding matrix and vector read

$\displaystyle\mathcal{M}_{1,1}$ $\displaystyle=1\ ,$ (190)
$\displaystyle\mathcal{M}_{1,2}$ $\displaystyle=0\ ,$ (191)
$\displaystyle\mathcal{M}_{1,3}$
$\displaystyle=\mathbb{E}[Y_{1,1}^{2}]+(n-1)\mathbb{E}[Y_{1,2}^{2}]\ ,$ (192)
$\displaystyle\mathcal{M}_{2,2}$ $\displaystyle=\mathbb{E}[Y_{1,1}^{2}]\ ,$
(193) $\displaystyle\mathcal{M}_{2,3}$
$\displaystyle=\mathbb{E}[Y_{1,1}^{3}]+(n-1)\mathbb{E}[Y_{1,1}Y_{1,2}^{2}]\ ,$
(194) $\displaystyle\mathcal{M}_{3,3}$
$\displaystyle=\mathbb{E}[Y_{1,1}^{4}]+2(n-1)\mathbb{E}[Y_{1,1}^{2}Y_{1,2}^{2}]+(n-1)\mathbb{E}[Y_{1,2}^{4}]+(n-1)(n-2)\mathbb{E}[Y_{1,2}^{2}Y_{2,3}^{2}]\
,$ (195) $\displaystyle\mathcal{R}_{1}$ $\displaystyle=0\ ,$ (196)
$\displaystyle\mathcal{R}_{2}$ $\displaystyle=\mathbb{E}[S_{1,1}Y_{1,1}]\ ,$
(197) $\displaystyle\mathcal{R}_{3}$
$\displaystyle=\mathbb{E}[S_{1,1}Y_{1,1}^{2}]+(n-1)\mathbb{E}[S_{1,1}Y_{1,2}^{2}]\
.$ (198)

The expressions of the moments of $Y$ stated above allow to give a completely
explicit form for the three-dimensional system of equations, with again an
expansion in the large $n$ limit obtained symbolically that yields Eq. (124).

##  Appendix E Additional details on the universality conjecture for arbitrary
$D$

We continue here the discussion of Sec. V.6 and explain more precisely the
concrete implementation of the computation. One missing ingredient was the
definition of the estimators $b_{G}(Y)$, a possible choice is constructed as
follows. For $H=(V,E)$ an (unmarked) graph (with possibly multiple edges and
self-loops), let us denote

$T_{H}(Y)=\sum_{\phi\in\mathcal{F}(V,[n])}\prod_{e=\\{a,b\\}\in
E}Y_{\phi(a),\phi(b)}\ ,$ (199)

where $\mathcal{F}(E,F)$ is the set of functions from $E$ to $F$. This
quantity is called the trace of the graph $H$, labelled by matrices $Y$, in
Ma20 ; CeDaMa16 . If $G=(V,E,v,w)$ is a marked graph we shall call ${\cal
C}_{\rm u}$ the set of its connected components that do not contain any marked
vertex, while $V_{\rm m}$ and $E_{\rm m}$ will denote the vertices and edges
of the connected components that contain at least one marked vertex. We then
define

$\displaystyle(b^{\prime}_{G}(Y))_{i,j}$
$\displaystyle=n^{x_{G}}\left(\sum_{\begin{subarray}{c}\phi\in\mathcal{F}(V_{\rm
m},[n])\\\ \phi(v)=i,\phi(w)=j\end{subarray}}\prod_{e=\\{a,b\\}\in E_{\rm
m}}Y_{\phi(a),\phi(b)}\right)\prod_{H\in{\cal C}_{\rm
u}}(T_{H}(Y)-\mathbb{E}[T_{H}(Y)])\ ,$ (200) $\displaystyle(b_{G}(Y))_{i,j}$
$\displaystyle=\frac{1}{2}((b^{\prime}_{G}(Y))_{i,j}+(b^{\prime}_{G}(Y))_{j,i})\
,$ (201)

with $x_{G}$ an exponent discussed later on. The linear combinations of such
$b_{G}$, with $G\in\mathcal{A}^{(D)}$, span all the permutation equivariant
polynomials of degree at most $D$. We call good graphs those which decompose
in terms of their connected components as a path of $p$ edges with the two
endvertices marked, and for $i=1,\dots,D$ $q_{i}$ simple cycles of $i$ edges,
and set $x_{G}=0$ for these good graphs. This definition coincides precisely
with the one of Eq. (79), thus the $\\{b_{G}\\}_{G\in\mathcal{A}^{(D,{\rm
g})}}$ span indeed the equivariant functions under the orthogonal group; all
other graphs are placed in the bad set $\mathcal{A}^{(D,{\rm b})}$. Note that
the $b_{G}$ associated to the good graphs were seen in Sec. IV to be properly
normalized, in the sense that the block $\mathcal{M}^{(n,D,{\rm gg})}$ admits
a finite invertible limit as $n\to\infty$.

Let us now discuss the choice of the exponent $x_{G}$ for the bad graphs, that
should yield a finite limit for the matrix elements of $\mathcal{M}^{(n,D)}$
of $\mathcal{R}^{(n,D)}$, and the invertibility of $\mathcal{M}^{(\infty,D)}$
(recall the objective set in Eq. (130)); a necessary (but not sufficient)
condition for the invertibility of a Gram matrix is that its diagonal elements
are strictly positive. We will argue shortly that there exists indeed a choice
of $x_{G}$ such that $\mathcal{M}^{(\infty,D)}_{G,G}>0$ (for the only bad
graph at $D=2$ one had $x_{G}=1/2$, see the definition of $b_{4}$ in (116));
note meanwhile that the Cauchy-Schwartz inequality gives

$|\mathcal{M}_{G,G^{\prime}}^{(n,D)}|\leq\sqrt{\mathcal{M}_{G,G}^{(n,D)}}\sqrt{\mathcal{M}_{G^{\prime},G^{\prime}}^{(n,D)}}\
\ ,\qquad|\mathcal{R}_{G}^{(n,D)}|\leq\sqrt{\mathbb{E}[{\rm
tr}((S^{(n)})^{2})]}\sqrt{\mathcal{M}_{G,G}^{(n,D)}}\ ,$ (202)

hence this choice of $x_{G}$ implies that the off-diagonal elements of
$\mathcal{M}^{(n,D)}$ and those of $\mathcal{R}^{(n,D)}$ are bounded (and the
forthcoming argument will show that they admit a finite limit as
$n\to\infty$).

A moment of thought reveals that $\mathcal{M}^{(n,D)}_{G,G^{\prime}}$ is of
the form

$n^{x_{G}+x_{G^{\prime}}-1}\sum_{i=1}^{s}C_{i}\prod_{j=1}^{t_{i}}\mathbb{E}[T_{H_{i,j}}(Y)]\
,$ (203)

where $s$, $\\{C_{i},t_{i}\\}$ and the unmarked graphs $\\{H_{i,j}\\}$ depends
on $(G,G^{\prime})$. Indeed when computing $\mathbb{E}[{\rm
tr}(b_{G}(Y)b_{G^{\prime}}(Y))]$ one has to join the two marked graphs $G$ and
$G^{\prime}$ by identifying their marked vertices, which yields an unmarked
graph; because of the centering of the contribution of the unmarked components
in (201) the expansion of the product gives rise to the form (203) for some
real coefficients $C_{i}$. We shall see that $\mathbb{E}[T_{H}(Y)]=n^{y_{H}}$
times a polynomial in $1/\sqrt{n}$, with $y_{H}$ a well-defined exponent,
hence $\mathcal{M}^{(n,D)}_{G,G}$ has an algebraic large $n$ behavior with an
exponent $2x_{G}-1+\max_{i}\sum_{j}y_{H_{i,j}}$. One can thus set $x_{G}$ in
such a way that the exponent of $\mathcal{M}^{(n,D)}_{G,G}$ vanishes, which
enforces the existence and the positivity of its limit,
$\mathcal{M}^{(\infty,D)}_{G,G}>0$.

In order to express $\mathbb{E}[T_{H}(Y)]$ in a more explicit form (see
Appendix D for some concrete examples of low orders) one should exploit the
independence of the random variables forming the Wishart and GOE matrices; it
is thus wise to introduce a variant of $T_{H}$ called the injective trace in
Ma20 ; CeDaMa16 ,

$T_{H}^{\rm inj}(Y)=\sum_{\phi\in\mathcal{I}(V,[n])}\prod_{e=\\{a,b\\}\in
E}Y_{\phi(a),\phi(b)}\ .$ (204)

By decomposing the sum over $\phi\in\mathcal{F}(E,F)$ as a sum over the
partitions $\pi$ of $E$ and a sum over the injections from the blocks of $\pi$
towards $F$, one can reduce the computation of $\mathbb{E}[T_{H}(Y)]$ to the
ones of $\mathbb{E}[T_{H^{\pi}}^{\rm inj}(Y)]$, with $H^{\pi}$ the graphs
obtained from $H$ by gluing together its vertices according to the various
partitions $\pi$. One advantage of this injective formulation is the
factorization over connected components it ensures, thanks to the independence
structure of GOE and Wishart matrices. Since in addition each edge of
$Y_{i,j}$ can be decomposed as $S_{i,j}+Z_{i,j}$, which are independent one of
the other, one further reduces the computation to the one of
$\mathbb{E}[T_{H}^{\rm inj}(S)]$ and $\mathbb{E}[T_{H}^{\rm inj}(Z)]$, where
$H$ is an arbitrary connected graph. For the GOE part one can write a very
explicit formula by exploiting the independence of the matrix elements:

$\mathbb{E}[T_{H}^{\rm
inj}(Z)]=\left(\frac{\Delta}{n}\right)^{|E|/2}n(n-1)\dots(n-|V|+1)\prod_{e\in
E_{\neq}}\mathbb{E}[B_{1,2}^{d_{e}}]\prod_{e\in
E_{=}}\mathbb{E}[B_{1,1}^{d_{e}}]\ ,$ (205)

where $|E|$ is the number of edges of $H$, $|V|$ its number of vertices, and
we called $E_{\neq}$ (resp. $E_{=}$) the set (without repetitions) of edges
between distinct variables (resp. self-loops), $d_{e}$ being the multiplicity
of the edge $e$ in $H$. By the centering of the $B_{i,j}$ this contribution
vanishes if there are edges with a multiplicity of 1, and one can use the
known values of the moments of Gaussian random variables to further simplify
this expression. The algebraic dependency on $n$ claimed above is apparent on
this expression. The Wishart counterpart, $\mathbb{E}[T_{H}^{\rm inj}(S)]$,
requires an additional step: from the expression (107) of the matrix elements
$S_{i,j}$ one has to introduce a summation over an index $\mu_{e}=1,\dots,r$
for each edge of $H$. This can be represented by a bipartite graph, with two
types of vertices representing the indices $i$ and $\mu$ respectively, and
edges corresponding to the random variables $X_{i,\mu}$ joining them.
Introducing once again a sum over the partitions of the $\mu$-vertices, and
identifying them along the blocks, one obtains a bipartite graph in which all
the $i$ indices and all the $\mu$ indices are distinct. The average over the
$X_{i,\mu}$ can then be performed very easily thanks to their independence,
and one obtains a formula similar to (205), with products of moments of
$X_{1,1}$ according to the degeneracy of each edge (with again a vanishing of
the corresponding contribution if some edges are simple). The combinatorics of
the sum over the $\mu$’s yields a polynomial factor in $r=n/\alpha$, which is
thus polynomial in $n$.

The conclusion of this series of observations is that $\mathbb{E}[T_{H}(Y)]$
has indeed the claimed dependency in $n$ of the form $n^{y_{H}}$ times a
polynomial in $1/\sqrt{n}$, with coefficients that depend on $\alpha,\Delta$,
and on the moments of $X$. The exponent $y_{H}$ depends on the graph $H$ in a
complicated way, through the various graphical transformations we just
exposed, but in principle allows to fix the exponents $x_{G}$ of the bad
graphs. Note that the computation of $\mathcal{R}^{(n,D)}_{G}$ follows from
the same graphical considerations, by adding to the graph $G$ one edge
labelled with $S$ between its marked vertices.

We ensured a necessary condition for the invertibility of
$\mathcal{M}^{(\infty,D)}$, it would remain to prove that this matrix is
indeed invertible, as well as the other conditions in (130), a seemingly
challenging task in combinatorics and graph theory. Note indeed that most
random matrix computations focus on the leading order behavior of
$\mathbb{E}[{\rm Tr}(Y^{p})]$, or covariances of the form $\mathbb{E}[{\rm
Tr}(Y^{p}){\rm Tr}(Y^{q})]-\mathbb{E}[{\rm Tr}(Y^{p})]\mathbb{E}[{\rm
Tr}(Y^{q})]$. In graphical terms this corresponds to $\mathbb{E}[T_{H}(Y)]$
for relatively simple graphs $H$ with one or two disjoint cycles, whose
injective trace is dominated respectively by cacti graphs CeDaMa16 ;
MaFoLaKrMeZd19  (a cactus graph being defined as a graph in which each edge
belongs exactly to one simple cycle), and annular non-crossing partitions
MaMiPeSp22 . On the contrary we require here a control of the trace of all
graphs, at their leading order; the theory of traffic distributions Ma20 ;
CeDaMa16  deals precisely with those objects, but with a different
normalization, that forces the limit of the injective trace of non-cacti
graphs to vanish CeDaMa16 . Upper bounds have been obtained for the traces of
arbitrary graphs in MiSp12 , but again we would need tighter results and
matching lower bounds to achieve our goal.

Note finally that (201) is by far not the only possible choice for the basic
functions $b_{G}$; its advantage was to make obvious the distinction between
good and bad graphs, yet replacing the summations there by their injective
versions, as well as replacing the monomials on the edges by the Hermite
polynomials, orthogonal with respect to the Gaussian measure (as done in
MoWe22  in the proof of their proposition 3.3) might instead simplify the
derivation of the properties (130).

## Acknowledgments

I warmly thank Guillaume Barraquand, Djalil Chafaï, Laura Foini, Florent
Krzakala, Antoine Maillard, Cris Moore, Soledad Villar and Lenka Zdeborová for
useful discussions, and in particular A. M. for providing some numerical data
from MaKrMeZd22  and for a careful reading of the manuscript, and C. M. for
sharing an unpublished draft and for pointing the reference Pr76 . I also
benefited from stimulating exchanges with Jean Barbier, Francesco Camilli and
Marc Mézard.

## References

  * (1) O. Ledoit and S. Péché.  Eigenvectors of some large sample covariance matrix ensembles.  Probab. Theory Relat. Fields, 151, 233–264 (2011).
  * (2) J. Bun, R. Allez, J.-P. Bouchaud, and M. Potters.  Rotational Invariant Estimator for General Noisy Matrices.  IEEE Transactions on Information Theory, 62(12), 7475–7490 (2016).
  * (3) J. Bun, J.-P. Bouchaud, and M. Potters.  Cleaning large correlation matrices: Tools from Random Matrix Theory.  Physics Reports, 666, 1–109 (2017).
  * (4) B. Olshausen and D. Field.  Emergence of simple-cell receptive field properties by learning a sparse code for natural images.  Nature, 381, 607–609 (1996).
  * (5) M. Zibulevsky and B. A. Pearlmutter.  Blind source separation by sparse decomposition in a signal dictionary.  Neural computation, 13(4), 863–882 (2001).
  * (6) K. Kreutz-Delgado, J. F. Murray, B. D. Rao, K. Engan, T.-W. Lee, and T. J. Sejnowski.  Dictionary learning algorithms for sparse representation.  Neural computation, 15(2), 349–396 (2003).
  * (7) D. A. Spielman, H. Wang, and J. Wright.  Exact Recovery of Sparsely-Used Dictionaries.  In Proceedings of the 25th Annual Conference on Learning Theory, volume 23, pages 37.1–37.18, 2012.
  * (8) L. Zdeborova and F. Krzakala.  Statistical physics of inference: thresholds and algorithms.  Advances in Physics, 65(5), 453–552 (2016).
  * (9) S. Rangan and A. K. Fletcher.  Iterative estimation of constrained rank-one matrices in noise.  In 2012 IEEE International Symposium on Information Theory Proceedings, pages 1246–1250, 2012.
  * (10) T. Lesieur, F. Krzakala, and L. Zdeborová.  Constrained low-rank matrix estimation: phase transitions, approximate message passing and applications.  Journal of Statistical Mechanics: Theory and Experiment, 2017(7), 073403 (2017).
  * (11) M. Lelarge and L. Miolane.  Fundamental limits of symmetric low-rank matrix estimation.  Probab. Theory Relat. Fields, 173, 859–929 (2019).
  * (12) J. Barbier and N. Macris.  The adaptive interpolation method: a simple scheme to prove replica formulas in Bayesian inference.  Probab. Theory Relat. Fields, 174, 1133–1185 (2019).
  * (13) A. Montanari and A. S. Wein.  Equivalence of Approximate Message Passing and Low-Degree Polynomials in Rank-One Matrix Estimation.  arXiv:2212.06996, (2022).
  * (14) Y. Kabashima, F. Krzakala, M. Mézard, A. Sakata, and L. Zdeborová.  Phase Transitions and Sample Complexity in Bayes-Optimal Matrix Factorization.  IEEE Transactions on Information Theory, 62(7), 4228–4265 (2016).
  * (15) H. C. Schmidt.  Statistical Physics of Sparse and Dense Models in Optimization and Inference.  In PhD thesis, 2018.
  * (16) A. Maillard, F. Krzakala, M. Mézard, and L. Zdeborová.  Perturbative construction of mean-field equations in extensive-rank matrix factorization and denoising.  Journal of Statistical Mechanics: Theory and Experiment, 2022(8), 083301 (2022).
  * (17) J. Barbier and N. Macris.  Statistical limits of dictionary learning: Random matrix theory and the spectral replica method.  Phys. Rev. E, 106, 024136 (2022).
  * (18) E. Troiani, V. Erba, F. Krzakala, A. Maillard, and L. Zdeborová.  Optimal denoising of rotationally invariant rectangular matrices.  Proceedings of Mathematical and Scientific Machine Learning (MSML), PMLR, 190, 97–112 (2022).
  * (19) F. Pourkamali, J. Barbier, and N. Macris.  Matrix Inference in Growing Rank Regimes.  arXiv:2306.01412, (2023).
  * (20) F. Pourkamali and N. Macris.  Bayesian Extensive-Rank Matrix Factorization with Rotational Invariant Priors.  arXiv:2306.04592, (2023).
  * (21) F. Pourkamali and N. Macris.  Rectangular Rotational Invariant Estimator for General Additive Noise Matrices.  arXiv:2304.12264, (2023).
  * (22) I. D. Landau, G. C. Mel, and S. Ganguli.  Singular vectors of sums of rectangular random matrices and optimal estimation of high-rank signals: The extensive spike model.  Phys. Rev. E, 108, 054129 (2023).
  * (23) F. Camilli and M. Mézard.  Matrix factorization with neural networks.  Phys. Rev. E, 107, 064308 (2023).
  * (24) F. Camilli and M. Mézard.  The Decimation Scheme for Symmetric Matrix Factorization.  arXiv:2307.16564, (2023).
  * (25) Harish-Chandra.  Differential Operators on a Semisimple Lie Algebra.  American Journal of Mathematics, 79(1), 87–120 (1957).
  * (26) C. Itzykson and J. Zuber.  The planar approximation. II.  Journal of Mathematical Physics, 21(3), 411–421 (1980).
  * (27) A. Matytsin.  On the large-$N$ limit of the Itzykson-Zuber integral.  Nuclear Physics B, 411(2), 805–820 (1994).
  * (28) A. Guionnet and O. Zeitouni.  Large Deviations Asymptotics for Spherical Integrals.  Journal of Functional Analysis, 188(2), 461–515 (2002).
  * (29) J. Bun, J. P. Bouchaud, S. N. Majumdar, and M. Potters.  Instanton Approach to Large $N$ Harish-Chandra-Itzykson-Zuber Integrals.  Phys. Rev. Lett., 113, 070201 (2014).
  * (30) S. B. Hopkins and D. Steurer.  Efficient Bayesian Estimation from Few Samples: Community Detection and Related Problems.  In 2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS), pages 379–390, 2017.
  * (31) D. Kunisky, A. S. Wein, and A. S. Bandeira.  Notes on Computational Hardness of Hypothesis Testing: Predictions Using the Low-Degree Likelihood Ratio.  In Mathematical Analysis, its Applications and Computation, pages 1–50. Springer International Publishing, 2022.
  * (32) G. Bresler and B. Huang.  The algorithmic phase transition of random $k$-sat for low degree polynomials.  In 2021 IEEE 62nd Annual Symposium on Foundations of Computer Science (FOCS), pages 298–309, 2022.
  * (33) T. Schramm and A. S. Wein.  Computational barriers to estimation from low-degree polynomials.  The Annals of Statistics, 50(3), 1833 – 1858 (2022).
  * (34) D. Jonsson.  Some limit theorems for the eigenvalues of a sample covariance matrix.  Journal of Multivariate Analysis, 12(1), 1–38 (1982).
  * (35) A. M. Khorunzhy, B. A. Khoruzhenko, and L. A. Pastur.  Asymptotic properties of large random matrices with independent entries.  J. Math. Phys., 37(10), 5033–5060 (1996).
  * (36) G. Anderson and O. Zeitouni.  A CLT for a band matrix model.  Probab. Theory Relat. Fields, 134, 283–338 (2006).
  * (37) G. W. Anderson, A. Guionnet, and O. Zeitouni.  An introduction to random matrices.  Cambridge University Press, 2010.
  * (38) Z. Bai and J. W. Silverstein.  Spectral Analysis of Large Dimensional Random Matrices.  Springer, 2010.
  * (39) J. A. Mingo and R. Speicher.  Second order freeness and fluctuations of random matrices: I. Gaussian and Wishart matrices and cyclic Fock spaces.  Journal of Functional Analysis, 235(1), 226–270 (2006).
  * (40) C. E. I. Redelmeier.  Real Second-Order Freeness and the Asymptotic Real Second-Order Freeness of Several Real Matrix Models.  International Mathematics Research Notices, 2014(12), 3353–3395 (2012).
  * (41) J. A. Mingo and M. Popa.  Real second order freeness and Haar orthogonal matrices.  J. Math. Phys., 54(5), 051701 (2013).
  * (42) C. Male, J. A. Mingo, S. Péché, and R. Speicher.  Joint global fluctuations of complex Wigner and deterministic matrices.  Random Matrices: Theory and Applications, 11(02), 2250015 (2022).
  * (43) D. Voiculescu.  Limit laws for Random matrices and free products.  Invent. math., 104, 201–220 (1991).
  * (44) M. Potters and J.-P. Bouchaud.  A first course in random matrix theory.  Cambridge University Press, 2021.
  * (45) J. A. Mingo and R. Speicher.  Free Probability and Random Matrices.  Springer, 2017.
  * (46) A. Nica and R. Speicher.  Lectures on the Combinatorics of Free Probability.  Cambridge University Press, 2006.
  * (47) P. Tarrago.  Spectral deconvolution of unitarily invariant matrix models.  arXiv:2006.09356, (2020).
  * (48) Z. D. Bai, B. Q. Miao, and G. M. Pan.  On asymptotics of eigenvectors of large sample covariance matrix.  The Annals of Probability, 35(4), 1532 – 1572 (2007).
  * (49) L. Erdős, B. Schlein, and H.-T. Yau.  Local Semicircle Law and Complete Delocalization for Wigner Random Matrices.  Commun. Math. Phys., 287, 641–655 (2009).
  * (50) F. Benaych-Georges.  Eigenvectors of Wigner matrices: universality of global fluctuations.  arXiv:1104.1219, (2011).
  * (51) A. Bouferroum.  Eigenvectors of Sample Covariance Matrices: Universality of global fluctuations.  arXiv:1306.4277, (2013).
  * (52) S. O’Rourke, V. Vu, and K. Wang.  Eigenvectors of random matrices: A survey.  Journal of Combinatorial Theory, Series A, 144, 361–442 (2016).
  * (53) R. Dudeja, S. Sen, and Y. M. Lu.  Spectral Universality of Regularized Linear Regression with Nearly Deterministic Sensing Matrices.  arXiv:2208.02753, (2022).
  * (54) T. Wang, X. Zhong, and Z. Fan.  Universality of Approximate Message Passing algorithms and tensor networks.  arXiv:2206.13037, (2022).
  * (55) G. W. Anderson and B. Farrell.  Asymptotically liberating sequences of random unitary matrices.  Advances in Mathematics, 255, 381–413 (2014).
  * (56) C. Procesi.  The invariant theory of $n\times n$ matrices.  Advances in Mathematics, 19(3), 306–381 (1976).
  * (57) R. Goodman and N. R. Wallach.  Symmetry, Representations, and Invariants.  Springer, 2009.
  * (58) B. Simon.  The Classical Moment Problem as a Self-Adjoint Finite Difference Operator.  Advances in Mathematics, 137(1), 82–203 (1998).
  * (59) P. Biane.  On the free convolution with a semi-circular distribution.  Indiana Univ. Math. J., 46, 705–718 (1997).
  * (60) F. Lehner.  Cumulants, lattice paths, and orthogonal polynomials.  Discrete mathematics, 270(1-3), 177–191 (2003).
  * (61) R. Stanley.  Enumerative Combinatorics: Volume 2.  Cambridge Studies in Advanced Mathematics. Cambridge University Press, 1997.
  * (62) R. Gustafson and S. Milne.  Schur functions, Good’s identity, and hypergeometric series well poised in $SU(n)$.  Advances in Mathematics, 48(2), 177–188 (1983).
  * (63) D. Voiculescu.  The analogues of entropy and of Fisher’s information measure in free probability theory, I.  Commun. Math. Phys., 155, 71–92 (1993).
  * (64) D. Guo, S. Shamai, and S. Verdu.  Mutual information and minimum mean-square error in Gaussian channels.  IEEE Transactions on Information Theory, 51(4), 1261–1282 (2005).
  * (65) D. Guo, S. Shamai, and S. Verdu.  Additive non-Gaussian noise channels: mutual information and conditional mean estimation.  In Proceedings. International Symposium on Information Theory, 2005\. ISIT 2005., pages 719–723, 2005.
  * (66) V. A. Marcenko and L. A. Pastur.  Distribution of eigenvalues for some sets of random matrices.  Math. USSR Sb., 1(4), 457 (1967).
  * (67) O. Puny, D. Lim, B. T. Kiani, H. Maron, and Y. Lipman.  Equivariant polynomials for graph neural networks.  In Proceedings of the 40th International Conference on Machine Learning, ICML’23. JMLR.org, 2023.
  * (68) C. Male.  Traffic Distributions and Independence: Permutation Invariant Random Matrices and the Three Notions of Independence.  Memoirs of the AMS, 267(1300) (2020).
  * (69) F. Benaych-Georges, J.-P. Bouchaud, and M. Potters.  Optimal cleaning for singular values of cross-covariance matrices.  The Annals of Applied Probability, 33(2), 1295–1326 (2023).
  * (70) P. Biane.  Free probability for probabilists.  Quantum Probability Communications, pages 55–71 (2003).
  * (71) B. Collins and P. Śniady.  Integration with Respect to the Haar Measure on Unitary, Orthogonal and Symplectic Group.  Commun. Math. Phys., 264, 773 (2006).
  * (72) B. Blum-Smith and S. Villar.  Machine learning and invariant theory.  arXiv:2209.14991, (2022).
  * (73) B. Collins, S. Matsumoto, and J. Novak.  The Weingarten Calculus.  Notices of the American Mathematical Society, 69(05), 734 (2022).
  * (74) G. Cébron, A. Dahlqvist, and C. Male.  Traffic Distributions and Independence II: Universal Constructions for Traffic Spaces.  arXiv:1601.00168, (2016).
  * (75) A. Maillard, L. Foini, A. L. Castellanos, F. Krzakala, M. Mézard, and L. Zdeborová.  High-temperature expansions and message passing algorithms.  J. Stat. Mech., 2019, 113301 (2019).
  * (76) J. A. Mingo and R. Speicher.  Sharp bounds for sums associated to graphs of matrices.  Journal of Functional Analysis, 262(5), 2272–2288 (2012).