Bayesian Extensive-Rank Matrix Factorization with Rotational Invariant Priors
Authors: Farzad Pourkamali
Nicolas Macris
Venue: NeurIPS 2023 spotlight (2023)
Keywords: Matrix factorization
Bayesian inference
rotation invariant estimators
random matrix theory
spherical integrals
replica method

Abstract: We consider a statistical model for matrix factorization in a regime where the rank of the two hidden matrix factors grows linearly with their dimension and their product is corrupted by additive noise. Despite various approaches, statistical and algorithmic limits of such problems have remained elusive. We study a Bayesian setting with the assumptions that (a) one of the matrix factors is symmetric, (b) both factors as well as the additive noise have rotational invariant priors, (c) the priors are known to the statistician. We derive analytical formulas for Rotation Invariant Estimators to reconstruct the two matrix factors, and conjecture that these are optimal in the large-dimension limit, in the sense that they minimize the average mean-square-error. We provide numerical checks which confirm the optimality conjecture when confronted to Oracle Estimators which are optimal by definition, but involve the ground-truth. Our derivation relies on a combination of tools, namely random matrix theory transforms, spherical integral formulas, and the replica method from statistical mechanics.

Bayesian Extensive-Rank Matrix Factorization with Rotational Invariant Priors
Farzad Pourkamali & Nicolas Macris
School of Computer and Communication Sciences,
Ecole Polytechnique Fédérale de Lausanne ffarzad. pourkamali, Marzaad.puuramaa... aanrraamliinncolass nicolas macris)@epfll. ch
Abstract
We consider statistical model for matrix factorization in a regime where the rank of the two hidden matrix factors grows linearly with their dimension and their product is corrupted by additive noise. Despite various approaches, statistical and algorithmic limits such problems have remained elusive. We study Bayesian setting with the assumptions that (a) one of the matrix factors is symmetric, (b) both factors as well as the additive noise have rotational invariant priors, (c) the priors are known to the statistician. We derive analytical formulas for Rotation Invariant Estimators to reconstruct the two matrix factors, and conjecture that these optimal in the large-dimension limit, in the sense that they minimize the average rrean-sqqare-ervee We provide numerical checks which confirm the optimality conjecture when confronted to Oracle Estimators which optimal by definition, but involve the ground- truth. Our derivation relies on a combination tools, namely random matrix theory transforms, spherical integral formulas, and the replica method from statistical mechanics.
Introduction
Matrix factorization (MF) is the problem of reconstructing matrices X and Y from the noisy observations of their product. Applications in signal processing and machine learning abound, such as
for example dimensionality reduction (12). sparse coding (355) representation learning [6) robust
principal components analysis 7 8 blind source separation or matrix completion [10
In this work approach the problem from Bayesian perspective and assume that observation data matrix S = VKXY= + W is given to a statistician who knows the prior distributions of X and Y well the prior of the additive noise matrix W and the signal-to-noiss ratio 0. The task of the statistician to construct estimators E Ex( (-). Ey(-) (-) for the matrix factors X, Y. that idcally, minimize the average mean-square-eroo (MSE) E X Ex Ex(S)) (S) and Y Ey) (I-IIr the Frobenius norm and the expectation w.r.t X, Y, W). We consider priors which are rotation invariant for all three matrices X, Y. W and for X furthermore impose that it is square and symmetric. These matrix ensembles are defined precisely in section 2. but the reader can keep in mind the examples of Wigner Wishart matrices for X, and general Gaussian Y and W with i.i.d elements. We look at the asymptotic regime where all matrix dimensions and ranks tend to infinity the same speed. We remark that the usual "rotation ambiguity" occuring MF not present because we impose that at least one of the two matrix factors is symmetric. We also remark that MF is different (and more difficult) than matrix denoising which would consist in constructing an estimator (Xy(S) for the signal
whole by minimizing E X (S)
The rotation invariance of the model implies that the estimators minimizing the MSE belong to the class of rotation invariant estimators (RIE). RIEs are matrix estimators which have the same singular vectors (or eigenvectors) as the observation or data matrix. These estimators have been
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

for covariance estimation, [17] for cross-covariance estimation, and (18). [9 20 for extensions to rectangular matrices). For the present MF model, we derive optimal estimators (minimizing the MSE) that belong to the RIE class and can be computed explicitly the large dimensional limit from the observation matrix and the knowledge of the priors. We propose:
proposed for matrix denoising problems (see references [12 [6] [
1. an explicit RIE to estimate X, which requires the knowledge of the priors of both X, Y and of the noise W. Moreover, under the assumption that X is positive-semi-defiice sub-optimal RIE can be derived which does not require any prior on X.
2. an explicit RIE to estimate Y, which requires the knowledge of the priors of the noise W and X only (the prior of Y is not required).
3. combined with the singular value decomposition (SVD) of the observation matrix, our explicit RIEs provide a spectral algorithm to reconstruct both factors X and Y.
The derivation of the proposed estimators relies on the replica method from statistical mechanics combined with techniques from random matrix theory and finite-rank spherical integrals Although the replica method is not rigorous and involves concentration assumptions, the derivation is entirely analytical and suggests that the estimators are optimal in the limit of large dimensions. This is corroborated by numerical calculations comparing our explicit RIEs with Oracle Estimators which are optimal by definition and involve the ground-trum -tunh matrices.
1.1 Related literature and discussion
When the matrices X and Y are assumed to have low-rank compared to their dimension, the mathematical theory of MF has enjoyed much progress under various settings (Bayesian, spectral, algorithmic) and fundamental information theoretical and algorithmic limits have been rigorously derived. The behavior of eigenvalues/singular values and eigenvector/singular vectors of finite-rank perturbations of a Gaussian matrix is studied in (23-25] which leads to spectral estimators when the noise matrix is Gaussian distributed. For the case of factorized prior, and Gaussian noise, closed form expressions have been established for the asymptotic Gayes-lptimal estimation error [26 30). and iterative algorithms based on approximate message passing has been proposed [3 [31,2]. 32). The low-rank matrix denoising problem has been addressed in various other settings, such as structured
noise matrix [4],3]] mismatched estimation problem [35 38) and estimation in the regime with
diverging aspect-ratio of matrices (39].
In extensive-rank regimes, when the rank grows like the matrix dimensions, despite various attempts there is no solid theory of MF. One approach is based on Approximate Message Passing (AMP) methods developed in [4). Despite acceptable performance in practical settings (43), as pointed
out in 44 the AMP al algorithms developed in these works are (theoretically) sub-optimal. Other approaches rooted in statistical physics have been considered in 47) but have not led to explicit reconstructions of matrix factors or algorithms. A practical probabilistic approach to MF problem is based on variational Bayesian approximations (48)50). in which one tries approximate the posterior distribution with proper distribution. In [5I] it is shown that under Gaussian priors, the solution to the MF problem a reweighted SVD of the observation matrix. We point out here that these estimators can be seen as a RIE and therefore there seems to be a rather close relation between the RIE studied here and the variational Bayesian approach. This also suggests that adapting RIEs real data an interesting direction for future research. Finally, let us also mention optimization approaches where one constructs estimators by following a gradient flow (or gradient descent) trajectory of a training loss of the type ||S XY il + reg. term (see [52) [33 for analysis in rotation invariant models). Benchmarking these various other algorithmic approaches against our explicit RIEs (conjectured to be optimal) is outside the scope of this work and is left for future work.
Constraints such as sparsity or non-negativity of the matrix entries which have important applications (54) are not covered by our theory. Despite this drawback, we believe that the proposed estimators are important both for theoretical and practical purposes. Even non-rotation invariant problems our explicit RIEs may serve as sub-optimal estimators, and as we show in an example they can be used as a "warmed-up" spectral initialization for more efficient algorithms (see for example (1555] for related ideas in other contexts). The methodology developed here may open up the way to further analysis in inference and learning problems perhaps also in the context of neural networks where extensive rank weight matrices must be estimated.

---

1.2 Organization and notations
2 we introduce the precise MF model, general class of RIEs, and the Oracle estimators. In section we present the explicit RIEs (and algorithm) to estimate X and Y. We provide the numerical examples and calculations in section 4 In section we sketch the derivation of RIE for X, while the one for is similar and deferred to the appendices.
In section
The following notations are used throughout. For vector E R. we denote by I E R NxM constructed as r I' ONN(M-N) with F' E R RN*N a diagonal matrix with diagonal 7. The same notations will also be used for the vector о and the corresponding matrix Σ and For sequence of non-symmetrie matrices of growing size, we denote the limiting empirical singular value distribution (ESD) by μA and the limiting empirical eigenvalue distribution of A by PA- For a sequence of symmetric matrices B of growing size, we denote the limiting empirical eigenvalue distribution by pB. and the limiting cigenvalue distribution of B by PBa.
a matrix
2 Matrix factorization model and rotation invariant estimators
2.1 Matrix factorization model
Let X -T NNNN
symmetric matrix distributed according to a rotationally invariant prior (x(X). i.e., for any orthogonal matrix O E RNNN we have Px Px(OXOT) Px (X). Let also Y E N be distributed according to a bi-rotationally invariant prior Py i.e. for any orthogonal matrices U E RNN, V E M we have Py (UUUVT) = PYY).. We observe the data matrix S E ERNNM
S VEXY W
where W NRM is also bi-rotationally invariant distributed, and K. R4 is proportional to the signal-to-noise-ratio (SNR). The goal is to recover both factors X and Y from the data matrix S. For definiteness, we consider the regime M 2 N with aspect ratio N/M + a E (0, 1] as N + op. The case of a > 1 can be analyzed in the same manner and is presented in appendix F Furthermore. we assume that the entries of X, Y and W are of the order O(//N). This scaling is such that the eigenvalues of X and singular values of Y. W and S are of the order O(1) as N cc.
The empirical eigenvalue distribution of converge weakly to measure px, and the ,yr ESD Y Y. W converge weakly to measures py Hw with bounded support on the real line. Moreover,
Assumption
these measures are known to the statistician. He can deduce (in principle) these measures from the priors on X.Y.W. W.
matrixfactorriin factorization the hidden matrices have dimensions X E N XH Y E RHxM and in the Bayesian framework with bi-rotational invariant priors for both factors, the optimal estimators are trivially the zero matrix. Indeed, from bi-rotational invariance fothfactors,
Remark 1. In a general formulation of matrix E
(Y), which implies that the Bayesian estimate is zero. Here, by imposing that E N is symmetric and Px (OXOT) = Px (X), we can break this symmetry and find non-trivial estimators. This is due the fact that the map X + -X cannot be realized as a (real) orthogonal transformation, so (x(-X) x)= Px (X) does not hold in general (various examples are given in section 4 and appendices). Of course, if the prior is even, e.g. Wigner ensemble, again the Bayesian posterior estimate is trivially zero for both factors. As we will see our RIES consistent with these observations.
have Px
(X). (-Y)
2.2 Rotation invariant estimators
To recover matrices X, Y from S, we consider two denoising problems. One is recovering X by treating both Y, W as "noise" matrices, and the other estimating Y by treating X, W as "noise" As will become clear the procedure is not iterative, and the two denoising problems are solved independently and simultaneously. In the following, for each of these two problems, we introduce two rotation invariant classes of estimators and discuss their optimum Oracle estimators. We then provide an explicit construction and algorithm for RIES which we conjecture have the optimum performance of Oracle estimators in the large N limit.
3

---

2.2.1 RIE class for X
Consider the SVD of S UsrV3, where Us E RNNN Vs RMM are orthogonal, and E RNX is diagonal matrix with singular values of S on its diagonal, (7) A rotational invaraint estimator for X is denoted Ex(S), and is constructed as: ( Us diag(f.1....... ..N) U3 where E.i..... ENN are the eigenvalues of the estimator.
(2)
First, we derive an Oracle estimator by minimizing the squared error N X -xx((S)]lf for a given instance, over the RIE class equivalently over the choice of the eigenvalues (Ema) ISIEN' Let the
eigen-decomposition of X be X Ei= A a with æ; RN eigenvectors of X. The error be expanded as:
N
N
N
N
X (S) =
Coi
Σ
X (u]x;))
N
N
N
where u,'s are columns of Ug. Minimizing over Ex, ss we find the optimum among the RIE class:
n x (S) Σ u,u e. - (uTa,) u Xu;
N
(3)
j=1
Expression B defines the Oracle estimator which requires the knowledge signal matrix X. Surprisingly, the large N limit, the optimal eigenvalues (Etx) 1SiSN can be computed from the observation matrix and knowledge of the measures px ; po,, In the next section, we show that this leads to explicitly computable (or algorithmic) RIE, which we conjecture to be optimal as N + 0o, the sense that its performance matches the one the Oracle estimator.
Now we remark that the Oracle estimator is not only optimal within the rotation invariant class but also Bayesian optimal. From the Bayesian estimation point of view, one wishes to minimize the average mean squared error (MSE) MSE x = E X (S) where the expectation is over
X. Y. W, and (S) is an estimator of X. The MSE minimized for X*(S) E[X [S] which is the posterior mean. Therefore, the posterior mean estimator has the minimum MSE (MMSE) among all possible estimators, particular MSF S MSFE; for any N. In appendix A.T we show that, for rotational invariant priors, the posterior mean estimator is inside the RIE class. Thus, since 3 (S) optimum among the RIE class MSFg, MSF Therefore, we conclude that the Oracle estimator (3 is Bayesian optimal in the sense that MSEE; MSE x. MMSE.
2.2.2 RIE class for Y
Estimators for Y from the rotation invariant class are denoted (S), and are constructed as:
Ey (S) Us| diag(Sy)...... .... Eyn)) ONx(M-N) VS
where Eyp...., EuN are the singular values of the estimator.
Let the SVD Y be Y = EN os ( (r) with ( e R y E RM the left and right singular vectors of Y. To derive an Oracle estimator, proceed as above. Expanding the error, have:
N
N
1 Y Ey N
2
(oTU,")
S
Su
N
N
T=1
j=1
where 's are columns of Vs. Minimizing Evy 's. find the optimum among the RIE class:
ES(S) = Σ C u.o E = Eo, ( uf Y
(5)
Expression (5) defines the Oracle estimator which requires the knowledge of signal matrix Y. Like for the case of X, in the large N limit we can derive the optimal singular values (5t) ISISN in terms

---

of the singular values of observation matrix and knowledge of the measures Px Pw. This leads to an explicitly computable (or algorithmic) RIE, which is conjectured be optimal as N + oc, in the sense that it has the same performance as the Oracle estimator. Note that unlike the estimator for X, we not need the knowledge of jay.
In appendix A.2 we show that for bi-rotationally invariant priors the posterior mean estimator (S) = S| belongs to the RIE class, which (by similar arguments to the case of X) implies that the Oracle estimator Bayesian optimal.
3 Algorithmic RIEs for the matrix factors
In this section, we present our explicit RIEs for X, Y and the corresponding algorithm. We conjecture that their performance matches the of Oracles estimators in the large N limit and they therefore Bayesian optimal in this limit. Let us first give a brief reminder on useful transforms in random matrix theory.
3.1 Preliminaries transforms in random matrix theory
For probability density function p(x) on IR, the Stieltjes (or Cauchy) transform defined as
9p(z) = [ p(xx)
pla) da for E C)supp(p)
By Plemelj formulae we have for x E R2
lim 9p(x - ic) = + Tip(#)
(6)
with H[p((z) = p.v. J p(t) dt the Hilbert transform of p (here p.v. stands for "principal value").
Denoting the inverse of (2) by 9p (z), the R-transform of p is defined as (56]:
R₂(2) 97
For probability density function μ with support contained in - KK with K > O, we define generating function of (even) moments M₁ : [0, K - + R+ as Mn(z) = (e)((() di 1. For
a E (0, 1]. define T( (z) = (az + += + 1), and H, (z) = 2T(1) (M,())) The rectangular
R-transform with aspect ratio is defined 57
(2)
H/a) (a) -1
3.2 Explicit RIE for X
The RIE for X is constructed as (S) = Em Siz u.u with eigenvalues
(fr IsiSN
Im
lim
{ 9ox
17 1ox
(V
1-ك + Gpx
(-
2 - C1
(7)
(7T
2- 10+
KC3
KCa
where %i the i-th singular value of S, as is the symmetrized limiting ESD of S, and
Si (2)
1
c(a) uw
(=)
+1-9) =
(8)
and (3 satisfies
[aG (2) + 1 1--1(( (1) ((== 11)
(z- -2) (x=
(9)
G
Remark 2. If px is symmetric measure, px px(-x)) then (-x(-2) -Gpx (2). This implies that the optimal eigenvalues (Ea) 7 are all zero, and (x = O, see figure
S1, Ca are the only parameters which appear in the final estimator. However, in derivation of the RIE, we have defined other parameters which do not appear in the final estimator and omit them here.
5

---

3.2.1 An estimator for X²
It is interesting to note that we can construct a RIE for X as Exx(S) = Ei= Lu with
eigenvalues (St))
Im lim 21 33
9ps 9ps (2)
(10)
K. kTA(()))
2-7-il
S3
with S3 in (8) Note that, S1, S3 can be evaluated using the observation matrix and the knowledge of jy, JAw, and therefore this time the statistician does not need to know the prior of X. Furthermore, assuming that X is positive semi-definite (PSD), we can construct sub-optimal RIE
for X by using V C22 for the eigenvalues of the estimator.
3.2.2 Case of Gaussian Y, W
Yf W have i..dd Gaussian entries with variance !/n, then Cis (2) C/ao (2) = 5/a. Consequently, S1, S3 can easily be computed to be Çi Ç3 Gas (2) + (a=a)/(az), thus the estimator 7 can be evaluated from observation matrix. particular, estimator (10) simplifies to:
(11)
a (Tps(t))) (TH(as](1s) + ))
3.3 Explicit RIE for Y
Our explicit RIE for Y is constructed as (S) = Ei= Sy, ucd with singular values
Im
lim
q4
(12)
Vk (T))(1))
2-
10+
where 7, the i-th singular value S. and 6a is the solution to the following system equations
81 eia) (4)
+ 2 91
Rox (6r + /9193)
Rox (01 v0143)
91
P1 RRpx (91 + /999s)
+ R,x (6r /919s)
(13)
91 Ops (2),
92 a
us
(2) + (1-a)
(1- Gns
2-8 0
3-91 (p)(2) 34
Q3
Similarly to the estimator derived for X, if px a symmetric measure then the optimal singular values for the estimator of Y are all zero, remark 5
X shifted Wigner matrix, i.e. X F + el with =T FT NNN having i.i.d. Gaussian entries with variance I/N and * 0 a real number, then Roy (2) + Moreover, if W is
Gaussian matrix with variance I/N, then the set of equations (13) simplifies great extent, and we can compute d1 analytically in terms of Gay (z), see appendix D.4
3.4 Algorithmic nature of the RIES
The explicit RIEs (7) and 122 proposed in this section, provide spectral algorithms estimate the matrix factors from the data matrix (and the priors). An essential ingredient that must be extracted from the data matrix is (2). This quantity can be approximated from the observation matrix using Cauchy kernel method introduced in (58) (see section 19.5.2), from which (.).) can be approximated using 6 Therefore, given an observation matrix S, the spectral algorithm proceeds as follows:
1. Compute the SVD of S.
2. Approximate Gps (z) from the singular values of S.
3. Construct the RIEs for X, Y as proposed in paragraphs 3.2 3.3
Like the for X. omit some of the parameters which do appear the final estimator.

---

OOadk: cximater (3(
eainaie Ej (S)
(5(s)
RIE. (7(
(Q-(8)
MNN
0.2
0.1
Figure I: MSE of estimating X. MSE normalized by the norm the signal, | X a Wishart matrix with aspect ratio 1/4, X = HH with H E R NNxN having i.i.d. Gaussian entries of variance M matrices with i.i.d.
Figure 2: MSE of estimating Y. MSE is normalized by the norm the signal, || li.. has uniform spec- tral density, 24 ([1.33 3)) X is a shifted Wigner matrix with e = 3. and W is a N M matrix with i..d. Gaussian entries of variance 1/N. RIE is applied to N 2000, MM 4000, and the results are averaged over 10 runs (crror bars are invisible).
1/N. Both Y and W are
Gaussian entries of variance I/N. RIE is applied to N 2000, M = 4000, and the results are averaged over 10 runs (error bars are invisible).
4 Numerical results
4.1 Performance of RIE for X
We consider the case where Y. W both have i.i.d. Gaussian entries of variance 1/N, and X is Wishart matrix, X HHT with H IR NxNN having i.i.d. Gaussian entries of variance 1/v. For various SNRs, we examine the performance of two proposed estimators, the RIE
77) and the square-root of the estimator (TO) (since X is PSD), which is sub-optimal. In figure I the MSEs of these algorithmic estimators are compared with the one of Oracle estimator (3) We see that the
average performance of the algorithmic RIE E (S) is very close to the (optimal) Oracle estimator (S) (relative errors are small and provided in the appendices) and we believe that the slight mismatch is due to the numerical approximations and finite-size effects. Note that, although the
estimator B x² (S) is sub-optimal, it does not use any prior knowledge of X. For more examples, details of the numerical experiments and the relative error of the estimators, we refer to appendix C.3
4.2 Performance of RIE for Y
We consider the case where W has i.i.d. Gaussian entries of variance I/N, and X is a shifted Wigner matrix with e 3. Matrix Y is constructed as Y Uy VV with Uy E NxN Vy E R M are Haar distributed, and the singular values are generated independently from the uniform distribution on [1, MSEs of the RIE 112 and the Oracle estimator (5) are illustrated figure 2
We see that
the performance of the algorithmic RIE E (S) ()) is very close to the optimal estimator Ey (S).
Non-rotational invariant prior In another example, which we omit here, with the same settings for X, W. we consider the case where Y is sparse matrix with entries distributed according to prior. The RIE is not optimal in this setting (since the prior is not bi-rotational invariant), however applying simple thresholding function on the matrix constructed by RIE yields an estimate with lower MSE. This observation suggests that for the case of general priors, the RIEs can provide a spectral initialization for more efficient estimators. For more details and examples, see appendix D.4
4.3 Comparing RIES of matrix factorization and denoising
The proposed RIEs, namely and 122 simplify greatly when the matrices W, Y are Gaussian, and X a shifted Wigner matrix. We perform experiments with these priors, where for given
observation matrix S, we look at the RIEs of X, Y for the MF problem, and simultaneously at the RIE of the product X as whole for the denoising problem with formulas introduced in [19 (which can also be obtained by taking X to be the identity matrix, see appendix D.3.I Figure
7

---

Oradle citinuoe Eq. RUBL E3:8)
03-
eximatee 8:18
V eitinnss (
AXIE, (8)
A (8)
0.
MSE
MIS
MSE
03-
0.1
0.2
(a) Estimating X
(b) Estimating Y
(c) Estimating XY
Figure 3: MSE of factorization problem. MSE is normalized by the norm the signal. a shifted Wigner matrix with c 1, and both and W are x M matrices with i.i.d. Gaussian entries of variance 1/x. RIE applicd to N = 2000, M = 4000. In each run, the ohservation matrix S is generated according to [ and the factors X, Y are estimated simultancously from S. Results are averaged over 10 runs (error bars invisible).
illustrates these experiments. In particular the MSE of the denoising RIE matches well the one of the associated Oracle estimator, and as expected is lower than the MSE of the product of MF-RIEs.
5 Derivation of the explicit RIES
In this section, we sketch the derivation of our explicit RIE for X. The RIE for Y is derived similarly. but requires more involved analysis and is presented in appendix D For simplicity, we take the SNR parameter in (T) to be I, so the model is XY + W. The optimal eigenvalues are constructed as 52 = Σj= j= A (u]ar)) We assume that in the large N limit, E can be approximated by its expectation and we introduce
N Σ (u(aa))
(14)
j=1
where the expectation is over the (left) singular vectors of the observation matrix S. Therefore, to compute these eigenvalues, we need to find the mean squared overlap E
(u. [{(uuq x² between
eigenvectors of X and singular vectors of S. In what follows, we will see that rescaling of) this quantity can be expressed in terms of i-th singular value of S and j-th eigenvector of X (and the limiting measures, indeed). Thus, we will use the notation Ox (7r, Ay) == NE (u]x,)² 2 following. In the next section, we discuss how the overlap can be computed from the resolvent of the "Hermitized" version of S.
in the
5.1 Relation between overlap and resolvent
Construct the matrix S E (NN+M))(N++M) from the observation matrix:
ONXN
S 0
s
ST
By Theorem 7.3.3 in [59]. S has the following eigen-decomposition:
Üs
Üs
N
0
0
ÚS
ÜS
T
s
0
0
(15)
(1)
(1)
N
v
(2)
0
with V33 (1)
v (2)
in which V E RMXN, And, V =
Ü
L
Us.
S
v2 S
V2
Eigenvalues of S are signed singular values of S. therefore the limiting eigenvalue distribution of S (ignoring zero eigenvalues) is the same as the limiting symmetrized singular value distribution of S. Define the resolvent S.
(2(2) = (zI - S)
We assume that as N -> oo and is not too close to the real axis, the matrix Gs(2) concentrates around its mean. Consequently, the value of Gs(z) becomes uncorrelated with the particular
8

---

realization of S. Specifically, as N + oc Gs(2) converges to a deterministic matrix for any fixed value E CR (independent of N). Denote the eigenvectors of S by s E RM+N 1...., M + N. For = ie with E IR and small we have:
2N
N+M Σ
+
Gc) - =
SESt +
S
x 7x)3
+
k=2N+1
where it are the eigenvalues S, which are in fact the (signed) singular values of S, 71 n,...., ,N = ,N, +N+1 --..... F2N -N
Define the vectors i; | 0ax]" for æ; eigenvectors of X. We have
N+M Σ (aTsk)² 2=
/ Im Gs(x - -
7x)² + (0.sk) 2
(16)
x2
Given the structure of S 's's in 15). second sum in (T6) is zero. We assume that in the limit of large N this quantity concentrates on Ox(7) A:) and depends only on the singular values and eigenvalue pairs (7g, Ai). We thus have:
(m[s)) = (=]
(+)s++) for ك N, and the
a! (Im Gs(x- - -ie)
N N-+0c
(x- Ox(1, 1x A.) dt
e
(17)
where the overlap function Ox (t, X;) is extended (continuously) to arbitrary values within the support of as (the symmetrized limiting singular value distribution of S) with the property that (t, Ox(ttA, for t supp(us) Sending 0, we find
a (Im Gs(r - ie)) + (xx()))x(x;Aa)
(18)
This is crucial relation as it allows us to study the overlap by means of the resolvent of S. In the next section, we establish connection between this resolvent and the signal X, which enables us to determine the optimal eigenvalues values En; in terms of the singular values of S.
5.2 Resolvent relation
To derive the resolvent relation between S and X, we fix the matrix X and consider the model
S XU + UN
with Y. W E N fixed matrices with limiting singular value distribution μy kw, and U₁, U2 e NNXN Vi, Va2 R xM independent random Haar matrices. Indeed, if we substitute the SVD of the matrices Y, W in model (T) we find the latter model. Now, the average over the singular vectors of S (with fixed X) is equivalent the average over matrices U1. U2. V1, Vz. appendix C. using the Replica trick, we derive the following relation in the limit N oo:
(Gs(=)) S3 Gxa
C
0
(19)
(2-
with <i, S2, Sa satisfying set of equations (4T). (.) is the expectation w.r.t. the singular vectors of S (or equivalently over U1, U2 V1, V2), and G xa the resolvent of X As stated earlier, assume that the resolvent G2(2) concentrates in the limit N + oc, therefore we drop the brackets in the following computation.
5.3
Overlaps and optimal eigenvalues
From
19
we find:
2 2- xi
Ox(1,Ai)
lim
x. T G
(π)
+0+
Go
(20)
1
Im
lim
Tps (7)
3 2 z C3>?
In Fig. we illustrate that the theoretical predictions (20) are in good agreement with numerical simulations for particular case of X Wigner matrix, and Y, W with i..dd Gaussian entries.

---

Simnlattir (6 200)
Stoddrtion - 500.
Figure 4: Comparison of the theoretical prediction 20 of the rescaled overlap with the numerical simulation. The rescaled overlap between 200-th and 800-th left singular vector of S and the eigen- vectors of X is illustrated. X = X E R has i.dd. Gaussian entries with variance 1/vW and is fixed. Both Y and Z are N M matrices with i.i.d. Gaussian entries of variance 1/x. The simu- lation results are average of 1000 experiments with fixed X, and N 1000, M = 2000. Some of the simulation points are dropped for clarity.
Tlaney (6 500)
Thess (K 500)
1.5
Nuqaa]
0.5
One that the overlap so the optimal eigenvalues E Zj=1 A (u. are all zero, as discussed remark 2
function of eigenvalues A..
-55
-0.5
0.5
1.5
Once we have the overlap, we can compute the optimal eigenvalues to be
N
1 N
X
Im lim
(21)
(P)(
2-
Cz>
j=1
j=1
With a bit of algebra, we find the estimator 7 in the limit - oo, see appendix
C.2
6 Conclusion
We studied the MF problem with extensive-rank hidden matrices, and proposed explicit (optimal) RIEs to recover the distinet factors. The model we considered, although limited, is the first analytically solvable model in this challenging regime, and we believe it paves the way for future investigations.. Extending the methodology developed here to the general case where both matrices can be non- symmetrie is interesting research direction (note that in this case factors can be recovered up to rotations). Morcover, adapting RIEs to incorporate additional structures/ constraints on the signals is problem with practical importance, that we leave for future investigations.
In general, the MF problem in the extensive-rank regime has not received an in in-depth exploration compared to its counterpart in the low-rank regime. While substantial progress has been made in understanding and devising algorithms for low-rank matrices, the challenges posed by extensive-rank matrices remain relatively uncharted. Specifically, the study of general non-rotation invariant priors with i.i.d entries (e.g. a prior supported on the positive real line as non-negative MF) stands out as an underdeveloped area from the theory side.
Acknowledgments
We are thankful Jean Barbier for interesting discussions on related problems. The work of has been supported by Swiss National Science Foundation grant no 200021-204119
References
[1] Ivana Tošić and Pascal Frossard. Dictionary learning. IEEE Signal Processing Magazine, 28(2):27-38, 2011.
[2] Julien Mairal, Jean Ponce, Guillermo Sapiro, Andrew Zisserman, and Francis Bach. Supervised dictionary learning. Advances in neural information processing systems, 21, 2008.
[3] Bruno Olshausen and David Field. Emergence of simple-cell receptive field properties by learning sparse code for natural images. Nature, 381 (6586)))600766666 1996.
[4] Bruno A Olshausen and David J Field. Sparse coding with an overcomplete basis A strategy employed by vI? Vision research, 37(23):3311-332 1997.
[5] Kenneth Kreutz-Delgado, Joseph F Murray, Bhaskar D Rao, Kjersti Engan, Te-Won Lee, and Terrence J Sejnowski. Dictionary learning algorithms for sparse representation. Neural computation, 15(2):349-396 2003.
10

---

[6] Yoshua Bengio, Aaron Courville, and Pascal Vincent. Representation learning: A review and new perspectives. IEEE transactions on pattern analysis and machine intelligence, 35(8):1798- 1828, 2013.
[7] Emmanuel J Candès, Xiaodong Li. Yi Ma, and John Wright. Robust principal component analysis? Journal of the ACM (JACM). 58(3):1-37, 2011.
[8] Amelia Perry, Alexander S Wein, Afonso S Bandeira, and Ankur Moitra. Optimality and sub- optimality of PCA Spiked random matrix models. The Annals of Statistics, 46(5):2416-244 2451, 2018.
[9] Adel Belouchrani, Karim Abed-Meraim, J-F Cardoso, and Eric Moulines, A blind source separation technique using second-order statistics. IEEE Transactions signal processing, 45(2):434 444, 1997.
[10] Emmanuel J Candès and Terence Tao. The power of convex relaxation: Near-optimal matrix completion. IEEE Transactions on Information Theory, 56(5): 2010.
[11] Emmanuel Candes and Benjamin Recht. Exact matrix completion via convex optimization. Communications of the ACM, 55(6):111-11,, 2012.
[12] Charles Stein. Estimation a covariance matrix. In 39th Annual Meeting IMS, Atlanta, GA, 1975, 1975.
[13] Akimichi Takemura. An orthogonally invariant minimax estimator of the covariance matrix of a multivariate normal population. Tsukuba journal of mathematics, 8(2):367-376, 1984.
[14] Joel Bun, Romain Allez, Jean-Philippe Bouchaud, and Mare Potters. Rotational invariant estimator for general noisy matrices. IEEE Transactions on Information Theory, 62(12):7447- 7490, 2016.
[15] Joël Bun, Jean--Phllippe Bouchaud, and Marc Potters. Cleaning large correlation matrices: tools from random matrix theory. Physics Reports, 666:1-109, 2017.
[16] Farzad Pourkamali, Jean Barbier, and Nicolas Macris. Matrix inference in growing rank regimes. arXiv preprint 2301::230060//1122 2023.
[17] Florent Benaych-Georges, Jean- Phiilppi Bouchaud, and Marc Potters. Optimal cleaning for
singular values of cross-covariance matrices. arXiv preprint arXiv:1901.055443 2019.
[18] Emanuele Troiani, Vittorio Erba, Florent Krzakala, Antoine Maillard. and Lenka Zdeborová. Op- timal denoising of rotationally invariant rectangular matrices. arXiv preprint arXiv:2203.07752 2022.
[19] Farzad Pourkamali and Nicolas Macris. Rectangular rotational invariant estimator for general additive noise matrices. In 2023 IEEE International Symposium on Information Theory (ISIT), pages 2081-2086, 2023.
[20] Itamar D Landau, Gabriel C Mel, and Surya Ganguli. Singular vectors of sums of rectangular random matrices and optimal estimators of high-rank signals: The extensive spike model. arXiv 2023.
preprint arXiv:2306.003:
[21] Alice Guionnet and Maïda. fourier view on the R-transform and related asymptotics of spherical integrals. Journal of functional analysis, 222(():435-49)) 2005.
[22] Florent
Rectangular R-transform as the limit of rectangular spherical integrals. Journal of Theoretical Probability, 24(4):969-99) 987, 2011.
[23] Jinho Baik, Gérard Ben Arous, and Sandrine Péché. Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices. Annals of Probability, pages 1643-1697, 2005.
[24] Florent Benaych-George and Raj Rao Nadakuditi. The eigenvalues and eigenvectors of finite, low rank perturbations of large random matrices. Advances in Mathematics, 227(1):444-521,
2011.
[25] Florent Benaych-Georges and Raj Rao Nadakuditi. The singular values and vectors of low rank perturbations of large rectangular random matrices. Journal of Mullivariate Analysis, 111:120-13: 135,0011,,111 2012.
[26] Thibault Lesieur, Florent Krzakala, and Lenka Zdeborová. Constrained low-rank matrix estimation: Phase transitions, approximate message passing and applications. Journal of Statistical Mechanics: Theory and Experiment, 2017(7: 2017.
11

---

[27] Mohamad Dia, Nicolas Macris, Florent Krzakala, Thibault Lesieur, Lenka Zdeborová, et al. Mutual information for symmetric rank-one matrix estimation: proof of the replica formula. Advances in Neural Information Processing Systems, 29, 2016.
[28] Mare Lelarge and Léo Miolane. Fundamental limits of symmetric low-rank matrix estimation. Probability Theory and Related Fields, 173:859-929, 2019.
[29] Léo Miolane. Fundamental limits of low-rank matrix estimation: the non-symmetric case. arXiv
preprint arXiv: / 707700000775 2017.
[30] Jean Barbier and Nicolas Macris. The adaptive interpolation method: a simple scheme to prove replica formulas in bayesian inference. Probability theory and related fields, 174(3): (74(3):1133-1133 133--1885 2019.
[31] Andrea Montanari and Ramji Venkataramanan, Estimation of low-rank matrices via approximate message passing. The Annals of Statistics, 49(1), 2021.
[32] Alyson K Fletcher and Sundeep Rangan. Iterative reconstruction of rank-one matrices in noise. Information and Inference: A Journal of the IMA, 7(3):531-562, 2018.
[33] Jean Barbier, Francesco Camilli, Marco Mondelli, and Manuel Sáenz. Fundamental limits in structured principal component analysis and how to reach them. Proceedings of the National Academy of Sciences, 220630) 12, 2023.
[34] Zhou Fan. Approximate message passing algorithms for rotationally invariant matrices. The Annals of Statistics, 50(1): 197-224, 2022.
[35] Farzad Pourkamali and Nicolas Macris. Mismatched estimation of symmetric rank-one matrices under gaussian noise. In International Zurich Seminar on Information and Communication (IZS 2022). Proceedings, pages 84 88. ETH Zurich, 2022.
[36] Jean Barbier, TianQi Hou, Marco Mondelli, and Manuel Sáenz. The price of ignorance: how much does it cost to forget noise structure in low-rank matrix estimation? arXiv preprint arXiv:2205. 10009, 2022.
[37] Farzad Pourkamali and Nicolas Macris. Mismatched estimation of non-symmetric rank-one matrices under gaussian noise. In 2022 IEEE International Symposium on Information Theory (ISIT), pages 1288-1293. IEEE, 2022.
[38] Alice Guionnet, Justin Ko, Florent Krzakala, and Lenka Zdeborová. Estimating rank-one matrices with mismatched prior and noise: universality and large deviations. arXiv preprint 2xXiv:2306.092833 2023.
[39] Andrea Montanari and Yuchen Wu. Fundamental limits of low-rank matrix estimation with diverging aspect ratios. arXiv preprini 12Xiv:2211.0048: 2022.
[40] Yoshiyuki Kabashima, Florent Krzakala, Mare Mézard, Ayaka Sakata, and Lenka Zdeborová. Phase transitions and sample complexity in Bayes-optimal matrix factorization. IEEE Transac- tions on information theory, 62(7):4228 4265, 2016.
[41] Jason T Parker, Philip Schniter, and Volkan Cevher. Bilinear generalized approximate message passing-pant i: Derivation. IEEE Transactions on Signal Processing, 62(22):5839-588) 2014.
[42] Qiuyun Zou, Haochuan Zhang, and Hongwen Yang. Multi-layer bilinear generalized approxi- mate message passing. IEEE Transactions Signal Processing, 69:4529 4543,, 2021.
[43] Jason T Parker, Philip Schniter, and Volkan Cevher. Bilinear generalized approximate message passing- ii: Applications. IEEE Trunsactions on Signal Processing, 62(22):5854-5867, 2014.
[44] Antoine Maillard, Florent Krzakala, Mare Mézard, and Lenka Zdeborová. Perturbative con- struction of mean-field equations in extensive-rank matrix factorization and denoising. Journal of Statistical Mechanics: Theory and Experiment, 2022(8):083301, 2022.
[45] Jean Barbier and Nicolas Macris. Statistical limits of dictionary learning: random matrix theory and the spectral replica method. Physical Review E. 106(2):024136, 2022.
[46] Francesco Camilli and Mare Mézard. Matrix factorization with neural networks. Physical Review E, 10766))064308,, 2023.
[47] Francesco Camilli and Marc Mézard. The decimation scheme for symmetric matrix factorization. arXiv preprint arXiv:2307.16564 2023.
12

---

[48] Christopher M Bishop. Variational principal components. 1999.
[49] Yew Jin Lim and Yee Whye Teh. Variational Bayesian approach to movie rating prediction. In Proceedings of KDD cup and workshop. volume 7, pages 15-21, 2007.
[50] Alexander Ilin and Tapani Raiko. Practical approaches to principal component analysis in the presence of missing values. The Journal of Machine Learning Research, 11:1957-2008 2010.
[51] Shinichi Nakajima, Masashi Sugiyama. Derin Babacan, and Ryota Tomioka. Global analytic solution of fully-observed variational bayesian matrix factorization. The Journal of Machine Learning Research, 14(1):1-37, 2013.
[52] Salma Tarmoun, Guilherme Franca, Benjamin D Haeffele, and Rene Vidal. Understanding the dynamics of gradient flow in overparameterized linear models. In International Conference on Machine Learning, pages 10153-10161. PMLR, 2021.
[53] Antoine Bodin and Nicolas Macris. Gradient flow extensive- rrenn positive semi- matrix denoising. arXiv preprint arXiv:2303.094 09474, 2023.
[54] Daniel Lee and H Sebastian Seung. Learning the parts objects by non-negative matrix factorization. Nature, 401 (0755))))788-779 1999.
[55] Marco Mondelli and Ramji Venkataramanan. Approximate message passing with spectral initialization for generalized linear models. International Conference on Artificcil Intelligence and Statistics, pages 397-405, PMLR, 2021.
[56] James A Mingo and Roland Speicher. Free probability and random matrices, volume 35. Springer, 2017.
[57] Florent Benaych-CGeorges. Rectangular random matrices, related convolution. Probability
Theory and Related Fields, 144:471-515, 2009.
[58] Mare Potters and Jean-Philippe Bouchaud. A First Course in Random Matrix Theory: For Physicists, Engineers and Data Scientists. Cambridge University Press, 2020.
[59] Roger A Horn and Charles R Johnson. Matrix analysis. Cambridge university press, 2012.
[60] Jean Zinn-Justin. Quantum field theory and critical phenomena, volume 171. Oxford university press, 2021.
[61] Dennis S Bernstein. Matrix mathematics. In Matrix Mathematics. Princeton university press, 2009.
[62] Harish-Chandre Differential operators semisimple lie algebra. American Journal of Mathematics, pages 87-120, 1957.
[63] Alice Guionnet and Jonathan Husson. Asymptotics of dimensional spherical integrals and applications. ALEA, 19:769-797, 2022.
13

---

A Posterior mean estimator is in the RIE class
In this section, we show that for rotational invariant priors, the posterior mean estimator is inside the RIE class. For each of the estimators of X, Y. we present an equivalent definition of the RIE, then we show that posterior mean estimator satisfies this definition.
A.1 Estimator
Lemma 1. Given the observation matrix S. let X (S) be an estimator of X. Then X(S) is RIE if and only for any orthogonal matrices E RRNX NNN E RMMM;
=
(22)
Proof. If (S) is RIE, then the property 22
clearly follows from the definition 2 Now we turn
to the converse.
Suppose that an estimator (S) satisfies 22 First, we show that if the observation matrix is diagonal, then the estimator also diagonal. Consider the observation matrix to be Sding diag(s),.... 8K) ONx(M-N) ]. Let I E RNNN J E I M MM be diagonal matrices with diagonal entries all one except the k-th entry which is -1. Note that for S < N. we have 5diag_ J Je Moreover, matrices I JE are indeed orthogonal. For any S k S N, from
the property we have:
X (Sdias)) (1, Soiea = 1 X (Srdiae)
(23)
This implies that all entries on the k-th row and k-th column X (Sddis)) are zero except the k-th entry on the diagonal. Since this holds for any h;, we conclude that (Sdiia)) is diagonal.
Now, for given general observation matrix with SVD S UsLV put U =UT, V = V in the
property 22 We have:
- -
From the argument above, the matrix the Ihs is diagonal. Consequently, the matrix (S) diagonal which implies that the columns of Us are eigenvectors of X (S). Therefore, X(S) is RIE.
Now, we prove that the posterior mean estimator (*)S) E X S] has the property (22). and therefore belongs to the RIE class. For simplicity, we drop the SNR factor Vk. For any orthogonal
matrices U E NNN V R RMM we have:
f JY dx X (Y XY)
E X USVT
f dY dx (Y) XY)
say dx UXUT UUS
(9 JX dx PXX)) -XXUUYV)
(0) [ JYY (ON
J dx (XX Jo
SE
where in (a). we changed variables X + UXUT, used det U| 1, and rotational invariance of Px. (X) = (U In (b), we changed variables Y + UYVT, used det U| det =1 and bi-rotational invariance of Py, (Y) (UY In (c), we used the bi-rotational invariance property of Pw. namely Pw(USVV - UXYVT) = -
14

---

A.2 Y Estimator
Lemma 2. Given the observation matrix S. let (()) be an estimator for Y. Then Y()) is a RIE if
and only for any orthogonal matrices U E RNNN V E RRMM.
(24)
Proof. If (S)) is RIE, then this property clearly follows from the definition 4 Let us now show the converse.
Suppose that an estimator (S)) satisfies (24 First, we show that if the observation matrix is diagonal, then the estimator is also diagonal. Consider the observation matrix to be Sding ONX(M-N) Let I E RNNN, J E MxM AM be diagonal matrices with diagonal entries all one except the k-th entry which is -1. Note that for I I N, have Sullag = I J Moreover, matrices : J are indeed orthogonal. For any S S N, from
diag(s),....
the property we have:
( (I, Sciiss - I (S JR
(25)
This implies that all entries on the k-th row and k-th column of Y (Stllas)) is zero except the k-th entry on the diagonal. Since this holds for any h we conclude that (Sella))) is diagonal.
Now, for given general observation matrix S UsIV put U UT, V V5 in the property
24). We have:
Y = USV
From the argument above, the matrix the Ihs is diagonal. Consequently, the matrix USY is diagonal which implies that the columns of Us, Vs are the left and right singular vectors of (S). Therefore, Y (S) is a RIE.
Now, we prove that the posterior mean estimator Y-)()) = B has the property 24 and it inside the RIE class. For simplicity, we drop the SNR factor vx. For any orthogonal matrices U E N V E MXM, we have:
dX dÝ Y (V -XX) XY)
E U
d dY PY)X) XX)
[ dX dÝ UÝ
f NÝ (V)
1 |XX DÝ
J dY
S dx DY Y (X
d AY (Y) xY) Ir
=
where in (a), we changed variables Y + UYVT, used det U det V = 1, and bi-rotational invariance of Py, (Ý) Py VU) In (b), we changed variables X UXUT, used det U = 1, and rotational invariance of Px. (x) = (XXUX In (c), we used the bi-rotational invariance property of Pw, namely (USVV UXÝVT) = Pw XY).
B The replica method for deriving the resolvent relation
In this section we present the replica method used to obtain the resolvent relation. For simplicity of notation we use G(z) Gz(z) for the resolvent of random matrix S.
15

