AppAppendix References

# Statistical properties of large data sets with linear latent features

Philipp Fleig  Department of Physics & Astronomy, University of Pennsylvania,
Philadelphia, PA 19104, USA  Ilya Nemenman  Department of Physics, Emory
University, Atlanta, GA 30322, USA  Department of Biology, Emory University,
Atlanta, GA 30322, USA  Initiative in Theory and Modeling of Living Systems,
Atlanta, GA 30322, USA

###### Abstract

Analytical understanding of how low-dimensional latent features reveal
themselves in large-dimensional data is still lacking. We study this by
defining a linear latent feature model with additive noise constructed from
probabilistic matrices, and analytically and numerically computing the
statistical distributions of pairwise correlations and eigenvalues of the
correlation matrix. This allows us to resolve the latent feature structure
across a wide range of data regimes set by the number of recorded variables,
observations, latent features and the signal-to-noise ratio. We find a
characteristic imprint of latent features in the distribution of correlations
and eigenvalues and provide an analytic estimate for the boundary between
signal and noise even in the absence of a clear spectral gap.

††preprint: APS/123-QED

Introduction. Massively parallel experiments are now standard in all fields of
science. They record the state of the system through a large number, $N$, of
variables $x_{i},\,i=1\dots N$. These variables could be positions of
particles, agents, or tracers in dusty plasmas [1], soft matter [2], insect
swarming [3], and dynamical systems [4]. Or they can be field values at
different spatial points in fluids [5], climate data [6], or activity of
“nodes” in gene expression networks [7], neural recordings [8], postures [9,
10], biodiversity [11], ecology [12], etc. Crucially, the number of recorded
variables is often much larger than the true number of (latent) degrees of
freedom in the system [13, 14, 15]. This allows to use correlations among the
measured variables to detect the latent state variables, which can then be
included in models of the system [16, 13, 17, 18]. Such approaches generally
go under the name of dimensionality reduction. While many dimensionality
reduction methods have been developed, the classic Principal Components
Analysis (PCA) still reigns supreme in its widespread use. PCA assumes the
recorded data to be Gaussian, and the latent features to be linear
combinations of the recorded variable activities. It then calculates the
variable-variable correlation matrix and identifies the eigendirections
(principal components) corresponding to the large eigenvalues of the
correlation matrix as the latent features. If there is a gap between $m\ll N$
largest eigenvalues and the rest, then $m$ principal components capture a
large part of the variance of all variables, achieving the dimensionality
reduction.

While the number of variables may be (much) larger than the number of latent
components, modern experiments are usually undersampled since nonstationarity
precludes making measurements for a long time. Thus the number of independent
measurements, $T$, is often of the same order as $N$. Because of this,
statistical fluctuations in the correlation matrix are large (in fact, the
matrix is degenerate for $T\leq N$), and the principal components calculated
from samples often only vaguely resemble their true values, with many
correlations emerging as statistical artifacts. To identify if findings of the
PCA analysis can be trusted, one refers to Random Matrix Theory [19], and in
particular the Marčenko-Pastur (MP) eigenvalue density of a pure noise
correlation matrix [20]. Specifically, one calculates the upper and the lower
bounds of eigenvalues expected by pure chance from $T$ measurements of $N$
independent variables, and eigenvalues outside this interval (and their
corresponding eigenvectors) are deemed to be statistically significant.

This straightforward approach assumes that signal-induced correlations among
the variables do not influence the spectrum of the noise-induced correlations.
This has never been proven and, as we will show, is, in fact, incorrect. More
generally, we are not aware of results extending the MP analysis to produce
correlation eigenvalue densities when the correlations come from the sampling
noise and from true, low-dimensional, latent signals, with the latter coming
from some known distribution in its own turn (though see work on spiked
covariance matrix models [21, 22, 23]). Even statistics of the entries of the
correlation matrix (rather than of its eigenvalues) have not been reported in
this case. In this Letter, we fill in these gaps and calculate – analytically
and numerically – various statistical properties of correlation matrices
emerging from data sets with low-dimensional latent feature structure. We show
that the distribution of pairwise correlations and the spectra of eigenvalues
of these correlation matrices carry signatures of the number of latent
features, allowing one not only to choose, rigorously, which of the principal
components are above the noise floor, but also to see if the overall model of
latent features plus noise is a good description for a particular data set.

We do our analysis in two limits. First is the classical statistics limit,
where the number of variables is constant (though possibly large) and the
number of latent features is small, while the number of observations grows to
infinity, such that $N/T\rightarrow 0$. Second is the intensive limit, where
both the number of the measured variables and the number of observations
grows, such that $N/T=\text{const}$, while the number of latent features stays
finite $m={\rm const}$. Both of these limits routinely happen in modern
datasets. We leave considerations of the extensive limit, where the number of
latent features grows with the number of the measured variables, to a future
publication. We believe that these results are an important step in the
development of analytical tools for understanding large datasets.

The model and its limits. We consider observations produced by a random matrix
model that combines latent signals and uncorrelated noise:

$\displaystyle\mathbf{X}=\mathbf{U}\mathbf{V}+\sigma\mathbf{R}\,.$ (1)

The component matrices of the latent signal term $\mathbf{U}$ and $\mathbf{V}$
have dimensions $T\times m$ and $m\times N$, respectively. Thus $m$ latent
features get randomly sampled $T$ times (matrix $\mathbf{U}$), and each of the
$N$ measured variables is a random linear combination of the latent features
(matrix $\mathbf{V}$). We assume $m\leq T,N$ throughout this work, such that
the rank of the signal matrix $\mathbf{U}\mathbf{V}$ is equal to the number of
latent features $m$. In other words, an estimate of the features can be
inferred from the samples uniquely.

The entries of $\mathbf{U}$ and $\mathbf{V}$ are Gaussian random variables
with zero mean and variances $\sigma_{U}^{2}$ and $\sigma_{V}^{2}$,
respectively:

$\displaystyle U_{t\mu}\sim\mathcal{N}(0,\sigma_{U}^{2})\,,\quad V_{\mu
n}\sim\mathcal{N}(0,\sigma_{V}^{2})\,,$ (2) $\displaystyle
t=1,\ldots,T,\;\mu=1,\ldots,m,\;n=1,\dots,N.$ (3)

We make this choice for analytic tractability; in applications to real data,
the variances of the entries may need to be matched to variances of each
measured variable. Finally, the elements of the noise matrix $\mathbf{R}$ are
i.i.d. unit variance Gaussian random variables, so that the noise in every
observation has variance $\sigma^{2}$.

From Eqs. (1, 2), the elements of the signal matrix $\mathbf{UV}$ are a sum of
$m$ products of two Gaussian variables with variances $\sigma_{U}^{2}$ and
$\sigma_{V}^{2}$. In Online Supplementary Materials A we use characteristic
functions to derive the probability density of these entries and show that
their variance is

$\displaystyle\sigma_{UV}^{2}=m\sigma_{U}^{2}\sigma_{V}^{2}\,.$ (4)

In other words, as expected from addition of independent random variables,
each of the latent components adds $\sigma^{2}_{U}\sigma^{2}_{V}$ to the
variance of the observations.

Further, for a large number of latent features $m$, the probability density of
${\mathbf{X}}$ approaches a Gaussian with zero mean, see Fig. S2. This allows
us to define a Gaussian signal-to-noise ratio,
$\textsf{SNR}\equiv\sigma_{UV}^{2}/\sigma^{2}=m\sigma_{U}^{2}\sigma_{V}^{2}/\sigma^{2}$.
Since our goal is to calculate properties of the data matrix independent of
the units chosen to measure each of the variables, we normalize the data
matrix

$\widetilde{\mathbf{X}}\equiv\mathbf{X}/\sigma_{X},\;\sigma_{X}^{2}\equiv\sigma_{UV}^{2}+\sigma^{2}=\sigma^{2}(1+\textsf{SNR}).$
(5)

Note that this normalization is by an expected standard deviation, and is
different from subtracting empirical means and standardizing by an empirical
standard deviation. However, we expect the difference to be of order
$T^{-1/2}$, and thus negligible in what follows.

In this work, we will be exploring the properties of the normalized empirical
covariance matrix (NECM)

${\mathbf{C}}=\frac{1}{T}\widetilde{\mathbf{X}}^{T}\widetilde{\mathbf{X}}=\frac{1}{T}(\widetilde{\mathbf{U}\mathbf{V}})^{T}(\widetilde{\mathbf{U}\mathbf{V}})+\tilde{\sigma}^{2}\mathbf{R}^{T}\mathbf{R}\\\
\quad+\tilde{\sigma}(\widetilde{\mathbf{U}\mathbf{V}})^{T}\mathbf{R}+\tilde{\sigma}\mathbf{R}^{T}\widetilde{\mathbf{U}\mathbf{V}}\,,$
(6)

as well as the matrix of correlation coefficients

$c_{pq}=\frac{C_{pq}}{\sqrt{C_{pp}}\sqrt{C_{qq}}}\,.$ (7)

To explore different regimes of the problem, we use the following parameters:

$\displaystyle q\equiv N/T,\;q_{U}\equiv m/T,\;\mbox{and}\;q_{V^{T}}\equiv
m/N,$ (8)

only two of which are independent. These parameters emerge naturally in our
theoretical analysis in Online Supplementary Materials C.1 as controlling
behavior of various observables in the model. Then the classical statistics
and the intensive limits, introduced above, become:

Classical stats.: $\displaystyle q\rightarrow 0,\;q_{U}\rightarrow
0,\;q_{V^{T}}=\text{const}\,,$ (9) Intensive: $\displaystyle
q=\text{const},\;q_{U}\rightarrow 0,\;q_{V^{T}}\rightarrow 0\,,$ (10)

together with $\textsf{SNR}=\text{const}$ in both limits. Notably,
$q_{U}^{-1}$, gives the number of observations available per latent feature to
be learned. Since the parameter $q_{U}$ is small in both limits, it means that
the latent features are well sampled in our model, even if the measured
variables (controlled by $q$) may not be.

_Density of pairwise correlations._

Figure 1: Distribution of pairwise correlations in the pure signal limit
$\textsf{SNR}\rightarrow\infty$, for $m=2,4$ and $30$ latent features.
Analytic form given by a symmetric Beta distribution (magenta) and simulated
data (gray). Each simulations is run with $N=80$ variables and $T=2048$
observations and constitutes $1000$ independent model realisations.

The first observable statistics we calculate is the probability density of
correlations $c_{pq}$ of the standardized variables $\widetilde{\mathbf{X}}$
in our latent features model. Our goal is to analyze the dependence of the
density of the matrix entries $c_{pq}$ on $m$, $T$, and the noise strength.
The numerator in the correlation matrix in Eq. (7) has three contributions:
$(\mathbf{UV})^{T}(\mathbf{UV})$ from the pure latent features signal,
$\mathbf{R}^{T}\mathbf{R}$ from the pure noise, and two cross terms between
the signal and the noise, e. g. $(\mathbf{UV})^{T}\mathbf{R}$. Each of these
terms is analyzed separately in Online Supplementary Materials B, and we
reduce each term to correlations between independent Gaussian vectors. Such
correlations are distributed according to the symmetric Beta distribution [24]

$\displaystyle\mathrm{pdf}(r)=\text{Beta}(r;\alpha,\alpha;\ell=-1;s=2)\,,$
(11)

where the location $\ell$ and scale $s$ of the Beta distribution are set such
that correlations fall on the interval $[-1,1]$. The shape parameter $\alpha$
is determined individually for each contribution. For the signal-signal,
noise-noise, and signal-noise contribution we find

$\displaystyle\alpha_{\text{s}}=\frac{m-1}{2}\,,\;\alpha_{\text{n}}=\frac{T-1}{2}\,,\;\alpha_{\text{sn}}=\frac{m^{1/2}T^{1/2}-1}{2}\,,$
(12)

respectively. The sum of Beta distributions can be well approximated by a
single Beta distribution [25]. This finally allows us to approximate the
distribution of the entries of the full correlation matrix by a single Beta
distribution of the form Eq. (11). In Online Supplementary Materials B, we
show that in the limit when contributions of $\mathcal{O}(T^{-1/2})$ and
$\mathcal{O}(m^{-1/2})$ can be neglected, the parameter of the approximating
distribution is

$\displaystyle\alpha\approx\frac{\left(\sqrt{m(1+1/\textsf{SNR})}^{-1}+\sqrt{T(1+\textsf{SNR})}^{-1}\right)^{-2}-1}{2}\,.$
(13)

Notably, $\alpha$ is a function of $m$, and hence the shape of the
distribution depends on $m$. Thus the number of latent dimensions in data can
be estimated from the empirical distribution of the correlation coefficients.

Numerical validation of this result in the pure signal limit
$\textsf{SNR}\rightarrow\infty$, is shown in Fig. 1. For small $m$, the
distribution distinctly changes shape as $m$ varies. However, when the number
of latent features becomes comparable to the number of variables, and thus
$q_{V^{T}}\rightarrow 1$, the distribution approaches a Gaussian, making it
increasingly difficult to infer the precise value of $m$ from its shape. The
quality of the analytic approximation increases the smaller $q_{U}$, i. e.,
when more observation per latent feature are available. The analytic
approximation is also exact in the large noise limit $\mathsf{SNR}\rightarrow
0$. However, deviations appear for finite $\mathsf{SNR}$ when $m$ is small,
see Fig. S3.

_Eigenvalue density._ We compute the eigenvalue density of NECM, $\mathbf{C}$,
cf. Eq. (6), from its Stieljtes transform,
$g_{N}(z)=N^{-1}\text{Tr}(z\mathbf{I}-\mathbf{C})^{-1}$, where $z$ is a
complex number and denote the large-$N$ limit of $g^{N}_{\mathbf{C}}$ by
$\mathfrak{g}_{\mathbf{C}}$ [19]. The density of eigenvalues is then obtained
from the Sokhotski–Plemelj formula

$\displaystyle\rho(\lambda)=\frac{1}{\pi}\lim_{\eta\rightarrow
0^{+}}\Im\mathfrak{g}_{\mathbf{C}}(z=\lambda-i\eta)\,,$ (14)

where $\Im$ denotes the imaginary part.

Figure 2: Eigenvalue density and bounds as a function of the SNR in the
classical statistics and intensive limit. (a) Classical statistics limit
($q=0.01$, $q_{V^{T}}=0.9$) at three levels of SNR. Top plot: zero noise
analytic density Eq. (20) (yellow); large noise limit given by the Marčenko-
Pastur distribution (black); and intermediate noise semi-analytic density
(blue) with the numerical simulation (gray) for comparison. Vertical dashed
blue line is the approximate boundary between noise (left) and signal (right).
Bottom: eigenvalue bounds $\lambda_{\pm}^{\textsf{SNR}}$, as a function of the
SNR. True bounds obtained as numerical solution of Eq. (15) (green) and
approximate bounds given by the analytic expression in Eq. (22) (magenta). The
approximate noise region is striped. Horizontal dotted line indicates the
$\mathsf{SNR}$ value of the blue density curve and the dashed line extending
across plots indicates the signal-noise boundary. (b) Intensive limit
($q=0.01$, $q_{V^{T}}=0.09$) with plots analogous to (a). Top: zero noise
density Eq. (23) (yellow); for intermediate noise (blue), the left bump
corresponds to noise and the right bump to latent feature signal. Bottom: for
$\textsf{SNR}\gtrapprox\textsf{SNR}_{\text{split}}$ (horizontal solid black
line), the density splits into two bumps. Simulations constitute $360$
independent realisations of the model with $N=300$ and
$\sigma_{U}^{2}\sigma_{V}^{2}=1$.

Full details of the computations are given in Supplementary Online Materials
C. Briefly, we consider the eigenvalue density in the classical and intensive
limits, Eqs. (9, 10). This allows us to simplify the calculations by
neglecting the cross terms between the signal $\mathbf{UV}$ and the noise
$\mathbf{R}$, cf. Supplementary Online Materials C.3. Then, in the classical
statistics limit, the Stieltjes transform satisfies the third order polynomial
equation

$\displaystyle
a\mathfrak{g}_{\mathbf{C}}^{3}+b\mathfrak{g}_{\mathbf{C}}^{2}+c\mathfrak{g}_{\mathbf{C}}+d=0\,,$
(15)

with coefficients

$\displaystyle a$ $\displaystyle=\frac{qz}{1+\textsf{SNR}}\,,$ (16)
$\displaystyle b$
$\displaystyle=-\frac{qq_{V^{T}}z}{\textsf{SNR}}+\frac{(q_{V^{T}}-1)q+1}{1+\textsf{SNR}}-z\,,$
(17) $\displaystyle c$
$\displaystyle=\frac{(q-1)q_{V^{T}}}{\textsf{SNR}}+q_{V^{T}}z\big{(}1+\textsf{SNR}^{-1}\big{)}-q_{V^{T}}+1\,,$
(18) $\displaystyle d$
$\displaystyle=-q_{V^{T}}\big{(}1+\textsf{SNR}^{-1}\big{)}\,.$ (19)

While one can solve this cubic equation analytically, the resulting equations
are unwieldy, allowing for little direct insights. Therefore we rely on
analytic approximations in the two noise limits as well as numerics. Taking
the zero noise limit, $\textsf{SNR}\rightarrow\infty$, the equation reduces to
a quadratic polynomial which we solve and evaluate Eq. (14), to find the exact
eigenvalue density

$\displaystyle\rho^{\infty}(\lambda)$
$\displaystyle=\frac{\sqrt{(\lambda-\lambda_{-}^{\infty})(\lambda_{+}^{\infty}-\lambda)}}{2\pi\lambda{\bar{\sigma}_{X}}^{-2}q_{V^{T}}^{-1}}+(1-q_{V^{T}})\delta(\lambda)\,,$
(20)

where
$\bar{\sigma}_{X}^{2}\equiv\sigma_{X}^{2}/m=\sigma_{U}^{2}\sigma_{V}^{2}$. The
delta function represents the $N-m$ eigenvalues of the NECM which are
trivially zero. The $m$ non-trivial eigenvalues lie in a finite interval with
bounds

$\displaystyle\lambda_{\pm}^{\infty}=\bar{\sigma}_{X}^{-2}\left(1\pm\sqrt{q_{V^{T}}}^{-1}\right)^{2}\,.$
(21)

The eigenvalue density vanishes everywhere else. For finite SNR we solve Eq.
(15) numerically. A comparison of eigenvalue densities in the different noise
regimes, including the MP density [20] representing the pure noise limit for
$\textsf{SNR}\rightarrow 0$, is shown in the top plot of Fig. 2(a).

In Online Supplementary materials C.4.1, we derive an analytic approximation
for the eigenvalue bounds, $\lambda^{\textsf{SNR}}_{\pm}$, at finite SNR,
given by a weighted average of the zero noise bounds,
$\lambda_{\pm}^{\infty}$, and those of the MP density
$\lambda_{\pm}^{\text{MP}}=(1\pm\sqrt{q})^{2}$ 111Marčenko-Pastur (MP)
distribution,
$\rho^{\text{MP}}(\lambda)=\sqrt{(\lambda-\lambda^{\text{MP}}_{-})(\lambda^{\text{MP}}_{+}-\lambda)}/(2\pi
q\lambda)$, with $\lambda^{\text{MP}}_{\pm}=(1\pm\sqrt{q})^{2}$, [19].:

$\displaystyle\lambda^{\textsf{SNR}}_{\pm}$
$\displaystyle\approx\big{(}1+\textsf{SNR}^{-1}\big{)}^{-1}\,\lambda_{\pm}^{\infty}+\big{(}1+\textsf{SNR}\big{)}^{-1}\,\lambda_{\pm}^{\text{MP}}\,.$
(22)

In the bottom plot of Fig. 2(a), we show a comparison between this
approximation and the true, numerically computed, bounds at different values
of the SNR. The approximation is good everywhere, with the largest deviation
at $\textsf{SNR}\sim 10^{-1}$. The magnitude of this value is a result of the
particular choice we have made to define the SNR. The part of the eigenvalue
spectrum associated with the pure latent feature signal lies outside of the
interval,
$(1+\textsf{SNR})^{-1}\times[\lambda_{-}^{\text{MP}},\lambda_{+}^{\text{MP}}]$.
Eigenvalues within this interval, shown as a striped band, correspond to
noise. Notice that as the SNR is increased the noise range is shifted to the
left compared to the MP range due to the presence of the latent features
signal renormalizing the NECM. Thus using the naïve MP bounds for rejection of
eigenvalues as noise – a common procedure in data analysis – may be a bad
practice.

A different picture emerges in the intensive limit. Here the polynomial
equation for $\mathfrak{g}_{\mathbf{C}}$ is of the sixth order with lengthy
expressions for the polynomial coefficients, cf. Eq. (185). For the
$\textsf{SNR}\rightarrow\infty$ limit, we find the following analytic
expression for the density

$\displaystyle\rho^{\infty}(\lambda)$
$\displaystyle=\frac{\sqrt{(\lambda-\lambda_{-}^{\infty})(\lambda_{+}^{\infty}-\lambda)}}{2\pi\lambda{\bar{\sigma}_{X}}^{-2}(1+q)q_{V^{T}}^{-1}}+(1-q_{V^{T}})\delta(\lambda)\,,$
(23)

with the eigenvalue bounds given by:

$\displaystyle\lambda_{\pm}^{\infty}=\bar{\sigma}_{X}^{-2}\left(\sqrt{1+q}\pm\sqrt{q_{V^{T}}}^{-1}\right)^{2}\,.$
(24)

The form of the eigenvalue density at different levels of noise is shown in
the top plot of Fig. 2(b). For large values of the SNR, there is a gap in the
density between the eigenvalues corresponding to noise, and those
corresponding to the signal. The gap closes at lower SNR, and the combined
density converges to the MP density for $\textsf{SNR}\rightarrow 0$. As for
the classical limit, the approximate expression Eq. (22) for the bounds of the
signal part of the density is good, showing the largest deviation for
$\textsf{SNR}\sim 10^{-1}$, as shown in the bottom plot of Fig. 2(b). In
Online Supplementary material C.4.2 we estimate the SNR value at which the
eigenvalue density splits to be
$\textsf{SNR}_{\text{split}}\approx(\lambda_{+}^{\text{MP}}-\lambda_{-}^{\text{MP}})/\lambda_{-}^{\infty}$.
The existence of the gap suggests that having more observed variables (and
hence more data to define the latent components) makes it easier to
distinguish the signal from noise. This suggests that high throughput data
sets, where individual variables cannot be well sampled ($T<N$), might still
be very valuable, as long as $m\ll N,T$.

_Discussion._ We calculated statistical properties of data with latent linear
features, including the density of pairwise correlations, and density of
eigenvalues of the NECM. Our analysis provides two important insights. First,
by looking at the distribution of the correlations and their eigenvalues, one
can understand whether the latent features model is a reasonable model for the
data at hand, and also ballpark the number of the latent components.
Importantly, even if the eigenvalue density does not have a prominent gap, one
can understand that the underlying model has a latent structure, which
manifests as a distortion of the MP sea of eigenvalues. This is because our
signal matrix in Eq. (1) is stochastic, in contrast to spiked covariance
models, where deterministic perturbations appear as delta functions in the
spectrum and are detectable only as true outliers [21, 22, 23]. Second, since
the spectrum of the noise correlations in the latent features model is shifted
compared to the MP model, one should not use a simple truncation at the right
edge of the MP density to distinguish which of the principal components are
statistically significant.

We note that a lot of ink has been expended to decide whether a few variables
measured well many times are more or less valuable than many variables
measured infrequently and with high noise. We find that, in the intensive
limit, the many measured variables lead to separation of the noise and the
signal eigenvalues, resulting in a potentially more accurate discrimination of
noise from signal. However, the number of observations, and their quality
contribute to the SNR, a high value of which is also required for the opening
of the signal-noise gap. Thus the quality and the quantity of measurements all
contribute to the value data in very specific ways, which we now understand.

Our analysis involved a few approximations. The strongest of these was in
neglecting the signal-noise contributions in the computation of the eigenvalue
density. Including these contributions would make the polynomial equation for
the Stieltjes transform substantially more complicated. However, we do not
expect significant qualitative changes to the structure of the eigenvalue
density in the limits considered. An additional limitation is that the
approximation for the bounds of the eigenvalue density, Eq. (22), is strictly
only valid in the extreme noise limits, $\textsf{SNR}\rightarrow 0$ and
$\textsf{SNR}\rightarrow\infty$. In deriving the analytic density of
correlations, we assumed $q_{U}\rightarrow 0$ in accordance with the limits in
Eq. (9) and Eq. $\eqref{eq:intensive}$, and worked to leading order in $T$. We
expect the quality of the analytic density to improve if these assumptions are
removed, in particular in the regime of finite $\mathsf{SNR}$, cf. Fig. S3.

Finally, to connect our results with the analysis of real data, additional
steps are required. To model data, the means and variances of rows and columns
of our matrix model have to be fit to the data. Similarly, methods to estimate
the SNR from the data and to determine whether the Gaussian assumption for the
distribution of the noise and the latent components is valid will need to be
developed. Some of such extension may remain analytically tractable. Further,
for our model and its extensions, it is also important to calculate the
expected overlap of empirical eigenvectors with their true values.

###### Acknowledgements.

We are grateful to Marc Potters for his insightful comments. IN thanks the
Aspen Center for Physics, partially funded by NSF Grant PHY-1607611, for
hospitality. PF thanks Mirna Kramar for continued discussions and support.
This work was supported in part by the Simons Foundation Grants 400425 (PF)
and 827661 (IN) and by NSF Grants BCS-1822677, PHY-2014173, and PHY-2010524
(IN).