---

assume First, we express the entries of the resolvent (22) using the Gaussian integral representation of an inverse matrix 60
Gi)(2) - V N2+) del (=I - -S) _S) /(
M+N
II dms) nhlli exp - 7
1 [5]]
M+N
II Ar
dma) nonl exp { n" (zI s)n}
(26)
M+N
( II
dml)
exp { n (=I - s)n}
For not close to the real axis, the resolvent is expected to exhibit self-averaging behavior in the limit of large N, meaning that it will not depend on the particular matrix realization. Thus, we can examine the resolvent Gz(z) by analyzing its ensemble average, denoted by (-) in the following.
/
M+N II
1 dmF)) nal exp { 5"d (=1- (21-5)m3)
(G](t))
(27)
z
where is the denominator in 26). Computing the average is, general, non-trivial. However, the replica method provides us with a technique to overcome this issue by employing the following identity:
M+N II
(Gig(z)) = lim
zn
/
dne))
7:nyy exp {
2 1 (=
n-0
k=1
(28)
M+N
lim
II
dm (T)
(1) (1) mi
mj exp
(r) (2I - (5)
So, the problem now is reduced computation of an average system (26). After computing the average value (the bracket) in (28), we can perform an analytical continuation of the result to real values of and then take the limit n 0. Throughout, as is common in the replica method, that the analytical continuation can be done with only n different sets of points. Of course, this is totally uncontrolled step that comes with no guarantees.
copies (or replicas) of the initial
C Derivation of the RIE for X
In this section, we consider estimating X, and treat both Y and W as noise. We consider to be fixed, and the observation model:
S= =UXUY U.WV
(29)
where Y. W N are fixed matrices with limiting singular value distribution Ny.Pw: and U1. U2 E RNN, VI. V2 E MMM are independent random Haar matrices.
Construct the hermitization E M from as
ONXN ST
S 0 x
S =
For simplicity of notation, we use T XU,YV, T (+++)) the hermitization of T, and W denotes the hermitization of the matrix U, W.T.
16

---

C.1 Resolvent relation
We want find a relation between G(2) = Gz((), and the signal matrix X. From (8 we have
W+M
(G))(z)) = lim
/(
II k=1 v=1
II dmk
(r) ) (1/))11)) (exp {
(1) (1)
1 Σ
(T))
-II- (5))
2 T=1
U....
/
N-A1 n
lim n-+00
II dmp ( Mi m)
(1)
exp Σ n
(7) n {
1=1
( Σ (- Tn/ (n(
1
(xxp( 1
cxp ( Σ (T) Wn)()))
exp
U...
U2.V3
(30)
Split each replica (7()) into two vectors (+)) E BA, (6) E RRM, ./// = (T) [[
The exponent in
b(T)
the first bracket 30 can written as:
(T) a!T) T XU YU,YV + + - VYYUU
TT XUGYVT
(31)
6 (T) YU
Using the formula for the rectangular spherical integral
[2
(see Theorem
H.T.. we find:
N
n
exp Trb()) XU,YVT))
O(a) By
NM X ( |
exp
U.
(32)
with 0 (z) Se (1) dt. In 32 we used that ((t) (+) X is rank-one matrix with non-zero singular value [ |||||| |.
Similarly, for the second bracket in 30 we can write:
1 (Tn
= a ( U,,VT6 + b V22 V.
2a (n U2WVVT
(33)
b U, WVV
which using the formula of rectangular spherical integrals, implics
N
exp Σ b
(n) U2WV
exp
ola)
la (1(7)19166))1111
(34)
U2.42
2
NM
From
30
we find:
N+M n /
(G))
lim
II dm
(7)) 1/ m
n 700
r=1
||a (11()112116777
112
x exp
Σ
2 1|
No pa
NM
(35)
Now, we introduce delta functions (p (r) - ||a)) ), (pb5²)
116(2112) and 56 ((p(5)
and using them, the integral in (35) can be written as (for brevity we drop the limit term):
N+M
(Gi(())) = / II II
dmk (7)
dp (T) dp (T) dpp))) m m
8 (pb)) |1a
|
|xa(t)113
x
8 (p(T))
8
P3
(36)
N
M
N
1=1
x expp
Σ
=|| (1
NOW (( ((p)))NQ0c
NQ(a)
Hw
1=1
17

---

In the next step, we replace each delta with its Fourier transform, 0(p *11a-113) x s dci exp
(ei +1|a=113) J. After rearranging, we find:
(G)( & / (II
dpi
dp2
dp3
dC1
dc2 dcs"))
T=1
N 2
+ 0 (p) pb
-
- - 15 Ds
(r)
exp
(p2 P3
N++M
x
dm
1i M2
k=
T=1
Σ
[| () 2
-
-
exp
-
(37)
The inner integral in 37 a Gaussian integral, and can be written as:
N+M
f
II dnl m ")
(1)
k=
2 (T)T
(1) (2- 1
IN IN
(r)
X2
17(1)
x exp
(a - 2 )
(7)
(38)
Denote the matrix in the exponent by ( Its determinant reads:
det C (z c (7) M II [---<1--- -5--
where A 's eigenvalues of X. So replacing the formula for the Gaussian integrals, 37 can be written as:
(Gig(z)) x / (II dol do (r) dps ac(T) ach act.)) (cl) ij
Nn F. PP,, S2; Ca) (39)
x exp { -
2
with
N
Fo (Pr.P, P3. S1, C2,
Σ
In(2 x1 (z) + In(= (5²))
((())
11
- o(a)) (p2 (p())
(T)
o (v) 1
(1) c (
20 + + ( (T)
(40)
In the large N limit, the integral in 339 can be computed using the saddle sooitts of the function F In the evaluation of this integral, we use the replicu symmetric ansatz that assumes saddle-point of the form:
pp = P1,
p5 = P2
ps = P3
Vr (1. (L...... :n.
- = Cı,
C = C2 C2,
C = C3
The saddle point solution of the set of equations:
cla my
(pip))
a (cf() (piPs) + cla) (pipt))
c(o) (p5p))
C
C
Pi
ps
P
(41)
p
Gpy2
p5
p3
3-
C3
xx
18

---

Now, since the relation
39 and the solutions 41 hold for arbitrary indices i.j, we can state the The inverse of Cxx and the block structure of (z(z) are computed in 112) we have (for sufficiently large N):
relation in matrix form.
sections H.2
From
(Gs(t)),,
IN + s ST
SG sG
g 1)
GsTS(3?)ST
(42)
xx2 G
0
IM
c-
With this relation, we proceed to simplify the equations The normalized trace of the upper-left blocks of (Gs(z)),, V..12
4Tp.
is:
N
N
+
N
22
7k
N
2
7
k=1
1
1
(43)
N
2/2 7
k=1
1
1
(p)
2N
2 - 7K
+ 7kk
k=1
and the normalized trace of the upper-left block in Cx is I ( (3 have pi Gas (z).
pi. Therefore, we
c3
The normalized trace of lower-right block of
reads:
1
+M N) 22 M (p) 1
1 N +
M - NN
09ps (1 a) (44)
M
=
2 22
M
k=1
and the normalized trace of the lower-right block in C -1 is 3- pz. Therefore, we have
(=)-1 p5 = a9ps (=) + (1 a) Moreover, we also have that C 0s9ag (+)+1-a'
Therefore, the saddle point equations 41 be rewritten in simplified form, which does not involve 8x2, as:
c(a) (pip))
:gg (2)-1
clo) By
(PIPS)
Ci
G = 0 asgps (3)+1-a
C
Pi
Pj
(45)
pi Gps(2)) (2),
p agps (z) (1 - a)
P =
S Gps
c3
Note that Çi. C2 can be computed from the observation matrix, and we only need to find C3 satisfying the following equation:
(2 - pi-Ci (2) 1- c(() 1 [aGps(z)
1- 1-1 [(= ((=-<i)gpss (i)
11)
(46)
C3
C.2 Overlaps and optimal eigenvalues
We restate the relation between the resolvent and the overlaps from the main text 18 For ã; [x:, TT with I, eigenvectors of X, we have:
a! (Im Gs(x ic)) n. (x)(( (s)
(47)
From
42
we find:
Ox(7,d:) ~ TT)s(>)
Im
lim
a.. 55 -
Gx
-
2 =-+y-10
(48)
1m
lim
TT)(())
2- 10
- Ci - C2
19

---

Once we have the overlap, we can compute the optimal eigenvalues from 14 in section 5 Note that, until now we had absorbed vi into X. Therefore, we should use with O x (7, ( This leads to:
N
E ~ N
Σ
==1
1 N
Xj
Im
lim
(T;)
10+
N j=1
2
2i
-
kgnd
1 (T)(
1 1
X1
Im
lim
Σ
2
10+
KC3
N
2-
X3
j=1
KC3
(49)
Im
lim
1
1/1
1 1
1
1 1
Σ
N
(T)I(())
s 7
10+
C3
2
3-
2
شー + X2
j=1
ACi
{
(V
2 2i
1
(V
2 Si
Im
lim
9px
9p-x
kT(is(%)
2- 10
12 53
KS3
2 53
kS3
Im
lim
1
9px
(V
(V=-(t)
+ 9px
(-1
(-v²
}}
(7)))((7))
G3
nL3
C.2.1 Estimating X²
The resolvent relation we have found in (42) is in terms of Gx.. Therefore, like other RIEs in other problems [14 19), we can express the estimator for X² without any knowledge about px or px2. One can see that, the optimal RIE for X is constructed in the same way as for X with eigenvalues denoted by ege- To compute the optimal eigenvalues, we absorb vk into X and we use the exact expression in In the end, we only need to divide by k to find an estimator for the true X².
S2 N
(7:
X2
Im
lim
(7:)
10+
N
2i
23
j=1
1
N
X3
Im
lim
(ps()))
*-+% 10+
C3
N
3- X3
j=1
N
3-6
23
3-62
Im
lim
11
C
(T))
2 10+
C3
N
X
(50)
j=1
1 N
2 - 72
Im
lim
Σ
(ps()))
2- -0
N
N
1- x
j=1
2 2
2-
Im
lim
Opxa
(T:(
2- 10+
C3
C32
C3
Im
lim
P3
(1:)
2- 10
1
Ci
1
Im
lim
gas (2)
(1s(1))
2
1
C
C
where in (a) we used
(1 and for (b) we used (45).
Thus, the optimal eigenvalues for X² read:
1
- Gps(2)
1
2
K (TA((T))
Im
lim
10+
C
(51)
20

---

Note that the parameters Çi. S can be computed from 45 Px2.
(6). without the knowledge of px or
Remark 3. The main barrier to find an estimator X is that the resolvent relation 42 is in terms of Gpga. Moreover, in the estimator for X, second equality in (49). we have the sum プド= -53
which cannot be written in terms of 9pxz
Remark 4. If we add the assumption that the matrix X is positive semi-definite, without any further knowledge on the prior, we can use v Sor the eigenvalues of (x)).. However, note that, this estimator sub-optimal for X V X3 (x] * ΣS= >3 (u?a;)
C.3 Numerical Examples
In this section, we will illustrate the derived formulas 42
48
and 49 with numerical experiments.
We consider matrices Y. W NRM to have i.i.d. Gaussian entries, so c/a) (z) = clo (=) = L which leads to simplification of saddle point equations (45
ci = P2.
=x) G = a2 02)
55 = Pp
(52)
p Gps (z),
p2 as (1)
a
P/ - 9ps (2)
1
C.3.1 Resolvent relation
We take x In model (29). without loss of generality we can consider X to be diagonal. In figures
5 and respectivel, we consider the X to be a diagonal matrix obtained by taking the eigenvalues of Wigner matrix and Wishart matrix respectively.
Note that μl and Gas (2) can be computed analytically using tools from random matrix theory, but the computation is highly involved. In our experiments,
instead a numerical estimation of Gas (2) obtained from the observation matrix with the help of a Cauchy kernel to compute the parameters ÇT. G3 (see section G and (58) for details the Cauchy kernel method).
Unlike the simpler models [5] for which the fluctuations are derived to be of the order 1//N, based on derivation cannot assess the order of fluctuations. However, from our numerics we observe that the fluctuations are of the order (NN). Moreover, fluctuations near the edge points density are larger (in particular for the last row in both figures 5 6). which is due to the fact that the limiting measures have higher fluctuations on their edge-points
Another observation, from comparison of figures 5 6 is that the fluctuations for the first example are relatively larger than the second one. One possible guess could be that this is due the symmetry of Px in the first example. However based on more extensive numerical observations (which we omit here) we speculate that this issue fact related to the existence of small eigenvalues of X. In other words, if X has eigenvalue 0 or small eigenvalues, we have higher fluctuations in the relation 47
C.3.2 Overlaps
To illustrate the formula for the overlap (48), we fix the matrix X and run experiments over various realization of the model (29). For each experiment, we record the overlap of k-th left singular vector of S and the eigenvectors of X. To compute the theoretical prediction, we find Gi = <3 for
7h- + where t is the average of k-th singular value of S in the experiments.
To find Çi
S3 use the set of equations
41 which for Y, W Gaussian can be written as:
Ci 1 p2:
G = pi P3
(53)
pi с
Opxe
2 1). p==
ps = z-Gi 572
Opxa ( 1)
Now we proceed to simplify the solution above:
C = +
Epxe -1
C1
7t
Si
21

---

TeC
055
3.
(a) Entry first diagonal entry in the upper-left block
6.8
0.6-
0.4
1.5
2.5
0.5
3.5
(b) Entry i 5000, last diagonal entry in the lower-right block
0.1
0.05
-0.05
-0.5
-0.1
-0.15
(c) Entry i 4999. j 5000, non-diagonal entry in the lower-right block
0.1
To.
In
0.5
33
2.0
5.5
(d) Entry 1, 5000, an entry in the upper-right block
0.5
-1
-1-
35
(c) Entry = 1, = 2, non-diagonal entry in the upper-left block
Figure Illustration
42 is diagonal matrix from eigenvalues Wigner matrix and Y,
are Gaussian matrices with N
2000, M 3000. The empirical estimate of Gs(z) (dashed blue line) is
computed for 2 7
2x for 1 < N. Theoretical estimate (solid orange linc) computed from the rhs of (42) with parameters obtained from the generated matrix. Note that, the theoretical estimate has also fluctuations because the parameters <i, G are given by the numerical estimate of Gpy (z).
22

---

03
-05
1.5
(a) Entry first diagonal entry in the upper-left block
0.6
6.4
0.3
0.2
0.1
0-
(b) Entry i 5000, last diagonal entry in the lower-right block
3e
-0.1
(c) Entry i 4999, j
5000, non-diagonal entry in the lower-right block
0.05
Be-
0.00
-0.15
0.2
(d) Entry 1,j 5000, an entry in the upper-right block
0.4
0.8
-0.2
-044
(e) Entry L,j 2, non-diagonal entry the upper- left block
Figure 6: Illustration of 42 X is diagonal matrix from the cigenvalues of Wishart matrix with aspect ratio 1/2 and Y, Z Gaussian matrices with N 2000, M 3000. The empirical estimate of Gs (dashed blue line) is computed for = computed from the rhs of 22 with parameters obtained from the generated matrix. Note that, the theoretical estimate has also fluctuations because the parameters Çī, C are given by the numerical estimate of Gag (z).
iv IN for 1 < i < N. The Theoretical estimate (solid orange line)
23

---

Sam. (k 500)
(k- 800)
Sum. (k 200)
200)
1.5
Nola:]
N (uf
0.5
22
115
0.5
0.5
1.5
Ag
(a) Wigner X
(b) Square root Wishart X with = 1/4
Figure Computation the rescaled overlap. Both Y and are matrices with i..d. Gaussian entries of variance 1/N, and aspect ratio N/M = 1/2. The simulation results are averaged over 1000 experiments with fixed X, and N = 1000, MM = 2000. Some of the simulation points are dropped for clarity.
2
C2
2i
9px
2 11 +
2 52
2 1) +
57 = 12 = _²
9xxa
- 1) = 212 + (1 (1- 7i
9xx2
(54)
C122 (1 -
1
57
px2
2
(si² +
+1 )
Si
Rxxx
51
Ci2 + (1
1)
a
Thus, Ç is the solution to (54 For each example, we solve this equation and compare the obtained theoretical overlap against the average over the experiments.
Wigner X. Let X E N be a Wigner matrix, then (2) = I-
1 Solving (54) we can
compute the overlap using (48) In Fig. 7a we compare the theoretical computation with simulations for N 1000, M 2000. As in previous cases ps( is approximated using Cauchy kernel [58
Square root Wishart X. Let X IR N be the square root of a Wishart matrix X
HH
with H e R N having i...d Gaussian entries. Then Rp (2) = a' - a' = N/v' Solving (54) we can compute the overlap using (48). In Fig. 7b we compare the theoretical computation with simulations for N 1000, N' 4000, M 2000.
C.3.3 RIE performance
In this section, we investigate the performance our proposed estimators for X. We compare performances of the optimal RIE (49 with the one of Oracle estimator (3. Moreover, we illustrate
the performance of the estimator for X (50). and the sub-optimal estimator of X derived from
remark 4
24

---

For Y. W with Gaussian i..d. entries,
511
simplifies to:
1
Si
1
Im
lim
Gps (z)
(T)(()))
2-
10
63
1
Im
lim
gps
-9ps ci
K. (Tūs(())
2-%--0
1
1
Im lim
Gps(2) (z) Gas(2) Gās (3) -
(Tjjs(())))
2- 10-
Ops (z) + 1-a
1 kTAs(())
Im
(*H(ps](71) + xitūs(?ii)
(+Hāūs((7)) + xi))((7) + l-a
- -(HH[ūs](7)) + xips(())) aH|js|(>)) + Tips (1r) + 1-a
1
mūs(T)) -1++ (-1+ a (TPp)(())2 (TH(ūs](<i)) + -1+1 2
K (Ts((T))
K
[))))))) (TH(ūs]((s)) + -1+)))
(55)
For our first example, we consider two priors for X:
Shifted Wigner X. We consider X F + cf where F = FT RNN has i..d. entries with variance I/N, and e * 0 a real number. Then, the spectrum of X is a shifted version of the Wigner law
(x(A) V (A- c)²
for < <
2T
and the Stieltjes transform reads:
9px(2) -- c-c- V 1 2- () + 22-0
Wishart X. Take X = A HHT with H E Nxx' having i.i.d. Gaussian entries, with N/v a' < 1. Then, the spectrum of X is the renowned Marchenko-Paue -Pattrr distribution:
M (o 1) (V2 +1)² A
xx (2)
for -- x
Va + +1)
2x
Va
and the Stieltjes transform reads:
( 1) V[= Val 1)
= - -va 1)
9px (2)
2z
In Figure & the MSE of Oracle estimator, RIE (49). and X2-RIE is illustrated for shifted Wigner X with 3, and Wishart with aspect-ratio a' 1/4. We see that the performance of RIE is close to
the one of Oracle estimator, which implies the optimality of the proposed estimator (49 Moreover,
we observe the sub-optimality of estimating X using
E (S). Note that, in the low-SNR regime,
the estimated eigenvalues E might be negative which makes the estimator
(S) undefined,
so the MSE is not depicted in this case.
In Figure the MSE of estimating X² shown. We see that in the high-SNR regimes the RIE 55 has the same performance as the Oracle estimator.
25

---

Oracie estimator, (q(S)
W Oracie estimator, (8(
0.08
RIE. (5
0.12
RIE. (5
Ex) (8)
Exx (8)
MSI
0.07
MSE
0.06
6000
(a) Shifted Wigner, c c= 3
(b) Wishart, a' = 1/4
Figure 8: Estimating X. The MSB is normalized by the norm of the signal, || Both Y and W are N x M matrices with i.i.d. Gaussian entries of variance 1/N, and aspect ratio N/M = 1/2. The RIE is applied to N = 2000, M = 4000, and the results are averaged over 10 runs (error bars are invisible). Average relative error between RIF (S) and Oracle estimator is also reported.
Oracle cstimator, (
V Oracle estimator, (S)
0.26
0.21
KIE, S ()
A KIE. E (S))
0.24
0.22-
MSB
0.22
0.18
0.18
0.10
1495
0.16
sire Mns
2
(a) Shifted Wigner, c c= 3
(b) Wishart, a' 1/4
Figure 9: Estimating x²_ The MSE is normalized by the norm of the signal, || Both Y and W are N x M matrices with i.i.d. Gaussian entries of variance 1/ I/N, and aspect ratio N /M = 1/2. The RIE is applied
to N = 2000, M 4000, and the results are averaged over 10 runs (error bars are invisible). Average relative error between RIE E (S) and Oracle estimator is also reported.
Bernoulli spectral distribution. In this case, the matrix X is constructed as x AU with Ux a N N orthogonal matrix distributed according to Haar measure on orthogonal matrices, and = diag(A A) where X has i. Bernoulli elements. Thus, Px = pōo + (1 - p)001 for p E (0. 1), and the Stieltjes transform is:
1
(xx = 2
+ p) -
For this prior, have that X = x2, so both estimators (S) and (S) should have the same performance. However, note that x (S) does not use any knowledge of px. In Figure 10 the MSE illustrated for these two estimators for two sparsity parameter, p = 0.5 and 0.9. We observe that, except for for the low-SNR regimes, both estimators have the same MSE. The poor performance of xra (S)) in the low-SNR regimes might be due to the fact that, some of the estimated eigenvalues S are negative although the true eigenvalue is O. This makes the estimation more difficult for the sparser prior, see Figure 10b However, this problem is resolved in E ES) by taking the knowledge of Gpx (z) into account.
Effect of aspect-ratio a. In Figure 11 we consider X to be shifted Wigner with 3, and the MSE is depicted for various values of the aspect-ratio a. As expected, as M increases (a decreases) and we have more observation or more data samples, the estimation error decreases.
26

---

U.S, Ornole (
U.S. Chacke (8)
0.8
0.5, RIB (8)
0.8
1.5
0.5. RIE E (S)
(1.,, (S)
Oracle (8)
(1., (S)
337%
p
RIE Ex
0.6
MSB
MSE
3200
0.
0.4
0.4
100
00.2
0022
(a) Estimating X
(b) Estimating X²
Figure 10: Estimating X and X² with Bernoulli spectral prior distribution. The MSE is normalized by the
norm of the signal, || = 2| Both Y and W are N x M matrices with i.dd Gaussian entries of variance 1/N, and aspect ratio N/M- 1/2. The RIE is applied to N 2000, M 4000, and the results averaged over 10 runs (error bars are invisible). Average relative error between RIE E (S) and Oracle estimator is also reported.
0.07
9:95
Oracie estimator. Ex (S) RIE. (S(S)
0.06
MSB
0.50
0.05
0225
0.04
000%
090
0.6
0.1
Figure 11: MSE estimating as function of aspect-ratio a. prior on is shifted Wigner with c = and K = 5. MSE is normalized by the norm of the signal, || Both Y and W are N x matrices with i.i.d. Gaussian entries of variance 1/N. The RIE is applicd to N = 2000, M = 1/aN, and the results are averaged over 10 rums (error bars are invisible). Average relative error between RIF (S) and Oracle estimator is also reported.
D Estimating Y
In this section, we present the derivation of the optimal RIE for Y. For simplicity, the SNR parameter in (T) absorbed into Y, so the model is S XY W. Therefore, the final estimator should be divided by 1/vR to give an estimate of the original Y.
The optimal singular values constructed E Σラ o, (uTv, ()) (vTu, We assume that, for large N, Eg can be approximated by its expectation:
N Σ oy E j=1
where the expectation is over the singular vectors of the observation matrix S. Therefore, compute the optimal singular values, need to find the mean overlap E (u. y/ (v7,) between singular vectors of Y and singular vectors of S. In the following we will see that (a rescaling of) this quantity can be expressed in terms of i-th singular value of S and j-th singular value of Y (and the limiting
measures, indeed). Thus, we will use the notation Oy (7;, oz)) = NE (u. y (1) (vTv, in what follows. In the nest section, we discuss how the overlap can be computed from the resolvent of the Hermitized matrix of S.
27

---

D.1 Relation between overlap and the resolvent
Construct the matrix S E (NN+M)XNN+M) from the observation matrix:
0 ONXN ST
S OMMM
s
By Theorem 7.3.3 in [59].
S has the following eigen-decomposition::
Üs
Üs s
Fx
0
Üs
Üs (
s
s
0
- N
0
(56)
(
V
(2)
0
0
with V3- = (1)
V (2)
in which VS" E RMXN, And, V
Ü== = +2 Us.
S
v2 $
Eigenvalues of S are signed singular values of S. therefore the limiting eigenvalue distribution of S (ignoring zero eigenvalues) is the same as the limiting symmetrized singular value distribution of S.
Define the resolvent of S
Gs(z) = ( =I-S) )
Denote the eigenvectors of S by s; E ERMN i= 1= 1,..., M + N. For 2 =
ic with IR and
6 >> N we have:
2N
N+
Gs(r - = x
SuSt
+
x + ie
Σ 2=2
S
)2 2
12
6
where 2 are the eigenvalues of S, which are in fact the (signed) singular values of S, f1 =
The IN = 7N, 7N: 7N+1 -...... 72N YN.
O2 Ui (r)
(1)
Define the vectors r, =
11 =
y
for " y; right/ left singular vectors of Y, we
00
have
2N
N+M
o. Ga(nGss--iit 16)]4=
(r.st)) (1T8k) +
x x ie
Σ
(rTsk) (1T
(ar t10 2 + 2
x2 +
k=1
2=
(57)
Given the structure of ss
(56). we have:
2 (uTyy ()) (v]y!")
for 1 k N for N W+1Sk$2N for 2N WW+IkkNNN
(rTss) (4Tsa) =
uh N (vE
In the limit of large N, the latter quantity is also self-averaging, due to the fact that as N co, these overlaps exhibit asymptotic independence, enabling the law of large numbers to be applied here. We can thus state that:
r. Gs(e - ic)) (e)) N-000
/ (x - 1)2 2 Or (1)ss)) di
(58)
where the overlap function (t, A1) is extended (continuously) to arbitrary values within the support of jás with the property that Oy (1 >x) (1, >s) for supp(us) Sending e - O, we find
r: (Im Gs(c i6) mūs Or (x.
(59)
In the next section, we establish a connection between the resolvent Gs(z) and the signal Y. which enables to determine the overlap and consequently the optimal singular values values Eg in terms of the singular values of the observation matrix S.
28

---

D.2 Resolvent relation for Y
In this section, we consider estimating Y, and treat both X and W noise. We consider the model to be:
S= OXOXY + +UWV
(60)
where X XT E W E RNM are fixed matrices with limiting eigenvaluc/singular value distribution PX, and O,U E R V E RRMx are independent random Haar matrices. For simplicity of notation, we use T = OXOTY, and T E MNN+ the hermitization of T.
And W denotes the hermitization of the matrix UWVT.
As in the case for X, we express the entries of G(z) = Gs(z) using Gaussian integral representation, and after applying the replica trick 28 we find:
(G)((z)
lim
/
II k=
N+-M
II
dmss
(T) (1) (1)
Ths "T1 exp Σ (T)
O.U.V (=l -8)((²1)
10000
N+M
lim
dr (7)
(1) (1)
f3 exp {
Σ
np (r) m7 (
n
1=
r=1
==1
exp ( Σ (T) Wn u/vv 17 (7) ( T (3) (
exp
2
(61)
a(r) b(r)
Split each replica v7 into two vectors a (+)
E R2,, b(T) RMM (,M,)(T) =
The exponent in
the first bracket in 61 can be written as
Tn (T)
(T) OXOTYb + r'XO0
TrOXOT Yo (T) 6 YT
(62)
where (T) is a symmetric N N matrix with two non-zero eigenvalues (+)) Y ||]))||| by lemma
1
Using the formula for the spherical integral
21
(see Theorem
we find:
1 T)OXOT(())
N
N (alt) Yoo) + ||a (11 (1))
exp 2
exp
Px
2
==1
+P
px
N
(al:))
yyb - |a 1|
(63)
By the same computation as previous section, for the second bracket we have:
exp
Σ Tr ( (WVTT)
x exp
N
Olo) pw
la 1191611111111111111
(64)
U.v 2
NM
From
61
63
64
we find:
N+M M
(Gt(())
= lim
dnj 113
(T)) (1) (1)
n n
k=1 rml
|aa (7) (11111111111 (lllsi
2||n 20)
No 1
1
X exp
2
NM
(69)
NPpp (alt) ||a -||a ( 
-| ||ا>|||| Soleeaan
( 1
((())
-NPpx
29

---

Now, we introduce delta functions (for brevity we drop the limit term):
/ N+M
(G(b(z))
II
dp (T)
dap (T) dap ( del ni 1 (1)
km
r=1
0(1 (T)
1 N
la 11)) 5(1 (1) M (6(2)112)
x
x x0 50 1 N
Yb (z) 199 8 ((01²) N
1 a (r) ( Yb(
(66)
2/|1 11p
N Now (qi %2
(T) )
x exp
NPpx (9) 1i
1 )() NP px (0l VV )
(7)
In the next step, we replace each delta with its Fourier transform. Note that for the parameters { T (qï
Q1 92, we use 608 11107113) x I dBf exp
* I| [11a_11)). and for @4 we use
5 5(1 (T) Tyb()) x f d&p exp {
NBT (11
(T) ( TYb()) After rearranging, we
find:
(G.)())) x /
n dap (T) daz dis (1) de asp a3g agp²) d3$"7)
(7)
N Σ
(T) (t) x) Oooi (gf (7) ((( + Ppx Pxx v (11gg /)7)) + Ppx (1/ q/
x exp
- -3/- 91
(7)
1
32
4/
// (7) - -28(- {1/)}
a
N+M
/
II dul
(1)
m fjj exp {
(1) (1)
1
21|m
2
[1)]
3/ la 11 - 36 1
2
8/471 |
a [Y66))
(67)
The inner integral is a Gaussian integral, and can be written as:
J(I N+M
II di
(r) mi m3
(1)
k=
(68)
(2- I
- Y 33- - YTY
n1)
X exp
Denote the matrix in the exponent by c(7)
Using the formula for determinant of block matrices (see
proposition 2.8.4 in 6T]). we have::
det c(7) = det - (==of()) IN - A(T) 2
((= ) IM
88 7) PYTY)
x det [(2 8 \M - -MM-- r
N
N
2 8/()
8//12
(= x - /// x (z - - M-N
k=1
z
T)
$ ok
k=1
N II
(z (/T) M
[ 7)7)((2
(2 - 1/57)
2 25 (7) Po
k=1
N
(2 - // M- N
(2 8f (= و6ا - (8(112 + + (2 - -
30

---

(7) where are the singular values of Y. So computing the Gaussian integrals, (67 can be written as:
(Gi)(2)) & / (II
R
dui
( d dp
da d3 dsl dsp d3(T)
==1
Nn Fo q1 92. q3 qı B1. P3. Ba)
exp {
2
(69)
with
(q1:
(q1: 92, 43, 12, B1)
1
Σ
(
1
1) In(2 -111)
n
a
N
Σ
(2
(/-h) ) (8/-)- + +13
(70)
k=1
O 92
(T)
(r)
Ppx (94 + V 9 (T)
Ppx (9 - V / 93²))
8
3 92
+ 8/7 43 + 2 94
(T)
We will evaluate the integral 67) using saddle-points of the function Fo From the replica symmetric ansatz at the saddle-point we have:
q 91,
9/ 92
43 = 93,
44 94
Vr E (1. m)
(TT B1,
BI = B2, 83 B3,
$I 34
Finally, we find the solution to be:
B = (io) (5)
23 + 28 V
RRxx (0a + /
Rpx (9i v9993)
B c(u) (5)
1 qi
R (9) + v (993) - RRx (2ā - (
2
V
98
B 2 Rpx (gi + / + Ryx (yi v
1
[1(2) = (2- (2-37) B5)
with
[z((2) = 312 + +,7(2 B7)
(x- D7
9pv (1(=)) + 87
J Zz(3)*
[2(3))
Z2(z]
ž-Pi Opy [2[z]]
[1(=)
+ 1-a
6 a 22(3)
2-3;
G5 (x- (z)
9pr (1)=)
2-37
22(z)²
Z2(=))
2z(z)
O = 22(3)= 9,y
(2)
(1(2)
Zz(3)
Za(z)
(71)
where py is the limiting eigenvalue distribution of YYT.
The relation (99 and the solutions (7 hold for arbitrary indices i, j. so we can state the relation in
the matrix form. Computing the inverse of Cy (see section H.2
we have:
([
I + SG ST Gsts(22)ST ST
SGSTs(22) 2Gsts(22)
(Gs(t))
IN +
YG
Z7(=) Zz(z)
YT
YGrry (1(t)
23))
2 8
Zg(z)
(1)=) YT
Gyry (1(=)
2x(z)
G Zz(z)
Za(z)
(2x(z)
(72)
With this relation, can further simplify the solution
We start with comparing the trace of upper-left block in (72). The normalized trace of the first block
in (es()))o.uu) is computed in (43 to be Gps (z). The normalized trace of the upper-left block in
31

---

Cy-1 is:
1 N
31²
(z) (<)(z)
Tr (= ---1% -11
IN+
YGYTY
YT
(= (2)
1
Σ
N
[1
812
o
N - 81
Z2 2
(8)
k=1
(2)
K
1
(2)
31
Σ
+1
N - B4
(2(2)
48
o
Zz(2)
(73)
Z1(+)
1 1 (1)
N
1
B3 (1)
Σ
N 57 Zz(2)
21(=)
57
Zz(z)
k=1
(z(#)
0%
[1(2)
88
9pv
+
Z2(2)²
Zz(2)
Zz(z))
= qi
Thus, of = -ppe(z)
The normalized trace of the lower-right block of (Gs(2)) is (s (2) + (1 a) (see (44)
The normalized trace of the lower-right block in Cy is:
= 83 GYty
21(2) ] Zz(z)
1 8
1
M - N 87 72(2)
Tr
+
M
Zz(2)
M Zz(2)
(1)=
M
Z2(z) Z1(2)
R=1 Z2(z]
ok
N 1 2 3
N
1
M N 2 - 89
M N Z2(2)
Σ
(1
2
+
M Z Z1(2)
(74)
k= Zz(=]
0k
19 9pv
71(2)
1-0
(z(
(z)
2 3
9/
So, 4/ - - (z) + (1 a
With a bit of algebra, we can express the parameters 93, 2 in terms of gī.
,7, 81:
% - 37)
2
B
- , 1
(75)
31
3
B
Therefore, the solution can be written without involving Spy as:
B ((o) (2)
+ 1 V o RRx (aā + /
Rpx (gt / 997%))
(gi
65 =
aw
85
Rox (gi + /oius) -R,x (4ī /
_ Rox (gs V + Rpx (at / v
(76)
qi GoS (z)
C = a Gas (2) + (1 a)
T 2712 Gng (=)
2-8-
a = Dj Gns (2)
Remark 5. The simplifications in (75) are derived with the assumption that 8i + O. However, in the initial set of equations (71), px is symmetric measure then 3 vi = 0 is a solution. pf is symmetric, then Rpx (-x====)---- - (z), and plugging o = the expression for 3 (1). we find that 37 0.
32

---

D.3
Overlaps and the optimal singular values
From
59
72
we find:
0r(7, =0
Im
lim
31 [2(2)
Vi
(T) Gyry
(z)
YTy:
TT)(())
2-+7-10
(2(2)
(77)
Im
lim
81 Zi(2) Z2(z)0
0,
Tūs(1)
2-by-10 10
From the overlap, we can compute the optimal singular values:
N of o
1
Oy (3)
N
N
of
Im
lim
31
N
(1(2)
Zz(2))
j=1
8
N
of
lim
(7)
10+
N (z(=)
(2)
(78)
j=1
(x(=)
J
N
21
lim
19
Σ
Z2(3)
xūs((1)
10
N Z2(z)
3=1
Zz(2)
(1)
0
1
(1 Z2(z)²
Z1(2)
B4
Im
lim
9pr
(ūs(
Zz(2)
Zz(z)
Im
lim
9i
(N))
-10+
where in the last equality we used the solution we have found in (71). Note that, based on (6)). we do not need to have any knowledge about py to compute 27. In the end, we need to divide the estimator by vK as we have absorbed it into Y.
D.3.1 Recovering the rectangular RIE for a denoising problem
Note that if in the model (60). we put X the model reduces to the additive denoising of Y, and recover the estimator recently proposed in [99] for the rectangular case.
For X
Rox (2) 1, so 76 reduces to:
Bi c(o) (9195)
c(o) Aay
(9] 9
8 =
33 0, 3
43
9ी Gps (2), O a (2) + (1 a)
(79)
g (z 2 (z)
-81), 07 (= Bi)9ps -
33

---

From
78 we have:
Im
lim
94
(1))
10+
Im
lim
2 -2) (2) -1
(P)(
2 10+
c(a)
(qi 42)
Im
lim
=9ps(3)
Juw
gps(2)
(T))(1))
10+
Im
lim
29ps( (2) -C(a) (9i92) -
(T))
2- 10+
Im
lim
29ps(2)
-
cla)
(=)
Im
79ps (7 100-i0+)) 100+)-C0 cla) Gp (7 10 (1- an ( -0+
Tūs (7:)
Im ((aa wa
1 1-a *H|ūs|(7)) + a ( (>x))
74
(T)((7))
7
+ +imps(1)) - + 20mH|(s](<x))
%
(80)
where in (a) we used the analyticity of rectangular R-transform [7] and in (b), we used Plemelj
formula 66. Note that, the final estimator should be divided by the vk.
D.4 Examples
Throughout the numerical experiments, we consider the matrix W have i..d. Gaussian entries with variance 1/N, so c(a) (2) 1 2. And, X = F + cI where F FT E N has i.id entries With these choices, the solution
with variance 1/N, and e j is real number, so R,x (z)
76 simplifies to:
3 1 45 + 43 = 33 = qÎ
B = qī +
a
oi 9ns (z), = a= (z) + a)
(81)
8 A1 (2)
8-8 a 32 Gps (2)
2
0
Note that in (11) 4T 2h are given in terms of the observation, so to find the solution we only need to find the parameters 95.92. q.. In (8T). one can see that we have the relation gī parameters 33,3 in terms of o T3 Th after a bit of algebra we have the following relation:
मं qt. Writing the
1 9
+ 44
(82)
2
In the expression for 9 in (8I. using (22 we rewrite 3 83 in terms of 92. 92. After some
manipulations we find that o is the solution to the following cubic equation:
2x +3cx2 + Gpp (z) (x) 1- 1----==-1-0
Gps 1 (2) 1
(n- (=
x-c
a: 9as (2) (z) =
(83)
Based on numerical simulations, pick the following root for g;:
12 30 + 6A
VB
1 2 B
(84)
with
A Gps (2)
Gns(z)
(1 - (ps(()))
-pss
B -21CAA 44/4(1) 362 + 6A)² + 54
34

---

Once we have gi, we can find Ø3 using 82 In the end, 5f ...... 34 can be evaluated. Note that, for the RIE, only qi is required. Other parameters are used evaluate the resolvent relation 722 and the overlap 777p.
D.4.1 Resolvent relation
We take λ 1. In model (60). without loss of generality we can consider Y to be diagonal.
In figure 12 Y is the diagonal matrix obtained from the singular values of a Gaussian matrix with i.id. entries of variance !/N. In figure [3] the non-zero entries (on main diagonal) Y are uniformly distributed in [1, 3]. As in previous cases, JAs. Gas (z) are estimated numerically using Cauchy kernel, from which the parameters 87 , are computed.
D.4.2 Overlap
To illustrate the formula for the overlap 777, fix the matrix Y and run experiments various realization of the model (60). For each experiment, we record the overlap of k-th singular vectors left and right) of S and singular vectors of Y. To compute the theoretical prediction, we evaluate the parameters 19, 19,99.88 34. for 2 = 7h 10 where "/k is the average of k-tth singular value S in the experiments.
In figure 14a the overlap is shown for with i.idd Gaussian entries variance N. so juy is the Marchenko- Pastur law with aspect-ratio a. In figure 14b] matrix Y is constructed Y Uy EVV. where Uy E NNx, Vy E MxM are Haar distributed orthogonal matrices, and singular values 01 ON are chosen independently uniformly from 3). so Joy L ((1,3))
D.4.3 RIE performance
In this section, we investigate the performance our proposed estimators for Y. To construct the RIE for Y. only need oi which we use (84). We compare performances of the optimal RIE (78) with the one oracle estimator (5).
In figures 15 16
the MSE of RIE and the oracle estimator is plotted for three cases of priors: Y with Gaussian entries, Y with uniform spectral density, and with Bernoulli spectral density. all cases, observe that the RIE has the same performance as the oracle estimator.
we take Y to have Gaussian entries (with variance ) ) and the MSE is depicted for various values of the aspect-ratio a. We see that as M increases (ax decreases) the estimation error (of Y) decreases.
Effect of aspect-ratio a. Figure 17
Sparse Y:
non-rotation invariant example. We consider Y to have i.i.d. entries from the -aooulll-aadmmaceer distribution,
1-p with probability
-
Yr
with probability p A ISiSN, 15jSM
with probability 1-P -
With normalization 1/VN, the spectrum of Y does not grow with the dimension and has a finite support, thus can apply our estimator to reconstruct Y. Note that the prior of is not rotationally invariant, and neither the oracle estimator nor the RIE are optimal. Therefore, taking the prior into account, we apply thresholding function on the entries of the matrix obtained from the RIE, (S). We apply the following function on each entry of the estimator:
if r
fx(x)
if |z|
for h E [0,
if
In figure 18 the MSE of the oracle estimator, RIE, and RIE+S((x) (with h = p) is plotted. A few remarks on this figure are in order. First, RIEs are not limited to rotationally invariant priors and can give non-trivial estimates for nontrotatiolally invariant priors, although they are sub-optimal. The
35

---

60-
15
25
(a) Entry i = j 1, first diagonal entry in the upper-left block
0.-
(b) Entry 1, 2, non-diagonal entry the upper-left block
02-
-0
155
0.5
(c) Entry = j 2001, first diagonal entry in the lower-right block
00-
0.5
(d) Entry i 2001. j 2002, non-diagonal entry in the lower-right block
(e) Entry = 1.j 2001, first entry the upper-right block
0.5
(f) Entry i 1.j = 2002, second entry in the upper-right block
Figure 12: Illustration of 72
R N is a diagonal matrix obtained from the singular values of a N x M matrix with i.i.d. entries of variance I/N, X = XT is shifted Wigner matrix with = 3, and Z is
Gaussian matrices with. The empirical estimate of (2 (dashed blue line) is computed
% - iv 2N 2000, .M 4000. Theoretical one (solid orange line) computed from the rhs of (2) with parameters computed from the generated matrix. Note that, the theoretical one has also fluctuations because the parameters 8f, 8 are computed from the numerical estimate of Gas (z).
for i & N.
36

---

15
(a) Entry i =
first diagonal entry in the upper-left block
-02
2.
(b) Entry 1, = 2, non-diagonal entry the upper-left block
Entry 2001, first diagonal entry in the lower-right block
02-
(d) Entry i 2001, 2002, non-diagonal entry the lower-right block
Entry 1,j 2001, first entry the upper-right block
35
( Entry 1 1,j 2002, second entry in the upper-right block
Figure 13: Illustration of distributed in [1,3], X
E N diagonal matrix with (main) diagonal entries uniformly X is shifted Wigner matrix with = 3, and Z a Gaussian matrices with.
The empirical estimate of Gz(z) (dashed blue line) computed for z Y - i + for N 2000, M 4000. Thcoretical one (solid orange line) computed from the rhs of (72) with parameters computed from the generated matrix. Note that, the thcoretical one has also fluctuations because the parameters .- Si are computed from the numerical estimate of Gas (2).
S N, for
37

---

Sim (k 800)
25-
1h. (K 800)
200)
Th. (k 200)
Niuly(") (viss)
1.5
1-
0.5
0.5
1.5
25
44
(a) Gaussian Y
(b) Y with Uniform spectrum, 24 (11,3))
Figure 14: Computation the rescaled overlap. X a shifted Wigner matrix with c 3, W has i.i.d. Gaussian entries of variance 1/N, and N/M = 1/2. The simulation results are average of 1000 experiments with fixed Y. and N = 1000, M = 2000. Some of the simulation points are dropped for clarity.
estimator, (S)
estimator, (5
0.5
RIE (S)
())
0.3
0.4
MSE
M886
0.3
0.2
045%
0293
Gaussian
(b) Uniform spectral, 2 ([1.3)
Figure 15: Estimating Y. MSE is normalized by the norm of the signal. | X is shifted Wigner matrix with e = and W has ii.d. Gaussian entries of variance 1/N, and N/M 1/2. The RIE applied to N = 2000, M = 4000, and the results are averaged over 10 runs (error bars are invisible). Average relative error between RIE Ey (S) and Oracle estimator is also reported.
LI24
0.5, Onciee (S)
p D.5, (7)
0.8
9.29%
P D.9. Oracle (S)
3.59
P D.9. RIE (7)
MSC
0.6
2.15
0.4
0445
0.2
0115
0.76
Figure 16: Estimating Y with Bernoulli spectral prior. MSE is normalized by the norm of the signal, |Y| Y has Bernoulli spectral distribution with parameter p. X a shifted Wigner matrix with 3, and W has i.idd. Gaussian entries of variance 1/N, and N/M = 1/2. The RIE is applied to N = 2000, M = 4000, and the results are averaged over 10 runs (error bars are Average relative error between RIE Ey (S) and Oracle estimator is also reported.
38

---

y Oracle estimator, 3;
à RIE, (S)
0.25%
0.00
MSE
0.08
0.07
0.8
0.6 0.4
Figure 17: MSE of estimating as function aspect-ratio a,
Gaussain entries of variance 1/x, and
K = 5. MSE normalized by the norm of the signal, ||Y fi. X is a shifted Wigner matrix with = 3, and W has i.i.d. Gaussian entries of variance 1/N. The RIE is applied to N 2000, M 1/a,, and the results are averaged over 10 runs (error bars are invisible). Average relative error between RIE E (8) and Oracle estimator is also reported.
Oraele estimator, Ei (8)
Orzzle estimator, E; (8)
RIE (5
RIE (): f,fx)
0.8
f,(x)
0.8
0.0
MSI
0.6
MSI
X
0.2
0-
5
(a) = 0.5
(b) 0.9
entries. MSE is normalized by the norm of signal, |Y| X is shifted Wigner matrix with c = 3, and W has i..d. Gaussian entries of variance 1/N, and N/M 1/2. The RIE is applied to N = 2000, M = 4000, and the results are averaged over 10 runs (error bars are invisible).
Figure 18: Estimating Y with
RIB's output can be refined, or used as a warmed-up initialization for other algorithms to get a better estimate.
In figure
19 for one experiment, the MSE is plotted for RIE and RIE+f(()) with the best h among {0.0.1. 0.1, 1.. We observe that for the particular case of Raadoull-Raaemmcheer prior, the thresh- olding stage can improve the MSE when SNR is greater than 1, however the parameter h should be chosen properly.
E Comparison of RIEs for MF and denoising
For estimating X, we have derived the estimator (4) for general priors px, /ow. This estimator simplifies greatly, with parameters in (2). when both μy, μw are Mauchenkk--aastee distribution, iec both Y. W having i.i.d. Gaussian entries of variance I/N. Similarly, although the RIE for Y in (78) is derived for the general priors, it reduces to a rather simple estimator if px, Pw are taken to be shifted Wigner, and Marchenko- Pastur distribution, respectively. Therefore, in our numerical examples on factorization problem, we consider X to be a shifted Wigner matrix, and Y, W to be Gaussian matrices.
In each experiment, the factors X, Y are estimated simultaneously using RIE from the observation matrix S. In addition to the MSE of estimating each factor, we also compute the MSE of estimating the product XY. We compare the MSE of the product with the MSE of the oracle estimator and the RIE introduced in [9 for the denoising problem. The oracle estimator for the denoising is
39

---

ES(S)
RIE, EV(S)
0.8
<RIE Au-1a))
RIE Ju-103
0.8
0.6
MSI
MSI
0.0
0.4
0.2
0-
0
2
(a) = 0.5
(b) 0.9
entries. MSE is normalized by the norm of signal. |Y p. X is shifted Wigner matrix with c = 3, and W has ii.d. Gaussian entries of variance 1/N, and N/M 1/2. The RIE is applied to N 2000, M = 4000, and thresholding function is applied with the best h among (0,0.1, 1}. Results averaged over 10 runs (error bars invisible).
Figure 19: Estimating Y with
Chaaie estimator, E, (8)
Ornale estimator, E; (8)
A RIF, E (8)
RIF. (T)
0.8
MSR
0.38
0.7
0.6
0.2
(a) Estimating X
(b) Estimating Y
Figure 20: MSE of factorization problem. MSE is normalized by the norm of the signal. X is a shifted Wigner matrix with = 1, and both Y and W are N M matrices with i.i.d. Gaussian entries of variance 1/N, and N/M 1/2. The RIE is applied to N 2000, M = 1000. In each run, the observation matrix S is generated according ( and the factors X, Y estimated simultaneously from Results averaged over 10 runs (error bars are invisible). Average relative error between RIEs and Oracle estimators also reported.
constructed as:
N
(S) = Σ
Says
I
Soy, = u XY U
(8:5)
where U, U,'s are left/right singular vectors of S. In the RIE proposed in 19]. the singular values estimated by (see section D.3.1
1
1
Im c(a)
(--a
TH]as](7)) + {)HH[[ss((x)) - (T))(
Eig;
74
V
(ps(7A)
%
+ imma((1)
+ 2amH)(ūs](ns))
7
(86)
Note that, in general the MSE of the denoising RIE E (S), is less than the MSE of the prdouct of the estimated factors E x (S).
In figures 20 2T the MSE of estimating the factors illustrated for c and c = 3 respectively.
The MSE of estimating the product is shown in figure 22
40

---

Oracle estimator, (8)
Ornale estimator, (8)
RIE. (S
0.5
RIE. EF:S)
0.08
MSI
MSR
0.3
0.08
0.2
140
0.1
0
(a) Estimating X
(b) Estimating Y
Figure 21: MSE factorization problem. MSE is normalized by norm signal. is shifted Wigner matrix with c = 3, and both Y and W are N x M matrices with i..d. Gaussian entries of variance 1/N, and N/M 1/2. The RIE is applied to N 2000, M = 4000. In each run. the observation matrix S is generated according @ and the factors X, Y estimated simultancously from S. Results averaged over 10 runs (error bars are invisible). Average relative error between RIEs and Oracle estimators also reported.
0.8
Oracke extmaatrr (8)
0.5
Oracle extiattrr (X)
ARIE, (S) * (S)))(()
A RIE (5)
0.4
(S)
MSB
MSB
0.3
0.4
0.2
0.2
0.1
0
(a)
(b)
Figure 22: MSE of the product of the factors. MSE is normalized by the norm of the signal |X || X is shifted Wigner matrix with e = c, 3, and both Y and W are N x M matrices with i.i.d. Gaussian entries of variance 1/N, and N/M 1/2. The RIE is applied N 2000, M 4000. Results are averaged over 10 rums (error bars are invisible).
F Case of a 2 1
In this section we consider the case where M S N and N/M - a aas N N
co. Throughout this section I' E N is (tall) matrix with I' M in its upper M x M block, and the rest zero entries. MM is diagonal matrix constructed from 7 E g. which are the singular values of
M(N+M)x(N+N) plays central role in deriving the RIES. For the case of M 2 N, with S = UsLV3, the matrix S has the following eigen-decompositton:
Similar to the case of < 1, resolvent of the matrix S
Ú
(1)
(1)
U
TM
ú(I)
-0(1)
U
((2)
]'
S
0
r M
S
s
(87)
Vs
0
0
0
Vs
0
with Us =
U
U
in which U E RRNXM
And, (
1 (1)
Vs = V2
Vs. The
V2
5
resolvent of S can be written as:
2M
M N
Gs(a - ic)
x
Syss
+
2 +ic
Σ
x 7x 2
2
x2
+
k=2
where io are the eigenvalues of S, which are in fact the (signed) singular values of S, 5 ......, TM 7MM714+ - F2M YM..
41

---

F1
Estimating X
The RIE for X is constructed in the same way in the case of S 1, 2 However, in the present case the observation matrix S has M (non-trivially zero) singular values and we need to estimate N cigenvalues for the RIE. As will be clear, the N .M eigenvalues are chosen to be equal.
F1111
Relation between overlap and the resolvent
Define the vectors 2 = | Our] for æ, eigenvectors of X. We have
M+N Σ
क (Im Gs(E ic) Σ (x 2 (8[Sk)
(88)
x2 +² k=2M+1
Given the structure of S& (87) have:
2 (o ux) 2
for S M
(aTs))2
2 (a. un uuuuu
for M + S S 2M
(I: ue]uu-))
for 2M + k M+N
We assume that in the limit of large N this quantity concentrates on Ox (7), A₂) and depends only the singular values and eigenvalue pairs (7g, A²). This assumption implies that the singular vectors associated with 0 singular values (u, for M + S s N) all have the same overlap with the eigenvectors of X, Ox (0, Az). We thus have:
a! ic))
N N-+00
a {
1
(x - 1)= 2x (x(ts))sss))d dt (1- x2 2 Ox (O,
(89)
where the overlap function Ox (t, is extended (continuously) to arbitrary values within the support of ïs (the symmetrized limiting singular value distribution of S) with the property that (1, As) = Ox (-1, A;) for L supp(us) Sending + O, we find
क (Im Gs(z - ie)) 1
A.) (1-)))(()00
(90)
F.1.2 Resolvent relation
We derive the resolvent relation for the same model as in 29 The derivation is similar to the procedure explained in section C.I. and we omit here. The final resolvent relation the same (42) with parameters satisfying:
((l Aa (pipi)
C (pirs)
Ci
S (ClV (pipt) +
Co/o) (pipt)),
S
(91)
pi 25 gpxz
p5
P
Again, with the same procedure as 43 (44) the saddle point equations 91 can be rewritten in simplified form, which does not involve pxe, as:
Ct ol (pTp))
((1/a))
C
C
Pi
Dng
Pg
p = gas (2) + (1 1 = (z), ps
2-S1 al3
Gos
+ 2 (
11- 1
(92)
with üs the limiting ESD of non-trivial singular values of S. Note that (t S2 can be computed from the observation matrix, and we only need to find çă satisfying the following equation:
(z (x- gms
(2) (2)
1 - -
c(/a)
Gas (=)(2 (2))==1) (2) ( (1 (1) (93)
a px
Note that both sets of equations
(90)
922 and
44
45
match for 1.
42

---

-10
estimater, (S)
0.17
Oracle estimator, (y(S)
RIB, (5)
(i(
9.5
0.10
MSE
MSE
1305
0.35
8.5
4599
0.14
1000
3
(a) Shifted Wigner, c
(b) Wishart, a' 1/4
Figure 23: Estimating X. The MSE is normalized by the norm of the signal. | P.. Both Y and W N x M matrices with i...d. Gaussian entries of variance 1/N, and aspect ratio N/M = The RIE is applied
to . = 2000, M 1000, and the results are averaged over 10 runs (error hars are invisible). Average relative error between RIE x (S) and Oracle estimator is also reported.
E.1.3 Overlaps and optimal eigenvalues
From 90
for non-trivially zero singular value of S we find:
Ox(7,A)
Im lim x-+--10
af 63 Gx2
2
Si
Tūs(7)
C3
(94)
Im lim
(1)
=-- 2 S7 - C3
And, in the case of M >> N, for zero singular values we have:
Ox (0, As) 2
Im
lim
a. 53 1
xxa
2 21
Di
(a 1)xT
S3
(95)
1
Im
lim
(a 1)TT
2--+0+
z Ci <3
Finally, the optimal eigenvalues can be derived in the same way as in 149
For 1 sis M, we have:
{ 1 [
(/>-55))
+ Spx (- (-1
/ ( }1}
E 2nT
Im lim
9px
(96)
x +0+
(33
KS3
nC8
And, for all M + +1sis
é 2x(a - 1)x
Im lim
9px
(V
x -
Gpx (-
}}
(97)
+
C
F.1.4 Numerical Examples
For matrices Y, W R XM with i..dd Gaussian entries of variance I/N and M >> N, we have that (V/o)
Cjr (z) (1/a) (z) which leads to simplification of equations 92
1i 20 C Dpg [z]+
G = 2
pi G (2) + (1 +
p3 Gps (z),
Ps =
(ps
+
x- (1
(111-)) + <
(98)
Therefore, Çi = G - (z).
In Figure 23 the MSE of the Oracle estimator and the RIE 96). 977) is illustrated for shifted Wigner X with c 3, and Wishart with aspect-ratio a' = 1/4.
Effect of aspect-ratio a. In Figure 24 the MSE is depicted for various values of the aspect-ratio a 1. We see that as M decreases (a increases) the estimation error (of Y) increases.
we take X to be a shifted Wigner matrix with = 3, and
43

---

0.1
A305
0.095
MSR
0.09
0.085
Orack estimator,
A (S)
1229
15
10
Figure 24: MSE estimating X function of aspect-ratio a 1, prior X is shifted Wigner with c
and x = 5. MSE normalized by the norm of the signal, || l Both Y and are Gaussian entries of variance 1/N. The RIE is applied to 2000, M 1/NN, and the results are averaged over 10 runs (error bars are invisible). Average relative eroor between RIF x (8) and Oracle estimator is also reported.
matrices with i.i.d.
F.2
Estimating Y
F2211
Relation between overlap and the resolvent
ON
y.
(1)
with 2 right/ left singular vectors of Y, we
For the vectors r;
1
y
0M
have
2M
M+N Σ
r. [Im Gs(x ((((xxxee) ic))4 = Σ
(r.Sk)
(1/5s) +
(r.sk) (1TsK)
(x 5x)² 2 + 2
x² 2² +
2M +1
(99)
Given the structure of s
87). we have:
1 (u (vTg!")
for 1 < k M
(rTss) (17
uT M (v
(r))
for M + 1 2M for 2M + + N + M
0
Therefore, in the limit N -> we have:
r. (Im Gs(x- (as(==ie)), ie)) N= 1 e
a / (c- + (r(t,ss dt
(100)
where the overlap function Oy Ax) is extended (continuously) arbitrary values within the support of As with the property that Oy X²) -Oy (t, A₂) for supp(μs) Sending e - 0, we find
iie) 1 ( r. (Im (amnGG(s-ii))
(101)
a
44

---

F.2.2 Resolvent relation
The resolvent relation for the model 60 with M < N is the same in 72 with parameters satisfying:
04 C (9i 92)
+ 2V
1i Rpx (9 + (9995)
Rxx (1) - /9993)
91
Clu) (0T
95 =
33 2 &
1 9i
Rpx (9t + 9993)) - Rpx (2x - 10--iiii (i73))
3 2 Rxx a + v /9/3))
Rpx (4i v 47
with
[Z](z) (z B1) - Z2(2z) = 31 + +gg(2 31)
J6 a Za(a)²
(e- 11-0 gpr (1(z))
(2(z)) a Z2(2)
+ a-1
a 3-97
Q 72(2) 9 Zi(t
-
72(=))
9 a (2)=
9pv
7z(=))
(1(z)
a (z(=))
1 8-81
94 Za(a)²
(2)
9pv
(1(=))
Z2(3))
a Z2(=)
8
(102)
With the same procedure form:
the saddle point equations (022) can be rewritten in simplified
(1)95)
81
"W
RRxx (91 + /9793) - Rpx (2i v 9i93)
oi
of
c(a) (gi%))
81
1
R (4ī + /9i53) - Rox ( - /
%
31 2 Rpx (4i + / (5)
+ RRx (9t v9793)
(103)
Pps (2) + {11 -1)
9 G (2)
4 (3-8 82
z-8
01 0
1 34 17 PT
Note that both sets of equations 101
103
and
59
76
match for a = 1.
F.2.3 Overlaps and optimal singular values
From
101.
we have:
3
(T) Gyry (z)
(r(7, 2e
Im
lim
y
YTy,
Tps(>) 2
Zz(z)
(z(2)
(104)
a
Im
lim 10
3
(n)(()
77(2)
Z2(2)0
Similar to 788. we can compute the optimal singular values to be:
S
Im
lim 4
(105)
Tps(7)
F.2.4 Numerical examples
We consider the matrix W to have i.idd Gaussian entries with variance I/N, so e(//2) (z) = 2. And, X F cf where F FT NxN has i.i.d. entries with variance I/N, and * 0 is real number, so R (z) = x + c. With these choices, the solution (03) simplifies to:
8i 2 + Ug
85 = qi
3 = 4i
3 = gi +
qî = g (2)
(1- 1) 2 = Gpe (z)
(106)
05 B1 1i
(3-81)²
2-8
o
2-9
1 3
31
45

---

1.91
Oraele estimator, E; (S)
vr Uracle estimator, (5-(8)
0.25
RIE. (S)
RIC. (5)
0.4
038
0.15
0.2
0.1
0.1
3
(a) Gaussian
(b) Uniform spectral, ([1,
Figure 25: Estimating Y. MSF is normalized by the norm of the signal, |Y X is a shifted Wigner matrix with = 3. and W has i..dd Gaussian entries of variance 1/N, and N/M = 2. The RIE is applied to N 2000, M 1000. and the results averaged over 10 runs (error bars invisible).
0.116
0.114
CURSS
MSE
0.112
0115
0.108
(S)(S)
0.106
RIE, (S)
Figure 26: MSE of estimating Y as function of aspect-ratio and 5. MSE normalized by the norm of the signal, || W has i.i.d. Gaussian entries of variance 1/N. The RIE is applied to N 2000, M 1/N,, and the results are averaged over 10 runs (error bars are invisible). Average relative error between RIE (S) and Oracle estimator is also reported.
Y has Gaussain entries of variance I/x, a shifted Wigner matrix with and
After a bit of algebra, we find that &i is the solution to the following qubic equation:
3r + 3 + 22 + - (z - Gps (z)) - Gas (2)
a a-1
(107)
- - - Gas (x)) (-)
a --1) - 11] =0
az
In figure 25 the MSE of RIE and the oracle estimator is plotted for two cases of priors: Y with Gaussian entries and Y with uniform spectral density.
Effeet of aspect-ratio a. In Figure 26 we take Y to have Gaussian entries (with variance // and the MSE is depicted for various values of the aspect-ratio a > 1. We see that as M decreases (a increases) the estimation error (of Y) increases.
G Details on numerical implementations
G.1 Numerical approximation of Gps (2)
The first step to construct the RIEs is to compute the Stieltjes transform of the observation matrix S. In section 19.5 of [58), several approaches have been proposed to approximate the Stieltjes transform of the spectral density of a given matrix. our implementations, we use the Cauchy kernel method in which for given matrix A with N singular values (or eigenvalues) (0:) ISiSN 9u, (2)
46

---

is approximated as:
N N
1
1
gua (z)
~
in
with 17; 's the "widths" of the kernel at each singular value (more precisely the imaginary part is sum of Lorentzians with width 77 around peaks at oj). The construction of the RIEs uses the Stieltjes transform of the limiting symmetrized measure of S. In the numerical experiments approximate this quantity as:
N
(2) ~
(108)
2N ==1
-
in
2 2+
7ss in)
with a fixed width 1j - v 1/2N. Note that for the case of a > (M < N). has M non-trivially zero singular values, and in the approximation above N should be replaced by M.
G.2 Construction of the RIEs
In the RIEs derived in 194, the final estimator for optimal singular values (eigenvalues) was rather simple and only required compute the Stieltjes transform the real line which can be easily and safely performed using the approximation above (see remark in section 19.5.2 in (58) ) However, in the RIEs of this work, need to solve system of equations in the limit
0 (2 close to the real
line). For this, to compute the optimal singular value Ege (or optimal eigenvalues ( we evaluate 9ps (z) for = = 7/ - i V2N" In this way, the other parameters (e.g. 4i) are evaluated for = very close to the real line, and the theoretical limit lim,-0 in (49)
78 can be estimated numerically, Moreover, as we considered fixed width in our numerical approximation of Sticltjes transform (T08). should be chosen to compensate the width for the cases where the support of as is wider. For example, for fixed N, as we increase SNR (from to 5) the support of pa grows, however we still have N singular values and the kernel's ' width in (108) is fixed, so E should be larger for higher SNRs to get a more accurate approximation of Gps(2)
G.3 Mismatch between RIEs and Oracle estimators
The RIEs are conjectured to have the same performance as the Oracle estimators in the limit N + co. Therefore, we believe that the mismatch between the proposed RIEs and the Oracle estimators is finite size effect. Moreover, this finiteness affects the accuracy of estimated parameters, since Gas (=) approximated numerically and we do not use random matrix theory to find its exact form.
Generically, the mismatch between the RIE and Oracle estimator is larger for the case of estimating X. We expect that this because of the extra approximation step in the derivation of the optimal cigenvalues. In the fifth line of (49), the sums are approximated by an integral which is the Stieltjes transform of px. This approximation does not appear in derivation of the optimal singular values for Y, see 78)
All in all, the small relative error (less than 1%) between RIEs and Oracle estimators in our numerical results validates our optimality conjecture and demonstrates that RIES can be successfully used in practice.
H Spherical integrals and matrix lemmas
H.1 Spherical Integrals
For two symmetric matrices A, B E R NxN the spherical integral is defined as:
N Tr AU
Iv(A,
exp 1
where the average is w.r.t. the Haar measure over group of (real) orthogonal N N matrices. The spherical integrals can also be defined w.r.t. the unitary symplectic group. These integrals are often referred to as Harish Chandra- Shandra-Ityykoon-P -uber (HCIZ) integrals in mathematical physics literature. The study of these objects dates back to the work of mathematician Harish Chandra 62
47

---

and they have since been extensively studied and developed in both physics and mathematics. In particular, (T] studied the limit of the integral in the case where one of the matrices, say A, has finite rank.
Theorem (Rank-one spherical integral, Guionnet and Maïda (2I)) Let 0 be the only non-zero eigenvalue of A (so it is rank one), and the empirical eigenvalue distribution of B converge weakly towards Pg. Then, for sufficiently small (see details in Theorem in (27)), have:
In (A. 2
S. Rp (1) dt = 2
lim N-000
(109)
When A has higher (but finite) rank, theorem 7 (2I) states that the limit is the sum over eigenvalues of the expression the rhs of (109.
Non-symmetric case. In the non-symmetric case the rectangular spherical integral is defined, for the matrices A E MM B E RNX. as:
JN(A, ( exp {VNM AUBV))
where U E -NXNN VE E RMM and the expectation is w.r.t. the Haar measure over orthogonal matrices of size N x N and M x M.
Theorem 2 (Rank-one rectangular spherical integral, Benaych-Georges 22 Let N/M -
(0. 1), and 6 be the only non-zero singular value of A. and the empirical singular value distribution of B converges weakly towards Un. Then, for 0 sufficiently small (see details in Theorem 2.2 in (22)
we have:
In JN (A. = 1 9 clo) (12)
0 cla) (1)
lim N N
di
1
di=
o(0) 2 Sua
(0
(110)
2
In our derivation, we use a generalization of this formula, namely when A has higher (but fixed) rank, the limit is the sum over singular values the expression on the rhs of (10). Although we are not aware if this generalization has been proved, we believe that the ideas found in 63 can be applied to show it holds.
Remark 6. It is known that additional terms may be present on the rhs of 109 and (10) when the parameter 0 is "large". This has been rigorously proved at least in the case of symmetric A and B (see theorem 6 in (21)). ). In the replica calculation the order of magnitude of this parameter is determined by the solutions of the saddle point equations, but it is difficult fully control its order of magnitude. However the numerics show very good agreement between our explicit RIEs and the
Oracle estimator, which strongly suggests sound to use (109 and 110p.
H.2 Matrix analysis tools
A C
Proposition (Inverse of block matrix, Bernstein [61].. For a block matrix F
with
D
A E R then,
N
B E NRN C MRN D E R M XM if A and D CA-1B, are non-singular,
A A-1 + CCA - - - CC --1
ID-CA--B)-CAA
(D-CA-
Block structure of Gs(z) The matrix Gs(z) is:
=IN -S ST
Gs(z) (II s)
ZIM
Using Proposition first we need to compute the inverse matrix (2Ins (2T))2))-(((SS) which simply reads:
(2IM 1
= ==(=IIM - STS) - =GTTs(#3)
Consequently, we find:
1 I + + SG
SGSts(22) =Gsts(22)
Gs(z)= =
(111)
Gs
48

---

Inverse of C For C3 since the blocks B, C' are zero, the inverse is simply:
Cx -1
[[2
IN cx
0
[( C [AS]
Si IN
x
(112)
=I< IM
Gxa -
IM
2-
Inverse of Cy
Let the block structure of Cy be as in Proposition
then
DD-CA -B)CAA'BB-__
(2 (2- { B
32
YTT)
-1
2 2- $9
(= B IM (3z +
31
YTV)
(1(2) IM - Zz(x) -1
= (= B1)
31 (z)
-1
IM - -YTT)
= Z2(z) Z2(z)
31
GYTy (1)
Z2(z)
Z2(z)
where Gyry the resolvent of the matrix YTY. So, we have
(2 B1)
IN
81
Y YGYYY Z1(2)
YT
91
Y YGrry
(1(=)
Cy
2)))
[2(=)
2z(=) 3 [x(=)
Zg(z)/
3: Z+(=)
GYTy
(1)) Yr
Gyry
(2)2)
Z,(=)
Zz(3)
Lemma 3. Consider two vectors x. y E RN The symmetric matrix xy + you has rank at most
two with non-zero eigenvalues x + ylll-
Proof. Construct the matrices A 2
E
Nx2
x2 as follows:
T
A
= x
y
Then, we have that zy + y = -A. Using the lemma we have that:
2² det (IIN - BA) = N det [zI2 AB)
So, the characteristic polynomial of ay + yæu N-2 det (=I, AB), which implies that the
Ty + yæu has eigenvalue O with multiplicity N 2, plus the eigenvalues of the 2 x matrix AB. The matrix AB is:
x
| xTy
AB
which has eigenvalues Iy + || |1yll-
Lemma 4. For matrices A E I 2.xx X B E NxM
we have:
M det (=IN - BA) N det (=IM - AB)
Proof. Construct the matrices C, D. E (MM+N) follows:
= 2I A
IM
OMXN ZIN
D=
B
IN
- B
We have:
det CD - N del (=IM
AB),
del DC = M del (=IN BA)
The result follows from the fact that det CD = det DC.
49


Reviews:
[{"review_id": "snjlYoJE6e", "reviewer": "Reviewer_LgE7", "summary": "For a matrix factorization model S = \\kappa XY + W, this paper proposes a method for estimating X and Y from S under the assumption that priors of X, Y, and W satisfy certain rotation invariance properties and their distributions of eigen/singular values are known. The proposed method is rather simple. First, using the singular value decomposition, we assess the left and right singular bases of S, which are eigen/singular bases of X and Y. Next, keeping the bases fixed, the eigen/singular values of X and Y are adjusted for minimizing the average mean squared errors, which can be analytically performed with the knowledge of limiting distributions of eigen/singular values of X, Y, and W. The Bayesian optimality is shown for the proposed estimators using the replica method from statistical mechanics and random matrix theory. ", "questions": "I am curious about what happens when c is set to zero for the shifted Wigner. Does it cause any singularity for the estimator? Or, does the estimator of X continuously converge to zero matrix as c -> 0?", "limitations": "Yes. ", "rating": 9, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "As long as I know, this is the first paper that shows a concrete practical method for constructing the Bayes optimal estimator for O(N) rank matrix factorization problem. ", "weaknesses": "The shown optimality holds under rather many assumptions (rotation invariance, knowledge of limiting eigen/singular value distributions) are necessary. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for positive comments. We address his/her question below:\n\n$\\bullet$ For $c=0$, the prior on spectral of $\\mathbf{X}$ is symmetric, $\\rho_X(x)=\\rho_X(-x)$, and using the definition of Stieltjes transform one can see that indeed from eq. (7) the estimator is 0 for all eigenvalues. We conjecture that the estimator is continuous at $c=0$.\nAnalyzing rigorously the continuity of the estimator as a function of $c$ is non-trivial, as one needs to consider the effect of $c$ on the limiting spectral measure of the observation $\\mu_S$ which enters into the parameters $\\zeta_1, \\zeta_3$. However, ignoring this technicality, the function $G(z) + G(-z)$ in the estimator converges continuously to 0 as c goes to 0. Moreover, our numerical checks support the continuity of the estimator. "}, {"review_id": "kXO7OHB7Fh", "reviewer": "Reviewer_LQDn", "summary": "This paper considers a matrix factorization problem in a setting where the rank of the factor matrices grow linearly with the ambient dimensions. They assume that the factors follow a prior distribution such that: (1) One of the matrix factor is symmetric, (2) Both factors and the noise are drawn from rotationally invariant distributions. (3) The priors are known to the statistician. \n\nThey propose to study a class of rotationally invariant estimators. They derive a closed form expression for the oracle estimator in this class, and show that it is Bayes optimal. They propose an estimator and conjecture that its performance matches the oracle. They provide evidence of their conjecture through experiments. ", "questions": "1. I am wondering how sensitive is the proposed estimators to misspecification of prior. \n2. If the prior is not presented, is there a way to estimate them? I feel the assumption that prior is known for both the factors and noise is a bit strong. Perhaps the authors should comment a little bit. For example, explain when it is reasonable to assume this information is given. \n3. Several literatures that might be useful to include: \n\nInformation-theoretic limits of MF: \n[1] Bayes-optimal limits in structured PCA, and how to reach them (arXiv:2210.01237)\n[2] Fundamental Limits of Low-Rank Matrix Estimation with Diverging Aspect Ratios (arXiv:2211.00488)\n\nAMP for rotationally invariant matrices: \n[1] Approximate Message Passing algorithms for rotationally invariant matrices (arXiv:2008.11892)\n", "limitations": "The limitations are clearly reflected in the model assumption, and societal impact not applicable. ", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "This paper seems to be the first one that explores MF in this challenging setting, and opens up a new research direction that might be of interest. They derive a closed form neat expression for the Bayes optimal estimator, although this cannot be implemented in practice. As an alternative, they propose an estimator that can be implemented in practice. Although their result is non-rigorous, simulation suggests that this is the correct thing to do. A sub-optimal estimator for one factor that does not require prior knowledge of $\\mu_X$ is also proposed. Their presentation is nice and clean.", "weaknesses": "Although this paper presents nice technical contributions, it is not clear why this prior structure should be considered in practice. It would be nice to give a few practical examples which past results can not cover but this work do. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the referee for the review and additional suggested references. We will add those (with possibly other references) in the final version.\n\n We agree that the required assumptions are restrictive and do not necessarily hold in practice, however as mentioned in the manuscript the proposed estimators can be used without the assumptions to get a sub-optimal estimate that can be processed further. Moreover, as  pointed out in the comments, this work is the first one to consider the problem in the high rank regime, and we believe it may open up a new way to study the problem in this regime. \n\nWe address his/her questions below:\n\n1. Analyzing the sensitivity of the estimators to the mismatched priors is an interesting object of study that requires independent investigation, and is beyond the scope of this manuscript. There exists a rich variety of scenarios that one can consider. Here, we give a brief speculation on the performance of RIE with mismatched priors. \n\n    $\\bullet$ Mismatched prior on the signal. The RIE for estimating $\\mathbf{Y}$ is independent of the prior of $\\mathbf{Y}$, and it is optimal in the case of mismatched prior on $\\mathbf{Y}$. \n\n    For the case of estimating $\\mathbf{X}$ with mismatched prior, the estimator $\\sqrt{\\widehat{\\mathbf{\\Xi}_{X^2}} (\\mathbf{S})}$ is applicable and is independent of the prior, although sub-optimal. However, the RIE for $\\mathbf{X}$ depends on the prior, and misspecification of the prior leads to the poor performance of RIE (comparing to the Oracle). But, we believe that the RIE can get a non-trivial estimate of $\\mathbf{X}$ (which is indeed sub-optimal), as the parameters $\\zeta_1, \\zeta_3$ in eq. (8,9) are correctly evaluated. We provide a numerical check of the performance of RIE in this case in the file uploaded as official comment above. \n\n    $\\bullet$ Mismatched prior on other factors. Misspecifying the priors for the other two factors (than the signal), leads to the incorrect evaluation of the parameters for RIE that can change the performance significantly. For example, consider estimating the matrix $\\mathbf{X}$, with $\\mathbf{Y}$ and $\\mathbf{W}$ both Gaussian matrices but with variances unknown to the statistician. Assuming mismatched variances (different than the real ones) for  $\\mathbf{Y}$ and  $\\mathbf{W}$, leads to rectangular R-transforms which are equal to the true ones times some constant. This results in a non-trivial mismatch in the parameters $\\zeta_1, \\zeta_3$ (eq. 8,9) used in the estimator, which will change the estimated optimal eigenvalue. We provide a numerical check of the performance of RIE in this case in the file uploaded as official comment above. \n\n\n2. We agree that the assumptions are strong, however note that these assumptions are required to show the optimality of the estimators.\nIn general, for $N$ large enough estimating the prior on one factor is possible if the prior of the other two factors are known. For example, knowing the prior of $\\mathbf{X}$ and the noise, we can estimate the spectral distribution of $\\mathbf{Y}$, using the spectral distribution of the observation. For this, one needs to go through the free additive/multiplicative convolutions: From spectral measure of the observation and the noise, one can find the rectangular R-transform of spectral distribution of the product $\\mathbf{XY}$, then again using the knowledge of prior of $\\mathbf{X}$ one can find an estimate of the spectral distribution of $\\mathbf{Y}$.\n\n     Moreover, even if we consider specific classes of priors for the factors we do not necessarily find unique estimates. For example, in the simpler setting of additive denoising $\\mathbf{X} + \\mathbf{Z}$, under the assumption that both $\\mathbf{X}$, $\\mathbf{Z}$ have i.i.d. Gaussian entries with unknown variances, it is not possible to estimate the variances uniquely for both $\\mathbf{X}$,$\\mathbf{Z}$.\n\n     Rotational invariance for the priors might not be very natural in practice, but our estimators could be used as an initialization for other algorithms (for example iterative ones - although for high rank these remain to be investigated). Additionally, assuming other priors than Gaussian does not lead to neat expressions for the estimators and the systems of equations must be solved numerically, which is impractical. On the other hand considering Gaussian priors, especially for the noise, is a common assumption in practice.\n \n"}, {"review_id": "arIYFB9ehK", "reviewer": "Reviewer_AAM7", "summary": "\nThis paper explores the Matrix Factorization problem, which involves estimating the matrices $X \\in \\mathbb{R}^{N \\times N}$ and $Y \\in \\mathbb{R}^{N \\times M}$ given the noisy matrix $S = \\sqrt{\\kappa} X Y + W$. \n\nThe focus is on the high-dimensional regime, where $N/M \\to \\alpha$, and the investigation includes bi-rotationally invariant $Y$ and $W$, as well as symmetric rotationally invariant $X$.\n\nThe authors examine rotationally invariant estimators, which are estimators that share the same singular vectors as the noisy matrix $S$.\nThe paper derives rotationally invariant estimators based on oracle knowledge of the target matrices and demonstrates that they are also Bayes optimal. By assuming concentration and utilizing replica methods, the paper derives explicitly computable estimators from the oracle estimators. The empirical performance of these derived estimators is investigated and shown to closely match that of the oracle estimators. This suggests that the estimators derived using non-rigorous methods from statistical physics are indeed optimal.", "questions": "Can the authors expand on the relevance of the methodology developed in the paper for the analysis of the weight matrices of neural networks?", "limitations": "Theoretical paper with no immediate negative societal impact.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The low-rank matrix factorization problem with finite-rank matrices is now a well-studied topic. Similarly, the low-rank matrix denoising problem with extensive (diverging) ranks has garnered recent interest. However, results on matrix factorization with extensive ranks have been relatively scarce. This paper aims to fill this gap in the literature by providing results for the matrix factorization problem with diverging ranks and under general rotationally-invariant priors.\n\nThe paper is well-written overall, and Section 5 provides a concise and easy-to-follow overview of the derivation of the results, which are otherwise quite complex.", "weaknesses": "The main results of the paper, which are the explicitly computed Rotationally Invariant Estimators, rely on non-rigorous methods from statistical physics. While the empirical results are compelling, it would be valuable in the future to establish a more solid theoretical foundation for these estimators.\n\nMoreover, the assumptions made on the matrices $X$ and $Y$ may be considered somewhat unnatural. It would be beneficial for the authors to provide additional motivation as to why the findings of this paper could be of interest to the NeurIPS community beyond the specific problem examined here. This could help clarify the broader significance and potential applications of the research.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the comments.\n\nWe agree that making the derivations mathematically rigorous is an interesting research direction. We want to bring your attention to the papers:\n \"Optimal cleaning for singular values of cross-covariance matrices\" (arXiv:1901.05543) , \"A short proof of Ledoit-Péché's RIE formula for covariance matrices ( arXiv:2201.05690),  that rigorously derive the optimal RIE for covariance estimation with Gaussian priors. The method developed in these works involves Gaussian integration by parts and Gaussian concentration. We believe that, using this technique we can establish the optimality of the proposed estimator at least for Gaussian priors, and this analysis would not involve replica method and spherical integrals. However such a rigorous mathematical analysis is at the moment beyond the scope of this paper.\n\nWe agree that the assumptions are strong and far from practice, however the settings considered in the manuscript is one of the first setting that the matrix factorization can be solved optimally in the high-rank regimes, and we believe that this work may open up a way to study this problem in more general settings. Moreover, some of the assumptions are required to show the optimality of the estimator, but the estimator can be used in practice under milder conditions to get a spectral estimate that can be refined.\n\nConcerning the specific question on neural networks:\n\nWe believe that our method can be applied to study the weight matrices in neural networks, at least when there is no non-linearity in the system. Non-linearities are prevalent in neural networks models which does not bode very well with rotation invariance, but it might be possible to consider \"linearizations\" such as an NTK approximation to circumvent this issue. These is an interesting open problem. \n\n"}, {"review_id": "aOR82kd9Oh", "reviewer": "Reviewer_PN68", "summary": "The authors consider the problem of matrix factorization of a noisy measurement in the setting of all matrices having an rotationally invariant prior. They provide a non-rigorous but comprehensive theoretical derivation of their results. They also provide a number of experiments validating there theoretical claims.", "questions": "- It would be nice to have a short summary and outlook at the end of the paper\n- In line 67 what do you mean by proper distribution?\n- Can you briefly elaborate on the relationship between your work and the works 34-36 mentioned in the introduction?", "limitations": "The limitations have been adequately addressed.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "- Explicit formulas for the reconstruction of the matrix factors\n- Strong experimental support for the theoretical claims", "weaknesses": "- The analysis is limited to the rotationally invariant setting", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the review. We address the questions below:\n\n1. We will surely include a conclusion section in the final version, also taking into account comments made by the referees.\n\n2. Depending on the nature of the problem, one can consider a general class of priors for the factors, and tries to find the parameters maximizing the posterior. For example, in ref [37] authors consider Gaussian priors and try to estimate the matrices. We will rephrase the wording \"proper distribution\" in the final version.\n\n3. Ref[34] is a general reference about variational inference. In [35,36], matrix factorization problem is solved when some of the data points (entries of the observation) are missing. In our setting, we have the full observation of data, therefore the problem setting are different.  Solving MF with full observed data based on variational inference is studied in ref[37]. In ref[37], authors showed that under Gaussian prior (which is rot. Inv.) the optimal estimate using variational inference approach is a reweighting of the SVD of the observation, which can be seen as a RIE.  "}]

Author rebuttal:
Numerical results on sensitivity of RIE to mismatched priors ( response to first question of Reviewer LQDn)

In figure 1, the spectral distribution of $\mathbf{X}$ is uniform on $[0,4]$, and both $\mathbf{Y}, \mathbf{W}$ are Gaussian matrices. We applied the RIE, assuming two different misspecified priors for $\mathbf{X}$:Shifted Wigner with $c=2$, Wishart with aspect-ratio $1$. Note that, the estimator $\sqrt{\widehat{{\mathbf{\Xi}_{X^2}^*}}(\mathbf{S})}$ does not require the knowledge of prior of $\mathbf{X}$, and we get the same performance for both cases, although sub-optimal. The RIE $\widehat{{\mathbf{\Xi}_X^*}}(\mathbf{S})$ performs worse than the Oracle estimator, but still it provides us with a non-trivial estimate. 

In figures 2, we consider estimating $\mathbf{X}$ with Wishart prior. $\mathbf{Y}$ has uniform spectral distribution on $[1,3]$, and the noise matrix $\mathbf{W}$ has Gaussian entries. The RIE is applied with the assumption that $\mathbf{Y}$ is Gaussian so that its spectral distribution is the square-root of Marchenko-Pastur law (whose support is $[1,3]$). As discussed in the comment, we see that the RIE performs poorly, comparing to the Oracle estimator. However, note that the normalized MSE (normalized by the norm of the signal) is below $1$, and we get a non-trivial estimate of the signal matrix $\mathbf{X}$.

Meta-review: The reviewers have good expertise and are extremely positive about the result. The manuscript considers an unusual matrix factorisation setting, one which seems to the AC to be better suited for a standard matrix analysis venue due in large part to limited attempt at motivating this question as being machine learning other than the second sentence of the introduction.  The authors and reviewers have had a useful discussion about the issues raised.
Decision: Accept (spotlight)