## References

  * Killer  _et al._ [2016] C. Killer, T. Bockwoldt, S. Schütt, M. Himpel, A. Melzer, and A. Piel, Phase separation of binary charged particle systems with small size disparities using a dusty plasma, Phys. Rev. Lett. 116, 115002 (2016).
  * Valentine  _et al._ [2001] M. T. Valentine, P. D. Kaplan, D. Thota, J. C. Crocker, T. Gisler, R. K. Prud’homme, M. Beck, and D. A. Weitz, Investigating the microenvironments of inhomogeneous soft materials with multiple particle tracking, Physical Review E 64, 061506 (2001).
  * Sinhuber  _et al._ [2019] M. Sinhuber, K. Van Der Vaart, R. Ni, J. G. Puckett, D. H. Kelley, and N. T. Ouellette, Three-dimensional time-resolved trajectories from laboratory insect swarms, Scientific Data 6, 1 (2019).
  * Lusch  _et al._ [2018] B. Lusch, J. N. Kutz, and S. L. Brunton, Deep learning for universal linear embeddings of nonlinear dynamics, Nature communications 9, 1 (2018).
  * Schanz  _et al._ [2016] D. Schanz, S. Gesemann, and A. Schröder, Shake-the-box: Lagrangian particle tracking at high particle image densities, Experiments in Fluids 57, 1 (2016).
  * [6] NOAA Physical Sciences Laboratory, Gridded climate data, https://psl.noaa.gov/data/gridded/, accessed: 2021-06-30.
  * Natale  _et al._ [2018] J. L. Natale, D. Hofmann, D. G. Hernández, and I. Nemenman, Reverse-Engineering Biological Networks From Large Data Sets, in  _Quantitative Biology: Theory, Computational Methods and Examples of Models_ , edited by B. Munsky, L. Tsimring, and W. S. Hlavacek (MIT Press, 2018).
  * Meshulam  _et al._ [2017] L. Meshulam, J. L. Gauthier, C. D. Brody, D. W. Tank, and W. Bialek, Collective Behavior of Place and Non-place Neurons in the Hippocampal Network., Neuron 96, 1178 (2017).
  * Stephens  _et al._ [2008] G. J. Stephens, B. Johnson-Kerner, W. Bialek, and W. S. Ryu, Dimensionality and dynamics in the behavior of c. elegans, PLoS Comput Biol 4, e1000028 (2008).
  * Berman  _et al._ [2014] G. J. Berman, D. M. Choi, W. Bialek, and J. W. Shaevitz, Mapping the stereotyped behaviour of freely moving fruit flies, Journal of The Royal Society Interface 11, 20140672 (2014).
  * Weisser  _et al._ [2017] W. Weisser, C. Roscher, S. Meyer, A. Ebeling, G. Luo, E. Allan, H. Bessler, R. Barnard, N. Buchmann, F. Buscot, _et al._ , Biodiversity effects on ecosystem functioning in a 15-year grassland experiment: Patterns, mechanisms, and open questions, Basic and Applied Ecology 23, 1 (2017).
  * Dell  _et al._ [2014] A. I. Dell, J. A. Bender, K. Branson, I. D. Couzin, G. G. de Polavieja, L. P. Noldus, A. Pérez-Escudero, P. Perona, A. D. Straw, M. Wikelski, and U. Brose, Automated image-based tracking and its application in ecology, Trends in Ecology & Evolution 29, 417 (2014).
  * Pandarinath  _et al._ [2018] C. Pandarinath, D. J. O’Shea, J. Collins, R. Jozefowicz, S. D. Stavisky, J. C. Kao, E. M. Trautmann, M. T. Kaufman, S. I. Ryu, L. R. Hochberg, _et al._ , Inferring single-trial neural population dynamics using sequential auto-encoders, Nature methods 15, 805 (2018).
  * Schwab  _et al._ [2014] D. J. Schwab, I. Nemenman, and P. Mehta, Zipf’s law and criticality in multivariate data without fine-tuning, Physical review letters 113, 068102 (2014).
  * Morrell  _et al._ [2021] M. C. Morrell, A. J. Sederberg, and I. Nemenman, Latent Dynamical Variables Produce Signatures of Spatiotemporal Criticality in Large Biological Systems, Physical review letters 126, 118302 (2021).
  * Gallego  _et al._ [2017] J. A. Gallego, M. G. Perich, L. E. Miller, and S. A. Solla, Neural manifolds for the control of movement, Neuron 94, 978 (2017).
  * Page  _et al._ [2020] J. Page, M. P. Brenner, and R. R. Kerswell, Revealing the state space of turbulence using machine learning, arXiv preprint arXiv:2008.07515 (2020).
  * Nieh  _et al._ [2021] E. H. Nieh, M. Schottdorf, N. W. Freeman, R. J. Low, S. Lewallen, S. A. Koay, L. Pinto, J. L. Gauthier, C. D. Brody, and D. W. Tank, Geometry of abstract learned knowledge in the hippocampus, Nature 595, 80 (2021).
  * Potters and Bouchaud [2020] M. Potters and J.-P. Bouchaud,  _A First Course in Random Matrix Theory: For Physicists, Engineers and Data Scientists_ (Cambridge University Press, 2020).
  * Marčenko and Pastur [1967] V. A. Marčenko and L. A. Pastur, DISTRIBUTION OF EIGENVALUES FOR SOME SETS OF RANDOM MATRICES, Mathematics of the USSR-Sbornik 1, 457 (1967).
  * Sengupta and Mitra [1999] A. M. Sengupta and P. P. Mitra, Distributions of singular values for some random matrices, Physical Review E 60, 3389 (1999).
  * Loubaton and Vallet [2011] P. Loubaton and P. Vallet, Almost Sure Localization of the Eigenvalues in a Gaussian Information Plus Noise Model. Application to the Spiked Models., Electronic Journal of Probability 16, 1934 (2011).
  * Capitaine and Donati-Martin [2016] M. Capitaine and C. Donati-Martin, Spectrum of deformed random matrices and free probability (2016), arXiv:1607.05560 [math.PR] .
  * Hotelling [1953] H. Hotelling, New light on the correlation coefficient and its transforms, Journal of the Royal Statistical Society. Series B (Methodological) 15, 193 (1953).
  * [25] H. (https://math.stackexchange.com/users/6460/henry), Sum of $n$ i.i.d Beta-distributed variables, Mathematics Stack Exchange, https://math.stackexchange.com/q/3096929 (version: 2019-02-02).
  * Note [1] Marčenko-Pastur (MP) distribution, $\rho^{\text{MP}}(\lambda)=\sqrt{(\lambda-\lambda^{\text{MP}}_{-})(\lambda^{\text{MP}}_{+}-\lambda)}/(2\pi q\lambda)$, with $\lambda^{\text{MP}}_{\pm}=(1\pm\sqrt{q})^{2}$, [19].

##  Appendix A Data distribution for the latent feature model with no noise,
its variance and large $m$ limit

Each entry $X_{ij}$ of the latent features data matrix $\mathbf{U}\mathbf{V}$
is given by the sum of $m$ products of two i.i.d. Gaussian random variables
$u\sim\mathcal{N}(0,\sigma_{U}^{2})$ and $v\sim\mathcal{N}(0,\sigma_{V}^{2})$:

$\displaystyle X_{ij}\sim\sum_{\mu=1}^{m}uv\,.$ (25)

The product, $x=uv$, is distributed according to the normal product
distribution [1]:

$\displaystyle
x\sim\frac{K_{0}\left(\frac{|x|}{\sigma_{U}\sigma_{V}}\right)}{\pi\sigma_{U}\sigma_{V}}\,,$
(26)

where $K_{\nu}$ is the modified Bessel function of the second kind:

$\displaystyle
K_{\nu}(x)=\frac{\Gamma\left(\nu+\frac{1}{2}\right)(2x)^{\nu}}{\sqrt{\pi}}\int_{0}^{\infty}dq\frac{\cos(q)}{(x^{2}+q^{2})^{\nu+1/2}}\,.$
(27)

To derive the probability density of the latent feature model entries
$X_{ij}$, we first compute the characteristic function $\varphi_{x}$ by taking
the Fourier transform of the normal product distribution. We then use the fact
that the characteristic function $\varphi_{X}$ of the sum of $m$ products $x$
is given by $\varphi_{X}=\left(\varphi_{x}\right)^{m}$. The inverse Fourier
transform of $\varphi_{X}$ then yields the sought after probability density.

Specifically, the characteristic function $\varphi_{x}$ of the normal product
distribution is

$\displaystyle\varphi_{x}(t)$
$\displaystyle=\mathbb{E}(e^{itx})=\int_{-\infty}^{\infty}dx\,\frac{K_{0}\left(\frac{|x|}{\sigma_{U}\sigma_{V}}\right)}{\pi\sigma_{U}\sigma_{V}}e^{itx}$
$\displaystyle=\int_{-\infty}^{\infty}dx\,\frac{K_{0}\left(|x|\right)}{\pi}e^{it\sigma_{U}\sigma_{V}x}$
$\displaystyle=\frac{1}{\pi}\int_{-\infty}^{\infty}dx\,\int_{0}^{\infty}dq\,\frac{\cos(q)}{\sqrt{|x|^{2}+q^{2}}}e^{it\sigma_{U}\sigma_{V}x}$
$\displaystyle=\frac{1}{\pi}\int_{-\infty}^{\infty}dx\,\int_{0}^{\infty}dq\,\frac{\cos(xq)}{\sqrt{1+q^{2}}}e^{it\sigma_{U}\sigma_{V}x}$
$\displaystyle=\frac{1}{\pi}\int_{0}^{\infty}dq\frac{1}{\sqrt{1+q^{2}}}$
$\displaystyle\quad\int_{-\infty}^{\infty}\frac{dx}{2\pi}\,\left(e^{ix(q+\sigma_{U}\sigma_{V}t)}+e^{ix(\sigma_{U}\sigma_{V}t-q)}\right)$
$\displaystyle=\int_{0}^{\infty}dq\frac{1}{\sqrt{1+q^{2}}}\left(\delta(\sigma_{U}\sigma_{V}t+q)+\delta(\sigma_{U}\sigma_{V}t-q)\right)$
$\displaystyle=\frac{1}{\sqrt{1+\sigma_{U}^{2}\sigma_{V}^{2}t^{2}}}$ (28)

for $t\in\mathbb{R}\setminus\\{0\\}$, and $\delta(\cdot)$ is the Dirac delta
function.

The characteristic function $\varphi_{X}$ of the sum of $m$ products $x$ is
given by

$\displaystyle\varphi_{X}=\left(\varphi_{x}\right)^{m}=\left(1+\sigma_{U}^{2}\sigma_{V}^{2}t^{2}\right)^{-\frac{m}{2}}$
(29)

Finally, performing the inverse transformation we obtain the probability
density function of the sum

$\displaystyle{\rm pdf}(X)$
$\displaystyle=\int_{-\infty}^{\infty}\frac{dt}{2\pi}\,e^{-itX}\frac{1}{\left(1+\sigma_{U}^{2}\sigma_{V}^{2}t^{2}\right)^{\frac{m}{2}}}$
$\displaystyle=\int_{-\infty}^{\infty}\frac{dt}{2\pi}\,e^{-itX}$
$\displaystyle\quad\int_{0}^{\infty}dq\,\frac{\delta(\sigma_{U}\sigma_{V}t+q)+\delta(\sigma_{U}\sigma_{V}t-q)}{\left(1+q^{2}\right)^{\frac{m}{2}}}$
$\displaystyle=\frac{1}{\sigma_{U}\sigma_{V}}\int_{0}^{\infty}dq$
$\displaystyle\quad\int_{-\infty}^{\infty}\frac{dt}{2\pi}\,e^{-\frac{itX}{\sigma_{U}\sigma_{V}}}\frac{\delta(t+q)+\delta(t-q)}{\left(1+q^{2}\right)^{\frac{m}{2}}}$
$\displaystyle=\frac{1}{\pi\sigma_{U}\sigma_{V}}\int_{0}^{\infty}dq\frac{\cos(\frac{q|X|}{\sigma_{U}\sigma_{V}})}{\left(1+q^{2}\right)^{\frac{m}{2}}}$
$\displaystyle=\left[\frac{|X|}{\sigma_{U}\sigma_{V}}\right]^{m-1}\frac{1}{\pi\sigma_{U}\sigma_{V}}\int_{0}^{\infty}dq\frac{\cos(q)}{\left(\frac{|X|^{2}}{\sigma_{U}^{2}\sigma_{V}^{2}}+q^{2}\right)^{\frac{m}{2}}}$
$\displaystyle=\left[\frac{|X|}{2}\right]^{\frac{m-1}{2}}\frac{K_{\frac{m-1}{2}}\left(\frac{|X|}{\sigma_{U}\sigma_{V}}\right)}{(\sigma_{U}\sigma_{V})^{\frac{m+1}{2}}\sqrt{\pi}\Gamma\left(\frac{m}{2}\right)}.$
(30)

Since the probability density function of $X$ is symmetric around zero, the
mean of the distribution vanishes:

$\displaystyle\mu_{{UV}}$ $\displaystyle=\int_{-\infty}^{\infty}dXX{\rm
pdf}(X)=0\,.$ (31)

Figure S1: Comparison between the variances computed from simulated data of
the latent feature model for $m=4,5,6,7,8$ (dashed gray lines) and the
numerically evaluated limit expression for $\sigma_{UV}^{2}$ in Eq. (37) as a
function of the limit value $Z$, for even $m$ values. We have chosen
$\sigma^{2}_{U}=\sigma^{2}_{V}=1$.

The variance is

$\displaystyle\sigma_{{UV}}^{2}$
$\displaystyle=\int_{-\infty}^{\infty}dXX^{2}{\rm pdf}(X)$
$\displaystyle=\frac{1}{2^{\frac{m-1}{2}}\sqrt{\pi}\Gamma\left(\frac{m}{2}\right)}\times$
$\displaystyle\quad\times\int_{-\infty}^{\infty}\frac{dX}{\sigma_{U}\sigma_{V}}\left[\frac{|X|}{\sigma_{U}\sigma_{V}}\right]^{\frac{m-1}{2}}|X|^{2}K_{\frac{m-1}{2}}\left(\frac{|X|}{\sigma_{U}\sigma_{V}}\right)$
$\displaystyle=\frac{2\sigma_{U}^{2}\sigma_{V}^{2}}{2^{\frac{m-1}{2}}\sqrt{\pi}\Gamma\left(\frac{m}{2}\right)}\int_{0}^{\infty}dX|X|^{\frac{m+3}{2}}K_{\frac{m-1}{2}}(|X|)\,.$
(32)

The integral above can be evaluated in terms of generalized hypergeometric
functions [2]. We present the calculation for when $m$ is even in detail:

$\displaystyle\int_{0}^{\infty}dXX^{\alpha-1}K_{\nu}(X)=\left[\Sigma(\nu,\alpha;Z)+\Sigma(-\nu,\alpha;Z)\right]_{0}^{\infty},$
(33)

where

$\displaystyle\Sigma(\nu,\alpha;Z)$ $\displaystyle\equiv-\frac{2^{\nu-1}\pi
Z^{\alpha-\nu}\text{csc}(\pi\nu)}{(\nu-\alpha)\Gamma(1-\nu)}\times$
$\phantom{}{}_{1}F_{2}\left(\frac{\alpha-\nu}{2};1-\nu,\frac{\alpha-\nu}{2}+1;\frac{Z^{2}}{4}\right)$
(34)

with parameters

$\displaystyle\alpha\equiv\frac{m+5}{2}\text{ and }\nu\equiv\frac{m-1}{2},$
(35)

and $\phantom{}{}_{1}F_{2}$ is the generalized hypergeometric function

$\phantom{}{}_{1}F_{2}(a_{1};b_{1},b_{2};z)=\sum_{k=0}^{\infty}\frac{(a_{1})z^{k}}{(b_{1})_{k}(b_{2})_{k}k!}\,.$
(36)

In the expression above, $(\cdot)_{k}$ is the Pochhammer symbol, and
$\text{csc}(\cdot)$ is the cosecant. Since $m$ is even, we also have
$\nu\notin\mathbb{Z}$. Putting everything together, we obtain the following
expression for the variance

$\displaystyle\sigma_{{UV}}^{2}$
$\displaystyle=\frac{\left[\Sigma(\nu,\alpha;Z)+\Sigma(-\nu,\alpha;Z)\right]_{0}^{\infty}}{2^{\frac{m-3}{2}}\sqrt{\pi}\Gamma\left(\frac{m}{2}\right)\sigma_{U}^{-2}\sigma_{V}^{-2}}$
$\displaystyle=\lim_{Z\rightarrow\infty}\frac{\Sigma(\nu,\alpha;Z)+\Sigma(-\nu,\alpha;Z)}{2^{\frac{m-3}{2}}\sqrt{\pi}\Gamma\left(\frac{m}{2}\right)\sigma_{U}^{-2}\sigma_{V}^{-2}}\,,$
(37)

where we have used the fact that the numerator after the first equality
vanishes at $Z=0$. We can evaluate the limit $Z\rightarrow\infty$, on the
right-hand side numerically as shown in Fig. S1 and find that the variance of
the latent feature data values is

$\displaystyle\sigma_{UV}^{2}=m\sigma_{U}^{2}\sigma_{V}^{2}\,.$ (38)

This is in agreement with the intuition that every latent dimension
contributes its own variance to the variance of the data.

We note that, for large values of the number of latent features $m$, the
distribution (30) becomes normal, in agreement with the law of large numbers:

$\displaystyle{\rm
pdf}(X_{ij})=\frac{1}{\sqrt{2\pi\sigma_{UV}^{2}}}e^{-\frac{X^{2}}{2\sigma_{UV}^{2}}}\,.$
(39)

Crucially, the variance of $X_{ij}$ remains $m$-dependent. Figure S2 for
compares exact analytical expression of the probability distribution and its
Gaussian approximation to numerical simulations.

As a final note, if we were interested in the distribution of data with noise,
we would need to convolve the density in Eq. (30) with the Gaussian density of
the noise.

Figure S2: Comparison of simulated data (gray) and the analytical distribution
(orange). In the limit of large $m$, the distribution approaches a Gaussian
form (blue). Simulated data constitutes a single realization of the model with
$\sigma_{U}^{2}=\sigma_{V}^{2}=1$.

##  Appendix B Probability density of the correlation coefficients

For our latent features model with noise, here we calculate the probability
distribution of entries in the empirical data correlation matrix. Before doing
this, a few notes are in order. First, the correlations depend on the basis,
in which variables are measured, becoming a diagonal matrix in the special
case when the measured variables are the principal axes of the data cloud.
Thus to make statements independent of the basis, we consider the distribution
of typical correlations, or correlations in the basis random w. r. t. the
principal axes of the data. For a given realization, the $N$-dimensional data
cloud is typically anisotropic, with $m<N$ long directions dominated by the
latent feature signal and $N-m$ short directions dominated by noise. When
$N\gg m$, principal axes of the data cloud do not align with the measured
variables for the vast majority of random rotations, and correlations between
any random pair of variables have contributions from all latent dimensions.
Thus we expect the number of latent dimensions to be imprinted in the
distribution of the elements of the correlation matrix, so that the statistics
of the elements carries information about the underlying structure of the
model.

###  B.1 Preliminaries: Density of the correlation coefficient of two random
Gaussian variables

The correlation coefficient of two independent zero-mean variables $x$ and $y$
sampled $T$ times is

$\displaystyle
r=\frac{1}{T}\sum_{t}\frac{x_{t}y_{t}}{\sigma_{x}\sigma_{y}}\,,$ (40)

where the vectors’ components are mutually independent, i.i.d. random
variables. The correlation coefficient is distributed according to [24]

$\displaystyle\mathrm{pdf}(r)=\frac{\Gamma\left(\frac{T}{2}\right)}{\Gamma\left(\frac{1}{2}\right)\Gamma\left(\frac{T-1}{2}\right)}(1-r^{2})^{\frac{T-3}{2}}\,.$
(41)

This can be rewritten in terms of a Beta distribution

$\displaystyle\text{Beta}(x;\alpha,\beta)=\frac{1}{\text{B}(\alpha,\beta)}x^{\alpha-1}(1-x)^{\beta-1}\,,$
(42)

where $x\in[0,1]$ and $\text{B}(\alpha,\beta)$ is the Beta function.
Specifically, the density of correlations is given by the symmetric Beta
distribution

$\displaystyle\mathrm{pdf}(r)=\text{Beta}\left(r;\alpha,\alpha;\ell=-1,s=2\right)\,,$
(43)

where the location $\ell$ and scale $s$ are set such that the density is
defined on the interval of correlation values [-1,1], and

$\displaystyle\alpha=\frac{T-1}{2}\,.$ (44)

We also note that the variance of a symmetric Beta distribution with the scale
$s=2$ is

$\displaystyle\text{var}=\frac{s^{2}}{4(2\alpha+1)}=\frac{1}{2\alpha+1}\,.$
(45)

###  B.2 Density of correlations in the latent feature model

Figure S3: Distribution of pairwise correlations in the regimes of finite and
small signal-to-noise ratio with $m=2,4$ and $30$ latent features. Analytic
form (magenta) and simulated data (gray). (a) $\mathsf{SNR}=20$ and (b)
$\mathsf{SNR}=10^{-5}$ (large noise limit). Each simulation is run with $N=80$
variables and $T=2048$ observations and constitutes $1000$ independent model
realisations.

There are multiple contributions to the correlations among the measured
variables. We compute them individually, and then combine the contributions.
We find that each contribution is distributed according to a symmetric Beta
distribution. To obtain the overall density, we approximate the sum of Beta
distributions by a single Beta distribution, the parameter of which is
obtained by matching the variance to the sum of the variances of the
individual components. To perform these analyses, we only keep terms to the
leading order in the $\textsf{SNR}\rightarrow 0$ or the
$\textsf{SNR}\rightarrow\infty$ limit. Further, we assume that $q_{U}$ is
small in accordance with the classical and intensive regimes limits.

We start with the pure noise contribution to the correlations

$\displaystyle(c_{R})_{pq}$
$\displaystyle=\frac{1}{T}\sum_{t}\frac{R^{T}_{pt}R_{tq}}{\sigma_{p}^{\mathrm{n}}\sigma_{q}^{\mathrm{n}}}\,,$
(46) $\displaystyle({\sigma_{q}^{\mathrm{n}}})^{2}$
$\displaystyle=\frac{1}{T}\sum_{t}R^{T}_{qt}R_{tq}\,.$ (47)

The expression on the right-hand side is the correlation coefficient between
two random Gaussian variables. Using Eq. (43), we arrive at

$\displaystyle\mathrm{pdf}((c_{R})_{pq})=\text{Beta}\left((c_{R})_{pq};\alpha_{\text{n}},\alpha_{\text{n}};-1,2\right),\;p\neq
q,$ (48)

with

$\displaystyle\alpha_{\text{n}}=\frac{T-1}{2}\,.$ (49)

and the variance of this density is

$\displaystyle\text{var}_{\text{n}}=T^{-1}\,.$ (50)

Next we compute the density of the pure signal contribution

$\displaystyle(c_{UV})_{pq}$
$\displaystyle=\frac{1}{T\sigma^{\mathrm{s}}_{p}\sigma^{\mathrm{s}}_{q}}\sum_{t}\left(\sum_{\mu}V_{p\mu}U_{\mu
t}\right)\left(\sum_{\nu}U_{t\nu}V_{\nu q}\right)\,,$ (51)
$\displaystyle(\sigma^{\mathrm{s}}_{p})^{2}$
$\displaystyle=\frac{1}{T}\sum_{t}\left(\sum_{\mu}V_{p\mu}U_{\mu
t}\right)\left(\sum_{\nu}U_{t\nu}V_{\nu p}\right)\,,$ (52)

and similarly for $\sigma_{q}$. Rearranging, we find

$\displaystyle(c_{UV})_{pq}$
$\displaystyle=\frac{1}{\sigma^{\mathrm{s}}_{p}\sigma^{\mathrm{s}}_{q}}\sum_{\mu\nu}V_{p\mu}V_{\nu
q}\left(\frac{1}{T}\sum_{t}U_{\mu t}U_{t\nu}\right)\,,$ (53)
$\displaystyle(\sigma^{\mathrm{s}}_{p})^{2}$
$\displaystyle=\sum_{\mu\nu}V_{p\mu}V_{\nu p}\left(\frac{1}{T}\sum_{t}U_{\mu
t}U_{t\nu}\right)\,.$ (54)

The expression in parentheses of both of the equations above is a
(co)-variance of Gaussian random numbers. For $\mu=\nu$, it follows the scaled
$\chi^{2}$-distribution with $T$ degrees of freedom. For $\mu\neq\nu$, it is
given by a rescaled version of the distribution in Eq. (30), with $T$ instead
of $m$. Crucially, the variance of either is $1/T$. Thus in the limit $q\to
0$, the terms in parentheses are
$\sigma^{2}_{U}\delta_{\mu\nu}+\mathcal{O}(T^{-1/2})$, where the correction
$\mathcal{O}(T^{-1/2})$ is probabilistic, but will be neglected in what
follows. We get

$\displaystyle(c_{UV})_{pq}$
$\displaystyle=\sigma_{U}^{2}\sum_{\mu\nu}\frac{V_{p\mu}V_{\nu
q}\delta_{\mu\nu}}{\sigma^{\mathrm{s}}_{p}\sigma^{\mathrm{s}}_{q}}=m\sigma_{U}^{2}\left(\frac{1}{m}\sum_{\mu}\frac{V_{p\mu}V_{\mu
q}}{\sigma^{\mathrm{s}}_{p}\sigma^{\mathrm{s}}_{q}}\right)\,,$ (55)
$\displaystyle(\sigma^{\mathrm{s}}_{p})^{2}$
$\displaystyle=m\sigma_{U}^{2}\left(\frac{1}{m}\sum_{\mu}V_{p\mu}^{2}\right)\,,$
(56)

We see that the sought after correlation is a correlation coefficient between
Gaussian variables, but with $m$ samples instead of $T$. Using again Eq. (43),
we write

$\displaystyle\mathrm{pdf}\left((c_{UV})_{pq}\right)=\text{Beta}\left((c_{UV})_{pq};\alpha_{\text{s}},\alpha_{\text{s}};-1,2\right),\,p\neq
q,$ (57)

with parameter

$\displaystyle\alpha_{\text{s}}=\frac{m-1}{2}\,.$ (58)

We remind the reader that Eq. (57) holds to $\mathcal{O}(T^{-1/2})$. The
variance of this density is

$\displaystyle\text{var}_{\text{s}}$ $\displaystyle=m^{-1}\,.$ (59)

This expression agrees with numerical simulations very well, cf. Fig. 1.

Finally, for the signal-noise cross terms in the correlation, we have

$\displaystyle(c_{(UV)^{T}R})_{pq}$
$\displaystyle=\frac{1}{T\sigma^{\mathrm{s}}_{p}\sigma_{q}^{\mathrm{n}}}\sum_{t}\sum_{\mu}V_{p\mu}U_{\mu
t}R_{tq}$
$\displaystyle=\frac{1}{\sigma^{\mathrm{s}}_{p}\sigma_{q}^{\mathrm{n}}}\sum_{\mu}V_{p\mu}\left(\frac{1}{T}\sum_{t}U_{\mu
t}R_{tq}\right)\,.$ (60)

For the quantity in parentheses in Eq. (60), we define

$\displaystyle r_{\mu p}\equiv\frac{1}{T}\sum_{t}U_{\mu t}R_{tq}\,.$ (61)

This is a covariance between two independent Gaussian random numbers and again
follows a rescaled form of the distribution in Eq. (30) with variance
$\sigma_{U}^{2}{\sigma_{q}^{\mathrm{n}}}^{2}T^{-1}$. Since $T$ is large, the
distribution approaches a Gaussian and we further define $r_{\mu
q}\equiv\sigma_{U}\sigma_{q}^{\mathrm{n}}T^{-1/2}r_{\mu q}^{\prime}$, such
that $r_{\mu q}^{\prime}$ is a unit Gaussian random variable. Thus we obtain

$(c_{(UV)^{T}R})_{pq}=\frac{\sigma_{U}m\sigma_{q}^{\mathrm{n}}T^{-1/2}}{\sigma^{\mathrm{s}}_{p}\sigma_{q}^{\mathrm{n}}}\left(\frac{1}{m}\sum_{\mu}V_{p\mu}r^{\prime}_{\mu
q}\right)\\\
=m^{1/2}T^{-1/2}\left(\frac{1}{m}\sum_{\mu}\frac{V_{p\mu}r^{\prime}_{\mu
q}}{\frac{1}{m}\sum_{\mu}{V_{p\mu}}^{2}}\right)\,,$ (62)

where we have extracted the factor of $m$ to highlight that the expression in
parenthesis is the correlation between Gaussian random numbers. From this,
using Eq. (43), we conclude that

$\displaystyle\mathrm{pdf}$ $\displaystyle\left((c_{(UV)^{T}R})_{pq}\right)=$
$\displaystyle\text{Beta}\left((c_{(UV)^{T}R})_{pq};\alpha_{\text{sn}},\alpha_{\text{sn}};-1,2\right),\,p\neq
q\,,$ (63)

with parameter

$\displaystyle\alpha_{\text{sn}}=\frac{m^{1/2}T^{1/2}-1}{2}\,.$ (64)

The variance of this density is

$\displaystyle\text{var}_{\text{sn}}$
$\displaystyle=m^{-1/2}T^{-1/2}=\sqrt{\text{var}_{\text{s}}\cdot\text{var}_{\text{n}}}\,.$
(65)

An analogous expression holds for the $R^{T}UV$ contribution.

The empirical correlation matrix is given by

$\displaystyle
c_{pq}=\frac{1}{T}\sum_{t}\frac{X_{pt}X_{tq}}{\sigma^{\mathrm{sn}}_{p}\sigma^{\mathrm{sn}}_{q}}\,,$
(66)

where

$\displaystyle(\sigma^{\mathrm{sn}}_{p})^{2}=(\sigma^{\mathrm{s}}_{p})^{2}+(\sigma^{\mathrm{n}}_{p})^{2}\,.$
(67)

Using Eqs (46), (51) and (60), the correlation matrix can be written as a
weighted sum of the three types of contributions

$\displaystyle c_{pq}$
$\displaystyle=\frac{\sigma^{\mathrm{s}}_{p}\sigma^{\mathrm{s}}_{q}}{\sigma^{\mathrm{sn}}_{p}\sigma^{\mathrm{sn}}_{q}}(c_{UV})_{pq}+\frac{\sigma^{\mathrm{s}}_{p}\sigma^{\mathrm{n}}_{q}}{\sigma^{\mathrm{sn}}_{p}\sigma^{\mathrm{sn}}_{q}}(c_{(UV)^{T}R})_{pq}$
$\displaystyle+\frac{\sigma^{\mathrm{n}}_{p}\sigma^{\mathrm{s}}_{q}}{\sigma^{\mathrm{sn}}_{p}\sigma^{\mathrm{sn}}_{q}}(c_{R^{T}UV})_{pq}+\frac{\sigma^{\mathrm{n}}_{p}\sigma^{\mathrm{n}}_{q}}{\sigma^{\mathrm{sn}}_{p}\sigma^{\mathrm{sn}}_{q}}(c_{R})_{pq}$
(68)

Each term on the right-hand side of this equation follows a Beta distribution
as computed above. However, the $\alpha$ parameter of each distribution is
modified by the corresponding weight in the above sum. Consequently, the
variance of each distribution is rescaled by the weight:

$\displaystyle\mathrm{var}^{\prime}_{\mathrm{s}}$
$\displaystyle=\frac{\sigma^{\mathrm{s}}_{p}\sigma^{\mathrm{s}}_{q}}{\sigma^{\mathrm{sn}}_{p}\sigma^{\mathrm{sn}}_{q}}\mathrm{var}_{\mathrm{s}}\,$
(69) $\displaystyle\mathrm{var}^{\prime}_{\mathrm{n}}$
$\displaystyle=\frac{\sigma^{\mathrm{n}}_{p}\sigma^{\mathrm{n}}_{q}}{\sigma^{\mathrm{sn}}_{p}\sigma^{\mathrm{sn}}_{q}}\mathrm{var}_{\mathrm{n}}\,$
(70) $\displaystyle\mathrm{var}^{\prime}_{\mathrm{sn}}$
$\displaystyle=\frac{\sigma^{\mathrm{s}}_{p}\sigma^{\mathrm{n}}_{q}}{\sigma^{\mathrm{sn}}_{p}\sigma^{\mathrm{sn}}_{q}}\mathrm{var}_{\mathrm{sn}}\,.$
(71)

To determine an expression for the combined distribution of signal and noise
correlations, we make use of the observation that the sum of Beta
distributions can be well approximated by a single Beta distribution [3]. We
determine the parameters of the Beta distribution by adding the means and
variances of the distributions in the sum and analytically match the parameter
of the single Beta distribution.

The means of the Beta distributions in Eq. (48), Eq. (57), and Eq. (B.2) are
zero and thus the mean of the density of the combined contributions is also
zero. Taking the sum of variances we obtain

$\displaystyle\mathrm{var}=\mathrm{var}^{\prime}_{\mathrm{s}}+\mathrm{var}^{\prime}_{\mathrm{n}}+\mathrm{var}^{\prime}_{\mathrm{sn}}+\mathrm{var}^{\prime}_{\mathrm{ns}}\,.$
(72)

In the limit when $T$ and $m$ are large enough such that contributions of
$\mathcal{O}(T^{-1/2})$ and $\mathcal{O}(m^{-1/2})$ can be neglected, we have
the following convergence of the empirical quantities

$\displaystyle(\sigma_{p}^{\mathrm{s}})^{2}$ $\displaystyle\rightarrow
m\sigma_{U}^{2}\sigma_{V}^{2}\,,$ (73)
$\displaystyle(\sigma_{p}^{\mathrm{n}})^{2}$
$\displaystyle\rightarrow\sigma^{2}\,,$ (74)
$\displaystyle(\sigma_{p}^{\mathrm{sn}})^{2}\,,(\sigma_{p}^{\mathrm{ns}})^{2}$
$\displaystyle\rightarrow m\sigma_{U}^{2}\sigma_{V}^{2}+\sigma^{2}\,.$ (75)

Consequently the variances of the contributions take the form

$\displaystyle\mathrm{var}^{\prime}_{\mathrm{s}}$
$\displaystyle\rightarrow\frac{m^{-1}}{1+\textsf{SNR}^{-1}}\,,$ (76)
$\displaystyle\mathrm{var}^{\prime}_{\mathrm{n}}$
$\displaystyle\rightarrow\frac{T^{-1}}{1+\textsf{SNR}}\,,$ (77)
$\displaystyle\mathrm{var}^{\prime}_{\mathrm{sn}}\,,\mathrm{var}^{\prime}_{\mathrm{ns}}$
$\displaystyle\rightarrow\frac{m^{-1/2}T^{-1/2}}{\sqrt{1+\textsf{SNR}}\sqrt{1+\textsf{SNR}^{-1}}}\,,$
(78)

Thus, in this limit, the variance of the Beta distribution, Eq. (72), is of
the form

var
$\displaystyle\approx\left(\frac{m^{-1/2}}{\sqrt{1+\textsf{SNR}^{-1}}}+\frac{T^{-1/2}}{\sqrt{1+\textsf{SNR}}}\right)^{2}\,.$
(80)

Finally, from the relation in Eq. (45), we obtain the parameter $\alpha$ of
the sought after Beta distribution.

$\displaystyle\alpha=\frac{\text{var}^{-1}-1}{2}\,.$ (81)

A comparison between the analytic form of the density and simulated data is
shown in Fig. 1 for $\mathsf{SNR}\rightarrow\infty$, and in Fig. S3 for finite
$\mathsf{SNR}$ and $\mathsf{SNR}\rightarrow 0$. In the extreme noise limits,
the analytic form closely matches the simulation. In the large noise limit of
$\mathsf{SNR}\rightarrow 0$, shown in Fig. S3 (b), the density is close to a
Gaussian, because the number of observations $T$ is large. In the regime of
finite $\mathsf{SNR}$, shown in Fig. S3 (a), deviations between the analytic
form and the simulation appear for small values of $m$. We expect that these
deviations will disappear by removing the various approximations made in the
above analytic derivation.

##  Appendix C Spectrum of the normalized empirical covariance Matrix

To compute the eigenvalue density of the NECM $\mathbf{C}$, we use methods of
Random Matrix Theory [4]. The standard approach is to compute the finite size
Stieljtes transform

$\displaystyle
g^{N}_{\mathbf{C}}(z)=\frac{1}{N}\text{Tr}(z\mathbf{I}-\mathbf{C})^{-1}\,,$
(82)

where $\mathbf{I}$ is the identity matrix, $z\in\mathbb{C}$ and
$g^{N}_{\mathbf{C}}$ is a complex function. In the limit of large matrices –
large $N$ or thermodynamic limit – the finite size Stieltjes transform
becomes, $\mathfrak{g}_{\mathbf{C}}(z)$. Then the eigenvalue density is
obtained as the imaginary part of the limit of the Stieltjes transform:

$\displaystyle\rho(\lambda)=\frac{1}{\pi}\lim_{\eta\rightarrow
0^{+}}\Im\mathfrak{g}(z=\lambda-i\eta)\,,$ (83)

where $\Im$ denotes the imaginary part.

We start with writing again the definition of the normalized empirical
covariance matrix (NECM), which differs from the correlation matrix only by
$\mathcal{O}(T^{-1/2})$:

$\displaystyle\mathbf{C}=$
$\displaystyle\frac{1}{T}\widetilde{\mathbf{X}}^{T}\widetilde{\mathbf{X}}=\frac{1}{T}(\widetilde{\mathbf{U}\mathbf{V}+\sigma\mathbf{R}})^{T}(\widetilde{\mathbf{U}\mathbf{V}+\sigma\mathbf{R}})$
$\displaystyle=$
$\displaystyle\frac{1}{T}\Big{(}(\widetilde{\mathbf{U}\mathbf{V}})^{T}(\widetilde{\mathbf{U}\mathbf{V}})+\tilde{\sigma}^{2}\mathbf{R}^{T}\mathbf{R}$
$\displaystyle+\tilde{\sigma}(\widetilde{\mathbf{U}\mathbf{V}})^{T}\mathbf{R}+\tilde{\sigma}\mathbf{R}^{T}\widetilde{\mathbf{U}\mathbf{V}}\Big{)}\,.$
(84)

The NECM contains three different contributions: the
$(\mathbf{UV})^{T}(\mathbf{UV})$ from the pure latent feature signal,
$\mathbf{R}^{T}\mathbf{R}$ from pure noise, and two terms of the type
$(\mathbf{UV})^{T}\mathbf{R}$, which are cross terms between the latent signal
and the noise. Each contribution is an $N\times N$ random matrix. Critical to
computing the eigenvalue density of random matrices is the concept of _matrix
freeness_ [5], which is the generalization of statistical independence to
matrices. The eigenvalue spectrum of sums and products of free matrices can be
computed from spectra of summands and factors using the $\mathcal{R}$\- and
the $\mathcal{S}$-transforms, which are related to the Stieltjes transform
$\mathfrak{g}$ and are additive and multiplicative, respectively. The signal-
signal and the noise-noise contributions in the NECM definition are certainly
free w. r. t. each other. We will argue in _Appendix C.3_ that, in our regimes
of interest (the zero-noise limit ($\mathsf{SNR}\rightarrow\infty$), the
classical statistics limit from Eq. (9), and intensive limit from Eq. (10)),
the cross-term contributions are negligible, so that we can drop them and
approximate the NECM as

$\displaystyle\mathbf{C}$
$\displaystyle\approx\frac{(\mathbf{U}\mathbf{V})^{T}(\mathbf{U}\mathbf{V})+\sigma^{2}\mathbf{R}^{T}\mathbf{R}}{\sigma_{X}^{2}T}:=\mathbf{C}_{\widetilde{\mathbf{U}\mathbf{V}}}+\mathbf{C}_{\widetilde{\sigma}\mathbf{R}}\,,$
(85)

so that free matrix theory applies.

###  C.1 Parameterizing the random matrix problem and the large matrix limit

To calculate the spectrum of the signal-signal contribution to the NECM,

$\displaystyle\mathbf{C}_{\widetilde{\mathbf{UV}}}=\frac{1}{\sigma_{X}^{2}T}(\mathbf{U}\mathbf{V})^{T}(\mathbf{U}\mathbf{V})\,,$
(86)

we note that, assuming $m<T,N$, this $N\times N$ matrix is of rank $m$. Thus
we can work in the basis, where

$\displaystyle\mathbf{C}_{\widetilde{\mathbf{UV}}}=\left(\begin{array}[]{cc}\mathbf{H}_{\widetilde{\mathbf{UV}}}&0\\\
0&0\end{array}\right)\,,$ (89)

and

$\displaystyle\mathbf{H}_{\widetilde{\mathbf{UV}}}=\frac{1}{\sigma_{X}^{2}T}\left(\mathbf{U}^{T}\mathbf{U}\right)\left(\mathbf{V}\mathbf{V}^{T}\right)\,.$
(90)

There are $m$ non-trivial eigenvalues associated with $\mathbf{H}$, while the
remaining $N-m$ eigenvalues are zero. The finite size Stieltjes transform,
$g_{\mathbf{C}}^{N}=N^{-1}\text{Tr}(z\mathbf{I}-{\mathbf{C}}_{\widetilde{\mathbf{UV}}})^{-1}$,
is then of the form

$\displaystyle g^{N}_{\mathbf{C}_{\widetilde{\mathbf{UV}}}}(z)$
$\displaystyle=\frac{1}{N}\left(m\,\frac{1}{m}\sum_{\mu=1}^{m}\frac{1}{z-\lambda_{\mu}}+\frac{N-m}{z}\right)$
$\displaystyle=\frac{1}{N}\left(m\,h^{m}_{\mathbf{H}_{\widetilde{\mathbf{UV}}}}(z)+\frac{N-m}{z}\right)\,,$
(91)

where $\lambda_{\mu}$ are the $m$ eigenvalues of
$\mathbf{H}_{\widetilde{\mathbf{UV}}}$ and
$h^{m}_{\mathbf{H}_{\widetilde{\mathbf{UV}}}}(z)$ is its finite size Stieltjes
transform.

Now we note that $\mathbf{H}_{\widetilde{\mathbf{UV}}}$ in Eq. (90) is the
product of two white Wishart matrices

$\displaystyle\mathbf{H}_{\widetilde{\mathbf{UV}}}=\frac{N}{\sigma_{X}^{2}}\,\mathbf{W}_{\mathbf{U}}\mathbf{W}_{\mathbf{V}^{T}}\,,$
(92)

where

$\displaystyle\mathbf{W}_{\mathbf{Y}}=\frac{1}{T}\mathbf{Y}^{T}\mathbf{Y}\,,$
(93)

is the Wishart matrix, and $\mathbf{Y}$ is a $T\times N$ matrix with i.i.d.
standard normal entries. The key parameter characterizing such standard
$\mathbf{W}_{\mathbf{Y}}$ is the ratio of the number of columns to that of
rows

$\displaystyle q\equiv\frac{N}{T}\,.$ (94)

Since $\mathbf{U}$ and $\mathbf{V}^{T}$ are $T\times m$ and $N\times m$
matrices, respectively, a natural characterisation of
$\mathbf{H}_{\widetilde{\mathbf{U}\mathbf{V}}}$ is then

$\displaystyle q\equiv\frac{N}{T},\quad q_{U}\equiv\frac{m}{T},\quad
q_{V^{T}}\equiv\frac{m}{N}\,,$ (95)

with $q\,q_{V^{T}}=q_{U}$, so that there are only two independent parameters.

It is now convenient to define

$\displaystyle\sigma_{X}^{2}=m\left(\sigma_{U}^{2}\sigma_{V}^{2}+\frac{\sigma^{2}}{m}\right)\equiv
m\bar{\sigma}_{X}^{2}\,,$ (96)

where we used Eq. (38), so that Eq. (92) becomes

$\displaystyle\mathbf{H}_{\widetilde{\mathbf{UV}}}=\frac{1}{q_{V^{T}}\bar{\sigma}_{X}^{2}}\,\mathbf{W}_{\mathbf{U}}\mathbf{W}_{\mathbf{V}^{T}}\,.$
(97)

In the following, we only consider the limit of large matrices. Here $T$, $N$,
$m$ and $\sigma^{2}$ go to infinity in such a way that $q$, $q_{V^{T}}$ and
SNR are all constant. Then in the thermodynamic limit the finite size
Stieltjes transform in Eq. (C.1) becomes

$\displaystyle\mathfrak{g}_{\mathbf{C_{\widetilde{\mathbf{UV}}}}}$
$\displaystyle=q_{V^{T}}\mathfrak{h}+\frac{1-q_{V^{T}}}{z}\,,$ (98)

where $\mathfrak{g}_{\mathbf{C_{\widetilde{\mathbf{UV}}}}}$ and $\mathfrak{h}$
are the large matrices limits of the Stieltjes transforms of
$g^{N}_{\mathbf{C}_{\widetilde{\mathbf{UV}}}}$ and
$h^{m}_{\mathbf{H}_{\widetilde{\mathbf{UV}}}}$, respectively.

###  C.2 The spectrum of $\mathbf{C}_{\widetilde{\mathbf{UV}}}$

We now compute the eigenvalue density of
$\mathbf{C}_{\widetilde{\mathbf{UV}}}$. The first step is to compute the
Stieljtes transform $\mathfrak{h}$. From Eq. (97), it is clear that this
reduces to the problem of computing the eigenvalue spectrum of a product of
two Wishart matrices.

The spectrum of a product of two free matrices can be computed with the help
of the $\mathcal{S}$-transform, which is defined for a random matrix
$\mathbf{A}$ as

$\displaystyle\mathcal{S}_{\mathbf{A}}(t)=\frac{t+1}{t{\mathcal{T}}_{\mathbf{A}}^{-1}(t)}\,,$
(99)

where ${\mathcal{T}}_{\mathbf{A}}^{-1}(t)$ is the functional inverse of the
$\mathcal{T}$-transform $\mathcal{T}_{\mathbf{A}}(z)$. In turn, the
$\mathcal{T}$-transform is related to the Stieltjes transform of $\mathbf{A}$
through the relation

$\displaystyle\mathcal{T}_{\mathbf{A}}(z)=z\mathfrak{g}_{\mathbf{A}}(z)-1\,.$
(100)

Crucially, for free matrices $\mathbf{A}$ and $\mathbf{B}$, the
$\mathcal{S}$-transform is multiplicative

$\displaystyle\mathcal{S}_{\mathbf{A}\mathbf{B}}(t)=\mathcal{S}_{\mathbf{A}}(t)\mathcal{S}_{\mathbf{B}}(t)\,,$
(101)

and, furthermore, for a scalar $a$,

$\displaystyle\mathcal{S}_{a\mathbf{A}}(t)=a^{-1}\mathcal{S}_{\mathbf{A}}(t)\,.$
(102)

For the white Wishart matrix, Eq. (93), the $\mathcal{S}$-transform is known
to be [4]

$\displaystyle\mathcal{S}_{\mathbf{W_{Y}}}(t)=\frac{1}{1+qt}\,.$ (103)

Thus we only need to use the multiplicative property of the
$\mathcal{S}$-transform to compute the signal-signal contributions to the
NECM. Specifically,

$\displaystyle\mathcal{S}_{\mathbf{H}_{\widetilde{\mathbf{UV}}}}(t)$
$\displaystyle=q_{V^{T}}\bar{\sigma}_{X}^{2}\mathcal{S}_{\mathbf{W}_{\mathbf{U}}}\mathcal{S}_{\mathbf{W}_{\mathbf{V}^{T}}}$
$\displaystyle=\frac{q_{V^{T}}\bar{\sigma}_{X}^{2}}{(1+q_{U}t)(1+q_{V^{T}}t)}\,.$
(104)

Equation (99) then yields

$\displaystyle{\mathcal{T}}^{-1}_{\mathbf{H}_{\widetilde{\mathbf{UV}}}}(t)=\frac{t+1}{t\,\mathcal{S}_{\mathbf{H}_{\widetilde{\mathbf{UV}}}}(t)}=\frac{t+1}{t}\frac{(1+q_{U}t)(1+q_{V^{T}}t)}{q_{V^{T}}\bar{\sigma}_{X}^{2}}\,.$
(105)

We now solve the equation for the functional inverse,
${\mathcal{T}}^{-1}(\mathcal{T}(z))=z$, using the definition of the
$\mathcal{T}$-transform, Eq. (100), and dividing by a common factor of $z$. We
obtain a cubic equation for the Stieltjes transform $\mathfrak{h}$:

$\displaystyle\mathfrak{h}^{3}$ $\displaystyle
z^{2}q_{U}q_{V^{T}}+\mathfrak{h}^{2}z\left(q_{V^{T}}(1-q_{U})+q_{U}(1-q_{V^{T}})\right)$
$\displaystyle+\mathfrak{h}\left((1-q_{U})(1-q_{V^{T}})-zq_{V^{T}}\bar{\sigma}_{X}^{2}\right)+q_{V^{T}}\bar{\sigma}_{X}^{2}=0\,.$
(106)

Finally, we divide by $q_{V^{T}}\bar{\sigma}_{X}^{2}$ to obtain

$\displaystyle\mathfrak{h}^{3}\frac{z^{2}q_{U}}{\bar{\sigma}_{X}^{2}}$
$\displaystyle+\mathfrak{h}^{2}\frac{z}{\bar{\sigma}_{X}^{2}}\left(1+q-2q_{U}\right)$
$\displaystyle+\mathfrak{h}\left(\frac{q_{V^{T}}^{-1}-q-1+q_{U}}{\bar{\sigma}_{X}^{2}}-z\right)+1=0\,.$
(107)

Similar equations for the Stieljtes transform of the product of two random
matrices have been stated in [6, 7, 8]. Their polynomials differ from Eq.
(C.2) in details, because we consider the Stieljtes transform of the
covariance matrix including a theoretical normalisation factor.

The next step is to solve Eq. (C.2) analytically in the classical statistics
limit and the intensive limit. We remind the reader that, for the pure signal
contribution, we work in the zero noise limit $\textsf{SNR}\rightarrow\infty$,
such that

$\displaystyle\bar{\sigma}_{X}^{2}\equiv\frac{\sigma_{X}^{2}}{m}=\sigma_{U}^{2}\sigma_{V}^{2}(1+\textsf{SNR}^{-1})=\sigma_{U}^{2}\sigma_{V}^{2}\,.$
(108)

####  C.2.1 Classical statistics limit

In the classical statistics limit, Eq. (9), the polynomial equation for the
Stieltjes transform, Eq. (C.2), becomes:

$\displaystyle\mathfrak{h}^{2}\frac{z}{\bar{\sigma}_{X}^{2}}+\mathfrak{h}\left(\frac{q_{V^{T}}^{-1}-1}{\bar{\sigma}_{X}^{2}}-z\right)+1=0\,.$
(109)

The discriminant is

$\displaystyle\Delta=z^{2}-2\frac{1+q_{V^{T}}^{-1}}{\bar{\sigma}_{X}^{2}}z+\left(\frac{q_{V^{T}}^{-1}-1}{\bar{\sigma}_{X}^{2}}\right)^{2},$
(110)

and the roots of the discriminant are

$\displaystyle\lambda_{\pm}^{\infty}=\bar{\sigma}_{X}^{-2}\left(1\pm\sqrt{q_{V^{T}}}^{-1}\right)^{2}.$
(111)

We thus obtain

$\displaystyle\mathfrak{h}_{\pm}=\frac{-\frac{q_{V^{T}}^{-1}-1}{\bar{\sigma}_{X}^{2}}+z\pm\sqrt{(z-\lambda_{-}^{\infty})(z-\lambda_{+}^{\infty})}}{2z{\bar{\sigma}_{X}}^{-2}}\,.$
(112)

To obtain $\mathfrak{g}_{\mathbf{C}_{\widetilde{\mathbf{UV}}}}$, we now need
to add the contribution of the zero eigenvalues:

$\displaystyle\mathfrak{g}_{\mathbf{C}_{\widetilde{\mathbf{UV}}}}$
$\displaystyle=q_{V^{T}}\mathfrak{h}_{\pm}+\frac{1-q_{V^{T}}}{z}$
$\displaystyle=\frac{1-q_{V^{T}}}{2z}+\frac{q_{V^{T}}}{2{\bar{\sigma}_{X}}^{-2}}\pm\frac{\sqrt{(z-\lambda_{-}^{\infty})(z-\lambda_{+}^{\infty})}}{2zq_{V^{T}}^{-1}{\bar{\sigma}_{X}}^{-2}}\,.$
(113)

We are now ready to obtain the eigenvalue density, as in Eq. (83). While this
is a standard calculation [9], we summarise it here for the reader’s benefit.
The second term on the right-hand side of Eq. (C.2.1) is real, does not
contribute to the imaginary part, and we ignore it. For the first and the
third terms, we multiply the numerators and the denominators by
$z^{*}=\lambda+i\eta$. The imaginary part of the first term is then

$\displaystyle\Im\left(\frac{1-q_{V^{T}}}{2z}\right)=\frac{(1-q_{V^{T}})\eta}{2(\eta^{2}+\lambda^{2})}=\frac{(1-q_{V^{T}})\pi}{2}\delta_{\eta}(\lambda)\,,$
(114)

where we have used the definition of the Lorentz curve,
$\delta_{\eta}(\lambda)=\pi^{-1}\eta/(\eta^{2}+\lambda^{2})$. For the third
term, the crucial step is to rewrite the square root using the relation

$\displaystyle\sqrt{a+ib}=P+iQ\,,$ (115)

where $a$ and $b$ are real, $b\neq 0$ and

$\displaystyle P$
$\displaystyle=\frac{1}{\sqrt{2}}\sqrt{\sqrt{a^{2}+b^{2}}+a},$ $\displaystyle
Q$
$\displaystyle=\frac{\text{sgn}(b)}{\sqrt{2}}\sqrt{\sqrt{a^{2}+b^{2}}-a}\,,$
(116)

where $\mathrm{sgn}(x)=1$ for $x>0$ and $-1$ for $x<0$ [10]. For the argument
of the square root in the third term of Eq. (C.2.1), we find

$\displaystyle a$
$\displaystyle=\lambda^{2}-\eta^{2}+\lambda_{+}^{\infty}\lambda_{-}^{\infty}-(\lambda_{+}^{\infty}+\lambda_{-}^{\infty})\lambda,$
$\displaystyle b$
$\displaystyle=(-2\lambda+\lambda_{+}^{\infty}+\lambda_{-}^{\infty})\eta\,.$
(117)

The imaginary part of the third term takes the form

$\displaystyle\Im{\left(\pm\frac{\sqrt{(z-\lambda_{-}^{\infty})(z-\lambda_{+}^{\infty})}}{2zq_{V^{T}}^{-1}{\bar{\sigma}_{X}}^{-2}}\right)}=$
$\displaystyle=\pm\frac{\Im\left(z^{*}\left[P+iQ\right]\right)}{2q_{V^{T}}^{-1}{\bar{\sigma}_{X}}^{-2}|z|^{2}}$
$\displaystyle=\pm\frac{1}{2q_{V^{T}}^{-1}{\bar{\sigma}_{X}}^{-2}}\left(\frac{\eta}{\eta^{2}+\lambda^{2}}P+\frac{\lambda}{\eta^{2}+\lambda^{2}}Q\right)$
$\displaystyle=\pm\frac{1}{2q_{V^{T}}^{-1}{\bar{\sigma}_{X}}^{-2}}\left(\pi\delta_{\eta}(\lambda)P+\frac{\lambda}{\eta^{2}+\lambda^{2}}Q\right).$
(118)

The final step to evaluate Eq. (83) and to obtain the eigenvalue density, is
to take the limit $\eta\rightarrow 0^{+}$. In this limit, the Lorentz curve in
Eq. (114) converges to the Dirac $\delta$-function. Combining Eqs. (114,
C.2.1) yields

$\displaystyle\rho^{\infty}(\lambda)$
$\displaystyle=\frac{1}{\pi}\lim_{\eta\rightarrow
0^{+}}\Im\,\mathfrak{g}_{\mathbf{C}_{\widetilde{\mathbf{UV}}}}$
$\displaystyle=\pm\frac{\lim_{\eta\rightarrow
0^{+}}P}{2q_{V^{T}}^{-1}\bar{\sigma}_{X}^{-2}}\delta(\lambda)\pm\frac{\lim_{\eta\rightarrow
0^{+}}Q}{2\pi\lambda\bar{\sigma}_{X}^{-2}q_{V^{T}}^{-1}}+\frac{1-q_{V^{T}}}{2}\delta(\lambda)$
(119)

with

$\displaystyle\lim_{\eta\rightarrow
0^{+}}P=\sqrt{\lambda_{+}^{\infty}\lambda_{-}^{\infty}}=\bar{\sigma}_{X}^{-2}(1-q_{V^{T}}^{-1})\,,$
(120)

where we have used the expression for the zero noise eigenvalue bounds in Eq.
(111), and

$\displaystyle\lim_{\eta\rightarrow 0^{+}}Q$
$\displaystyle=\frac{\text{sgn}(b)}{\sqrt{2}}\sqrt{2|(\lambda-\lambda_{-}^{\infty})(\lambda_{+}^{\infty}-\lambda)|}$
$\displaystyle=\text{sgn}(b)\sqrt{(\lambda-\lambda_{-}^{\infty})(\lambda_{+}^{\infty}-\lambda)},$
(121)

when $\lambda\in[\lambda_{-}^{\infty},\lambda_{+}^{\infty}]$, and the
expression vanishes elsewhere. The $\pm$ signs in Eq. (C.2.1) are chosen such
as to obtain a physically meaningful eigenvalue density. Finally, we find the
following form of the eigenvalue density

$\displaystyle\rho^{\infty}(\lambda)$
$\displaystyle=\frac{\sqrt{(\lambda-\lambda_{-}^{\infty})(\lambda_{+}^{\infty}-\lambda)}}{2\pi\lambda{\bar{\sigma}_{X}}^{-2}q_{V^{T}}^{-1}}+(1-q_{V^{T}})\delta(\lambda)\,,$
(122)

with
$\bar{\sigma}_{X}^{2}\equiv\sigma_{X}^{2}/m=\sigma_{U}^{2}\sigma_{V}^{2}$. We
note that in [11] an expression for an eigenvalue density was given in the
special case when $T=N=m$ and not including our theoretical normalisation
factor.

####  C.2.2 Intensive limit

For the intensive limit, Eq. (10), the polynomial equation for the Stieltjes
transform, Eq. (C.2), becomes

$\displaystyle\mathfrak{h}^{2}\frac{z}{\bar{\sigma}_{X}^{2}}\left(1+q\right)+\mathfrak{h}\left(\frac{q_{V^{T}}^{-1}-q-1}{\bar{\sigma}_{X}^{2}}-z\right)+1=0\,.$
(123)

The discriminant is

$\displaystyle\Delta=z^{2}-2\frac{q_{V^{T}}^{-1}+q+1}{\bar{\sigma}_{X}^{2}}z+\left(\frac{q_{V^{T}}^{-1}-q-1}{\bar{\sigma}_{X}^{2}}\right)^{2}\,.$
(124)

The roots of the discriminant are

$\displaystyle\lambda_{\pm}^{\infty}=\bar{\sigma}_{X}^{-2}\left(\sqrt{1+q}\pm\sqrt{q_{V^{T}}}^{-1}\right)^{2}\,.$
(125)

Then the solution of the quadratic equation is

$\displaystyle\mathfrak{h}_{\pm}=\frac{-\frac{q_{V^{T}}^{-1}-q-1}{\bar{\sigma}_{X}^{2}}+z\pm\sqrt{(z-\lambda_{-}^{\infty})(z-\lambda_{+}^{\infty})}}{2z{\bar{\sigma}_{X}}^{-2}(1+q)}.$
(126)

Following a calculation analogous to the classical limit, we now add the
contribution of the zero eigenvalues and then determine the density of the
eigenvalues. We find:

$\displaystyle\rho^{\infty}(\lambda)$
$\displaystyle=\frac{\sqrt{(\lambda-\lambda_{-}^{\infty})(\lambda_{+}^{\infty}-\lambda)}}{2\pi\lambda{\bar{\sigma}_{X}}^{-2}(1+q)q_{V^{T}}^{-1}}+(1-q_{V^{T}})\delta(\lambda)\,,$
(127)

with
$\bar{\sigma}_{X}^{2}\equiv\sigma_{X}^{2}/m=\sigma_{U}^{2}\sigma_{V}^{2}$.

###  C.3 Approximation to neglect the signal-noise cross terms

Now we explore when the contribution of the signal-noise cross terms to the
NECM can be neglected. Specifically, we will show that it can be done if
$q_{U}\rightarrow 0$ (that is, the number of measurements is much larger than
the number of latent features), which we always assume. To show this, we
compute the eigenvalue bounds, $\lambda_{\pm}^{\text{signal-noise}}$, of the
signal-noise contribution and compare their scaling with $T$ to the scaling of
the pure signal and the pure noise eigenvalue bounds.

For the pure signal contribution, the previous section shows that the
eigenvalue bounds $\lambda^{\infty}_{\pm}$ are
$\bar{\sigma}_{X}^{-2}\sim\mathcal{O}(T^{0})$. The pure noise eigenvalue
bounds, given by the Marčenko-Pastur bounds, scale as

$\displaystyle\lambda_{\pm}^{\mathrm{MP}}\sim 1\pm T^{-1/2}\,,$ (128)

where $1$ is due to self correlations. On the other hand, the signal-noise
cross terms do not have self-correlations, and thus we expect their bounds to
scale as

$\displaystyle\lambda^{\mathrm{signal-noise}}_{\pm}\sim T^{-1/2}\,,$ (129)

becoming negligible for $T\to\infty$. In Appendix C.3.1, we show this
analytically in the classical statistics limit. We have not been able to
achieve similar results more generally. However, since $q_{U}\rightarrow 0$
also in the intensive limit, we expect similar results to hold there too. To
show this, we resort to numerical simulations.

Figure S4: Difference between the eigenvalue density of the NECM spectrum with
and without the signal-noise cross terms quantified by $d^{\prime}$. (a)
Classical limit (orange), intensive limit (green) and neither of the two
limits (blue). (b) Magnified view of $d^{\prime}$ in the limits of interest:
classical limit (blue, orange and green) and intensive limit (red and purple).
Eigenvalue densities are computed from $120$ realisations of the random matrix
model.

Specifically, we numerically estimate the Jensen-Shannon divergence between
the numerically evaluated eigenvalue densities of the NECM, computed with and
without the signal-noise cross terms. To obtain a perceptually intuitive
measure of the difference between these distributions, we convert the Jensen-
Shannon divergence to the effective sensitivity index $d^{\prime}$ – the
distance between the means of two unit variance normal distribution with the
same Jensen-Shannon divergence as the two analyzed eigenvalue spectra. We
investigate the dependence of $d^{\prime}$ on various choices of our model
parameters. The comparison between the spectra of the full and the approximate
NECM is shown in Fig. S4. We observe that the sensitivity index reaches a
maximum for $\mathsf{SNR}\sim 10^{-1}$, and falls off in the limits of small
or large $\mathsf{SNR}$, where the noise or the signal dominate, respectively.
Crucially, the maximum value of $d^{\prime}$ is small for $q_{U}\to 0$. Thus
neglecting the cross-term contributions to the NECM spectrum in our limits of
interest is warranted.

####  C.3.1 Scaling behavior of the signal-noise eigenvalue bounds in the
classical limit

We now derive the scaling of the signal-noise eigenvalue spectrum bounds in
Eq. (129) in the classical statistics limit. From the NECM in Eq. (C), the
signal-noise cross terms are of the form

$\displaystyle\mathbf{M}\equiv\frac{\tilde{\sigma}}{T}(\widetilde{\mathbf{U}\mathbf{V}})^{T}{\mathbf{R}}\,.$
(130)

To compute the spectrum of this matrix, we use the following trick. The
singular values of $\mathbf{M}$ are equal to the square roots of the non-zero
eigenvalues of its square,

$\displaystyle\mathbf{M}^{2}$
$\displaystyle=\mathbf{M}\mathbf{M}^{T}=\frac{\tilde{\sigma}^{2}}{T^{2}}(\widetilde{\mathbf{U}\mathbf{V}})^{T}{\mathbf{R}}{\mathbf{R}}^{T}\widetilde{\mathbf{U}\mathbf{V}}\,.$
(131)

In turn, the non-zero eigenvalues of this $N\times N$ matrix, are equal to the
eigenvalues of the $T\times T$ matrix

$\displaystyle\widehat{\mathbf{M}}^{2}\equiv
q^{2}\tilde{\sigma}^{2}\frac{{\mathbf{R}}{\mathbf{R}^{T}}}{N}\frac{(\widetilde{\mathbf{U}\mathbf{V}})(\widetilde{\mathbf{U}\mathbf{V}})^{T}}{N}\,.$
(132)

We note that the r.h.s. of the above equation is a product of the $T\times T$
dual correlation matrices $\mathbf{C}_{\tilde{\sigma}\mathbf{R}^{T}}$ and
$\mathbf{C}_{(\widetilde{\mathbf{UV}})^{T}}$. To compute the spectrum of the
product, we employ the $\mathcal{S}$-transform formalism, as explained above.
The first step is to obtain the $\mathcal{S}$-transforms of the dual
correlation matrices, which we compute from the Stieltjes transform [4]. For
the noise part, we have

$\displaystyle\mathfrak{g}_{\mathbf{C}_{\mathbf{R}^{T}}}(z)=q^{2}\mathfrak{g}_{\mathbf{C}_{\mathbf{R}}}(qz)+\frac{1-q}{z}.$
(133)

For the signal part, we have

$\displaystyle\mathfrak{g}_{\mathbf{C}_{(\widetilde{\mathbf{UV}})^{T}}}(z)=q^{2}\mathfrak{g}_{\mathbf{C}_{\widetilde{\mathbf{UV}}}}(qz)+\frac{1-q}{z}\,.$
(134)

For the noise Wishart matrix, the Stieljtes transform is (cf. Eq. (103)):

$\displaystyle\mathcal{S}_{\mathbf{C}_{{{\mathbf{R}}}^{T}}}=\frac{1}{1+q^{-1}t}\,.$
(135)

Including the renormalized noise strength, $\tilde{\sigma}$, and the
additional factor of $q$ from Eq. (132), by using the scaling relation Eq.
(102), the $\mathcal{S}$-transform is

$\displaystyle\mathcal{S}_{\mathbf{C}_{{{q\tilde{\sigma}\mathbf{R}}}^{T}}}=\frac{q^{-2}\tilde{\sigma}^{-2}}{1+q^{-1}t}\,.$
(136)

Next, we write the $\mathcal{S}$-transform of the pure signal part. Evaluating
Eq. (134) using Eq. (C.2.1), we find

$\mathfrak{g}_{\mathbf{C}_{(\widetilde{\mathbf{UV}})^{T}}}=\\\
\frac{\frac{2-q-qq_{V^{T}}}{\bar{\sigma}_{X}^{2}}+q^{2}q_{V^{T}}z\pm
qq_{V^{T}}\sqrt{(qz-\lambda^{\infty}_{+})(qz-\lambda^{\infty}_{-})}}{2z\bar{\sigma}_{X}^{-2}}\,,$
(137)

from which we obtain the following equation for the Stieltjes transform:

$\displaystyle\left(2z\bar{\sigma}_{X}^{-2}\mathfrak{g}_{\mathbf{C}_{(\widetilde{\mathbf{UV}})^{T}}}-\frac{2-q-qq_{V^{T}}}{\bar{\sigma}_{X}^{2}}-q^{2}q_{V^{T}}z\right)^{2}$
$\displaystyle-q^{2}q_{V^{T}}^{2}(qz-\lambda_{+}^{\infty})(qz-\lambda_{-}^{\infty})=0\,.$
(138)

Using the relation, $\mathcal{T}=z\mathfrak{g}-1$, we find the equation for
the $\mathcal{T}$-transform

$\mathcal{T}^{2}_{\mathbf{C}_{(\widetilde{\mathbf{UV}})^{T}}}+\\\
\mathcal{T}_{\mathbf{C}_{(\widetilde{\mathbf{UV}})^{T}}}(-zq^{2}q_{V^{T}}\bar{\sigma}_{X}^{2}+qq_{V^{T}}+q)+q^{2}q_{V^{T}}=0\,.$
(139)

Interpreted as an equation for the functional inverse transform
$\mathcal{T}^{-1}$, this becomes

$\displaystyle
t^{2}+t(-\mathcal{T}^{-1}_{\mathbf{C}_{(\widetilde{\mathbf{UV}})^{T}}}q^{2}q_{V^{T}}\bar{\sigma}_{X}^{2}+qq_{V^{T}}+q)+q^{2}q_{V^{T}}=0\,.$
(140)

Now solving for the functional inverse transform, we find

$\displaystyle\mathcal{T}^{-1}_{\mathbf{C}_{(\widetilde{\mathbf{UV}})^{T}}}=\frac{(t+q)(t+q_{U})}{tqq_{U}\bar{\sigma}_{X}^{2}},$
(141)

from which we determine the $\mathcal{S}$-transform

$\displaystyle\mathcal{S}_{\mathbf{C}_{(\widetilde{\mathbf{UV}})^{T}}}=\frac{t+1}{t\mathcal{T}^{-1}_{\mathbf{C}_{(\widetilde{\mathbf{UV}})^{T}}}}=\frac{qq_{U}\bar{\sigma}_{X}^{2}(t+1)}{(t+q)(t+q_{U})}\,.$
(142)

The $\mathcal{S}$-transform of the product now reads

$\displaystyle\mathcal{S}_{\widehat{\mathbf{M}}^{2}}$
$\displaystyle=\mathcal{S}_{\mathbf{C}_{(\widetilde{\mathbf{UV}})^{T}}}\mathcal{S}_{\mathbf{C}_{{q\tilde{\sigma}{\mathbf{R}}}^{T}}}$
$\displaystyle=\frac{t+1}{t}\frac{tqq_{U}\bar{\sigma}_{X}^{2}}{(t+q)(t+q_{U})}\frac{q^{-2}\tilde{\sigma}^{-2}}{1+q^{-1}t}.$
(143)

From this we read off the inverse transform

$\displaystyle\mathcal{T}^{-1}_{\widehat{\mathbf{M}}^{2}}$
$\displaystyle=\frac{(t+q)(t+q_{U})(1+q^{-1}t)}{tq^{-1}q_{U}\bar{\sigma}_{X}^{2}\tilde{\sigma}^{-2}}$
$\displaystyle=\frac{(t+q)(t+q_{U})(1+q^{-1}t)}{tq^{-1}q_{U}\mathsf{SNR}^{-1}}.$
(144)

The equation for the $\mathcal{T}$-transform is now

$\displaystyle(\mathcal{T}_{\widehat{\mathbf{M}}^{2}}+q)(\mathcal{T}_{\widehat{\mathbf{M}}^{2}}+q_{U})(1+q^{-1}\mathcal{T}_{\widehat{\mathbf{M}}^{2}})=\dfrac{zq_{U}\mathcal{T}_{\widehat{\mathbf{M}}^{2}}}{q\,\mathsf{SNR}}.$
(145)

Using $\mathcal{T}=z\mathfrak{g}-1$ we write down the cubic polynomial
equation for $\mathfrak{g}$:

$\displaystyle
a\mathfrak{g}_{\widehat{\mathbf{M}}^{2}}^{3}+b\mathfrak{g}_{\widehat{\mathbf{M}}^{2}}^{2}+c\mathfrak{g}_{\widehat{\mathbf{M}}^{2}}+d=0\,,$
(146)

with coefficients

$\displaystyle a$ $\displaystyle=z^{3},$ (147) $\displaystyle b$
$\displaystyle=2qz^{2}+q_{U}z^{2}-3z^{2},$ (148) $\displaystyle c$
$\displaystyle=q^{2}z+2qq_{U}z-4qz-2q_{U}z+3z-\frac{q_{U}z^{2}}{\textsf{SNR}},$
(149) $\displaystyle d$
$\displaystyle=q^{2}q_{U}-q^{2}-2qq_{U}+2q+q_{U}-1+\frac{q_{U}z}{\textsf{SNR}}.$
(150)

The eigenvalue density is non-zero for complex solutions of the equation. The
equation admits complex solutions when the discriminant $\Delta$ is negative:

$\displaystyle\Delta=4P^{3}+27Q^{2}\,,$ (151)

where

$\displaystyle P$ $\displaystyle=\frac{3ac-b^{2}}{3a^{2}},$ (152)
$\displaystyle Q$ $\displaystyle=\frac{2b^{2}-9abc+27a^{2}d}{27a^{3}}.$ (153)

Written out explicitly, the determinant takes the form

$\displaystyle\Delta=\frac{4(3ac-b^{2})^{3}+(2b^{3}-9abc+27a^{2}d)^{2}}{27a^{6}}\,.$
(154)

The equation $\Delta=0$ yields a quadratic equation in $z$, giving the bounds
on the eigenvalue density of $\widehat{\mathbf{M}}^{2}$

$\displaystyle
z_{\pm}=\frac{8q^{2}+20qq_{U}-q_{U}^{2}\pm\sqrt{q_{U}(8q+q_{U})^{3}}}{8q_{U}\textsf{SNR}^{-1}}.$
(155)

From the definitions of $q=N/T$ and $q_{U}=m/T$, we see that these bounds
scale as

$\displaystyle z_{\pm}\sim T^{-1}\,.$ (156)

Since the singular values of $\mathbf{M}$ are equal to the square root of the
eigenvalues of $\widehat{\mathbf{M}}^{2}$, the eigenvalue bounds of the
signal-noise cross terms thus scale as

$\displaystyle\lambda^{\mathrm{signal-noise}}_{\pm}\sim T^{-1/2}\,.$ (157)

Thus the contribution of the cross-terms can be neglected in the classical
statistics limit.

###  C.4 Adding the noise contribution
$\mathbf{C}_{\mathbf{\tilde{\sigma}\mathbf{R}}}$

In Appendix C.2 we computed the spectrum of the pure signal contribution
$\mathbf{C}_{\widetilde{\mathbf{U}\mathbf{V}}}$ to the NECM in the classical
statistics and the intensive limits. Now we will add the pure noise
contribution $\mathbf{C}_{\tilde{\sigma}\mathbf{R}}$ to obtain the spectrum of
the approximate NECM. Since the noise and the signal contributions are free
matrices with respect to each other, the spectrum of their sum can be computed
using the $\mathcal{R}$-transform. The $\mathcal{R}$-transform of a random
matrix $\mathbf{A}$ is

$\displaystyle\mathcal{R}_{\mathbf{A}}(z)=\mathcal{B}_{\mathbf{A}}(z)-1/z\,,$
(158)

where the $\mathcal{B}$-transform is the functional inverse of the Stieltjes
transform

$\displaystyle\mathcal{B}_{\mathbf{A}}[\mathfrak{g}_{\mathbf{A}}]=z.$ (159)

The $\mathcal{R}$-transform is additive for free matrices:

$\displaystyle\mathcal{R}_{\mathbf{A}+\mathbf{B}}(z)=\mathcal{R}_{\mathbf{A}}(z)+\mathcal{R}_{\mathbf{B}}(z)\,.$
(160)

It scales according to

$\displaystyle\mathcal{R}_{a\mathbf{A}}(z)=a\mathcal{R}_{\mathbf{A}}(az)\,,$
(161)

where $a$ is a real number. For a white Wishart matrix, Eq. (93), the
$\mathcal{R}$-transform is known to be [4]

$\displaystyle\mathcal{R}_{\mathbf{W}_{\mathbf{Y}}}(z)=\frac{1}{1-qz}\,.$
(162)

For the pure noise contribution to the NECM,
$\mathbf{C}_{\tilde{\sigma}\mathbf{R}}\equiv\tilde{\sigma}^{2}\mathbf{R}^{T}\mathbf{R}/T$,
this results in

$\displaystyle\mathcal{R}_{\mathbf{C}_{\widetilde{\sigma}\mathbf{R}}}(z)=\frac{\widetilde{\sigma}^{2}}{1-qz\widetilde{\sigma}^{2}}\,.$
(163)

Our goal is to first obtain the $\mathcal{R}$-transform of the sum of the
signal and the noise contributions

$\displaystyle\mathcal{R}_{\mathbf{C}}(z)=\mathcal{R}_{\mathbf{C}_{\widetilde{\mathbf{U}\mathbf{V}}}}(z)+\mathcal{R}_{\mathbf{C}_{\widetilde{\sigma}\mathbf{R}}}(z),$
(164)

and from this to compute the Stieltjes transform to extract the eigenvalue
density. Computing the $\mathcal{R}$-transform of the pure signal contribution
$\mathbf{C}_{\widetilde{\mathbf{UV}}}$ in the classical statistics and the
intensive limit requires additional steps.

####  C.4.1 Classical statistics limit

First, we compute the $\mathcal{R}$-transform of
$\mathfrak{g}_{\mathbf{C}_{\widetilde{\mathbf{UV}}}}$. This is done by solving
the functional inverse equation,
$\mathfrak{g}_{\mathbf{C}_{\widetilde{\mathbf{UV}}}}[\mathcal{B}_{\mathbf{C}_{\widetilde{\mathbf{U}\mathbf{V}}}}]=z$,
which gives us the $\mathcal{B}$-transform, from which we compute the
$\mathcal{R}$-transform. From Eq. (C.2.1), we see that $\mathcal{B}$-transform
satisfies

$\displaystyle\mathcal{B}_{\mathbf{C}_{\widetilde{\mathbf{U}\mathbf{V}}}}\left(\mathcal{B}_{\mathbf{C}_{\widetilde{\mathbf{U}\mathbf{V}}}}z(-q_{V^{T}}\bar{\sigma}_{X}^{2}+z)+q_{V^{T}}\bar{\sigma}_{X}^{2}+q_{V^{T}}z-z\right)=0.$
(165)

A non-trivial solution of this equation is

$\displaystyle\mathcal{B}_{\mathbf{C}_{\widetilde{\mathbf{U}\mathbf{V}}}}=\frac{q_{V^{T}}\bar{\sigma}_{X}^{2}+q_{V^{T}}z-z}{z(q_{V^{T}}\bar{\sigma}_{X}^{2}-z)}\,.$
(166)

Using Eq. (158), this gives us the $\mathcal{R}$-transform,
$\mathcal{R}_{\mathbf{C}_{\widetilde{\mathbf{U}\mathbf{V}}}}$. To it we add
the $\mathcal{R}$-transform of the noise, Eq. (163), and subtract $-1/z$ to
get the $\mathcal{B}$-transform of the approximate NECM:

$\mathcal{B}_{\mathbf{C}}=\\\
\frac{\widetilde{\sigma}^{2}z(q_{V^{T}}\bar{\sigma}_{X}^{2}-z)+(-q\widetilde{\sigma}^{2}z+1)(q_{V^{T}}\bar{\sigma}_{X}^{2}+q_{V^{T}}z-z)}{z(q_{V^{T}}\bar{\sigma}_{X}^{2}-z)(-q\widetilde{\sigma}^{2}z+1)}\,.$
(167)

The final step is to write down and solve the inverse function equation
$\mathcal{B}_{\mathbf{C}}[\mathfrak{g}_{\mathbf{C}}]=z$. This is now
equivalent to solving the third order polynomial equation

$\displaystyle
a\mathfrak{g}_{\mathbf{C}}^{3}+b\mathfrak{g}_{\mathbf{C}}^{2}+c\mathfrak{g}_{\mathbf{C}}+d=0\,,\;\mbox{with}$
(168) $\displaystyle a=qz\widetilde{\sigma}^{2}\,,$ (169) $\displaystyle
b=-qq_{V^{T}}z\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{2}+\big{(}(q_{V^{T}}-1)q+1\big{)}\widetilde{\sigma}^{2}-z\,,$
(170) $\displaystyle
c=(q-1)q_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{2}+q_{V^{T}}z\bar{\sigma}_{X}^{2}-q_{V^{T}}+1\,,$
(171) $\displaystyle d=-q_{V^{T}}\bar{\sigma}_{X}^{2}\,.$ (172)

Written in terms of the signal-to-noise ratio, SNR, the coefficients take the
form

$\displaystyle a$ $\displaystyle=\frac{qz}{1+\textsf{SNR}}\,,$ (173)
$\displaystyle b$
$\displaystyle=-\frac{qq_{V^{T}}z}{\textsf{SNR}}+\frac{(q_{V^{T}}-1)q+1}{1+\textsf{SNR}}-z\,,$
(174) $\displaystyle c$
$\displaystyle=\frac{(q-1)q_{V^{T}}}{\textsf{SNR}}+q_{V^{T}}z\big{(}1+\textsf{SNR}^{-1}\big{)}-q_{V^{T}}+1\,,$
(175) $\displaystyle d$
$\displaystyle=-q_{V^{T}}\big{(}1+\textsf{SNR}^{-1}\big{)}\,.$ (176)

It is possible to solve this cubic equation analytically. However, the
expressions become lengthy and provide little insight. Therefore, we rely on
the numerical solution of the equation as shown in Fig. 2, as well as on the
following analyses in the limits of small and large noise.

First, in the limit of the pure signal, $\textsf{SNR}\rightarrow\infty$, we
recover Eq. (20). Similarly, by truncating the polynomial coefficients in the
pure noise limit, $\textsf{SNR}\rightarrow 0$, at order
$\mathcal{O}(\textsf{SNR}^{-1})$, the MP density is recovered.

We also derive an approximate analytic expression for the bounds of the
eigenvalue density, $\lambda_{\pm}^{\textsf{SNR}}$, which is valid around both
noise limits, $\textsf{SNR}\rightarrow 0$ and $\textsf{SNR}\rightarrow\infty$.
For this, we approximate the coefficients of the polynomial, noting that the
smallest contribution to the coefficients common to both limits comes from
terms of order $\mathcal{O}(q/(1+\textsf{SNR}))$ (recall that $q\rightarrow 0$
in the classical limit). Neglecting these terms leads to a quadratic
polynomial equation for the Stieltjes transform

$\displaystyle
r\mathfrak{g}_{\mathbf{C}}^{2}+s\mathfrak{g}_{\mathbf{C}}+t\approx
0\,,\;\mbox{with}$ (177) $\displaystyle
r=-\frac{qq_{V^{T}}z}{\textsf{SNR}}+\frac{1}{1+\textsf{SNR}}-z\,,$ (178)
$\displaystyle
s=\frac{(q-1)q_{V^{T}}}{\textsf{SNR}}+q_{V^{T}}z\big{(}1+\textsf{SNR}^{-1}\big{)}-q_{V^{T}}+1\,,$
(179) $\displaystyle t=-q_{V^{T}}\big{(}1+\textsf{SNR}^{-1}\big{)}\,.$ (180)

The approximate bounds of the eigenvalue density then are given by the roots
of the discriminant, $\Delta\mathfrak{g}_{\mathbf{C}}\approx s^{2}-4rt$, which
gives the following bounds for the nonzero range of the eigenvalue density:

$\displaystyle\lambda_{\pm}^{\textsf{SNR}}$
$\displaystyle\approx\frac{1+q_{V^{T}}^{-1}}{1+\textsf{SNR}^{-1}}+\frac{1+q}{1+\textsf{SNR}}\pm
2\sqrt{\frac{q_{V^{T}}^{-1}}{\left(1+\textsf{SNR}^{-1}\right)^{2}}+\frac{q}{\left(1+\textsf{SNR}\right)^{2}}+\frac{q}{\left(\sqrt{\textsf{SNR}}+\sqrt{\textsf{SNR}}^{-1}\right)^{2}}}$
$\displaystyle\approx\frac{1+q_{V^{T}}^{-1}}{1+\textsf{SNR}^{-1}}+\frac{1+q}{1+\textsf{SNR}}\pm
2\sqrt{\frac{q_{V^{T}}^{-1}}{\left(1+\textsf{SNR}^{-1}\right)^{2}}+\frac{q}{\left(1+\textsf{SNR}\right)^{2}}}$
$\displaystyle=\frac{1+q_{V^{T}}^{-1}}{1+\textsf{SNR}^{-1}}+\frac{1+q}{1+\textsf{SNR}}\pm
2\sqrt{\left[\sqrt{\frac{q_{V^{T}}^{-1}}{\left(1+\textsf{SNR}^{-1}\right)^{2}}}+\sqrt{\frac{q}{\left(1+\textsf{SNR}\right)^{2}}}\right]^{2}-\frac{2\sqrt{qq_{V^{T}}^{-1}}}{(1+\textsf{SNR}^{-1})(1+\textsf{SNR})}}$
$\displaystyle\approx\frac{1+q_{V^{T}}^{-1}}{1+\textsf{SNR}^{-1}}+\frac{1+q}{1+\textsf{SNR}}\pm
2\left[\sqrt{\frac{q_{V^{T}}^{-1}}{\left(1+\textsf{SNR}^{-1}\right)^{2}}}+\sqrt{\frac{q}{\left(1+\textsf{SNR}\right)^{2}}}\right]$
$\displaystyle=\frac{1}{1+\textsf{SNR}^{-1}}\lambda_{\pm}^{\infty}+\frac{1}{1+\textsf{SNR}}\lambda_{\pm}^{\text{MP}}\,.$
(181)

In the second line, we drop the third term under the square root, since it is
small in either of the two noise limits. In the third line, we have used
$(\sqrt{a}+\sqrt{b})^{2}=a+b+2\sqrt{ab}$, and in the fourth line, we dropped
the last term under the square root since it is also small in either of the
two noise limits. In the final line we recognize that the terms form the
weighted average of $\lambda^{\infty}_{\pm}$, Eq. (21), and the Marčenko-
Pastur bounds $\lambda_{\pm}^{\text{MP}}$.

####  C.4.2 Intensive limit

First we compute the $\mathcal{R}$-transform of
$\mathfrak{g}_{\mathbf{C}_{\widetilde{\mathbf{UV}}}}$. This is obtained by
solving the functional inverse equation
$\mathfrak{g}_{\mathbf{C}_{\widetilde{\mathbf{UV}}}}[\mathcal{B}_{\mathbf{C}_{\widetilde{\mathbf{UV}}}}]=z$,
which gives us the $\mathcal{B}$-transform, from which we compute the
$\mathcal{R}$-transform. To do this, we employed the symbolic algebra Python
library SymPy v1.6.2. The $\mathcal{B}$-transform satisfies the quadratic
equation

$\mathcal{B}_{\mathbf{C}_{\widetilde{\mathbf{UV}}}}^{2}z(q^{2}z-qq_{V^{T}}\bar{\sigma}_{X}^{2}+2qz-
q_{V^{T}}\bar{\sigma}_{X}^{2}+z)\\\
+\mathcal{B}_{\mathbf{C}_{\widetilde{\mathbf{UV}}}}(q^{2}q_{V^{T}}z-2q^{2}z+qq_{V^{T}}\bar{\sigma}_{X}^{2}+2qq_{V^{T}}z-3qz+q_{V^{T}}\bar{\sigma}_{X}^{2}+q_{V^{T}}z-z)-q^{2}q_{V^{T}}+q^{2}-qq_{V^{T}}+q=0\,,$
(182)

for which we find the solution

$\displaystyle\mathcal{B}_{\mathbf{C}_{\widetilde{\mathbf{UV}}}}$
$\displaystyle=\Big{(}-qq_{V^{T}}z+2qz-
q_{V^{T}}\bar{\sigma}_{X}^{2}-q_{V^{T}}z+z$
$\displaystyle-\sqrt{q^{2}q_{V^{T}}^{2}z^{2}-2qq_{V^{T}}^{2}\bar{\sigma}_{X}^{2}z+2qq_{V^{T}}^{2}z^{2}-2qq_{V^{T}}z^{2}+q_{V^{T}}^{2}\bar{\sigma}_{X}^{4}+2q_{V^{T}}^{2}\bar{\sigma}_{X}^{2}z+q_{V^{T}}^{2}z^{2}-2q_{V^{T}}\bar{\sigma}_{X}^{2}z-2q_{V^{T}}z^{2}+z^{2}}\Big{)}$
$\displaystyle\Big{/}\big{(}2z(qz-q_{V^{T}}\bar{\sigma}_{X}^{2}+z)\big{)}\,.$
(183)

Using Eq. (158), this gives us
$\mathcal{R}_{\mathbf{C}_{\widetilde{{\mathbf{U}\mathbf{V}}}}}$, to which we
add the $\mathcal{R}$-transform of the noise, Eq. (163), to obtain the
$\mathcal{R}$-transform, $\mathcal{R}_{\mathbf{C}}$, of the NECM. Subtracting,
$-1/z$, gives us the corresponding form of the $\mathcal{B}$-transform:

$\displaystyle\mathcal{B}_{\mathbf{C}}$
$\displaystyle=\Big{(}2\widetilde{\sigma}^{2}z(qz-
q_{V^{T}}\bar{\sigma}_{X}^{2}+z)+(-q\widetilde{\sigma}^{2}z+1)(-qq_{V^{T}}z+2qz-
q_{V^{T}}\bar{\sigma}_{X}^{2}-q_{V^{T}}z+z$
$\displaystyle-\sqrt{q^{2}q_{V^{T}}^{2}z^{2}-2qq_{V^{T}}^{2}\bar{\sigma}_{X}^{2}z+2qq_{V^{T}}^{2}z^{2}-2qq_{V^{T}}z^{2}+q_{V^{T}}^{2}\bar{\sigma}_{X}^{4}+2q_{V^{T}}^{2}\bar{\sigma}_{X}^{2}z+q_{V^{T}}^{2}z^{2}-2q_{V^{T}}\bar{\sigma}_{X}^{2}z-2q_{V^{T}}z^{2}+z^{2}})\Big{)}$
$\displaystyle\Big{/}\big{(}2z(-q\widetilde{\sigma}^{2}z+1)(qz-
q_{V^{T}}\bar{\sigma}_{X}^{2}+z)\big{)}\,.$ (184)

The final step is to write down and solve the inverse functional equation
$\mathcal{B}_{\mathbf{C}}[\mathfrak{g}_{\mathbf{C}}]=z$. The sixth order
polynomial equation that we need to solve is of the form

$\displaystyle
a\mathfrak{g}_{\mathbf{C}}^{6}+b\mathfrak{g}_{\mathbf{C}}^{5}+c\mathfrak{g}_{\mathbf{C}}^{4}+d\mathfrak{g}_{\mathbf{C}}^{3}+e\mathfrak{g}_{\mathbf{C}}^{2}+f\mathfrak{g}_{\mathbf{C}}+g=0\,,$
(185)

with coefficients

$\displaystyle a=$ $\displaystyle
q^{5}\widetilde{\sigma}^{6}z^{2}+2q^{4}\widetilde{\sigma}^{6}z^{2}+q^{3}\widetilde{\sigma}^{6}z^{2}\,,$
(186) $\displaystyle b=$ $\displaystyle
q^{5}q_{V^{T}}\widetilde{\sigma}^{6}z-2q^{5}\widetilde{\sigma}^{6}z-2q^{4}q_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{6}z^{2}+2q^{4}q_{V^{T}}\widetilde{\sigma}^{6}z-q^{4}\widetilde{\sigma}^{6}z-3q^{4}\widetilde{\sigma}^{4}z^{2}-2q^{3}q_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{6}z^{2}+q^{3}q_{V^{T}}\widetilde{\sigma}^{6}z$
$\displaystyle+3q^{3}\widetilde{\sigma}^{6}z-6q^{3}\widetilde{\sigma}^{4}z^{2}+2q^{2}\widetilde{\sigma}^{6}z-3q^{2}\widetilde{\sigma}^{4}z^{2}\,,$
(187) $\displaystyle c=$
$\displaystyle-q^{5}q_{V^{T}}\widetilde{\sigma}^{6}+q^{5}\widetilde{\sigma}^{6}-q^{4}q_{V^{T}}^{2}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{6}z+3q^{4}q_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{6}z-3q^{4}q_{V^{T}}\widetilde{\sigma}^{4}z-q^{4}\widetilde{\sigma}^{6}+6q^{4}\widetilde{\sigma}^{4}z+q^{3}q_{V^{T}}^{2}\bar{\sigma}_{X}^{4}\widetilde{\sigma}^{6}z^{2}$
$\displaystyle-q^{3}q_{V^{T}}^{2}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{6}z-2q^{3}q_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{6}z+6q^{3}q_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{4}z^{2}+2q^{3}q_{V^{T}}\widetilde{\sigma}^{6}-6q^{3}q_{V^{T}}\widetilde{\sigma}^{4}z-2q^{3}\widetilde{\sigma}^{6}+5q^{3}\widetilde{\sigma}^{4}z+3q^{3}\widetilde{\sigma}^{2}z^{2}$
$\displaystyle-4q^{2}q_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{6}z+6q^{2}q_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{4}z^{2}+q^{2}q_{V^{T}}\widetilde{\sigma}^{6}-3q^{2}q_{V^{T}}\widetilde{\sigma}^{4}z+q^{2}\widetilde{\sigma}^{6}-5q^{2}\widetilde{\sigma}^{4}z+6q^{2}\widetilde{\sigma}^{2}z^{2}+q\widetilde{\sigma}^{6}-4q\widetilde{\sigma}^{4}z+3q\widetilde{\sigma}^{2}z^{2}\,,$
(188) $\displaystyle d=$ $\displaystyle
q^{4}q_{V^{T}}^{2}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{6}-q^{4}q_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{6}+3q^{4}q_{V^{T}}\widetilde{\sigma}^{4}-3q^{4}\widetilde{\sigma}^{4}-q^{3}q_{V^{T}}^{2}\bar{\sigma}_{X}^{4}\widetilde{\sigma}^{6}z-q^{3}q_{V^{T}}^{2}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{6}$
$\displaystyle+3q^{3}q_{V^{T}}^{2}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{4}z+3q^{3}q_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{6}-9q^{3}q_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{4}z+q^{3}q_{V^{T}}\widetilde{\sigma}^{4}+3q^{3}q_{V^{T}}\widetilde{\sigma}^{2}z+q^{3}\widetilde{\sigma}^{4}-6q^{3}\widetilde{\sigma}^{2}z+2q^{2}q_{V^{T}}^{2}\bar{\sigma}_{X}^{4}\widetilde{\sigma}^{6}z$
$\displaystyle-3q^{2}q_{V^{T}}^{2}\bar{\sigma}_{X}^{4}\widetilde{\sigma}^{4}z^{2}-q^{2}q_{V^{T}}^{2}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{6}+3q^{2}q_{V^{T}}^{2}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{4}z+2q^{2}q_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{4}z-6q^{2}q_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{2}z^{2}-4q^{2}q_{V^{T}}\widetilde{\sigma}^{4}$
$\displaystyle+6q^{2}q_{V^{T}}\widetilde{\sigma}^{2}z+5q^{2}\widetilde{\sigma}^{4}-7q^{2}\widetilde{\sigma}^{2}z-q^{2}z^{2}-2qq_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{6}+8qq_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{4}z-6qq_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{2}z^{2}-2qq_{V^{T}}\widetilde{\sigma}^{4}$
$\displaystyle+3qq_{V^{T}}\widetilde{\sigma}^{2}z+q\widetilde{\sigma}^{2}z-2qz^{2}-\widetilde{\sigma}^{4}+2\widetilde{\sigma}^{2}z-z^{2}\,,$
(189) $\displaystyle e=$
$\displaystyle-3q^{3}q_{V^{T}}^{2}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{4}+3q^{3}q_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{4}-3q^{3}q_{V^{T}}\widetilde{\sigma}^{2}+3q^{3}\widetilde{\sigma}^{2}-q^{2}q_{V^{T}}^{2}\bar{\sigma}_{X}^{4}\widetilde{\sigma}^{6}+3q^{2}q_{V^{T}}^{2}\bar{\sigma}_{X}^{4}\widetilde{\sigma}^{4}z+2q^{2}q_{V^{T}}^{2}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{4}$
$\displaystyle-3q^{2}q_{V^{T}}^{2}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{2}z-6q^{2}q_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{4}+9q^{2}q_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{2}z-2q^{2}q_{V^{T}}\widetilde{\sigma}^{2}-q^{2}q_{V^{T}}z+q^{2}\widetilde{\sigma}^{2}+2q^{2}z+qq_{V^{T}}^{2}\bar{\sigma}_{X}^{4}\widetilde{\sigma}^{6}$
$\displaystyle-4qq_{V^{T}}^{2}\bar{\sigma}_{X}^{4}\widetilde{\sigma}^{4}z+3qq_{V^{T}}^{2}\bar{\sigma}_{X}^{4}\widetilde{\sigma}^{2}z^{2}+2qq_{V^{T}}^{2}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{4}-3qq_{V^{T}}^{2}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{2}z-2qq_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{4}+2qq_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{2}z$
$\displaystyle+2qq_{V^{T}}\bar{\sigma}_{X}^{2}z^{2}+2qq_{V^{T}}\widetilde{\sigma}^{2}-2qq_{V^{T}}z-3q\widetilde{\sigma}^{2}+3qz+2q_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{4}-4q_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{2}z+2q_{V^{T}}\bar{\sigma}_{X}^{2}z^{2}+q_{V^{T}}\widetilde{\sigma}^{2}$
$\displaystyle-q_{V^{T}}z-\widetilde{\sigma}^{2}+z\,,$ (190) $\displaystyle
f=$ $\displaystyle
3q^{2}q_{V^{T}}^{2}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{2}-3q^{2}q_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{2}+q^{2}q_{V^{T}}-q^{2}+2qq_{V^{T}}^{2}\bar{\sigma}_{X}^{4}\widetilde{\sigma}^{4}-3qq_{V^{T}}^{2}\bar{\sigma}_{X}^{4}\widetilde{\sigma}^{2}z-qq_{V^{T}}^{2}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{2}+qq_{V^{T}}^{2}\bar{\sigma}_{X}^{2}z$
$\displaystyle+3qq_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{2}-3qq_{V^{T}}\bar{\sigma}_{X}^{2}z+qq_{V^{T}}-q-q_{V^{T}}^{2}\bar{\sigma}_{X}^{4}\widetilde{\sigma}^{4}+2q_{V^{T}}^{2}\bar{\sigma}_{X}^{4}\widetilde{\sigma}^{2}z-q_{V^{T}}^{2}\bar{\sigma}_{X}^{4}z^{2}-q_{V^{T}}^{2}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{2}+q_{V^{T}}^{2}\bar{\sigma}_{X}^{2}z$
$\displaystyle+2q_{V^{T}}\bar{\sigma}_{X}^{2}\widetilde{\sigma}^{2}-2q_{V^{T}}\bar{\sigma}_{X}^{2}z\,,$
(191) $\displaystyle g=$ $\displaystyle-
qq_{V^{T}}^{2}\bar{\sigma}_{X}^{2}+qq_{V^{T}}\bar{\sigma}_{X}^{2}-q_{V^{T}}^{2}\bar{\sigma}_{X}^{4}\widetilde{\sigma}^{2}+q_{V^{T}}^{2}\bar{\sigma}_{X}^{4}z\,.$
(192)

We solve this polynomial equation numerically, looking for complex roots which
yield non-zero values of the eigenvalue density. For large signal-to-noise
ratio, we encounter numerical instabilities trying to determine the eigenvalue
density bounds. We run into these instabilities in the determination of the
true density bounds in Fig. 2(d). To fix this, we start at the peak of the
density and determine the values of $\lambda$, for which it hits zero for the
first time, to either side of the peak. All other zero density crossings are
assumed to be due to numerical instabilities.

The ranges of nonzero density are shown in Fig. 2(d). We see that, as the
signal-to-noise ratio increases, there is a bifurcation point where the
density splits into two bumps. The left bump is associated with the noise, and
the right bump is associated with pure latent feature signal. From our
approximate expression for the eigenvalue bounds, Eq. (22), we can estimate
the value of the SNR at which the splitting occurs. The defining equation for
this is given by the intersection between the right boundary of the noise
region and the left boundary of the signal part of the densit in Eq. (22),
resulting in:

$\displaystyle(1+\textsf{SNR}^{-1})\lambda_{+}^{\text{MP}}=\lambda_{-}^{\textsf{SNR}}\,.$
(193)

Solving for SNR, we obtain the following estimation for the splitting point:

$\displaystyle\textsf{SNR}_{\text{split}}\approx\frac{\lambda_{+}^{\text{MP}}-\lambda_{-}^{\text{MP}}}{\lambda_{-}^{\infty}}\,.$
(194)

## References

  * Wishart and Bartlett [1932] J. Wishart and M. S. Bartlett, The distribution of second order moment statistics in a normal system, Mathematical Proceedings of the Cambridge Philosophical Society 28, 455–459 (1932).
  * [2] Wolfram Research, Inc., functions.wolfram.com, http://functions.wolfram.com/03.04.21.0008.01, accessed: 2021-06-30.
  * [3] H. (https://math.stackexchange.com/users/6460/henry), Sum of $n$ i.i.d Beta-distributed variables, Mathematics Stack Exchange, https://math.stackexchange.com/q/3096929 (version: 2019-02-02).
  * Potters and Bouchaud [2020] M. Potters and J.-P. Bouchaud,  _A First Course in Random Matrix Theory: For Physicists, Engineers and Data Scientists_ (Cambridge University Press, 2020).
  * Voiculescu  _et al._ [1992] D. V. Voiculescu, K. J. Dykema, and A. Nica,  _Free random variables_ , 1 (American Mathematical Soc., 1992).
  * Müller [2002] R. R. Müller, A random matrix model of communication via antenna arrays, IEEE Transactions on information theory 48, 2495 (2002).
  * Burda  _et al._ [2010] Z. Burda, A. Jarosz, G. Livan, M. A. Nowak, and A. Swiech, Eigenvalues and singular values of products of rectangular gaussian random matrices, Phys. Rev. E 82, 061114 (2010).
  * Dupic and Castillo [2014] T. Dupic and I. P. Castillo, Spectral density of products of Wishart dilute random matrices. Part I: the dense case (2014), arXiv:1401.7802 [cond-mat.dis-nn] .
  * Livan  _et al._ [2018] G. Livan, M. Novaes, and P. Vivo, Introduction to random matrices, SpringerBriefs in Mathematical Physics 10.1007/978-3-319-70885-0 (2018).
  * Rabinowitz [1993] S. Rabinowitz, How to find the square root of a complex number, Mathematics and Informatics Quarterly 3, 54 (1993).
  * Cui  _et al._ [2020] W. Cui, J. W. Rocks, and P. Mehta, The perturbative resolvent method: spectral densities of random matrix ensembles via perturbation theory (2020), arXiv:2012.00663 [cond-mat.dis-nn] .