# Optimal Estimation and Computational Limit of Low-rank Gaussian Mixtures

Zhongyuan Lyu and Dong Xia111Dong Xia’s research was partially supported by
Hong Kong RGC Grant ECS 26302019, GRF 16303320 and GRF 16300121.
Hong Kong University of Science and Technology

(())

###### Abstract

Structural matrix-variate observations routinely arise in diverse fields such
as multi-layer network analysis and brain image clustering. While data of this
type have been extensively investigated with fruitful outcomes being
delivered, the fundamental questions like its statistical optimality and
computational limit are largely under-explored. In this paper, we propose a
low-rank Gaussian mixture model (LrMM) assuming each matrix-valued observation
has a planted low-rank structure. Minimax lower bounds for estimating the
underlying low-rank matrix are established allowing a whole range of sample
sizes and signal strength. Under a minimal condition on signal strength,
referred to as the information-theoretical limit or statistical limit, we
prove the minimax optimality of a maximum likelihood estimator which, in
general, is computationally infeasible. If the signal is stronger than a
certain threshold, called the computational limit, we design a computationally
fast estimator based on spectral aggregation and demonstrate its minimax
optimality. Moreover, when the signal strength is smaller than the
computational limit, we provide evidences based on the low-degree likelihood
ratio framework to claim that no polynomial-time algorithm can consistently
recover the underlying low-rank matrix. Our results reveal multiple phase
transitions in the minimax error rates and the statistical-to-computational
gap. Numerical experiments confirm our theoretical findings. We further
showcase the merit of our spectral aggregation method on the worldwide food
trading dataset.

##  1 Introduction

The recent decade has witnessed a burgeoning demand in processing and
analyzing large-scale matrix-variate data which routinely arise in diverse
fields. In gene expression analysis, e.g., the BHL (brain, heart and lung)
dataset (BHL, ; Mai et al., 2021), the measurement of gene expression on
different types of tissues is often repeated for multiple times. The resultant
observation for each tissue becomes a matrix and thus the cluster analysis is
operated on matrix-valued observations. A multi-layer network (Le et al.,
2018; Jing et al., 2021; Lyu et al., 2021; Paul and Chen, 2020) usually
consists of multiple networks on the same set of vertices. Since each observed
layer is equivalently represented as an adjacent matrix, problems such as
community detection (Paul and Chen, 2020), layer clustering (Jing et al.,
2021) and common probability matrix estimation (Le et al., 2018) are generally
attacked by statistical analysis on a collection of adjacency matrices. Other
notable examples include brain image clustering (Sun and Li, 2019; Wang et
al., 2017), EEG data analysis (Hu et al., 2020; Gao et al., 2021), etc.
Oftentimes, the dimensions of observed matrices are ultra-large or the number
of matrix-valued observations is relatively small, which has motivated the
exploration of hidden low-dimensional structures, e.g. sparsity and low-
rankness, in matrix-valued observations. All the aforementioned works assumed,
among others, certain types of low-rank structures for the underlying
parameters of interest and have delivered fruitful outcomes in real-world
applications.

Inspired by those foregoing works, throughout this paper, we assume that each
matrix-valued observation has a low-rank expectation which might vary for
different observations. Towards that end, several specific low-rank
statistical models, tailored for concrete applications, and respective
estimating procedures have been proposed. For instance, a mixture multi-layer
stochastic block model (SBM) was introduced in Jing et al. (2021) for
uncovering the global and local communities in multi-layer networks. At the
core of this model is the assumption that every layer has a low-rank expected
adjacency matrix. Their estimator was based on the (regularized) low-rank
tensor decomposition. A special multi-layer SBM was proposed by Paul and Chen
(2020) and estimated by a spectral method. In order to analyze the brain fMRI
data, Sun and Li (2019) proposed a tensor Gaussian mixture model and designed
an estimator via (fusedly-)truncated low-rank tensor decomposition. Despite
these prior efforts, usually motivated by particular applications, on the low-
rank estimates from a mixture of matrix-valued observations, many fundamental
questions remain unanswered. What is the role and benefit of low-rankness? How
do the sample size and signal strength (see the definition after eq.(2))
characterize the intrinsic difficulty, i.e., are there any phase transitions?
What is the statistically optimal rate, which estimator can achieve the rate
and is this estimator computationally feasible? What is the fastest error rate
achievable by estimators requiring only polynomial-time complexity? This paper
aims to answer all these questions and provides a complete picture for the
statistical and computational limits in the low-rank estimation from a mixture
of matrix-valued observations.

We now introduce the low-rank Gaussian mixture model (LrMM) to formalize the
questions. For simplicity, we focus on the mixture of two components and will
briefly discuss the case of multiple components in Section 7. The $d_{1}\times
d_{2}$ matrix $\mathbf{X}$ is said to follow an isotropic matrix normal (Gupta
and Nagar, 2018) distribution ${\cal
N}(\mathbf{M},\mathbf{I}_{d_{1}}\otimes\mathbf{I}_{d_{2}})$ if ${\it
vec}(\mathbf{X})\sim{\cal N}\big{(}{\it
vec}(\mathbf{M}),\mathbf{I}_{d_{1}d_{2}}\big{)}$, where $\mathbf{I}_{d}$
represents the $d\times d$ identity matrix and $\mathbf{M}$ is a deterministic
matrix. Clearly, this implies that $\mathbf{X}$ is equal to
$\mathbf{M}+\mathbf{Z}$ in distribution where $\mathbf{Z}$ has i.i.d. standard
normal entries. Denote222With a slight abuse of notation, we also denote
$p_{\mathbf{M}}$ the associated probability density function.

$p_{\mathbf{M}}=\frac{1}{2}{\cal
N}(\mathbf{M},\mathbf{I}_{d_{1}}\otimes\mathbf{I}_{d_{2}})+\frac{1}{2}{\cal
N}(-\mathbf{M},\mathbf{I}_{d_{1}}\otimes\mathbf{I}_{d_{2}})$ (1)

the symmetric mixture of two-component Gaussian mixture model. Then
$\mathbf{X}\sim p_{\mathbf{M}}$ means that $\mathbf{X}$ is sampled from ${\cal
N}(\mathbf{M},\mathbf{I}_{d_{1}}\otimes\mathbf{I}_{d_{2}})$ and ${\cal
N}(-\mathbf{M},\mathbf{I}_{d_{1}}\otimes\mathbf{I}_{d_{2}})$ with probability
both $1/2$, respectively. Put it differently, $\mathbf{X}$ equals
$s\mathbf{M}+\mathbf{Z}$ in distribution with $s$ being a Rademacher random
variable, called the label of $\mathbf{X}$, satisfying $\mathbb{P}(s=\pm
1)=1/2$. Throughout the paper, we assume that $\mathbf{M}$ has a small rank,
i.e., $r={\rm rank}(\mathbf{M})\ll\min\\{d_{1},d_{2}\\}$. Note that, under
model (1), the marginal expectation of $\mathbf{X}$ is actually zero. The
former claim of low-rank expectation in the last paragraph actually refers to
the conditional expectation $\mathbb{E}(\mathbf{X}|s)=s\mathbf{M}$ which is
low-rank. We remark that the condition of equal prior probabilities is not
essential and can be slightly relaxed. The assumption of symmetry of the two
components is only for ease of exposition. If the two components have distinct
mean matrices, say $\mathbf{M}_{1}$ and $\mathbf{M}_{2}$, respectively, one
can first estimate the average $(\mathbf{M}_{1}+\mathbf{M}_{2})/2$, subtract
it from all observations and reduce the problem to the symmetric case.
Similarly, the assumption of isotropic noise is relaxable as long as the
covariance tensor is known. The case of unknown covariance is much more
challenging (Davis et al., 2021; Bakshi et al., 2020; Belkin and Sinha, 2010;
Cai et al., 2019; Ge et al., 2015; Moitra and Valiant, 2010) even in the
vector case and is beyond the scope of the current paper.

Given i.i.d. observations $\mathbf{X},\mathbf{X}_{1},\cdots,\mathbf{X}_{n}$
sampled from the mixture distribution $p_{\mathbf{M}}$ in (1), our goals are
to estimate the latent low-rank matrix $\mathbf{M}$, establish the minimax
error rates and design computationally efficient estimators. We assume
$d_{1}\asymp d_{2}\asymp d$ meaning that there exist absolute constants
$c_{0},C_{0}>0$ satisfying
$c_{0}d\leq\min\\{d_{1},d_{2}\\}\leq\max\\{d_{1},d_{2}\\}\leq C_{0}d$. The
parameter space of interest is, for any $\lambda>0$,

${\cal
M}_{d_{1},d_{2}}(r,\lambda):=\left\\{\mathbf{M}\in\mathbb{R}^{d_{1}\times
d_{2}}:{\rm rank}(\mathbf{M})=r,\
\lambda\asymp\sigma_{r}(\mathbf{M})\leq\cdots\leq\sigma_{1}(\mathbf{M})\asymp\lambda\right\\}$
(2)

where $\sigma_{k}(\cdot)$ denotes the $k$-th largest singular value of a
matrix. For notational brevity, we shall write ${\cal M}(r,\lambda)$ for
short. The signal strength of low-rank models is usually determined by the
smallest non-zero singular value (Koltchinskii and Xia, 2016; Zhang and Xia,
2018; Xia, 2021; Cheng et al., 2021; Gavish and Donoho, 2014). The set ${\cal
M}_{d_{1},d_{2}}(r,\lambda)$ is the collection of all $d_{1}\times d_{2}$
rank-$r$ matrices whose signal strength is of order $\lambda$. For simplicity,
we only focus on the well-conditioned matrices, i.e., with a bounded condition
number. The minimax error rate of estimating $\mathbf{M}$ is defined by
$\inf_{\widehat{}\mathbf{M}}\sup_{\mathbf{M}\in{\cal
M}(r,\lambda)}\mathbb{E}\ell(\widehat{}\mathbf{M},\mathbf{M})$ , where the
infimum is taken over all possible estimator $\widehat{}\mathbf{M}$
constructed from the i.i.d. observations
$\mathbf{X}_{1},\cdots,\mathbf{X}_{n}$ and the loss function is
$\ell(\widehat{}\mathbf{M},\mathbf{M}):=\min_{\eta=\pm
1}\|\widehat{}\mathbf{M}-\eta\mathbf{M}\|_{\rm F}$ with $\|\cdot\|_{\rm F}$
standing for the Frobenius norm. Note that, due to the symmetry of model (1),
$\mathbf{M}$ is estimable up to a sign flip.

If $\mathbf{M}$ has a full rank with $r=\min\\{d_{1},d_{2}\\}$, model (1)
reduces to the canonical two component isotropic Gaussian mixture model (GMM)
in the dimension $d_{1}d_{2}\asymp d^{2}$, which has been extensively
investigated in the literature. See Balakrishnan et al. (2017); Chen (1995);
Ho and Nguyen (2016a); Xu et al. (2016); Wu and Yang (2020) and references
therein. For instance, Wu and Zhou (2019) proved that the minimax rate 333Note
that there is an additional term $d^{2}(\theta n)^{-1}$ derived in Wu and Zhou
(2019) which is actually negligible if inspecting all other terms carefully.
is

$\displaystyle\inf_{\widehat{}\mathbf{M}}\sup_{\|\mathbf{M}\|_{\rm
F}=\theta}\mathbb{E}\ell(\widehat{}\mathbf{M},\mathbf{M})\asymp\min\left\\{\frac{1}{\theta}\frac{d}{n^{1/2}}+\frac{d}{n^{1/2}},\
\theta\right\\}$ (3)

, and showed that a simple spectral method, together with a trivial estimate
for the case of small $\theta$, is minimax optimal. This rate implies
intriguing phenomenons of phase transitions concerning the sample size $n$ and
signal strength $\theta$. For instance, if the sample size $n\geq d^{2}$,
their result reveals three different minimax rates: $\theta$ for $\theta\leq
d^{1/2}n^{-1/4}$, $\theta^{-1}dn^{-1/2}$ for $d^{1/2}n^{-1/4}\leq\theta\leq 1$
and $dn^{-1/2}$ for $\theta\geq 1$. Interestingly, it also implies that non-
trivial estimate is impossible, i.e., information-theoretically impossible, if
the signal strength is smaller than $d^{1/2}n^{-1/4}$. Undoubtedly, if
$\mathbf{M}$ is low-rank with $r\ll d$, one can naturally foresee the
existence of multiple phase transitions for the minimax error rates.
Establishing these rates becomes more challenging for several reasons. On the
methodological front, a naive spectral method cannot attain the minimax
optimal rate and thus additional procedures are necessary. On the theoretical
front, the low-rank structure dictates a smaller intrinsic dimension and
brings about new behaviors to the phase transitions of the minimax error
rates. See, e.g. Koltchinskii and Xia (2015); Ma and Wu (2015) and references
therein. On the computational front, it is well recognized that the low-
rankness sometimes bears a so-called statistical-to-computational gap (Barak
and Moitra, 2016; Zhang and Xia, 2018) in the sense that there exist regimes
where statistically optimal estimators can be computationally infeasible,
e.g., requiring an exponential-time complexity.

The summary of our contributions is as follows. We establish the minimax rate
of estimating the rank-$r$ matrix $\mathbf{M}$ for the LrMM model that reads
as

$\displaystyle\inf_{\widehat{}\mathbf{M}}\ \sup_{\mathbf{M}\in{\cal
M}_{d_{1},d_{2}}(r,\lambda)}\mathbb{E}\ell(\widehat{}\mathbf{M},\mathbf{M})\
\asymp\
\min\left\\{\frac{1}{\lambda}\Big{(}\frac{d}{n}\Big{)}^{1/2}+\Big{(}\frac{dr}{n}\Big{)}^{1/2},\
\lambda r^{1/2}\right\\}$ (4)

where the infimum is taken over all possible estimators, regardless of their
computational feasibility. This rate implies that, when the sample size $n\geq
dr$, it is information-theoretically impossible to estimate $\mathbf{M}$ if
the signal strength $\lambda$ is smaller than
$d^{1/4}(rn)^{-1/4}+(d/n)^{1/2}$. Under minimal conditions, we prove that the
maximum likelihood estimator (MLE) can achieve the rate (4) up to a
logarithmic factor. Unfortunately, there are no known polynomial-time
algorithms with guaranteed performance to solve MLE. Earlier works (Tosh and
Dasgupta, 2017; Sanjeev and Kannan, 2001) show that solving MLE is generally
NP-hard. We then propose a computationally fast estimator based on spectral
aggregation. This approach can be viewed as a modified method of second moment
(Pearson, 1894; Wu and Yang, 2020) adapted with a spectral projection to
leverage the low-rank structure. We prove that this computationally efficient
estimator can achieve the minimax rate (4) as long as the signal strength
$\lambda$ is larger than $d^{1/2}n^{-1/4}$, which is much stronger than the
information-theoretical requirement for the minimal signal strength. This
difference unveils the statistical-to-computational gap in LrMM. Lastly, we
adopt the low-degree likelihood ratio framework (Kunisky et al., 2019) to
conjecture that no polynomial-time estimator is consistent if $\lambda$ is
smaller than $d^{1/2}n^{-1/4}$. The minimax rates, phase transitions and
statistical-to-computational gaps are illustrated in Figure 1.

Figure 1: The minimax rates, phase transitions and statistical-to-
computational gaps of LrMM, model (1). Here $r$ is the rank, the matrix
dimension $d_{1}\asymp d_{2}\asymp d$, $n$ is the sample size and $\lambda$
denotes the smallest non-zero singular value. There exist three regimes
concerning the sample size which are colored in blue. The minimax error rates
(up to logarithmic factor) of estimating $\mathbf{M}\in{\cal M}(r,\lambda)$ in
different regimes are colored in red. Here information-theoretically
impossible means that non-trivial estimates are impossible because of weak
signal strength. Within the low-degree likelihood ratio framework (Kunisky et
al., 2019), we provide evidence showing that no polynomial-time algorithms can
consistently estimate $\mathbf{M}$ if $\lambda$ is smaller than
$d^{1/2}n^{-1/4}$.

Our results are closely related yet crucially different from several existing
works. In Chen et al. (2021), a low-rank mixture model was proposed for linear
regression which is generally more challenging than our model (1). They
designed a computationally efficient estimator but provided no results
respecting the statistical optimality or computational limits. A multi-graph
network model was introduced by Wang et al. (2019) which allows heterogeneous
structure on each matrix-valued observation. However, their model has no
mixture nature and there is no guarantee on minimax optimality. More recently,
Jing et al. (2021) proposes a mixture multi-layer SBM and establishes the
minimax error rate of spectral estimate only for the special regime when the
sample size $n$ is smaller than $d$ and the signal strength, reflected by the
network sparsity, is strong enough. In addition, our LrMM is directly related
to low-rank tensor literature. By stacking the matrix observations slice by
slice, we end up with a tensor of size $n\times d_{1}\times d_{2}$ whose
expectation, under model (1), has a low Tucker rank $(1,r,r)$. See, e.g.,
Zhang and Xia (2018); Jing et al. (2021) and references therein. Minimax rates
for low-rank tensor denoising and noisy tensor completion have been
investigated by Zhang and Xia (2018) and Xia et al. (2021), respectively.
However, they both require the sample size $n$ to be of the same order of $d$,
which becomes unrealistic in the low-rank mixture model. Finally, it worths to
remark that our bound (4) reduces to the minimax bound of GMM (3) if we let
$\mathbf{M}$ be full-rank. To see this, one can just replace $\lambda$ and $r$
in our bound (4) by $\theta d^{-1/2}$ and $d$, respectively.

The rest of paper is organized as follows. We establish the minimax lower
bound in Section 2 and prove that the maximum likelihood estimator, albeit
computationally infeasible in general, achieves the minimax optimal rates. A
computationally fast estimator based on spectral aggregation is proposed in
Section 3 which attains minimax optimal rates as long as the signal strength
is strong. Section 4 justifies the statistical-to-computational gap by showing
that there exists some regime where the MLE can attains minimax optimal rates
but no-polynomial time algorithms can consistently recover the underlying low-
rank matrix. We then showcase results of numerical simulations in Section 5,
present a real-world data experiment in Section 6, and discuss open questions
and potential directions in Section 7.

##  2 Maximum likelihood estimator and minimax optimality

We slightly abuse the notation and denote $p_{\mathbf{M}}(\cdot)$ the
probability density function of $\mathbf{X}\in\mathbb{R}^{d_{1}\times d_{2}}$
under the LrMM model (1). The family of density functions parameterized by
${\cal M}_{d_{1},d_{2}}(r,\lambda)$ is written as (note that we assume
$d_{1}\asymp d_{2}\asymp d$)

${\cal P}_{d_{1},d_{2}}(r,\lambda):=\Big{\\{}p_{\mathbf{M}}:\mathbf{M}\in{\cal
M}_{d_{1},d_{2}}(r,\lambda)\Big{\\}}$

which is indexed by rank-$r$ matrices with the signal strength $\lambda$.
Given $i.i.d.$ observations $\mathbf{X}_{1},\cdots,\mathbf{X}_{n}$ sampled
from $p_{\mathbf{M}}$, the maximum likelihood estimator (not necessarily
unique) is defined by

$\displaystyle p_{\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}}}:=\underset{p_{\mathbf{M}}\in{\cal
P}_{d_{1},d_{2}}(r,\lambda)}{\arg\max}\
\sum_{i=1}^{n}\log\big{(}p_{\mathbf{M}}(\mathbf{X}_{i})\big{)}$ (5)

While the MLE estimator (5) is generally NP-hard to compute, it often serves
as a benchmark for understanding the information-theoretical limit of a
statistical model.

We begin with the regime $n=\tilde{\Omega}(dr)$444Here, $\tilde{\Omega}$
stands for the standard big-$\Omega$ notation up to a logarithmic factor.,
which falls into the typical low-dimensional setting555The low-dimensional
setting here refers to the case that dimension $d$ is allowed to grow with
sample size $n$, while the order of $n$ still dominates.. The convergence rate
of MLE in this regime has been thoroughly investigated for Gaussian mixture
model. See, for instance, Leroux (1992); Van de Geer (1993); Chen (1995);
Genovese and Wasserman (2000); Ghosal and Van Der Vaart (2001). The standard
tool, e.g. Van de Geer (1993) and (Van de Geer, 2000, Theorem 7.4),
establishes the convergence rate of MLE in the Hellinger distance defined by
$d_{H}(p_{\mathbf{M}_{1}},p_{\mathbf{M}_{2}}):=1-\int
p_{\mathbf{M}_{1}}^{1/2}(\mathbf{X})p_{\mathbf{M}_{2}}^{1/2}(\mathbf{X})d\mathbf{X}$
for two density functions $p_{\mathbf{M}_{1}}(\cdot)$ and
$p_{\mathbf{M}_{2}}(\cdot)$. According to this tool, it suffices to bound the
bracketing entropy number of a class of square root density functions around
the truth $p_{\mathbf{M}}^{1/2}$. While existing literature (Ho and Nguyen,
2016a, b; Maugis and Michel, 2011) have developed respective bracketing
entropy bounds for Gaussian mixture model, they only focus on the fixed
dimension $d$ and their method is inapplicable to matrix-variate observations
with a planted low-rank structure. By a covering argument and the construction
of bracket functions, we establish such a bracketing entropy bounds for LrMM
and derive the upper bound in Hellinger distance for
$d_{\textsf{H}}(\widehat{p}_{\mathbf{M}_{\textsf{\tiny
MLE}}},p_{\mathbf{M}})$.

To bridge the density estimation and parameter estimation, we resort to a
sharp characterization for the total variation distance (similarly, the
Hellinger distance) between Gaussian mixture densities established recently by
Davies et al. (2021).

######  Lemma 1.

(Lower bound of Hellinger distance) Let $\mathbf{M}_{1}$ and $\mathbf{M}$ be
two matrices, and denote $p_{\mathbf{M}_{1}}$ and $p_{\mathbf{M}}$ the two
density functions defined by (1). There exists absolute constants
$c_{0},c_{1},c_{2}>0$ such that, if $\|\mathbf{M}\|_{\rm
F}+\|\mathbf{M}_{1}\|_{\rm F}\leq c_{0}$ then

$d_{H}(p_{\mathbf{M}_{1}},p_{\mathbf{M}})\geq c_{1}\big{(}\|\mathbf{M}\|_{\rm
F}+\|\mathbf{M}_{1}\|_{\rm F}\big{)}\ell(\mathbf{M}_{1},\mathbf{M})$

Otherwise

$d_{H}(p_{\mathbf{M}_{1}},p_{\mathbf{M}})\geq
c_{2}\min\big{\\{}1,\ell(\mathbf{M}_{1},\mathbf{M})\big{\\}}$

where
$\ell(\mathbf{M}_{1},\mathbf{M}):=\min\\{\|\mathbf{M}_{1}-\mathbf{M}\|_{\rm
F},\ \|\mathbf{M}_{1}+\mathbf{M}\|_{\rm F}\\}$.

Together with the upper bound of Hellinger distance
$d_{\textsf{H}}(\widehat{p}_{\mathbf{M}_{\textsf{\tiny MLE}}},p_{\mathbf{M}})$
and Lemma 1, we obtain the error rate of the maximum likelihood estimator when
$n=\tilde{\Omega}(dr)$, namely the first part of Theorem 1.

However, the above argument fails when it comes to the regime
$n=\tilde{O}(dr)$666Again, $\tilde{O}$ stands for the standard big-$O$
notation up to a logarithmic factor., corresponding to an ultra high-
dimensional setting. The reason is that the minimax lower bound, as we will
see later in Theorem 2, suggests that the optimal error rate should be of
order $(dr/n)^{1/2}$, which can be larger than 1. Consequently, the Hellinger
distance is no longer an appropriate metric777The error rate of other bounded
metric, say, the Wasserstein distance considered in Doss et al. (2020), also
becomes trivial when $d>n$., for instance, the lower bound in Lemma 1 becomes
trivial. To this end, we turn to Kullback-Leibler (KL) divergence defined by
$D_{\textsf{KL}}(p_{\mathbf{M}_{1}}\|p_{\mathbf{M}_{2}}):=\int
p_{\mathbf{M}_{1}}(\mathbf{X})\log(p_{\mathbf{M}_{1}}(\mathbf{X})/p_{\mathbf{M}_{2}}(\mathbf{X}))d\mathbf{X}$.
Though KL divergence is not a metric itself, in many cases its convergence
also implies consistency of parameter estimate in some metric of interest (Van
de Geer, 2000). Moreover, the KL divergence in its form is closely related to
MLE and its unboundedness property is beneficial for our purpose since
$(dr/n)^{1/^{\prime}2}$ possibly diverges. By carefully characterizing the
distribution of
$\log(p_{\mathbf{M}_{1}}(\mathbf{X})/p_{\mathbf{M}_{2}}(\mathbf{X}))$ and
exploiting the concentration inequality of suprema of an empirical process
(Adamczak, 2008, Theorem 4), we are able to derive an upper bound for the KL
divergence
$D_{\textsf{KL}}(p_{\mathbf{M}}\|p_{\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}}})$. We also establish the following lower bound relating KL divergence
to the distance in the parameter space. Combining Lemma 2 with the upper bound
of $D_{\textsf{KL}}(p_{\mathbf{M}}\|p_{\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}}})$ leads to the desired error rate in the regime $n=\tilde{O}(dr)$,
i.e., the second part of Theorem 1.

######  Lemma 2.

(Lower bound of KL divergence) Let $\mathbf{M}_{1}$ and $\mathbf{M}$ be two
matrices, and denote $p_{\mathbf{M}_{1}}$ and $p_{\mathbf{M}}$ the two density
functions defined by (1). There exists absolute constants
$C_{0},C_{1}>1,c_{0}>0$ such that if $\|\mathbf{M}\|_{\rm F}\geq C_{0}$ and
$\|\mathbf{M}-\mathbf{M}_{1}\|_{\rm F}\geq C_{1}$, then

$D_{\textsf{KL}}(p_{\mathbf{M}}\|p_{\mathbf{M}_{1}})\geq
c_{0}\cdot\ell^{2}(\mathbf{M},\mathbf{M}_{1})$

Collecting two pieces, the error rate of the maximum likelihood estimator is
summarized in the following theorem.

######  Theorem 1.

Suppose $\mathbf{M}\in{\cal M}(r,\lambda)$ and let
$\widehat{}\mathbf{M}_{\textsf{\tiny MLE}}$ denote the maximum likelihood
estimator by (5).

  1. (1)

If $dr\log nd<n$, then there exist absolute constants
$c_{1},c_{2},C_{1},C_{2},C_{3}>0$ such that the following bound holds with
probability at least $1-\exp(-c_{1}d\log^{2}(nd))$,

$\displaystyle\ell(\widehat{}\mathbf{M}_{\textsf{\tiny MLE}},\mathbf{M})\leq
C_{1}\left(\sqrt{\frac{dr\log(nd)}{n}}+\frac{1}{\lambda}\sqrt{\frac{d\log(nd)}{n}}\right)$
(6)

If further assume $\lambda\leq C_{2}\exp(c_{2}d\log^{2}(nd))$, then

$\mathbb{E}\ell(\widehat{}\mathbf{M}_{\textsf{\tiny MLE}},\mathbf{M})\leq
C_{3}\left(\sqrt{\frac{dr\log(nd)}{n}}+\frac{1}{\lambda}\sqrt{\frac{d\log(nd)}{n}}\right)$

  2. (2)

If $dr\log nd\geq n$, then there exist absolute constants
$C_{4},C_{5},C_{6},C_{7}>0$ such that if
$C_{4}r^{-1/2}\leq\lambda\leq C_{5}d^{1/2}$, then the following bound holds
with probability at least $1-(nd)^{-4}$,

$\displaystyle\ell(\widehat{}\mathbf{M}_{\textsf{\tiny MLE}},\mathbf{M})\leq
C_{6}\sqrt{\frac{dr\log(nd)}{n}}$ (7)

And the following bound in expectation holds,

$\mathbb{E}\ell(\widehat{}\mathbf{M}_{\textsf{\tiny MLE}},\mathbf{M})\leq
C_{7}\sqrt{\frac{dr\log(nd)}{n}}$

We note that the logarithmic factor in (6) emerges from the bracketing entropy
bound and that in (7) arises from the tail inequality for suprema of empirical
processes of unbounded functions. The high probability bound in the first part
of Theorem 1 is proved without conditions on the sample size $n$, the rank $r$
or on the signal strength $\lambda$. It suggests intriguing phase transitions
in the regime $n=\tilde{\Omega}(dr)$. When $\lambda>r^{-1/2}$, the MLE attains
the rate $\tilde{O}\big{(}(rd/n)^{1/2}\big{)}$, growing with respect to the
rank $r$, which is the best achievable rate even if the labels of observations
are all known, namely in the oracle scenario. On the other hand, if
$\lambda<r^{-1/2}$, the MLE attains the rate
$\tilde{O}\big{(}\lambda^{-1}(d/n)^{1/2}\big{)}$ that is free of the
underlying rank $r$. Moreover, a trivial estimate by
$\widehat{}\mathbf{M}={\bf 0}$ attains the error rate $r^{1/2}\lambda$.
Therefore, the MLE becomes pointless if $\lambda$ is smaller than
$d^{1/4}(rn)^{-1/4}+(d/n)^{1/2}$, which is referred to as the information-
theoretically impossible regime. In the second statement of Theorem 1, a more
stringent condition is imposed on signal strength ($\lambda=O(d^{1/2})$) for
technical difficulty, though we believe that MLE could attain the optimal rate
$\tilde{O}((rd/n)^{1/2})$ in a wider range of $\lambda$ via more sophisticated
analysis. On the other hand, as long as $\lambda=\Omega(d^{1/2}n^{-1/4})$, a
computationally efficient estimator (see Section 3) is already able to attain
the optimal rate. As we intend to reveal the optimal estimation rate under
different signal strength, we only appeal to MLE when the signal strength is
not strong enough. Therefore, the technical condition of $\lambda$ for MLE is
not essential.

The next theorem demonstrates the minimax optimality of the MLE by
establishing a matching minimax lower bound up to the logarithmic factor. We
note that the minimax lower bound (8) is a statistical lower bound because it
takes no considerations of the computational feasibility. In Section 3, we
introduce a computationally fast estimator that achieves these lower bounds
but requires much more stringent conditions.

######  Theorem 2.

There exists an absolute constant $c_{1}>0$ such that

$\inf_{\widehat{}\mathbf{M}}\
\sup_{\mathbf{M}\in\mathcal{M}(r,\lambda)}\mathbb{E}\ell(\widehat{}\mathbf{M},\mathbf{M})\geq
c_{1}\left(\sqrt{\frac{dr}{n}}+\frac{1}{\lambda}\sqrt{\frac{d}{n}}\right)\bigwedge\lambda\sqrt{r},$
(8)

where the infimum is taken over all possible estimators and $a\wedge
b=\min\\{a,b\\}$.

##  3 Computationally efficient estimator by spectral aggregation

Since the MLE (5) is generally computationally infeasible, it is of crucial
importance to design an estimator which is polynomial-time computable. While
existing works have demonstrated the optimality of spectral method for both
estimation (Wu and Zhou, 2019) and clustering (Löffler et al., 2019) under the
GMM, it turns out that a naive spectral estimate is statistically sub-optimal
for our LrMM and additional subsequent treatments are necessary.

For technical simplicity, we adopt the sample splitting in our estimating
procedure. It will inevitably affect the constant factor in the error rate,
e.g., the $C_{1},C_{3}$ as in Theorem 1. Since our main interest concerns only
the convergence rate in terms of the model parameters, we spare no efforts to
improve the constant factor.

Without loss of generality, assume the sample size $n=4n_{0}$. We randomly
split the original sample $\mathbf{X}_{1},\cdots,\mathbf{X}_{n}$ into four
disjoint subsets of equal size, denoted by
$\\{\mathbf{X}_{i}^{(k)}\\}_{i=1}^{n_{0}}$ for $k=1,2,3,4$. Our estimating
procedure consists of three major steps:

  * -

Step 1 (Spectral initialization). Stack the observations column by column into
a $d_{1}\times(n_{0}d_{2})$ matrix
$[\mathbf{X}_{1}^{(1)},\cdots,\mathbf{X}_{n_{0}}^{(1)}]$, extract its leading
left singular vector, denoted by $\widehat{}\mathbf{u}_{1}$. Then, construct
the $d_{2}\times n_{0}$ matrix
$[\mathbf{X}_{1}^{(2)\top}\widehat{}\mathbf{u}_{1},\cdots,\mathbf{X}_{n_{0}}^{(2)\top}\widehat{}\mathbf{u}_{1}]$
and extract its left singular vector, denoted by $\widehat{}\mathbf{v}_{1}$.

  * -

Step 2 (Spectral refinement). Extract the top-$r$ left and right singular
vectors of

$\displaystyle\tilde{}\mathbf{U},\tilde{}\mathbf{V}\ \
\stackrel{{\scriptstyle{\rm SVD}_{r}}}{{\parbox{42.67912pt}{\leftarrowfill}}}\
\
\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}(\widehat{}\mathbf{u}_{1}^{\top}\mathbf{X}_{i}^{(3)}\widehat{}\mathbf{v}_{1})\mathbf{X}_{i}^{(3)}-\widehat{}\mathbf{u}_{1}\widehat{}\mathbf{v}_{1}^{\top}$
(9)

  * -

Step 3 (Aggregation). Denote $\check{}\mathbf{M}$ the best rank-$r$
approximation of

$\displaystyle\check{\mathbf{M}}\ \ \stackrel{{\scriptstyle{\rm rank}-r\ {\rm
approx.}}}{{\parbox{56.9055pt}{\leftarrowfill}}}\ \
\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}{\rm
Tr}(\tilde{}\mathbf{U}^{\top}\mathbf{X}_{i}^{(4)}\tilde{}\mathbf{V})\mathbf{X}_{i}^{(4)}-\tilde{}\mathbf{U}\tilde{}\mathbf{V}^{\top}$
(10)

Compute the scaling factor by

$\widehat{\Lambda}\ \ \longleftarrow\ \
\left[\max\left\\{\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}{\rm
Tr}^{2}(\tilde{}\mathbf{U}^{\top}\mathbf{X}_{i}^{(4)}\tilde{}\mathbf{V})-r,\
\frac{dr^{2}}{\sqrt{n}}\right\\}\right]^{1/2}$

The final estimator is defined by
$\widehat{}\mathbf{M}=\widehat{\Lambda}^{-1}\check{}\mathbf{M}$.

Due to eq. (10), we refer to this procedure as the spectral aggregation. Note
that (10) is, in spirit, similar to the method of second moment as in Gaussian
mixture model (Wu and Yang, 2020). The additional projection onto
$\tilde{}\mathbf{U}$ and $\tilde{}\mathbf{V}$ serves the purpose of denoising
to leverage the low-rank structure. In this regard, our estimating procedure
can also be viewed as a method of projected moments. The spectral
initialization in Step 1 is very similar to the tensor literature. See, for
instance, Montanari and Richard (2014); Zhang and Xia (2018); Xia and Zhou
(2019) and references therein. A crucial difference here is that the estimate
$\widehat{}\mathbf{v}_{1}$ relies on the estimate $\widehat{}\mathbf{u}_{1}$
to ensure that they are properly correlated in the sense that
$\widehat{}\mathbf{u}_{1}^{\top}\mathbf{M}\widehat{\bm{v}}_{1}$ is bounded
away from zero, which is a critical requirement for the refinement step (9).

Note that the expectation of the RHS of eq. (10), with respect to the
randomness of $\\{\mathbf{X}_{i}^{(4)}\\}_{i=1}^{n_{0}}$, is ${\rm
Tr}(\tilde{}\mathbf{U}^{\top}\mathbf{M}\tilde{}\mathbf{V})\mathbf{M}$.
Therefore, its best rank-$r$ approximation needs to scaled to serve as a valid
estimator for $\mathbf{M}$. The quantity $\widehat{\Lambda}$ is an estimate of
this scaling factor. The performance of the final estimator
$\widehat{}\mathbf{M}$ is guaranteed by the following theorem where we assume
$\mathbf{M}\in{\cal M}(r,d)$ defined in (2) and $d_{1}\asymp d_{2}\asymp d$.

######  Theorem 3.

There exist absolute constants $c_{0},C_{0},C_{1},C_{2},C_{3},C_{4}>0$ such
that if the signal strength $\lambda\geq C_{0}d^{1/2}n^{-1/4}$ and
$\min\\{d,n\\}\geq C_{1}r\log r$, then with probability at least
$1-\exp(-c_{0}(n\wedge r^{-1}d))$,

$\ell(\widehat{}\mathbf{M},\mathbf{M})\leq
C_{2}\left(\sqrt{\frac{dr}{n}}+\frac{1}{\lambda}\sqrt{\frac{d}{n}}\right)$

, if we further assume $\lambda\leq C_{3}\exp(c_{0}(n\wedge r^{-1}d)-\log n)$,
then

$\mathbb{E}\ell(\widehat{}\mathbf{M},\mathbf{M})\leq
C_{4}\left(\sqrt{\frac{dr}{n}}+\frac{1}{\lambda}\sqrt{\frac{d}{n}}\right)$

By Theorem 2 and Theorem 3, we conclude that the estimator
$\widehat{}\mathbf{M}$ can attain the minimax optimal error rate as long as
the signal strength is larger than $d^{1/2}n^{-1/4}$, which we refer to as the
strong signal phase. This is much more stringent than the information-
theoretical limit $d^{1/4}(rn)^{-1/4}+(d/n)^{1/2}$ suggested by the maximum
likelihood estimator and minimax lower bound in Section 2.

##  4 Statistical and computational tradeoffs

Section 2 and Section 3 indicate the existence of a gap in the signal
strengths required by the, in general, computationally infeasible maximum
likelihood estimator and the computationally fast spectral-based estimator.
Gap of this type is usually called the statistical-to-computational gap. In
this section, we provide evidences claiming that no polynomial-time algorithm
can consistently estimate $\mathbf{M}$ if the signal strength is smaller than
$d^{1/2}n^{-1/4}$. Our evidence is built on the low-degree likelihood ratio
framework for hypothesis testing (Kunisky et al., 2019; Löffler et al., 2020;
Hopkins, 2018). This framework delivered convincing evidences justifying the
statistical-to-computational gap for sparse Gaussian mixture model (Löffler et
al., 2020) and tensor PCA model, and demonstrated the sharp phase transitions
for the spiked Wigner matrix model (Kunisky et al., 2019).

The low-degree likelihood ratio framework aims to test two sequences of
hypothesis. For our purpose, consider the following hypothesis testing:

$\displaystyle H_{0}^{(n)}:\mathbf{M}={\bf 0}\quad{\rm versus}\quad
H_{1}^{(n)}:\mathbf{M}\in{\cal M}(1,\lambda)$ (11)

where $n$ denotes the sample size. By observing i.i.d. matrices
$\mathbf{X}_{1},\cdots,\mathbf{X}_{n}$ sampled from the mixture model (1), the
interest is to test whether the data is pure noise or there is a planted low-
rank matrix. Without loss of generality, it suffices to focus on the rank-one
case since the “information” strength $\|\mathbf{M}\|_{\rm F}$ increases if
the rank is larger and, as a result, the hypothesis testing becomes easier for
larger ranks.

Classical textbook results, say, by Neyman-Pearson Lemma, dictate that the
likelihood ratio test has preferable power and is uniformly most powerful
under some scenarios. Direct computation of the likelihood ratio for testing
(11) is rather involved due to the composite hypothesis in $H_{1}^{(n)}$. For
simplicity, under the alternative hypothesis, we impose a prior distribution
on $\mathbf{M}$ assuming that $\mathbf{M}=\lambda\mathbf{u}\mathbf{v}^{\top}$
with a fixed $\lambda$ and the entries of $\mathbf{u}$ and $\mathbf{v}$
independently taking the values $\pm d_{1}^{-1/2}$ and $\pm d_{2}^{-1/2}$,
respectively, with probability $1/2$. Denote $\mathbb{P}_{n}$, treated as the
alternative hypothesis, the distribution of
$(\mathbf{X}_{1},\cdots,\mathbf{X}_{n})$ under LrMM (1) with $\mathbf{M}$
sampled from the aforementioned prior distribution. Note that, for brevity, we
suppress the dependence of $\mathbb{P}_{n}$ on $\lambda$. Let $\mathbb{Q}_{n}$
be the distribution of $(\mathbf{X}_{1},\cdots,\mathbf{X}_{n})$ under the null
hypothesis, i.e., each $\mathbf{X}_{i}$ is sampled from LrMM (1) with
$\mathbf{M}={\bf 0}$. Instead of (11), we consider the following hypothesis
testing

$\displaystyle H_{0}^{(n)}:\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\
\stackrel{{\scriptstyle{\rm i.i.d.}}}{{\sim}}\ \mathbb{Q}_{n}\quad{\rm
versus}\quad H_{1}^{(n)}:\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\
\stackrel{{\scriptstyle{\rm i.i.d.}}}{{\sim}}\ \mathbb{P}_{n}$ (12)

Denote
$L_{n}({\bm{\mathcal{X}}}):=d\mathbb{P}_{n}/d\mathbb{Q}_{n}(\mathbf{X}_{1},\cdots,\mathbf{X}_{n})$
the likelihood ratio, where ${\bm{\mathcal{X}}}\in\mathbb{R}^{d_{1}\times
d_{2}\times n}$ is constructed by stacking $n$ data matrices. A well-
recognized fact is that the two distributions $\mathbb{P}_{n}$ and
$\mathbb{Q}_{n}$ are statistically indistinguishable if
$\|L_{n}\|^{2}:=\mathbb{E}_{\mathbb{Q}_{n}}[L_{n}({\bm{\mathcal{X}}})^{2}]$
remains bounded as $n\to\infty$. Here statistically indistinguishable means
that no test can have both type I and type II error probabilities vanishing
asymptotically.

Let $L_{n}^{\leq D}({\bm{\mathcal{X}}})$ denote the orthogonal projection of
$L_{n}({\bm{\mathcal{X}}})$ onto the linear subspace of polynomials
$\mathbb{R}^{d_{1}\times d_{2}\times n}\mapsto\mathbb{R}$ of degree at most
$D$. Similarly, define $\|L_{n}^{\leq
D}\|^{2}:=\mathbb{E}_{\mathbb{Q}_{n}}[L_{n}^{\leq
D}({\bm{\mathcal{X}}})^{2}]$. At the core of low-degree likelihood ratio
framework is the following conjecture888We note that a recent work Zadik et
al. (2021) introduces a very special counter-example to Conjecture 1. However,
our LrMM is more closely related to the spiked matrix and tensor model where
Conjecture 1 has contributed convincing evidences to the computational
hardness. Therefore, we still postulate the correctness of Conjecture 1 for
our LrMM., adapted to the matrix-variate case for our purpose. Here, a test
$\phi_{n}(\cdot)$ taking value $1$ means rejecting the null hypothesis and
takes value $0$ if the null hypothesis is not rejected.

######  Conjecture 1.

Consider $\mathbb{P}_{n}$ and $\mathbb{Q}_{n}$ defined in (12). If there
exists $\varepsilon>0$ and $D=D_{n}\geq(\log nd)^{1+\varepsilon}$ for which
$\|L_{n}^{\leq D}\|=1+o(1)$, then there is no polynomial-time test
$\phi_{n}:\mathbb{R}^{d_{1}\times d_{2}\times n}\mapsto\\{0,1\\}$ such that
the sum of type-I error and type-II error probabilities

$\mathbb{E}_{\mathbb{Q}_{n}}[\phi_{n}({\bm{\mathcal{X}}})]+\mathbb{E}_{\mathbb{P}_{n}}[1-\phi_{n}({\bm{\mathcal{X}}})]\to
0\quad{\rm as}\quad n\to\infty.$

Basically, Conjecture 1 means that the two distributions $\mathbb{P}_{n}$ and
$\mathbb{Q}_{n}$ are indistinguishable by polynomial-time algorithms if
$\|L_{n}^{\leq D}\|=1+o(1)$. Under the low-degree framework, we now state the
computational lower bound of our signal strength for testing (12).

######  Theorem 4.

Consider $\mathbb{P}_{n}$ and $\mathbb{Q}_{n}$ defined in (12). If
$\lambda=o(d^{1/2}n^{-1/4})$, then $\|L_{n}^{\leq D}\|^{2}=1+o(1)$.

By Theorem 4, conditioned on Conjecture 1, detecting the signal matrix in LrMM
as in (12) becomes computationally hard as long as the signal strength
$\lambda$ is at a smaller order of $d^{1/2}n^{-1/4}$. In principle, the
estimation of signal matrix is at least as hard (computationally) as detection
as in (12), as the latter one only concerns the mere existence thereof, and
hence we would expect, at least, the same lower bound also holds for
estimation problem in LrMM. Notably, if $n=1$, LrMM reduces to the typical
matrix perturbation model (Cai and Zhang, 2018; Xia, 2021) where there exists
no statistical-to-computational gap and the signal strength requirement
$O(d^{1/2})$ is both the statistical and computational limit. Interestingly,
if $n$ is at the order of $d$, the computational hardness occurs at the signal
strength $O(d^{1/4})$ which coincides with the prior literature on spiked
tensor model. See Zhang and Xia (2018); Kunisky et al. (2019) and references
therein.

##  5 Numerical simulations

In this section, we present numerical experiments to confirm our theoretical
findings in the strong signal phase and showcase the performance of our
algorithm. Particularly, we apply the spectral aggregate algorithm on $n$
independent data matrices generated from LrMM model in (1), with a signal
matrix $\mathbf{M}\in\mathbb{R}^{d\times d}$ of rank $r$ constructed as
follows. We first generate two uniformly random $d\times r$ orthonormal
matrices $\mathbf{U}$ and $\mathbf{V}$, say, by computing the column span
(i.e. the image) of a random $d\times r$ Gaussian random matrix with i.i.d.
$\mathcal{N}(0,1)$ entries. Then we fix the smallest and the largest singular
value to be $\lambda_{r}=\lambda$ and $\lambda_{1}=1.5\lambda$, respectively,
and form a diagonal matrix
$\mathbf{\Lambda}=\text{diag}(\lambda_{1},\cdots,\lambda_{r})$ where the
values of the diagonal terms are equally spaced in a decreasing order. Finally
we get our signal matrix
$\mathbf{M}=\mathbf{U}\mathbf{\Lambda}\mathbf{V}^{\top}$. We study the effect
of parameters $(n,d,r,\lambda)$ on the error
$\ell(\widehat{}\mathbf{M},\mathbf{M})$ via varying one/two parameters while
fixing the rest of them. In each experiment (for a given parameter group
$(n,d,r,\lambda)$), the value of the error is the average based on 100
independent simulations with the same signal matrix $\mathbf{M}$. As the aim
of sample splitting step is to facilitate the theoretical analysis, we apply
the spectral aggregation algorithm on all samples without sample splitting in
all numerical experiments. For brevity, Regime 1 is referred to the case when
$\bm{n\leq dr}$, Regime 2 is referred to the case when $\bm{dr\leq
n\leq(dr)^{2}}$ and Regime 3 is referred to the case when
$\bm{n\geq(dr)^{2}}$. The information are summarized as follows:

  * •

Experiment 1: $n=300$, $d=250$, $r=2$ (Regime 1). $\lambda$ is varying from
$3\sqrt{d}n^{-1/4}$ to $10\sqrt{d}n^{-1/4}$.

  * •

Experiment 2: $n=500$, $d=100$, $r=2$ (Regime 2). $\lambda$ is varying from
$3\sqrt{d}n^{-1/4}$ to $10\sqrt{d}n^{-1/4}$.

  * •

Experiment 3: $n=3000$, $d=20$, $r=2$ (Regime 3). $\lambda$ is varying from
$3\sqrt{d}n^{-1/4}$ to $10\sqrt{d}n^{-1/4}$.

  * •

Experiment 4: $d\in\\{100,200\\}$, $r=2$. $n$ is varying from $100$ to $1000$
with $\lambda=3\sqrt{d}n^{-1/4}$.

  * •

Experiment 5: $n\in\\{100,200\\}$, $r=2$. $d$ is varying from $100$ to $500$
with $\lambda=3\sqrt{d}n^{-1/4}$.

  * •

Experiment 6: $n=10000$, $d=10$, $\lambda\in\\{\sqrt{d}n^{-1/4},5\\}$ (Regime
3). $r$ is varying from $2$ to $10$.

In Experiment 1 & 2 (Regime 1 & 2), the error stays almost constant as
$\lambda$ increases. Both cases fall into the strong signal phase and an
optimal rate of $O((dr/n)^{1/2})$ can be attained, suggested by Theorem 3.
While in Experiment 3 (Regime 3) with the same range of $\lambda$, the phase
transition effect is clearly demonstrated in the bottom panel of Figure 2:
when $\lambda$ varies from $C_{1}d^{1/2}n^{-1/4}$ to $C_{2}r^{-1/2}$, the
optimal rate $O(\lambda^{-1}(d/n)^{1/2})$ is linear in $\lambda^{-1}$; when
$\lambda\geq C_{2}r^{-1/2}$, the optimal rate $O((dr/n)^{1/2})$ is again
independent of $\lambda$.
In Experiment 4 & 5, we screen the effect of varying $n$ and $d$, respectively
in Figure 3. As expected, the error becomes smaller as $n$ grows (or $d$
decreases). The linearity between the error rate and $n^{-1/2}$ (or $d^{1/2}$)
can be verified in the right panels, which is in accordance with Theorem 3.
In Experiment 6, we let $r$ vary with other parameters fixed and focus on
Regime 3, which is the most interesting case due to the phase transition
effect in terms of rank $r$. As shown in Figure 5, the error rate
$O(\lambda^{-1}(d/n)^{1/2})$ is constant in $r$ with
$\lambda\in(C_{1}d^{1/2}n^{-1/4},C_{2}r^{-1/2})$ and when $\lambda\geq
C_{2}r^{-1/2}$, the error rate increases with $r$.

Figure 2: Experiments with $\lambda$ varying. Top-left panel: Regime 1; Top-
right panel: Regime 2; Bottom panel: Regime 3.

Figure 3: Experiments with $n$ varying. Left panel:
$\ell(\widehat{}\mathbf{M},\mathbf{M})$ against $n$; Right panel:
$\ell(\widehat{}\mathbf{M},\mathbf{M})$ against $n^{-1/2}$.

Figure 4: Experiments with $d$ varying. Left panel:
$\ell(\widehat{}\mathbf{M},\mathbf{M})$ against $d$; Right panel:
$\ell(\widehat{}\mathbf{M},\mathbf{M})$ against $d^{1/2}$.  Figure 5:
Experiments with $r$ varying. Blue curve $\lambda=d^{1/2}n^{-1/4}$ corresponds
to the case where error rate is of order $\lambda^{-1}(d/n)^{1/2}$; Red curve
$\lambda=5(\geq\max_{r}r^{-1/2})$ corresponds to case where the error rate is
of order $(dr/n)^{1/2}$.

##  6 Real data experiment

We present an application of our algorithm on a real-world dataset, which is a
collection of multiple layers of worldwide food trading networks (De Domenico
et al. (2015)), recording the trade flows of 30 food products between 99
countries. We pre-process the data the same as in Jing et al. (2021) and end
up with a 3-rd order binary tensor ${\bm{\mathcal{X}}}$ of dimension $99\times
99\times 30$. Each layer of this tensor $[{\bm{\mathcal{X}}}]_{\cdot\cdot
i}=\mathbf{X}_{i}$ represents the adjacency matrix of one specific type of
food product $i$, and nodes are different countries/regions which are common
across all layers. As shown in Jing et al. (2021), the layers could be
clustered into two groups, one of which mainly consists of raw or unprocessed
food and another is made of processed food. We adopt this clustering result as
ground truth and assume all layers are generated independently according to
two expected adjacency matrices
$\mathbf{M}_{1},\mathbf{M}_{2}\in\mathbb{R}^{99\times 99}$. Note that though
throughout the paper the noise matrix $\mathbf{Z}_{i}$ is assumed to be
Gaussian , we believe the spectral aggregation can be applied to more general
setting (for instance, observations with sub-gaussian noise). Our goal is to
recover $\mathbf{M}_{1}$ and $\mathbf{M}_{2}$. To make it adapted to our
framework (as mentioned in Section 1), we first construct centered
observations $\widetilde{}\mathbf{X}_{i}=\mathbf{X}_{i}-\bar{}\mathbf{X}$,
where $\bar{}\mathbf{X}=n^{-1}\sum_{i=1}^{n}\mathbf{X}_{i}$ is the sample
average of adjacency matrices over all layers. Here, $\bar{}\mathbf{X}$ serves
as an estimate of $(\mathbf{M}_{1}+\mathbf{M}_{2})/2$. Then we apply the
spectral aggregation algorithm with rank $r=10$ to
$\\{\widetilde{}\mathbf{X}_{i}\\}_{i=1}^{n}$ to get $\widehat{\mathbf{M}}$. It
turns out that the final result is not sensitive to choice of rank $r$.
Finally we can construct
$\widehat{}\mathbf{M}_{1}=\bar{}\mathbf{X}+\widehat{\mathbf{M}}$ and
$\widehat{}\mathbf{M}_{2}=\bar{}\mathbf{X}-\widehat{\mathbf{M}}$. To
appropriately visualize our result, we rearrange the order of columns and rows
of $\widehat{}\mathbf{M}_{1}$ and $\widehat{}\mathbf{M}_{2}$ in the same way
as in Jing et al. (2021), which is based on the community labels estimated by
tensor method therein, in order to have a glance of community structures. In
Figure 6, the mean matrix in the left panel demonstrates a strong trend of
global trading, while the other one shows the dominance of regional trading.
These findings coincides with results in Jing et al. (2021), whereas we are
estimating the difference of two center matrices instead of clustering all
observations. Note that the results in Jing et al. (2021) require layer
clustering before producing
$\widehat{}\mathbf{M}_{1},\widehat{}\mathbf{M}_{2}$ but our method does not.

Figure 6: Heatmaps for $\widehat{}\mathbf{M}_{1}$ and
$\widehat{}\mathbf{M}_{2}$.

##  7 Discussion

Our main focus in this paper is on the optimal estimation and computational
limits for the two-component low-rank Gaussian mixtures. It is of great
interest to investigate the minimax optimal estimation when the number of
components is greater than two. Unfortunately, our spectral aggregation method
is inapplicable and we cannot immediately see an easy generalization of the
maximum likelihood estimator to the multi-component case. There are several
possibilities. For instance, unlike the two-component case, it might be
necessary to, at least partially, recover the latent labels before estimating
the underlying low-rank components. Indeed, the linear regression low-rank
mixture model (Chen et al., 2021) was treated by this way. However, it is well
recognized that consistent clustering often requires a much stronger condition
on the signal strength. See, for instance, Löffler et al. (2019); Wu and Zhou
(2019) and references therein. For the two-component symmetric case as in
model (1), consistent clustering requires a signal strength at least999To see
this, one can simply assume the singular vectors $\mathbf{U}$ and $\mathbf{V}$
are available before hand. in the order of $\Omega(1)$ when $r$ is a constant,
which can be much more stringent than the condition required by the spectral
aggregation method in Regime 3. It therefore indicates another possibility:
there might exist some method that can reliably estimate the multiple low-rank
components without the prerequisite of meaningful clustering. We leave this
for future works.

## References

  * (1) https://www.ncbi.nlm.nih.gov/sites/GDSbrowser?acc=GDS1083.
  * Abramowitz and Stegun (1948) Milton Abramowitz and Irene A Stegun.  _Handbook of mathematical functions with formulas, graphs, and mathematical tables_ , volume 55.  US Government printing office, 1948.
  * Adamczak (2008) Radoslaw Adamczak.  A tail inequality for suprema of unbounded empirical processes with applications to markov chains.  _Electronic Journal of Probability_ , 13:1000–1034, 2008\.
  * Bakshi et al. (2020) Ainesh Bakshi, Ilias Diakonikolas, He Jia, Daniel M Kane, Pravesh K Kothari, and Santosh S Vempala.  Robustly learning mixtures of $k$ arbitrary gaussians.  _arXiv preprint arXiv:2012.02119_ , 2020.
  * Balakrishnan et al. (2017) Sivaraman Balakrishnan, Martin J Wainwright, and Bin Yu.  Statistical guarantees for the em algorithm: From population to sample-based analysis.  _The Annals of Statistics_ , 45(1):77–120, 2017\.
  * Barak and Moitra (2016) Boaz Barak and Ankur Moitra.  Noisy tensor completion via the sum-of-squares hierarchy.  In _Conference on Learning Theory_ , pages 417–445. PMLR, 2016.
  * Belkin and Sinha (2010) Mikhail Belkin and Kaushik Sinha.  Toward learning gaussian mixtures with arbitrary separation.  In _COLT_ , pages 407–419, 2010.
  * Cai and Zhang (2018) T Tony Cai and Anru Zhang.  Rate-optimal perturbation bounds for singular subspaces with applications to high-dimensional statistics.  _The Annals of Statistics_ , 46(1):60–89, 2018\.
  * Cai et al. (2013) T Tony Cai, Zongming Ma, and Yihong Wu.  Sparse pca: Optimal rates and adaptive estimation.  _The Annals of Statistics_ , 41(6):3074–3110, 2013.
  * Cai et al. (2019) T Tony Cai, Jing Ma, and Linjun Zhang.  Chime: Clustering of high-dimensional gaussian mixtures with em algorithm and its optimality.  _The Annals of Statistics_ , 47(3):1234–1267, 2019.
  * Chen (1995) Jiahua Chen.  Optimal rate of convergence for finite mixture models.  _The Annals of Statistics_ , pages 221–233, 1995.
  * Chen et al. (2021) Yanxi Chen, Cong Ma, H Vincent Poor, and Yuxin Chena.  Learning mixtures of low-rank models.  _IEEE Transactions on Information Theory_ , 2021.
  * Cheng et al. (2021) Chen Cheng, Yuting Wei, and Yuxin Chen.  Tackling small eigen-gaps: Fine-grained eigenvector estimation and inference under heteroscedastic noise.  _IEEE Transactions on Information Theory_ , 67(11):7380–7419, 2021.
  * Davies et al. (2021) Sami Davies, Arya Mazumdar, Soumyabrata Pal, and Cyrus Rashtchian.  Lower bounds on the total variation distance between mixtures of two gaussians.  _arXiv preprint arXiv:2109.01064_ , 2021.
  * Davis et al. (2021) Damek Davis, Mateo Diaz, and Kaizheng Wang.  Clustering a mixture of gaussians with unknown covariance.  _arXiv preprint arXiv:2110.01602_ , 2021.
  * De Domenico et al. (2015) Manlio De Domenico, Vincenzo Nicosia, Alexandre Arenas, and Vito Latora.  Structural reducibility of multilayer networks.  _Nature communications_ , 6(1):1–9, 2015.
  * Do (2003) Minh N Do.  Fast approximation of kullback-leibler distance for dependence trees and hidden markov models.  _IEEE signal processing letters_ , 10(4):115–118, 2003.
  * Doss et al. (2020) Natalie Doss, Yihong Wu, Pengkun Yang, and Harrison H Zhou.  Optimal estimation of high-dimensional location gaussian mixtures.  _arXiv preprint arXiv:2002.05818_ , 2020.
  * Gao et al. (2021) Xu Gao, Weining Shen, Liwen Zhang, Jianhua Hu, Norbert J Fortin, Ron D Frostig, and Hernando Ombao.  Regularized matrix data clustering and its application to image analysis.  _Biometrics_ , 77(3):890–902, 2021.
  * Gavish and Donoho (2014) Matan Gavish and David L Donoho.  The optimal hard threshold for singular values is $4/\sqrt{3}$.  _IEEE Transactions on Information Theory_ , 60(8):5040–5053, 2014.
  * Ge et al. (2015) Rong Ge, Qingqing Huang, and Sham M Kakade.  Learning mixtures of gaussians in high dimensions.  In _Proceedings of the forty-seventh annual ACM symposium on Theory of computing_ , pages 761–770, 2015.
  * Genovese and Wasserman (2000) Christopher R Genovese and Larry Wasserman.  Rates of convergence for the gaussian mixture sieve.  _The Annals of Statistics_ , 28(4):1105–1127, 2000.
  * Ghosal and Van Der Vaart (2001) Subhashis Ghosal and Aad W Van Der Vaart.  Entropies and rates of convergence for maximum likelihood and bayes estimation for mixtures of normal densities.  _Annals of Statistics_ , pages 1233–1263, 2001.
  * Golub and van Loan (2013) Gene H. Golub and Charles F. van Loan.  _Matrix Computations_.  JHU Press, fourth edition, 2013.  ISBN 1421407949 9781421407944.  URL http://www.cs.cornell.edu/cv/GVL4/golubandvanloan.htm.
  * Gupta and Nagar (2018) Arjun K Gupta and Daya K Nagar.  _Matrix variate distributions_ , volume 104.  CRC Press, 2018.
  * Ho and Nguyen (2016a) Nhat Ho and XuanLong Nguyen.  Convergence rates of parameter estimation for some weakly identifiable finite mixtures.  _The Annals of Statistics_ , 44(6):2726–2755, 2016a.
  * Ho and Nguyen (2016b) Nhat Ho and XuanLong Nguyen.  On strong identifiability and convergence rates of parameter estimation in finite mixtures.  _Electronic Journal of Statistics_ , 10(1):271–307, 2016b.
  * Hopkins (2018) Samuel Hopkins.  _Statistical inference and the sum of squares method_.  PhD thesis, Cornell University, 2018.
  * Hu et al. (2020) Wei Hu, Weining Shen, Hua Zhou, and Dehan Kong.  Matrix linear discriminant analysis.  _Technometrics_ , 62(2):196–205, 2020.
  * Jing et al. (2021) Bing-Yi Jing, Ting Li, Zhongyuan Lyu, and Dong Xia.  Community detection on mixture multilayer networks via regularized tensor decomposition.  _The Annals of Statistics_ , 49(6):3181–3205, 2021.
  * Koltchinskii and Lounici (2017) Vladimir Koltchinskii and Karim Lounici.  Concentration inequalities and moment bounds for sample covariance operators.  _Bernoulli_ , 23(1):110–133, 2017.
  * Koltchinskii and Xia (2015) Vladimir Koltchinskii and Dong Xia.  Optimal estimation of low rank density matrices.  _J. Mach. Learn. Res._ , 16(53):1757–1792, 2015\.
  * Koltchinskii and Xia (2016) Vladimir Koltchinskii and Dong Xia.  Perturbation of linear forms of singular vectors under gaussian noise.  In _High Dimensional Probability VII_ , pages 397–423. Springer, 2016\.
  * Koltchinskii et al. (2011) Vladimir Koltchinskii, Karim Lounici, and Alexandre B Tsybakov.  Nuclear-norm penalization and optimal rates for noisy low-rank matrix completion.  _The Annals of Statistics_ , 39(5):2302–2329, 2011.
  * Kunisky et al. (2019) Dmitriy Kunisky, Alexander S Wein, and Afonso S Bandeira.  Notes on computational hardness of hypothesis testing: Predictions using the low-degree likelihood ratio.  _arXiv preprint arXiv:1907.11636_ , 2019.
  * Laurent and Massart (2000) Beatrice Laurent and Pascal Massart.  Adaptive estimation of a quadratic functional by model selection.  _Annals of Statistics_ , pages 1302–1338, 2000.
  * Le et al. (2018) Can M Le, Keith Levin, and Elizaveta Levina.  Estimating a network from multiple noisy realizations.  _Electronic Journal of Statistics_ , 12(2):4697–4740, 2018.
  * Ledoux and Talagrand (1991) Michel Ledoux and Michel Talagrand.  _Probability in Banach Spaces: isoperimetry and processes_ , volume 23.  Springer Science & Business Media, 1991.
  * Leroux (1992) Brian G Leroux.  Consistent estimation of a mixing distribution.  _The Annals of Statistics_ , pages 1350–1360, 1992.
  * Löffler et al. (2019) Matthias Löffler, Anderson Y Zhang, and Harrison H Zhou.  Optimality of spectral clustering in the gaussian mixture model.  _arXiv preprint arXiv:1911.00538_ , 2019.
  * Löffler et al. (2020) Matthias Löffler, Alexander S Wein, and Afonso S Bandeira.  Computationally efficient sparse clustering.  _arXiv preprint arXiv:2005.10817_ , 2020.
  * Lyu et al. (2021) Zhongyuan Lyu, Dong Xia, and Yuan Zhang.  Latent space model for higher-order networks and generalized tensor decomposition.  _arXiv preprint arXiv:2106.16042_ , 2021.
  * Ma and Wu (2015) Zongming Ma and Yihong Wu.  Volume ratio, sparsity, and minimaxity under unitarily invariant norms.  _IEEE Transactions on Information Theory_ , 61(12):6939–6956, 2015.
  * Mai et al. (2021) Qing Mai, Xin Zhang, Yuqing Pan, and Kai Deng.  A doubly enhanced em algorithm for model-based tensor clustering.  _Journal of the American Statistical Association_ , pages 1–15, 2021\.
  * Maugis and Michel (2011) Cathy Maugis and Bertrand Michel.  A non asymptotic penalized criterion for gaussian mixture model selection.  _ESAIM: Probability and Statistics_ , 15:41–68, 2011.
  * Moitra and Valiant (2010) Ankur Moitra and Gregory Valiant.  Settling the polynomial learnability of mixtures of gaussians.  In _2010 IEEE 51st Annual Symposium on Foundations of Computer Science_ , pages 93–102. IEEE, 2010.
  * Montanari and Richard (2014) Andrea Montanari and Emile Richard.  A statistical model for tensor pca.  _arXiv preprint arXiv:1411.1076_ , 2014.
  * Paul and Chen (2020) Subhadeep Paul and Yuguo Chen.  Spectral and matrix factorization methods for consistent community detection in multi-layer networks.  _The Annals of Statistics_ , 48(1):230–250, 2020\.
  * Pearson (1894) Karl Pearson.  Contributions to the mathematical theory of evolution.  _Philosophical Transactions of the Royal Society of London. A_ , 185:71–110, 1894.
  * Sanjeev and Kannan (2001) Arora Sanjeev and Ravi Kannan.  Learning mixtures of arbitrary gaussians.  In _Proceedings of the thirty-third annual ACM symposium on Theory of computing_ , pages 247–257, 2001.
  * Sun and Li (2019) Will Wei Sun and Lexin Li.  Dynamic tensor clustering.  _Journal of the American Statistical Association_ , 114(528):1894–1907, 2019.
  * Tosh and Dasgupta (2017) Christopher Tosh and Sanjoy Dasgupta.  Maximum likelihood estimation for mixtures of spherical gaussians is np-hard.  _J. Mach. Learn. Res._ , 18:175–1, 2017.
  * Van de Geer (1993) Sara Van de Geer.  Hellinger-consistency of certain nonparametric maximum likelihood estimators.  _The Annals of Statistics_ , pages 14–44, 1993.
  * Van de Geer (2000) Sara Van de Geer.  _Empirical Processes in M-estimation_ , volume 6.  Cambridge university press, 2000.
  * Van Der Vaart et al. (1996) Aad W Van Der Vaart, Adrianus Willem van der Vaart, Aad van der Vaart, and Jon Wellner.  _Weak convergence and empirical processes: with applications to statistics_.  Springer Science & Business Media, 1996.
  * Vershynin (2010) Roman Vershynin.  Introduction to the non-asymptotic analysis of random matrices.  _arXiv preprint arXiv:1011.3027_ , 2010.
  * Wainwright (2019) Martin J Wainwright.  _High-dimensional statistics: A non-asymptotic viewpoint_ , volume 48.  Cambridge University Press, 2019.
  * Wang et al. (2019) Lu Wang, Zhengwu Zhang, and David Dunson.  Common and individual structure of brain networks.  _The Annals of Applied Statistics_ , 13(1):85–112, 2019.
  * Wang et al. (2017) Xiao Wang, Hongtu Zhu, and Alzheimer’s Disease Neuroimaging Initiative.  Generalized scalar-on-image regression models via total variation.  _Journal of the American Statistical Association_ , 112(519):1156–1168, 2017.
  * Wu and Yang (2020) Yihong Wu and Pengkun Yang.  Optimal estimation of gaussian mixtures via denoised method of moments.  _The Annals of Statistics_ , 48(4):1981–2007, 2020.
  * Wu and Zhou (2019) Yihong Wu and Harrison H Zhou.  Randomly initialized em algorithm for two-component gaussian mixture achieves near optimality in ${O}(\sqrt{n})$ iterations.  _arXiv preprint arXiv:1908.10935_ , 2019.
  * Xia (2021) Dong Xia.  Normal approximation and confidence region of singular subspaces.  _Electronic Journal of Statistics_ , 15(2):3798–3851, 2021.
  * Xia and Zhou (2019) Dong Xia and Fan Zhou.  The sup-norm perturbation of hosvd and low rank tensor denoising.  _The Journal of Machine Learning Research_ , 20(1):2206–2247, 2019.
  * Xia et al. (2021) Dong Xia, Ming Yuan, and Cun-Hui Zhang.  Statistically optimal and computationally efficient low rank tensor completion from noisy entries.  _The Annals of Statistics_ , 49(1):76–99, 2021\.
  * Xu et al. (2016) Ji Xu, Daniel Hsu, and Arian Maleki.  Global analysis of expectation maximization for mixtures of two gaussians.  _arXiv preprint arXiv:1608.07630_ , 2016.
  * Zadik et al. (2021) Ilias Zadik, Min Jae Song, Alexander S Wein, and Joan Bruna.  Lattice-based methods surpass sum-of-squares in clustering.  _arXiv preprint arXiv:2112.03898_ , 2021.
  * Zhang and Xia (2018) Anru Zhang and Dong Xia.  Tensor svd: Statistical and computational limits.  _IEEE Transactions on Information Theory_ , 64(11):7311–7338, 2018.

##  Appendix A Proofs for main results

###  A.1 Proof of Theorem 1

For technical reasons discussed in Section 2, we split our proof into two
cases, corresponding to the first and second statement in Theorem 1.

#### Case 1: $dr\log(nd)<n$

In this regime, the standard tool to establish the convergence rate of MLE is
applicable. To this end, we need to introduce the following notations. Define

$\bar{\mathcal{P}}_{d_{1},d_{2}}(r,\lambda):=\left\\{\frac{p_{\mathbf{M}}+p_{\mathbf{M}^{\prime}}}{2}:\mathbf{M}^{\prime}\in\mathcal{M}_{d_{1},d_{2}}(r,\lambda)\right\\},\quad\bar{\mathcal{P}}^{1/2}_{d_{1},d_{2}}(r,\lambda):=\left\\{p^{\frac{1}{2}}:p\in\bar{\mathcal{P}}_{d_{1},d_{2}}(r,\lambda)\right\\}$

and for any small $\delta>0$, define a Hellinger ball centered at
$p_{\mathbf{M}}$ with radius $\delta$ by

$\bar{\mathcal{P}}^{1/2}_{d_{1},d_{2}}(r,\lambda,\delta):=\left\\{\bar{p}^{\frac{1}{2}}\in\bar{\mathcal{P}}^{1/2}_{d}(\lambda):d_{\textsf{H}}\left(\bar{p},p_{\mathbf{M}}\right)\leq\delta\right\\}$

We refer to
$H_{B}(\epsilon,\bar{\mathcal{P}}^{1/2}_{d_{1},d_{2}}(r,\lambda,\delta),L_{2}(\mu))$
as the $\epsilon$-bracketing entropy of
$\bar{\mathcal{P}}^{1/2}_{d_{1},d_{2}}(r,\lambda,\delta)$ under $L_{2}(\mu)$
metric with Lebesgue measure $\mu$ and view
$\mathcal{J}_{B}(\delta,\bar{\mathcal{P}}^{1/2}_{d_{1},d_{2}}(r,\lambda,\delta),L_{2}(\mu))$
as the entropy integral of
$\bar{\mathcal{P}}^{1/2}_{d_{1},d_{2}}(r,\lambda,\delta)$, which is defined as

$\mathcal{J}_{B}(\delta,\bar{\mathcal{P}}^{1/2}_{d_{1},d_{2}}(r,\lambda,\delta),L_{2}(\mu)):=\int_{\delta^{2}/2^{13}}^{\delta}H_{B}^{1/2}(\epsilon,\bar{\mathcal{P}}^{1/2}_{d_{1},d_{2}}(r,\lambda,\delta),L_{2}(\mu))d\epsilon\vee\delta$

Now we state Theorem 7.4 in Van de Geer (2000) (adapted to our notation),
which establishes the rate of convergence of MLE.

######  Lemma 3 (Van de Geer (2000)).

Take
$\Psi(\delta)\geq\mathcal{J}_{B}(\delta,\bar{\mathcal{P}}^{1/2}_{d_{1},d_{2}}(r,\lambda,\delta),L_{2}(\mu))$
in such a way that $\Psi(\delta)/\delta^{2}$ is a non-increasing function of
$\delta$. Then for a universal constant $c$, and for

$\sqrt{n}\delta_{n}^{2}\geq c\Psi(\delta_{n})$

we have for all $\delta\geq\delta_{n}$

$\mathbb{P}(d_{\textsf{H}}(p_{\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}}},p_{\mathbf{M}})>\delta)\leq
c\exp\left(-\frac{n\delta^{2}}{c^{2}}\right)$

A combination of Lemma 1 and Lemma 3 implies that the convergence rate of
$\widehat{}\mathbf{M}_{\textsf{\tiny MLE}}$ would entail an upper bound on the
$\epsilon$-bracketing entropy
$H_{B}(\epsilon,\bar{\mathcal{P}}^{1/2}_{d_{1},d_{2}}(r,\lambda,\delta),L_{2}(\mu))$.
Notice that for any $\delta>0$,

$\displaystyle
H_{B}(\epsilon,\bar{\mathcal{P}}^{1/2}_{d_{1},d_{2}}(r,\lambda,\delta),L_{2}(\mu))$
$\displaystyle\overset{(a)}{\leq}H_{B}(\epsilon,\bar{\mathcal{P}}^{1/2}_{d_{1},d_{2}}(r,\lambda),L_{2}(\mu))\overset{(b)}{=}H_{B}({\epsilon}/{\sqrt{2}},\bar{\mathcal{P}}_{d_{1},d_{2}}(r,\lambda),d_{\textsf{H}})$
$\displaystyle\overset{(c)}{\leq}H_{B}({\epsilon},{\mathcal{P}}_{d_{1},d_{2}}(r,\lambda),d_{\textsf{H}})$
(13)

where (a) is due to
$\bar{\mathcal{P}}^{1/2}_{d_{1},d_{2}}(r,\lambda,\delta)\subset\bar{\mathcal{P}}^{1/2}_{d_{1},d_{2}}(r,\lambda)$,
(b) follows from the definition of Hellinger distance $d_{\textsf{H}}$ and (c)
is due to the following fact (cf. Lemma 4.2 in Van de Geer (2000), Ho and
Nguyen (2016a)): for any
$\bar{p}_{1}=\frac{1}{2}(p_{\mathbf{M}_{1}}+p_{\mathbf{M}})\in\bar{\mathcal{P}}_{d_{1},d_{2}}(r,\lambda),\bar{p}_{2}=\frac{1}{2}(p_{\mathbf{M}_{2}}+p_{\mathbf{M}})\in\bar{\mathcal{P}}_{d_{1},d_{2}}(r,\lambda)$

$d_{\textsf{H}}^{2}(\bar{p}_{1},\bar{p}_{2})\leq\frac{1}{2}d_{\textsf{H}}^{2}(p_{\mathbf{M}_{1}},p_{\mathbf{M}_{2}})$

In view of (A.1), it suffices to bound
$H_{B}(\epsilon,\mathcal{P}_{d_{1},d_{2}}(r,\lambda),d_{\textsf{H}})$. The
following lemma characterize the size of bracketing entropy of
$\mathcal{P}_{d_{1},d_{2}}(r,\lambda)$.

######  Lemma 4.

Assume $d_{1}\asymp d_{2}\asymp d$ then we have

$H_{B}(\epsilon,\mathcal{P}_{d_{1},d_{2}}(r,\lambda),d_{\textsf{H}})\lesssim
dr\log\left(\frac{d}{\epsilon}\right)$

Using relation (A.1) and Lemma 4 we can arrive at

$\mathcal{J}_{B}(\delta,\bar{\mathcal{P}}^{1/2}_{d_{1},d_{2}}(r,\lambda,\delta),L_{2}(\mu))\lesssim\int_{\delta^{2}/2^{13}}^{\delta}\sqrt{dr\log\left(\frac{d}{\epsilon}\right)}d\epsilon\vee\delta\lesssim\delta\sqrt{dr\log\left(\frac{d}{\delta}\right)}$

Now we can take
$\Psi(\delta)=C\delta\sqrt{dr\log\left(\frac{d}{\delta}\right)}$ for some
absolute constant $C>0$ and $\delta=\delta_{n}=\sqrt{\frac{dr}{n}\log(nd)}$,
then we have
$\Psi(\delta)/\delta^{2}=C\frac{1}{\delta}\sqrt{dr\log\left(\frac{d}{\delta}\right)}$
is a non-increasing function of $\delta$ and that

$\sqrt{n}\delta_{n}^{2}=\frac{dr}{\sqrt{n}}\log(nd)\geq
c\frac{dr}{\sqrt{n}}\sqrt{\log(nd)}\sqrt{\log\left(\frac{\sqrt{nd}}{\sqrt{r\log(nd)}}\right)}=c\Psi(\delta_{n})$

By Lemma 3, with probability at least $1-\exp(-cd\log^{2}(nd))$ we have

$d_{\textsf{H}}(p_{\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}}},p_{\mathbf{M}})\leq C\sqrt{\frac{dr\log(nd)}{n}}$

It suffices to use Lemma 1 to connect the density estimation and parameter
estimation. Notice that $\|\widehat{}\mathbf{M}_{\textsf{\tiny MLE}}\|_{\rm
F}+\|\mathbf{M}\|_{\rm F}\asymp\lambda\sqrt{r}$. By Lemma 1, if
$\lambda\sqrt{r}\lesssim 1$, with probability at least
$1-\exp(-cd\log^{2}(nd))$:

$\ell(\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}},\mathbf{M})\lesssim(\lambda\sqrt{r})^{-1}\cdot
d_{\textsf{H}}(p_{\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}}},p_{\mathbf{M}})\leq\frac{1}{\lambda}\sqrt{\frac{d\log(nd)}{n}}$

If $\lambda\sqrt{r}\gtrsim 1$, note that in this case ($dr<n\log(nd)$), we
have with probability at least $1-\exp(-cd\log^{2}(nd))$:

$\min\\{1,\ell(\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}},\mathbf{M})\\}\lesssim
d_{\textsf{H}}(p_{\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}}},p_{\mathbf{M}})\leq\sqrt{\frac{dr\log(nd)}{n}}<1$

implying that $\ell(\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}},\mathbf{M})\lesssim\sqrt{{dr\log(nd)}/{n}}$. Combining two pieces we
conclude that with probability at least $1-\exp(-cd\log^{2}(nd))$:

$\ell(\widehat{}\mathbf{M}_{\textsf{\tiny MLE}},\mathbf{M})\leq
C\left(\sqrt{\frac{dr\log(nd)}{n}}\vee\frac{1}{\lambda}\sqrt{\frac{d\log(nd)}{n}}\right)$

We can further have a bound in expectation:

$\mathbb{E}\ell(\widehat{}\mathbf{M}_{\textsf{\tiny MLE}},\mathbf{M})\leq
C\left(\sqrt{\frac{dr\log(nd)}{n}}\vee\frac{1}{\lambda}\sqrt{\frac{d\log(nd)}{n}}\right)$

provided that $\lambda\leq\exp(cd\log^{2}(nd))$.

#### Case 2: $dr\log(nd)\geq n$

In this regime, our ultimate goal is to have
$\|\widehat{}\mathbf{M}_{\textsf{\tiny MLE}}-\mathbf{M}\|_{\rm
F}\lesssim\sqrt{dr\log(nd)/n}$ with high probability and in expectation and
hence we can assume $\|\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}}-\mathbf{M}\|_{\rm F}\geq c_{0}\sqrt{dr\log(nd)/n}$ for some absolute
constant $c_{0}>0$ (otherwise we have the desired result). Without loss of
generality, we assume $\|\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}}-\mathbf{M}\|_{\rm F}\leq\|\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}}+\mathbf{M}\|_{\rm F}$. Unlike Case 1, we resort to KL divergence instead
of Hellinger distance to establish the convergence rate. Let $P_{\mathbf{M}}$
denote the distribution of (1) and recall the definition of KL divergence, for
any $\mathbf{M},\mathbf{M}^{\prime}\in\mathcal{M}_{d_{1},d_{2}}(r,\lambda)$ we
have

$\displaystyle
D_{\textsf{KL}}\left(p_{\mathbf{M}}\big{\|}p_{\mathbf{M}^{\prime}}\right)=\int\left(\log\frac{p_{\mathbf{M}}}{p_{\mathbf{M}^{\prime}}}\right)dP_{\mathbf{M}}$

Note that for fixed $\mathbf{M}$ and $\mathbf{M}^{\prime}$, we simply have
$D_{\textsf{KL}}\left(p_{\mathbf{M}}\big{\|}p_{\mathbf{M}^{\prime}}\right)=\mathbb{E}\log\left({p_{\mathbf{M}}}(\mathbf{X})/{p_{\mathbf{M}^{\prime}}(\mathbf{X})}\right)$
for $\mathbf{X}\sim p_{\mathbf{M}}$. On the other hand, by the definition of
the maximum likelihood estimator $\widehat{}\mathbf{M}_{\textsf{\tiny MLE}}$,
we have

$\frac{1}{n}\sum_{i=1}^{n}\log\frac{p_{\mathbf{M}}(\mathbf{X}_{i})}{p_{\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}}}(\mathbf{X}_{i})}\leq 0$

Therefore, we can have that

$\displaystyle
D_{\textsf{KL}}\left(p_{\mathbf{M}}\big{\|}p_{\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}}}\right)\leq-\frac{1}{n}\sum_{i=1}^{n}\log\frac{p_{\mathbf{M}}(\mathbf{X}_{i})}{p_{\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}}}(\mathbf{X}_{i})}+D_{\textsf{KL}}\left(p_{\mathbf{M}}\big{\|}p_{\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}}}\right)$ (14)

Now we give an upper bound of RHS of (14). To this end, we consider a ball in
$\mathcal{M}_{d_{1},d_{2}}(r,\lambda)$ with radius $\delta$, i.e.,
$\mathcal{M}(\delta):=\\{\mathbf{M}^{\prime}\in\mathcal{M}_{d_{1},d_{2}}(r,\lambda):\|\mathbf{M}^{\prime}-\mathbf{M}\|_{\rm
F}\leq\delta\\}$. Our aim is to bound the following quantity:

$\theta_{n}(\delta):=\sup_{\mathbf{M}^{\prime}\in\mathcal{M}(\delta)}\left|\frac{1}{n}\sum_{i=1}^{n}\log\frac{p_{\mathbf{M}}(\mathbf{X}_{i})}{p_{\mathbf{M}^{\prime}}(\mathbf{X}_{i})}-D_{\textsf{KL}}\left(p_{\mathbf{M}}\big{\|}p_{\mathbf{M}^{\prime}}\right)\right|$

Observe that

$\displaystyle\log\frac{p_{\mathbf{M}}(\mathbf{X})}{p_{\mathbf{M}^{\prime}}(\mathbf{X})}$
$\displaystyle=\log\left(\frac{e^{-\frac{1}{2}\|\mathbf{X}-\mathbf{M}\|_{\rm
F}^{2}}+e^{-\frac{1}{2}\|\mathbf{X}+\mathbf{M}\|_{\rm
F}^{2}}}{e^{-\frac{1}{2}\|\mathbf{X}-\mathbf{M}^{\prime}\|_{\rm
F}^{2}}+e^{-\frac{1}{2}\|\mathbf{X}+\mathbf{M}^{\prime}\|_{\rm
F}^{2}}}\right)=\frac{1}{2}\|\mathbf{M}^{\prime}\|_{\rm
F}^{2}-\frac{1}{2}\|\mathbf{M}\|_{\rm
F}^{2}+\log\left(\frac{e^{\langle\mathbf{X},\mathbf{M}\rangle}+e^{-\langle\mathbf{X},\mathbf{M}\rangle}}{e^{\langle\mathbf{X},\mathbf{M}^{\prime}\rangle}+e^{-\langle\mathbf{X},\mathbf{M}^{\prime}\rangle}}\right)$

By log-sum-exp inequality, we have

$\displaystyle|\langle\mathbf{X},\mathbf{M}\rangle|-|\langle\mathbf{X},\mathbf{M}^{\prime}\rangle|-\log
2\leq\log\left(\frac{e^{\langle\mathbf{X},\mathbf{M}\rangle}+e^{-\langle\mathbf{X},\mathbf{M}\rangle}}{e^{\langle\mathbf{X},\mathbf{M}^{\prime}\rangle}+e^{-\langle\mathbf{X},\mathbf{M}^{\prime}\rangle}}\right)\leq\log
2+|\langle\mathbf{X},\mathbf{M}\rangle|-|\langle\mathbf{X},\mathbf{M}^{\prime}\rangle|$

Hence we have

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\log\frac{p_{\mathbf{M}}(\mathbf{X}_{i})}{p_{\mathbf{M}^{\prime}}(\mathbf{X}_{i})}-D_{\textsf{KL}}\left(p_{\mathbf{M}}\big{\|}p_{\mathbf{M}^{\prime}}\right)\leq
2\log
2+\frac{1}{n}\sum_{i=1}^{n}\left[|\langle\mathbf{X}_{i},\mathbf{M}\rangle|-|\langle\mathbf{X}_{i},\mathbf{M}^{\prime}\rangle|\right]-\mathbb{E}\left[|\langle\mathbf{X},\mathbf{M}\rangle|-|\langle\mathbf{X},\mathbf{M}^{\prime}\rangle|\right]$

which implies $\theta_{n}(\delta)\leq 2\log 2+\tilde{\theta}_{n}(\delta)$,
where

$\displaystyle\tilde{\theta}_{n}(\delta):=\sup_{\mathbf{M}^{\prime}\in\mathcal{M}(\delta)}\left|\frac{1}{n}\sum_{i=1}^{n}\left[|\langle\mathbf{X}_{i},\mathbf{M}\rangle|-|\langle\mathbf{X}_{i},\mathbf{M}^{\prime}\rangle|\right]-\mathbb{E}\left[|\langle\mathbf{X},\mathbf{M}\rangle|-|\langle\mathbf{X},\mathbf{M}^{\prime}\rangle|\right]\right|$

To get a high probability bound for $\tilde{\theta}_{n}(\delta)$, we first
upper bound its expectation. By symmetrization (see, e.g., in (Van Der Vaart
et al., 1996, Lemma 2.3.1)), we have

$\displaystyle\mathbb{E}\tilde{\theta}_{n}(\delta)\leq
2\mathbb{E}\left(\sup_{\mathbf{M}^{\prime}\in\mathcal{M}(\delta)}\left|\frac{1}{n}\sum_{i=1}^{n}\varepsilon_{i}\left[|\langle\mathbf{X}_{i},\mathbf{M}\rangle|-|\langle\mathbf{X}_{i},\mathbf{M}^{\prime}\rangle|\right]\right|\right)$

where $\\{\varepsilon_{i}\\}_{i=1}^{n}$ are independent Rademacher random
variables, which is independent of $\\{\mathbf{X}_{i}\\}_{i=1}^{n}$. Denote
$\phi_{i}(\mathbf{M}^{\prime})=|\langle\mathbf{X}_{i},\mathbf{M}\rangle|-|\langle\mathbf{X}_{i},\mathbf{M}^{\prime}\rangle|$,
for any $\mathbf{M}_{1},\mathbf{M}_{2}\in\mathcal{M}(\delta)$ we have

$|\phi_{i}(\mathbf{M}_{1})-\phi_{i}(\mathbf{M}_{2})|\leq|\langle\mathbf{X}_{i},\mathbf{M}_{1}-\mathbf{M}_{2}\rangle|=|\langle\mathbf{X}_{i},\mathbf{M}_{1}-\mathbf{M}\rangle-\langle\mathbf{X}_{i},\mathbf{M}_{2}-\mathbf{M}\rangle|$

which means $\phi_{i}(\mathbf{M}^{\prime})$ is $1$-Lipschitz in
$\langle\mathbf{X}_{i},\mathbf{M}^{\prime}-\mathbf{M}\rangle$. By comparison
theorem ((Ledoux and Talagrand, 1991, Theorem 4.12)), we deduce that

$\displaystyle\mathbb{E}\left(\sup_{\mathbf{M}^{\prime}\in\mathcal{M}(\delta)}\left|\frac{1}{n}\sum_{i=1}^{n}\varepsilon_{i}\left[|\langle\mathbf{X}_{i},\mathbf{M}\rangle|-|\langle\mathbf{X}_{i},\mathbf{M}^{\prime}\rangle|\right]\right|\right)\leq\mathbb{E}\left(\sup_{\mathbf{M}^{\prime}\in\mathcal{M}(\delta)}\left|\frac{1}{n}\sum_{i=1}^{n}\varepsilon_{i}\langle\mathbf{X}_{i},\mathbf{M}^{\prime}-\mathbf{M}\rangle\right|\right)$

Hence we proceed as

$\displaystyle\mathbb{E}\tilde{\theta}_{n}(\delta)$ $\displaystyle\leq
2\mathbb{E}\left(\sup_{\mathbf{M}^{\prime}\in\mathcal{M}(\delta)}\left|\frac{1}{n}\sum_{i=1}^{n}\varepsilon_{i}\langle\mathbf{X}_{i},\mathbf{M}^{\prime}-\mathbf{M}\rangle\right|\right)=2\mathbb{E}\left(\sup_{\mathbf{M}^{\prime}\in\mathcal{M}(\delta)}\left|\frac{1}{n}\sum_{i=1}^{n}\varepsilon_{i}\langle
s_{i}\mathbf{M}+\mathbf{Z}_{i},\mathbf{M}^{\prime}-\mathbf{M}\rangle\right|\right)$
$\displaystyle\leq
2\mathbb{E}\left(\sup_{\mathbf{M}^{\prime}\in\mathcal{M}(\delta)}\left|\left\langle\frac{1}{n}\sum_{i=1}^{n}\varepsilon_{i}\mathbf{M},\mathbf{M}^{\prime}-\mathbf{M}\right\rangle\right|\right)+2\mathbb{E}\left(\sup_{\mathbf{M}^{\prime}\in\mathcal{M}(\delta)}\left|\left\langle\frac{1}{n}\sum_{i=1}^{n}\varepsilon_{i}\mathbf{Z}_{i},\mathbf{M}^{\prime}-\mathbf{M}\right\rangle\right|\right)$
$\displaystyle\overset{(a)}{\leq}2\delta\|\mathbf{M}\|_{\rm
F}\mathbb{E}\left|\frac{1}{n}\sum_{i=1}^{n}\varepsilon_{i}\right|+2\sqrt{2r}\delta\mathbb{E}\left\|\frac{1}{n}\sum_{i=1}^{n}\mathbf{Z}_{i}\right\|$
$\displaystyle\overset{(b)}{\leq}2\delta\lambda\sqrt{\frac{r}{n}}+2\sqrt{2}\delta\sqrt{\frac{dr}{n}}\overset{(c)}{\lesssim}\delta\sqrt{\frac{dr}{n}}$

where in (a) we’ve used
$\|\mathbf{M}^{\prime}-\mathbf{M}\|_{*}\leq\text{rank}(\mathbf{M}^{\prime}-\mathbf{M})\cdot\|\mathbf{M}^{\prime}-\mathbf{M}\|_{\rm
F}\leq\sqrt{2r}\|\mathbf{M}^{\prime}-\mathbf{M}\|_{\rm F}$, in (b) we have a
simple bound for
$\mathbb{E}\left|\frac{1}{n}\sum_{i=1}^{n}\varepsilon_{i}\right|\leq
1/\sqrt{n}$ by Jensen’s inequality, and (c) is due to the assumption
$\lambda\lesssim\sqrt{d}$. Define
$\sigma^{2}:=\sup_{\mathbf{M}^{\prime}\in\mathcal{M}(\delta)}\sum_{i=1}^{n}\mathbb{E}\left[|\langle\mathbf{X}_{i},\mathbf{M}\rangle|-|\langle\mathbf{X}_{i},\mathbf{M}^{\prime}\rangle|\right]^{2}$,
notice that

$\displaystyle\left[|\langle\mathbf{X}_{i},\mathbf{M}\rangle|-|\langle\mathbf{X}_{i},\mathbf{M}^{\prime}\rangle|\right]^{2}$
$\displaystyle\leq|\langle\mathbf{X}_{i},\mathbf{M}-\mathbf{M}^{\prime}\rangle|^{2}=|\langle
s_{i}\mathbf{M}+\mathbf{Z}_{i},\mathbf{M}-\mathbf{M}^{\prime}\rangle|^{2}$
$\displaystyle=\langle\mathbf{M},\mathbf{M}-\mathbf{M}^{\prime}\rangle^{2}+\langle\mathbf{Z}_{i},\mathbf{M}-\mathbf{M}^{\prime}\rangle^{2}+2s_{i}\langle\mathbf{M},\mathbf{M}-\mathbf{M}^{\prime}\rangle\langle\mathbf{Z}_{i},\mathbf{M}-\mathbf{M}^{\prime}\rangle$

Observe that
$\langle\mathbf{Z}_{i},\mathbf{M}-\mathbf{M}^{\prime}\rangle\sim\mathcal{N}(0,\|\mathbf{M}-\mathbf{M}^{\prime}\|_{\rm
F}^{2})$, we have $\sigma^{2}\leq n\|\mathbf{M}\|_{\rm
F}^{2}\delta^{2}+n\delta^{2}\lesssim n\lambda^{2}r\delta^{2}$, the last
inequality is due to $\|\mathbf{M}\|_{\rm F}\geq\lambda\sqrt{r}\gtrsim 1$ in
this regime. Moreover, by Lemma 2.2.2 in Van Der Vaart et al. (1996), we have

$\displaystyle\left\|\max_{i}\sup_{\mathbf{M}^{\prime}\in\mathcal{M}(\delta)}\left|\phi_{i}(\mathbf{M}^{\prime})-\mathbb{E}\phi_{i}(\mathbf{M}^{\prime})\right|\right\|_{\psi_{2}}\lesssim\sqrt{\log
n}\max_{i}\left\|\sup_{\mathbf{M}^{\prime}\in\mathcal{M}(\delta)}\left|\phi_{i}(\mathbf{M}^{\prime})-\mathbb{E}\phi_{i}(\mathbf{M}^{\prime})\right|\right\|_{\psi_{2}}$

It suffices to note that for each $i\in[n]$,

$\displaystyle\sup_{\mathbf{M}^{\prime}\in\mathcal{M}(\delta)}\left|\phi_{i}(\mathbf{M}^{\prime})-\mathbb{E}\phi_{i}(\mathbf{M}^{\prime})\right|$
$\displaystyle\leq\sup_{\mathbf{M}^{\prime}\in\mathcal{M}(\delta)}\left|\langle\mathbf{X}_{i},\mathbf{M}^{\prime}-\mathbf{M}\rangle\right|+\sup_{\mathbf{M}^{\prime}\in\mathcal{M}(\delta)}\mathbb{E}\left|\langle\mathbf{X}_{i},\mathbf{M}^{\prime}-\mathbf{M}\rangle\right|$
$\displaystyle\leq\delta\|\mathbf{M}\|_{\rm
F}+\delta\sqrt{r}\|\mathbf{Z}_{i}\|+\mathbb{E}\left(\delta\|\mathbf{M}\|_{\rm
F}+\delta\sqrt{r}\|\mathbf{Z}_{i}\|\right)$
$\displaystyle\lesssim\delta\sqrt{dr}+\delta\sqrt{r}\|\mathbf{Z}_{i}\|$

where in the last inequality we’ve used
$\mathbb{E}\|\mathbf{Z}_{i}\|\lesssim\sqrt{d}$ and $\lambda\lesssim\sqrt{d}$.
By random matrix theory, we know $\mathbb{E}\|\mathbf{Z}_{i}\|\asymp\sqrt{d}$
and $\|\mathbf{Z}_{i}\|-\mathbb{E}\|\mathbf{Z}_{i}\|$ is sub-gaussian, then
$\|\|\mathbf{Z}_{i}\|\|_{\psi_{2}}\lesssim\sqrt{d}$. Hence

$\displaystyle\left\|\max_{i}\sup_{\mathbf{M}^{\prime}\in\mathcal{M}(\delta)}\left|\phi_{i}(\mathbf{M}^{\prime})-\mathbb{E}\phi_{i}(\mathbf{M}^{\prime})\right|\right\|_{\psi_{2}}\lesssim\delta\sqrt{dr\log
n}+\delta\sqrt{r\log
n}\max_{i}\left\|\|\mathbf{Z}_{i}\|\right\|_{\psi_{2}}\lesssim\delta\sqrt{dr\log
n}$

Now we can invoke concentration inequality for suprema of empirical processes
of unbounded functions (Adamczak, 2008, Theorem 4), we have for any $t\geq 0$:

$\displaystyle\mathbb{P}\left(\tilde{\theta}_{n}(\delta)\geq
C\delta\sqrt{\frac{dr\log
n}{n}}+\delta\sqrt{\frac{dr}{n}}t\right)\leq\exp\left(-ct^{2}\right)$ (15)

Note that (15) only holds for any given $\delta>0$. Consider any
$\delta\in[\sqrt{dr/n},2\sqrt{dr}]$, let $\delta_{j}=2^{j}\sqrt{dr/n}$ for
$j=0,1,\cdots,k^{*}+1$ with $k^{*}:=\lfloor\log_{2}(2\sqrt{n})\rfloor$, then
$\delta\in\bigcup_{j=1}^{k^{*}}[\delta_{j},\delta_{j+1}]$. By construction,
for any $\delta\in[\delta_{j},\delta_{j+1}]$, we have
$\delta\asymp\delta_{j}\asymp\delta_{j+1}$. Hence for fixed $j$, (15) holds
for any $\delta\in[\delta_{j},\delta_{j+1}]$ up to change in constants $c>0$
and $C>0$. Take a union bound over all $j=0,1,\cdots,k^{*}+1$, we have (15)
holds for any $\delta\in[\sqrt{dr/n},2\sqrt{dr}]$ and any
$t\geq\log(k^{*}+2)\gtrsim\log\log n$. Combined with (14), we conclude that

$\displaystyle
D_{\textsf{KL}}\left(p_{\mathbf{M}}\big{\|}p_{\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}}}\right)$
$\displaystyle\leq\left|\frac{1}{n}\sum_{i=1}^{n}\log\frac{p_{\mathbf{M}}(\mathbf{X}_{i})}{p_{\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}}}(\mathbf{X}_{i})}-D_{\textsf{KL}}\left(p_{\mathbf{M}}\big{\|}p_{\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}}}\right)\right|\leq\tilde{\theta}_{n}(\|\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}}-\mathbf{M}\|_{\rm F})+2\log 2$ $\displaystyle\leq
C\|\widehat{}\mathbf{M}_{\textsf{\tiny MLE}}-\mathbf{M}\|_{\rm
F}\sqrt{\frac{dr\log(nd)}{n}}$ (16)

where the last inequality holds with probability at least $1-(nd)^{-4}$, due
to the facts that $\|\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}}-\mathbf{M}\|_{\rm F}\gtrsim\sqrt{dr\log(nd)/n}\geq 1$. Since
$\|\widehat{}\mathbf{M}_{\textsf{\tiny MLE}}-\mathbf{M}\|_{\rm F}\gtrsim 1$
and $\|\mathbf{M}\|_{\rm F}\gtrsim 1$ in this regime, it turns out that we can
apply Lemma 2 to get a lower bound of
$D_{\textsf{KL}}\left(p_{\mathbf{M}}\big{\|}p_{\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}}}\right)$, hence we have with probability at least $1-(nd)^{-4}$ that

$\displaystyle\|\widehat{}\mathbf{M}_{\textsf{\tiny MLE}}-\mathbf{M}\|_{\rm
F}\leq C\sqrt{\frac{dr\log(nd)}{n}}$

Finally, we can have a bound in expectation

$\displaystyle\mathbb{E}\|\widehat{}\mathbf{M}_{\textsf{\tiny
MLE}}-\mathbf{M}\|_{\rm F}\leq C\sqrt{\frac{dr\log(nd)}{n}}$

given that $\lambda\lesssim\sqrt{d}$. ∎

###  A.2 Proof of Theorem 3

In the proof, we consider the conditional model with sample splitting, i.e.,
$\mathbf{X}_{i}^{(k)}\overset{d}{=}s_{i}^{(k)}\mathbf{M}+\mathbf{Z}_{i}^{(k)}$
for $k=1,2,3,4$ and $i\in[n_{0}]$. Let
$\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$ denote the thin SVD of the signal
matrix $\mathbf{M}$ and recall that $d_{1}\asymp d_{2}\asymp d$.

#### Step 1:

Denote
$\mathbf{X}_{f}=[\mathbf{X}_{1}^{(1)},\cdots,\mathbf{X}_{n_{0}}^{(1)}]\in\mathbb{R}^{d\times{n_{0}}d}$.
A key observation is that $\widehat{}\mathbf{u}_{1}$ is also the leading
eigenvector of
$\frac{1}{n_{0}}\mathbf{X}_{f}\mathbf{X}_{f}^{\top}-d\mathbf{I}_{d}$. Then we
have

$\displaystyle\frac{1}{n_{0}}\mathbf{X}_{f}\mathbf{X}_{f}^{\top}-d\mathbf{I}_{d}$
$\displaystyle=\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\mathbf{X}_{i}^{(1)}\mathbf{X}_{i}^{(1)\top}-d\mathbf{I}_{d}=\mathbf{M}\mathbf{M}^{\top}+\Delta$
(17)

where

$\Delta:=\mathbf{M}\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}s_{i}^{(1)}\mathbf{Z}_{i}^{(1)\top}\right)+\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}s_{i}^{(1)}\mathbf{Z}_{i}^{(1)}\right)\mathbf{M}^{\top}+\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\mathbf{Z}_{i}^{(1)}\mathbf{Z}_{i}^{(1)\top}-d\mathbf{I}_{d_{1}}$

Note that $\sum_{i=1}^{n_{0}}s_{i}^{(1)}\mathbf{Z}_{i}^{(1)}$
($\sum_{i=1}^{n_{0}}s_{i}^{(1)}\mathbf{Z}_{i}^{(1)\top}$) is a $d_{1}\times
d_{2}$ ($d_{2}\times d_{1}$) matrix of independent centered Gaussian entries
with variance $n_{0}$, then by random matrix theory (e.g. Vershynin (2010))
with probability at least $1-\exp(-cd)$, we have:

$\left\|\mathbf{M}\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}s_{i}^{(1)}\mathbf{Z}_{i}^{(1)\top}\right)\right\|\lesssim\lambda_{1}\sqrt{\frac{d}{n}},\quad\left\|\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}s_{i}^{(1)}\mathbf{Z}_{i}^{(1)}\right)\mathbf{M}^{\top}\right\|\lesssim\lambda_{1}\sqrt{\frac{d}{n}}$

Furthermore, since

$\displaystyle\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\mathbf{Z}_{i}^{(1)}\mathbf{Z}_{i}^{(1)\top}-d\mathbf{I}_{d_{1}}=d\left(\frac{1}{n_{0}d}\sum_{i=1}^{n_{0}}\sum_{j=1}^{d}[\mathbf{Z}_{i}^{(1)}]_{:j}[\mathbf{Z}_{i}^{(1)}]^{\top}_{:j}-\mathbf{I}_{d_{1}}\right)$

where $[\mathbf{Z}_{i}^{(1)}]_{:j}$ is the $j$-th column of
$\mathbf{Z}_{i}^{(1)}$. By concentration of sample covariance operator (e.g.
Koltchinskii and Lounici (2017)), we have with probability at least
$1-\exp(-cd)$:

$\displaystyle\left\|d\left(\frac{1}{n_{0}d}\sum_{i=1}^{n_{0}}\sum_{j=1}^{d}[\mathbf{Z}_{i}^{(1)}]_{:j}[\mathbf{Z}_{i}^{(1)}]^{\top}_{:j}-\mathbf{I}_{d_{1}}\right)\right\|\lesssim\frac{d}{\sqrt{n}}$

By (17) and eigenvalue perturbation theory (e.g. Corollary 8.1.6 in Golub and
van Loan (2013)), we have101010With slight abuse of notation, we use
$\lambda_{j}(\cdot)$ to denote the $j$-th largest eigenvalue of a given
matrix, while $\lambda_{j}$’s themselves are singular values of $\mathbf{M}$.

$\displaystyle\lambda_{1}^{2}-\|\Delta\|\leq\lambda_{1}\left(\frac{1}{n_{0}}\mathbf{X}_{f}\mathbf{X}_{f}^{\top}-d\mathbf{I}_{d_{1}}\right)\leq\lambda_{1}^{2}+\|\Delta\|$

Therefore, we obtain

$\displaystyle\lambda_{1}^{2}-2\|\Delta\|\leq\widehat{}\mathbf{u}_{1}^{\top}\mathbf{M}\mathbf{M}^{\top}\widehat{}\mathbf{u}_{1}\leq\lambda_{1}^{2}+2\|\Delta\|$
(18)

Hence with probability at least $1-\exp(-cd)$ we have

$\displaystyle\widehat{}\mathbf{u}_{1}^{\top}\mathbf{M}\mathbf{M}^{\top}\widehat{}\mathbf{u}_{1}\geq\lambda_{1}^{2}-2\|\Delta\|\geq\lambda_{1}^{2}-C\left(\lambda_{1}\sqrt{\frac{d}{n}}+\frac{d}{\sqrt{n}}\right)\gtrsim\lambda_{1}^{2}$
(19)

where the last inequality holds provided that $\lambda^{2}\geq
C_{0}\frac{d}{\sqrt{n}}$ for some large absolute constant $C_{0}>0$.

#### Step 2:

Observe that $\widehat{}\mathbf{v}_{1}$ is the leading eigenvector of
$\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\mathbf{X}_{i}^{(2)\top}\widehat{}\mathbf{u}_{1}\widehat{}\mathbf{u}_{1}^{\top}\mathbf{X}_{i}^{(2)}-\mathbf{I}_{d_{2}}$
and we have the following decomposition:

$\displaystyle\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\mathbf{X}_{i}^{(2)\top}\widehat{}\mathbf{u}_{1}\widehat{}\mathbf{u}_{1}^{\top}\mathbf{X}_{i}^{(2)}-\mathbf{I}_{d_{2}}=\mathbf{M}^{\top}\widehat{}\mathbf{u}_{1}\widehat{}\mathbf{u}_{1}^{\top}\mathbf{M}+\Delta^{\prime}$
(20)

where

$\Delta^{\prime}:=\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}s_{i}^{(2)}\mathbf{Z}_{i}^{(2)\top}\widehat{}\mathbf{u}_{1}\right)\widehat{}\mathbf{u}_{1}^{\top}\mathbf{M}+\mathbf{M}^{\top}\widehat{}\mathbf{u}_{1}\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\widehat{}\mathbf{u}_{1}^{\top}s_{i}^{(2)}\mathbf{Z}_{i}^{(2)}\right)+\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\mathbf{Z}_{i}^{(2)\top}\widehat{}\mathbf{u}_{1}\widehat{}\mathbf{u}_{1}^{\top}\mathbf{Z}_{i}^{(2)}-\mathbf{I}_{d_{2}}$

Due to the independence of $\widehat{}\mathbf{u}_{1}$ and
$\\{\mathbf{Z}_{i}^{(2)}\\}_{i=1}^{n_{0}}$, we conclude that
$\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}s_{i}^{(2)}\mathbf{Z}_{i}^{(2)\top}\widehat{}\mathbf{u}_{1}\sim{\cal
N}(0,\frac{1}{n_{0}}\mathbf{I}_{d_{2}})$, hence with probability at least
$1-\exp(-cd)$:

$\displaystyle\left\|\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}s_{i}^{(2)}\mathbf{Z}_{i}^{(2)\top}\widehat{}\mathbf{u}_{1}\right)\widehat{}\mathbf{u}_{1}^{\top}\mathbf{M}\right\|\lesssim\lambda_{1}\sqrt{\frac{d}{n}},\quad\left\|\mathbf{M}^{\top}\widehat{}\mathbf{u}_{1}\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\widehat{}\mathbf{u}_{1}^{\top}s_{i}^{(2)}\mathbf{Z}_{i}^{(2)}\right)\right\|\lesssim\lambda_{1}\sqrt{\frac{d}{n}}$

Notice that
$\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\mathbf{Z}_{i}^{(2)\top}\widehat{}\mathbf{u}_{1}\widehat{}\mathbf{u}_{1}^{\top}\mathbf{Z}_{i}^{(2)}\overset{d}{=}\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\mathbf{z}_{i}\mathbf{z}_{i}^{\top}$,
where $\mathbf{z}_{i}\sim{\cal N}(0,\mathbf{I}_{d_{2}})$ and
$\mathbf{z}_{i}$’s are independent. Again, by concentration of sample
covariance operator, with probability at least $1-\exp(-cd)$:

$\displaystyle\left\|\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\mathbf{Z}_{i}^{(2)\top}\widehat{}\mathbf{u}_{1}\widehat{}\mathbf{u}_{1}^{\top}\mathbf{Z}_{i}^{(2)}-\mathbf{I}_{d_{2}}\right\|\lesssim\sqrt{\frac{d}{n}}\vee\frac{d}{n}$

Therefore, (20) and eigenvalue perturbation theory imply that

$\displaystyle\lambda_{1}\left(\mathbf{M}^{\top}\widehat{}\mathbf{u}_{1}\widehat{}\mathbf{u}_{1}^{\top}\mathbf{M}\right)-\|\Delta^{\prime}\|\leq\lambda_{1}\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\mathbf{X}_{i}^{(2)\top}\widehat{}\mathbf{u}_{1}\widehat{}\mathbf{u}_{1}^{\top}\mathbf{X}_{i}^{(2)}-\mathbf{I}_{d_{2}}\right)\leq\lambda_{1}\left(\mathbf{M}^{\top}\widehat{}\mathbf{u}_{1}\widehat{}\mathbf{u}_{1}^{\top}\mathbf{M}\right)+\|\Delta^{\prime}\|$

Combined with the decomposition (20), we can arrive at

$\displaystyle\lambda_{1}^{2}-2\left(\|\Delta\|+\|\Delta^{\prime}\|\right)\leq\widehat{}\mathbf{v}_{1}^{\top}\mathbf{M}^{\top}\widehat{}\mathbf{u}_{1}\widehat{}\mathbf{u}_{1}^{\top}\mathbf{M}\widehat{}\mathbf{v}_{1}\leq\lambda_{1}^{2}+2\left(\|\Delta\|+\|\Delta^{\prime}\|\right)$
(21)

Thus, we get

$\displaystyle|\widehat{}\mathbf{u}_{1}^{\top}\mathbf{M}\widehat{}\mathbf{v}_{1}|\asymp\lambda_{1}$
(22)

with probability at least $1-\exp(-cd)$, provided that $\lambda^{2}\geq
C_{0}\frac{d}{\sqrt{n}}$.

#### Step 3:

We proceed our analysis by conditioning on the event
$\\{\eqref{step2-lowerbound}\text{~{}holds}\\}$. Observe that

$\displaystyle\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}(\widehat{}\mathbf{u}_{1}^{\top}\mathbf{X}_{i}^{(3)}\widehat{}\mathbf{v}_{1})\mathbf{X}_{i}^{(3)}-\widehat{}\mathbf{u}_{1}\widehat{}\mathbf{v}_{1}^{\top}$
$\displaystyle=\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}(\widehat{}\mathbf{u}_{1}^{\top}(s_{i}^{(3)}\mathbf{M}+\mathbf{Z}_{i}^{(3)})\widehat{}\mathbf{v}_{1})(s_{i}^{(3)}\mathbf{M}+\mathbf{Z}_{i}^{(3)})$
$\displaystyle=:(\widehat{}\mathbf{u}_{1}^{\top}\mathbf{M}\widehat{}\mathbf{v}_{1})\mathbf{M}+\Upsilon$
(23)

where

$\Upsilon:=\widehat{}\mathbf{u}_{1}^{\top}\mathbf{M}\widehat{}\mathbf{v}_{1}\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}s_{i}^{(3)}\mathbf{Z}_{i}^{(3)}\right)+\mathbf{M}\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}s_{i}^{(3)}(\widehat{}\mathbf{u}_{1}^{\top}\mathbf{Z}_{i}^{(3)}\widehat{}\mathbf{v}_{1})\right)+\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}(\widehat{}\mathbf{u}_{1}^{\top}\mathbf{Z}_{i}^{(3)}\widehat{}\mathbf{v}_{1})\mathbf{Z}_{i}^{(3)}-\widehat{}\mathbf{u}_{1}\widehat{}\mathbf{v}_{1}^{\top}$

Now we give an upper bound for $\|\Upsilon\|$. By random matrix theory we know
with probability at least $1-\exp(-cd)$:

$\left\|\widehat{}\mathbf{u}_{1}^{\top}\mathbf{M}\widehat{}\mathbf{v}_{1}\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}s_{i}^{(3)}\mathbf{Z}_{i}^{(3)}\right)\right\|\lesssim\lambda_{1}\sqrt{\frac{d}{n}}$

Next, notice that
$\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}s_{i}^{(3)}(\widehat{}\mathbf{u}_{1}^{\top}\mathbf{Z}_{i}^{(3)}\widehat{}\mathbf{v}_{1})\sim{\cal
N}(0,\frac{1}{n_{0}})$, we have with probability at least $1-\exp(-cd)$:

$\left\|\mathbf{M}\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}s_{i}^{(3)}(\widehat{}\mathbf{u}_{1}^{\top}\mathbf{Z}_{i}^{(3)}\widehat{}\mathbf{v}_{1})\right)\right\|\lesssim\lambda_{1}\sqrt{\frac{d}{n}}$

It remains to bound
$\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}(\widehat{}\mathbf{u}_{1}^{\top}\mathbf{Z}_{i}^{(3)}\widehat{}\mathbf{v}_{1})\mathbf{Z}_{i}^{(3)}-\widehat{}\mathbf{u}_{1}\widehat{}\mathbf{v}_{1}^{\top}$.
Notice the following decomposition:

$\displaystyle\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}(\widehat{}\mathbf{u}_{1}^{\top}\mathbf{Z}_{i}^{(3)}\widehat{}\mathbf{v}_{1})\mathbf{Z}_{i}^{(3)}=\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}(\widehat{}\mathbf{u}_{1}^{\top}\mathbf{Z}_{i}^{(3)}\widehat{}\mathbf{v}_{1})\left[\mathcal{P}_{\widehat{}\mathbf{u}_{1}}\mathbf{Z}_{i}^{(3)}\mathcal{P}_{\widehat{}\mathbf{v}_{1}}+\mathcal{P}^{\perp}_{\widehat{}\mathbf{u}_{1}}\mathbf{Z}_{i}^{(3)}\mathcal{P}_{\widehat{}\mathbf{v}_{1}}+\mathcal{P}^{\perp}_{\widehat{}\mathbf{u}_{1}}\mathbf{Z}_{i}^{(3)}\mathcal{P}^{\perp}_{\widehat{}\mathbf{v}_{1}}+\mathcal{P}_{\widehat{}\mathbf{u}_{1}}\mathbf{Z}_{i}^{(3)}\mathcal{P}^{\perp}_{\widehat{}\mathbf{v}_{1}}\right]$

where $\mathcal{P}_{\mathbf{u}}$ is the projection matrix onto the column
space of $\mathbf{u}$ and $\mathcal{P}_{\mathbf{u}}^{\perp}$ is the projection
matrix onto orthogonal complement of the column space of $\mathbf{u}$. Since
$\sum_{i=1}^{n_{0}}(\widehat{}\mathbf{u}_{1}^{\top}\mathbf{Z}_{i}^{(3)}\widehat{}\mathbf{v}_{1})^{2}\sim\chi^{2}_{n_{0}}$,
by concentration for chi-square random variable with $n_{0}$ degrees of
freedom (see Laurent and Massart (2000)), we have with probability at least
$1-\exp(-c\sqrt{d(d\wedge n)})$:

$\left\|\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}(\widehat{}\mathbf{u}_{1}^{\top}\mathbf{Z}_{i}^{(3)}\widehat{}\mathbf{v}_{1})\mathcal{P}_{\widehat{}\mathbf{u}_{1}}\mathbf{Z}_{i}^{(3)}\mathcal{P}_{\widehat{}\mathbf{v}_{1}}-\widehat{}\mathbf{u}_{1}\widehat{}\mathbf{v}_{1}^{\top}\right\|=\left\|\widehat{}\mathbf{u}_{1}\widehat{}\mathbf{v}_{1}^{\top}\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}(\widehat{}\mathbf{u}_{1}^{\top}\mathbf{Z}_{i}^{(3)}\widehat{}\mathbf{v}_{1})^{2}-1\right)\right\|\lesssim{\sqrt{\frac{d}{n}}}$

By property of Gaussian matrices we have
$(\widehat{}\mathbf{u}_{1}^{\top}\mathbf{Z}_{i}^{(3)}\widehat{}\mathbf{v}_{1})\mathcal{P}^{\perp}_{\widehat{}\mathbf{u}_{1}}\mathbf{Z}_{i}^{(3)}\mathcal{P}_{\widehat{}\mathbf{v}_{1}}\stackrel{{\scriptstyle
d}}{{=}}g_{i}\mathbf{Z}_{i}^{(3)}$, where $g_{i}\overset{i.i.d}{\sim}{\cal
N}(0,1)$ and $\\{g_{i}\\}_{i=1}^{n_{0}}$ is independent of
$\\{\mathbf{Z}_{i}^{(3)}\\}_{i=1}^{n_{0}}$. Hence we have

$\mathbb{P}\left(\left\|\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}(\widehat{}\mathbf{u}_{1}^{\top}\mathbf{Z}_{i}^{(3)}\widehat{}\mathbf{v}_{1})\mathcal{P}^{\perp}_{\widehat{}\mathbf{u}_{1}}\mathbf{Z}_{i}^{(3)}\mathcal{P}_{\widehat{}\mathbf{v}_{1}}\right\|\leq\frac{\sqrt{d}}{n}\sqrt{\sum_{i=1}^{n_{0}}g_{i}^{2}}\Bigg{|}\\{g_{i}\\}_{i=1}^{n_{0}}\right)\geq
1-\exp(-cd)$ (24)

In addition, by concentration for chi-square random variable we have
$\sqrt{\sum_{i=1}^{n_{0}}g_{i}^{2}}\lesssim\sqrt{n}$ with probability at least
$1-\exp(-cn)$. Combined with (24), we arrive at

$\left\|\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}(\widehat{}\mathbf{u}_{1}^{\top}\mathbf{Z}_{i}^{(3)}\widehat{}\mathbf{v}_{1})\mathcal{P}^{\perp}_{\widehat{}\mathbf{u}_{1}}\mathbf{Z}_{i}^{(3)}\mathcal{P}_{\widehat{}\mathbf{v}_{1}}\right\|\lesssim\sqrt{\frac{d}{n}}$

with probability at least $1-\exp(-c(d\wedge n))$. Similar arguments can be
applied to
$(\widehat{}\mathbf{u}_{1}^{\top}\mathbf{Z}_{i}^{(3)}\widehat{}\mathbf{v}_{1})\mathcal{P}^{\perp}_{\widehat{}\mathbf{u}_{1}}\mathbf{Z}_{i}^{(3)}\mathcal{P}^{\perp}_{\widehat{}\mathbf{v}_{1}}$
and
$(\widehat{}\mathbf{u}_{1}^{\top}\mathbf{Z}_{i}^{(3)}\widehat{}\mathbf{v}_{1})\mathcal{P}_{\widehat{}\mathbf{u}_{1}}\mathbf{Z}_{i}^{(3)}\mathcal{P}^{\perp}_{\widehat{}\mathbf{v}_{1}}$.
Collecting four parts we can bound the last term of $\Upsilon$ as

$\displaystyle\left\|\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}(\widehat{}\mathbf{u}_{1}^{\top}\mathbf{Z}_{i}^{(3)}\widehat{}\mathbf{v}_{1})\mathbf{Z}_{i}^{(3)}-\widehat{}\mathbf{u}_{1}\widehat{}\mathbf{v}_{1}^{\top}\right\|\lesssim\sqrt{\frac{d}{n}}$
(25)

with probability at least $1-\exp(-c(d\wedge n))$. Hence we have the following
bound for $\|\Upsilon\|$ with probability at least $1-\exp(-c(d\wedge n))$:

$\displaystyle\|\Upsilon\|\lesssim\lambda_{1}\sqrt{\frac{d}{n}}+\sqrt{\frac{d}{n}}$
(26)

For any $j\in[r]$, denote $\tilde{\lambda}_{j}$ is the $j$-th largest singular
value of
$\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}(\widehat{}\mathbf{u}_{1}^{\top}\mathbf{X}_{i}^{(3)}\tilde{}\mathbf{V}_{1})\mathbf{X}_{i}^{(3)}-\widehat{}\mathbf{u}_{1}\widehat{}\mathbf{v}_{1}^{\top}$
and $\tilde{}\mathbf{u}_{j},\tilde{}\mathbf{v}_{j}$ the corresponding left and
right singular vectors. By (A.2) and perturbation theory for singular values
we have

$\displaystyle\sigma_{j}\left((\widehat{}\mathbf{u}_{1}^{\top}\mathbf{M}\widehat{}\mathbf{v}_{1})\mathbf{M}\right)-\|\Upsilon\|\leq\tilde{\lambda}_{j}\leq\sigma_{j}\left((\widehat{}\mathbf{u}_{1}^{\top}\mathbf{M}\widehat{}\mathbf{v}_{1})\mathbf{M}\right)+\|\Upsilon\|$
(27)

By definition of singular value and singular vectors, we have that

$\displaystyle\tilde{\lambda}_{j}=\tilde{}\mathbf{u}_{j}^{\top}\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}(\widehat{}\mathbf{u}_{1}^{\top}\mathbf{X}_{i}^{(3)}\widehat{}\mathbf{v}_{1})\mathbf{X}_{i}^{(3)}-\widehat{}\mathbf{u}_{1}\widehat{}\mathbf{v}_{1}^{\top}\right)\tilde{}\mathbf{v}_{j}=(\widehat{}\mathbf{u}_{1}^{\top}\mathbf{M}\widehat{}\mathbf{v}_{1})(\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j})+\tilde{}\mathbf{u}_{j}^{\top}\Upsilon\tilde{}\mathbf{v}_{j}$

which implies

$\displaystyle\sigma_{j}\left((\widehat{}\mathbf{u}_{1}^{\top}\mathbf{M}\widehat{}\mathbf{v}_{1})\mathbf{M}\right)-2\|\Upsilon\|\leq(\widehat{}\mathbf{u}_{1}^{\top}\mathbf{M}\widehat{}\mathbf{v}_{1})(\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j})\leq\sigma_{j}\left((\widehat{}\mathbf{u}_{1}^{\top}\mathbf{M}\widehat{}\mathbf{v}_{1})\mathbf{M}\right)+2\|\Upsilon\|$
(28)

Using (22), it follows that with probability at least $1-\exp(-c(d\wedge n))$
such that for all $j\in[r]$:

$\displaystyle\left|\left|\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}\right|-\lambda_{j}\right|\lesssim\frac{\|\Upsilon\|}{\lambda_{1}}=o(\lambda)$
(29)

given that $\lambda^{2}\gtrsim\frac{d}{\sqrt{n}}$. Notice that this implies
with overwhelming probability, we have

  1. 1.

$|\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}|\asymp\lambda_{j}$

  2. 2.

$\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}$ share the same
sign with $\widehat{}\mathbf{u}_{1}^{\top}\mathbf{M}\widehat{}\mathbf{v}_{1}$

As a consequence,
$|\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}|\asymp\sum_{j=1}^{r}\lambda_{j}$
with probability at least $1-\exp(-c(d\wedge n))$. These facts will be used in
the following derivations.

#### Step 4:

We proceed by conditioning on the event
$\\{\eqref{step3-bound},\eqref{step3-absbound}\text{~{}holds}\\}$. Consider
the rank-$r$ approximation of
$\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{X}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)\mathbf{X}_{i}^{(4)}-\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}\tilde{}\mathbf{v}_{j}^{\top}$,
which admits the following decomposition:

$\displaystyle\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{X}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)\mathbf{X}_{i}^{(4)}-\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}\tilde{}\mathbf{v}_{j}^{\top}$
$\displaystyle=\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}(\tilde{}\mathbf{u}_{j}^{\top}(s_{i}^{(4)}\mathbf{M}+\mathbf{Z}_{i}^{(4)})\tilde{}\mathbf{v}_{j}\right)(s_{i}^{(4)}\mathbf{M}+\mathbf{Z}_{i}^{(4)})$
$\displaystyle=:\sum_{j=1}^{r}\left(\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}\right)\mathbf{M}+\Upsilon^{\prime}$
(30)

where

$\displaystyle\Upsilon^{\prime}$
$\displaystyle:=\sum_{j=1}^{r}\left(\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}\right)\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}s_{i}^{(4)}\mathbf{Z}_{i}^{(4)}\right)+\mathbf{M}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}s_{i}^{(4)}\mathbf{Z}_{i}^{(4)}\right)\tilde{}\mathbf{v}_{j}\right)$
$\displaystyle+\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)\mathbf{Z}_{i}^{(4)}-\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}\tilde{}\mathbf{v}_{j}^{\top}$

Similar to that in step 3, we need to upper bound $\|\Upsilon^{\prime}\|$, the
spectral norm of the perturbation term. The following bound is clear, which
holds with probability at least $1-\exp(-cd)$:

$\displaystyle\left\|\sum_{j=1}^{r}\left(\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}\right)\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}s_{i}^{(4)}\mathbf{Z}_{i}^{(4)}\right)\right\|\lesssim\sum_{j=1}^{r}\left(\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}\right)\sqrt{\frac{d}{n}}\lesssim\left(\sum_{j=1}^{r}\lambda_{j}\right)\sqrt{\frac{d}{n}}$

Due to the rotation invariance of Gaussian and the orthogonality of
$\tilde{}\mathbf{u}_{j}$’s and $\tilde{}\mathbf{v}_{j}$’s, we have with
probability at least $1-\exp(-cd/r)$:

$\displaystyle\left\|\mathbf{M}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}s_{i}^{(4)}\mathbf{Z}_{i}^{(4)}\right)\tilde{}\mathbf{v}_{j}\right)\right\|\lesssim\lambda_{1}\sqrt{\frac{d}{n}}$

The following decomposition is similar to that in step 3:

$\displaystyle\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)\mathbf{Z}_{i}^{(4)}=\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)\left[\mathcal{P}_{\tilde{}\mathbf{U}}\mathbf{Z}_{i}^{(4)}\mathcal{P}_{\tilde{}\mathbf{V}}+\mathcal{P}^{\perp}_{\tilde{}\mathbf{U}}\mathbf{Z}_{i}^{(4)}\mathcal{P}_{\tilde{}\mathbf{V}}+\mathcal{P}^{\perp}_{\tilde{}\mathbf{U}}\mathbf{Z}_{i}^{(4)}\mathcal{P}^{\perp}_{\tilde{}\mathbf{V}}+\mathcal{P}_{\tilde{}\mathbf{U}}\mathbf{Z}_{i}^{(4)}\mathcal{P}^{\perp}_{\tilde{}\mathbf{V}}\right]$
(31)

By the property of Gaussian matrices, we have

$\displaystyle\left\|\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)\mathcal{P}_{\tilde{}\mathbf{U}}\mathbf{Z}_{i}^{(4)}\mathcal{P}_{\tilde{}\mathbf{V}}\right\|{=}\left\|\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\text{Tr}\left(\tilde{}\mathbf{U}^{\top}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{V}\right){\tilde{}\mathbf{U}^{\top}}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{V}\right\|\overset{d}{=}\left\|\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\text{Tr}\left(\mathbf{Z}_{r,i}\right)\mathbf{Z}_{r,i}\right\|$

where $\\{\mathbf{Z}_{r,i}\\}_{i=1}^{n_{0}}$ are independent matrices of
dimension $r\times r$ with i.i.d standard normal entries. Hence

$\displaystyle\left\|\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)\mathcal{P}_{\tilde{}\mathbf{U}}\mathbf{Z}_{i}^{(4)}\mathcal{P}_{\tilde{}\mathbf{V}}-\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}\tilde{}\mathbf{v}_{j}^{\top}\right\|\overset{d}{=}\left\|\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\text{Tr}\left(\mathbf{Z}_{r,i}\right)\mathbf{Z}_{r,i}-\mathbf{I}_{r}\right\|$

The following lemma gives the concentration inequality of the above term.

######  Lemma 5.

Let $\mathbf{Z},\mathbf{Z}_{1},\cdots,\mathbf{Z}_{n}$ be $r\times r$
independent matrices with i.i.d standard normal entries. Then there exists a
constant $C>0$ such that, for all $t>0$, with probability at least $1-e^{-t}$:

$\left\|\frac{1}{n}\sum_{i=1}^{n}\text{Tr}\left(\mathbf{Z}_{i}\right)\mathbf{Z}_{i}-\mathbf{I}_{r}\right\|\leq
C\left(r\sqrt{\frac{t+\log(2r)}{n}}+r\frac{t+\log(2r)}{n}\right)$

By Lemma 5, if $d<nr$, we take $t=d/r$, then we have
$\left\|\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\text{Tr}\left(\mathbf{Z}_{i,r}\right)\mathbf{Z}_{i,r}-\mathbf{I}_{r}\right\|\lesssim\sqrt{\frac{dr}{n}}$
with probability at least $1-\exp(-d/r)$, provided that $d\gtrsim r\log r$. If
$d>nr$, we can take $t=\sqrt{nd/r}$, then we have
$\left\|\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\text{Tr}\left(\mathbf{Z}_{i,r}\right)\mathbf{Z}_{i,r}-\mathbf{I}_{r}\right\|\lesssim\sqrt{\frac{dr}{n}}$
with probability at least $1-\exp(-n)$, provided that $nd\gtrsim
r^{2}\log^{2}r$. In summary, we have with probability at least
$1-\exp(-c(d/r\wedge n))$:

$\displaystyle\left\|\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\text{Tr}\left(\mathbf{Z}_{i,r}\right)\mathbf{Z}_{i,r}-\mathbf{I}_{r}\right\|\lesssim\sqrt{\frac{dr}{n}}$

In addition, we have
$\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)\mathcal{P}^{\perp}_{\tilde{}\mathbf{U}}\mathbf{Z}_{i}^{(4)}\mathcal{P}_{\tilde{}\mathbf{V}}\overset{d}{=}\sqrt{r}g_{i}\mathbf{Z}_{i}^{(4)}$,
where $g_{i}\overset{i.i.d}{\sim}N(0,1)$ and $\\{g_{i}\\}_{i=1}^{n_{0}}$ is
independent of $\\{\mathbf{Z}_{i}^{(4)}\\}_{i=1}^{n_{0}}$. Hence we have

$\displaystyle\mathbb{P}\left(\left\|\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)\mathcal{P}^{\perp}_{\tilde{}\mathbf{U}}\mathbf{Z}_{i}^{(4)}\mathcal{P}_{\tilde{}\mathbf{V}}\right\|\geq\frac{\sqrt{dr}}{n}\sqrt{\sum_{i=1}^{n_{0}}g_{i}^{2}}\Bigg{|}\\{g_{i}\\}_{i=1}^{n_{0}}\right)\leq\exp(-cd)$

Note that by concentration for chi-square random variable
$\sqrt{\sum_{i=1}^{n_{0}}g_{i}^{2}}\lesssim\sqrt{n}$ with probability at least
$1-\exp(-cn)$. Then we can conclude that with probability at least
$1-\exp(-c(d\wedge n))$:

$\left\|\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)\mathcal{P}^{\perp}_{\tilde{}\mathbf{U}}\mathbf{Z}_{i}^{(4)}\mathcal{P}_{\tilde{}\mathbf{V}}\right\|\lesssim\sqrt{\frac{dr}{n}}$

The bounds for
$\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)\mathcal{P}^{\perp}_{\tilde{}\mathbf{U}}\mathbf{Z}_{i}^{(4)}\mathcal{P}^{\perp}_{\tilde{}\mathbf{V}}$
and
$\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)\mathcal{P}_{\tilde{}\mathbf{U}}\mathbf{Z}_{i}^{(4)}\mathcal{P}^{\perp}_{\tilde{}\mathbf{V}}$
can be obtained similarly. We have with probability at least
$1-\exp(-c(d/r\wedge n))$:

$\displaystyle\|\Upsilon^{\prime}\|\lesssim\left(\sum_{j=1}^{r}\lambda_{j}\right)\sqrt{\frac{d}{n}}+\sqrt{\frac{dr}{n}}$
(32)

#### Step 5:

Again, we continue on the event
$\\{\eqref{step4-upsilonprimebound}\text{~{}holds}\\}$. In this step, we first
construct $\widehat{\Lambda}$, which is an estimator for the pre-factor
$\Lambda^{*}:=|\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}|$
of the signal part in (A.2). Notice that

$\displaystyle\frac{1}{n_{0}}$
$\displaystyle\left[\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{X}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)^{2}-r\right]$
(33)
$\displaystyle=\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}\right)^{2}+\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}\right)\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)+\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)^{2}-r$
(34)

The second term
$\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}\right)\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)\overset{d}{=}\left|\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}\right|\sqrt{\frac{r}{n_{0}}}g$,
with $g$ being standard normal. Therefore, with probability at least
$1-\exp(-cd)$:

$\displaystyle\left|\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}\right)\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)\right|\lesssim\left(\sum_{j=1}^{r}\lambda_{j}\right)\sqrt{\frac{dr}{n}}$

The third term
$\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)^{2}-r\overset{d}{=}\frac{r}{n_{0}}(b-n_{0})$,
where $b{\sim}\chi_{n_{0}}^{2}$, hence we have with probability at least
$1-\exp(-c\sqrt{d(d\wedge n)})$:

$\displaystyle\left|\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)^{2}-r\right|\lesssim
r\sqrt{\frac{d}{n}}$

Therefore, we have with probability at least $1-\exp(-c\sqrt{d(d\wedge n)})$:

$\displaystyle\frac{1}{n_{0}}\left[\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{X}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)^{2}-r\right]\gtrsim\left(\sum_{j=1}^{r}\lambda_{j}\right)^{2}-\left(\sum_{j=1}^{r}\lambda_{j}\right)\sqrt{\frac{dr}{n}}-r\sqrt{\frac{d}{n}}\gtrsim\frac{dr^{2}}{\sqrt{n}}$

where we’ve used the fact
$\sum_{j=1}^{r}\lambda_{j}\gtrsim\frac{\sqrt{d}r}{n^{1/4}}$. Hence by
definition of $\widehat{\Lambda}$, with probability at least
$1-\exp(-c\sqrt{d(d\vee n)})$:

$\displaystyle\widehat{\Lambda}^{2}=\max\left\\{\frac{1}{n_{0}}\left[\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{X}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)^{2}-r\right],\frac{dr^{2}}{\sqrt{n}}\right\\}=\frac{1}{n_{0}}\left[\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{X}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)^{2}-r\right]$
(35)

The concentration inequalities of second and third term of (33) also imply
that with probability at least $1-\exp(-c\sqrt{d(d\vee n)})$:

$\displaystyle|\widehat{\Lambda}-\Lambda^{*}|=\frac{|\widehat{\Lambda}^{2}-\Lambda^{*2}|}{\widehat{\Lambda}+\Lambda^{*}}\lesssim\frac{\left(\sum_{j=1}^{r}\lambda_{j}\right)\sqrt{\frac{dr}{n}}+r\sqrt{\frac{d}{n}}}{\sum_{j=1}^{r}\lambda_{j}}\leq\sqrt{\frac{dr}{n}}+\frac{r}{\sum_{j=1}^{r}\lambda_{j}}\sqrt{\frac{d}{n}}$
(36)

Since the RHS of (36) is of order $o\left(\sum_{j=1}^{r}\lambda_{j}\right)$,
we have $\widehat{\Lambda}\asymp\Lambda^{*}\asymp\sum_{j=1}^{r}\lambda_{j}$
with probability at least $1-\exp(-c\sqrt{d(d\vee n)})$. Next, denote
$\check{}\mathbf{U}$, $\check{}\mathbf{V}$ the left and right leading $r$
singular vectors
of$\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{X}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)\mathbf{X}_{i}^{(4)}-\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}\tilde{}\mathbf{v}_{j}^{\top}$,
then the best rank-$r$ approximation is given by

$\displaystyle\check{}\mathbf{M}$
$\displaystyle=\check{}\mathbf{U}\check{}\mathbf{U}^{\top}\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{X}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)\mathbf{X}_{i}^{(4)}-\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}\tilde{}\mathbf{v}_{j}^{\top}\right)\check{}\mathbf{V}\check{}\mathbf{V}^{\top}$
$\displaystyle=\sum_{j=1}^{r}\left(\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}\right)\check{}\mathbf{U}\check{}\mathbf{U}^{\top}\mathbf{M}\check{}\mathbf{V}\check{}\mathbf{V}^{\top}+\check{}\mathbf{U}\check{}\mathbf{U}^{\top}\Upsilon^{\prime}\check{}\mathbf{V}\check{}\mathbf{V}^{\top}$

The error of low-rank approximation is characterized by the perturbation term,
given by the following lemma.

######  Lemma 6.

Consider a rank-$r$ matrix $\mathbf{M}\in\mathbb{R}^{d_{1}\times d_{2}}$ with
its thin-SVD form $\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$, where
$\mathbf{U}\in\mathbb{O}_{d_{1},r},\mathbf{V}\in\mathbb{O}_{d_{2},r}$ and
$\mathbf{\Sigma}=\text{diag}(\sigma_{1},\cdots,\sigma_{r})$,
$\sigma_{1}\geq\sigma_{2}\geq\cdots\sigma_{r}>0$, let $E$ be a $d_{1}\times
d_{2}$ perturbation matrix and $\widehat{}\mathbf{M}=\mathbf{M}+\mathbf{E}$.
Denote $\widehat{}\mathbf{M}_{r}$ the best rank-$r$ approximation of
$\widehat{}\mathbf{M}$. Suppose that $\sigma_{r}\geq 3\|\mathbf{E}\|$, then
there exists some absolute constant $C_{0}>0$ such that

$\displaystyle\|\widehat{}\mathbf{M}_{r}-\mathbf{M}\|_{\rm F}\leq
C_{0}\min\\{\|\mathbf{E}\|_{\rm F},\sqrt{r}\|\mathbf{E}\|\\}$

By Lemma 6, we have

$\displaystyle\left\|\check{}\mathbf{M}-\sum_{j=1}^{r}\left(\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}\right)\mathbf{M}\right\|_{\text{F}}\lesssim\sqrt{r}\|\Upsilon^{\prime}\|$

Recall that $\widehat{}\mathbf{M}=\check{}\mathbf{M}/\widehat{\Lambda}$.
Denote
$\eta^{*}=\text{sign}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}\right)$,
hence we have the following bound

$\displaystyle\|\widehat{}\mathbf{M}-\eta^{*}\mathbf{M}\|_{\rm
F}\leq\left\|\widehat{\Lambda}^{-1}\check{}\mathbf{M}-\widehat{\Lambda}^{-1}\sum_{j=1}^{r}\left(\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}\right)\mathbf{M}\right\|_{\text{F}}+\left\|\widehat{\Lambda}^{-1}\sum_{j=1}^{r}\left(\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}\right)\mathbf{M}-\eta^{*}\mathbf{M}\right\|_{\text{F}}$
(37)

Using (32) and (36), the first term can be bounded with probability at least
$1-\exp(-c(d/r\wedge n))$:

$\displaystyle\left\|\widehat{\Lambda}^{-1}\check{}\mathbf{M}-\widehat{\Lambda}^{-1}\sum_{j=1}^{r}\left(\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}\right)\mathbf{M}\right\|_{\text{F}}\lesssim\frac{\sqrt{r}\|\Upsilon^{\prime}\|}{\widehat{\Lambda}}\lesssim\sqrt{\frac{dr}{n}}+\frac{r}{\sum_{j=1}^{r}\lambda_{j}}\sqrt{\frac{d}{n}}$

Using (36), the second term can be bounded with probability at least
$1-\exp(-c\sqrt{d(d\wedge n)})$:

$\displaystyle\left\|\widehat{\Lambda}^{-1}\sum_{j=1}^{r}\left(\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}\right)\mathbf{M}-\eta^{*}\mathbf{M}\right\|_{\text{F}}=\left|\frac{\Lambda^{*}-\widehat{\Lambda}}{\widehat{\Lambda}}\right|\|\mathbf{M}\|_{\text{F}}\lesssim\sqrt{\frac{dr}{n}}+\frac{r}{\sum_{j=1}^{r}\lambda_{j}}\sqrt{\frac{d}{n}}$

where we’ve used the fact $\|\mathbf{M}\|_{\rm
F}\leq\sum_{j=1}^{r}\lambda_{j}$. Take union bound over all events that we’ve
conditioned on in previous steps, we conclude that with probability at least
$1-\exp(-c(d/r\wedge n))$:

$\displaystyle\min_{\eta\in\\{\pm
1\\}}\|\widehat{}\mathbf{M}-\eta\mathbf{M}\|_{\rm
F}\lesssim\sqrt{\frac{dr}{n}}+\frac{r}{\sum_{j=1}^{r}\lambda_{j}}\sqrt{\frac{d}{n}}$
(38)

#### Bound in expectation:

Denote the event $Q:=\\{\eqref{step5-hpbound}\text{~{}holds}\\}$. Then we have
the following bound in expectation:

$\displaystyle\mathbb{E}\min_{\eta\in\\{\pm
1\\}}\|\widehat{}\mathbf{M}-\eta\mathbf{M}\|_{\rm
F}=\mathbb{E}\min_{\eta\in\\{\pm
1\\}}\|\widehat{}\mathbf{M}-\eta\mathbf{M}\|_{\rm
F}\mathbb{I}_{Q}+\mathbb{E}\min_{\eta\in\\{\pm
1\\}}\|\widehat{}\mathbf{M}-\eta\mathbf{M}\|_{\rm F}\mathbb{I}_{Q^{c}}$ (39)

Note that by Von Neumann’s trace inequality, we have

$\displaystyle\Lambda^{*}$
$\displaystyle=\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}=\text{Tr}(\tilde{}\mathbf{U}^{\top}\mathbf{M}\tilde{}\mathbf{V})=\text{Tr}(\tilde{}\mathbf{U}^{\top}\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}\tilde{}\mathbf{V})\leq\sum_{j=1}^{r}\sigma_{j}(\tilde{}\mathbf{U}^{\top}\mathbf{U})\sigma_{j}(\mathbf{\Sigma}\mathbf{V}^{\top}\tilde{}\mathbf{V})\leq\sum_{j=1}^{r}\lambda_{j}$
(40)

Since $\check{}\mathbf{M}$ is a rank-$r$ projection of
$\Lambda^{*}\mathbf{M}+\Upsilon^{\prime}$ and $\widehat{\Lambda}$ is lower
bounded by $\frac{\sqrt{d}r}{n^{1/4}}$, we have the following upper bound
using (40):

$\displaystyle\|\widehat{}\mathbf{M}\|_{\rm
F}=\|\widehat{\Lambda}^{-1}\check{}\mathbf{M}\|_{\rm
F}\leq\frac{\sqrt{r}}{\widehat{\Lambda}}\left(\Lambda^{*}\|\mathbf{M}\|+\|\Upsilon^{\prime}\|\right)\leq\frac{n^{1/4}}{\sqrt{dr}}\left(\lambda_{1}\sum_{j=1}^{r}\lambda_{j}+\|\Upsilon^{\prime}\|\right)$
(41)

Now we turn to bound $\mathbb{E}\|\Upsilon^{\prime}\|$, note that by
definition we have

$\displaystyle\mathbb{E}\|\Upsilon^{\prime}\|$
$\displaystyle\leq\mathbb{E}\left\|\sum_{j=1}^{r}\left(\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}\right)\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}s_{i}^{(4)}\mathbf{Z}_{i}^{(4)}\right)\right\|+\mathbb{E}\left\|\mathbf{M}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}s_{i}^{(4)}\mathbf{Z}_{i}^{(4)}\right)\tilde{}\mathbf{v}_{j}\right)\right\|$
$\displaystyle+\mathbb{E}\left\|\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)\mathbf{Z}_{i}^{(4)}-\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}\tilde{}\mathbf{v}_{j}^{\top}\right\|$
(42)

The first term of (A.2) can be bounded as

$\displaystyle\mathbb{E}\left\|\sum_{j=1}^{r}\left(\tilde{}\mathbf{u}_{j}^{\top}\mathbf{M}\tilde{}\mathbf{v}_{j}\right)\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}s_{i}^{(4)}\mathbf{Z}_{i}^{(4)}\right)\right\|\lesssim\left(\sum_{j=1}^{r}\lambda_{j}\right)\sqrt{\frac{d}{n}}$

The second term of (A.2) can be bounded as

$\displaystyle\mathbb{E}\left\|\mathbf{M}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\left(\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}s_{i}^{(4)}\mathbf{Z}_{i}^{(4)}\right)\tilde{}\mathbf{v}_{j}\right)\right\|\lesssim\lambda_{1}\sqrt{\frac{r}{n}}$

For the last term of (A.2), recall the decomposition (31), we have

$\displaystyle\mathbb{E}\left\|\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)\mathcal{P}^{\perp}_{\tilde{}\mathbf{U}}\mathbf{Z}_{i}^{(4)}\mathcal{P}_{\tilde{}\mathbf{V}}\right\|$
$\displaystyle{=}\mathbb{E}\left[\mathbb{E}\left[\left\|\frac{\sqrt{r}}{n_{0}}\sum_{i=1}^{n_{0}}g_{i}Z_{i}^{(4)}\right\|\Bigg{|}\\{g_{i}\\}_{i=1}^{n}\right]\right]\lesssim\mathbb{E}\left[\frac{\sqrt{dr}}{n_{0}}\sqrt{\sum_{i=1}^{n_{0}}g_{i}^{2}}\right]\lesssim\sqrt{\frac{dr}{n}}$

Similar bounds hold for
$\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)\mathcal{P}^{\perp}_{\tilde{}\mathbf{U}}\mathbf{Z}_{i}^{(4)}\mathcal{P}^{\perp}_{\tilde{}\mathbf{V}}$
and
$\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)\mathcal{P}_{\tilde{}\mathbf{U}}\mathbf{Z}_{i}^{(4)}\mathcal{P}^{\perp}_{\tilde{}\mathbf{V}}$.
It remains to find
$\mathbb{E}\left\|\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)\mathcal{P}_{\tilde{}\mathbf{U}}\mathbf{Z}_{i}^{(4)}\mathcal{P}_{\tilde{}\mathbf{V}}-\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}\tilde{}\mathbf{v}_{j}^{\top}\right\|$.
For simplicity, denote
$\Gamma:=\left\|\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}^{\top}\mathbf{Z}_{i}^{(4)}\tilde{}\mathbf{v}_{j}\right)\mathcal{P}_{\tilde{}\mathbf{U}}\mathbf{Z}_{i}^{(4)}\mathcal{P}_{\tilde{}\mathbf{V}}-\sum_{j=1}^{r}\tilde{}\mathbf{u}_{j}\tilde{}\mathbf{v}_{j}^{\top}\right\|$,
then using Lemma 5 we can get

$\displaystyle\mathbb{E}\Gamma$
$\displaystyle=\int_{0}^{\infty}\mathbb{P}\left(\Gamma\geq
t\right)dt=\int_{0}^{2r\sqrt{\frac{\log(2r)}{n_{0}}}}\mathbb{P}\left(\Gamma\geq
t\right)dt+\int_{2r\sqrt{\frac{\log(2r)}{n_{0}}}}^{\infty}\mathbb{P}\left(\Gamma\geq
t\right)dt$ $\displaystyle\leq
2r\sqrt{\frac{\log(2r)}{n_{0}}}+\int_{2r\sqrt{\frac{\log(2r)}{n_{0}}}}^{\infty}\frac{r}{2n_{0}}\left(\sqrt{\frac{n_{0}}{u+\log(2r)}}+2\right)\mathbb{P}\left(\Gamma\geq
r\sqrt{\frac{u+\log(2r)}{n_{0}}}+r\frac{u+\log(2r)}{n_{0}}\right)du$
$\displaystyle\lesssim r\sqrt{\frac{\log
r}{n}}+\frac{r}{n}\int_{2r\sqrt{\frac{\log(2r)}{n_{0}}}}^{\infty}\left(\sqrt{\frac{n_{0}}{u+\log(2r)}}+2\right)\exp(-u)du\lesssim
r\sqrt{\frac{\log r}{n}}$

Hence we can conclude that

$\displaystyle\mathbb{E}\|\Upsilon^{\prime}\|\leq\left(\sum_{j=1}^{r}\lambda_{j}\right)\sqrt{\frac{d}{n}}+\sqrt{\frac{dr}{n}}$
(43)

provided that $d\gtrsim r\log r$. By (41), we have

$\displaystyle\mathbb{E}\|\widehat{}\mathbf{M}\|_{\rm
F}\lesssim\frac{n^{1/4}}{\sqrt{dr}}\left(\lambda_{1}\sum_{j=1}^{r}\lambda_{j}+\left(\sum_{j=1}^{r}\lambda_{j}\right)\sqrt{\frac{d}{n}}+\sqrt{\frac{dr}{n}}\right)\lesssim\frac{n^{1/4}}{\sqrt{dr}}\left(\lambda_{1}\sum_{j=1}^{r}\lambda_{j}\right)$

Hence

$\displaystyle\mathbb{E}\min_{\eta\in\\{\pm
1\\}}\|\widehat{}\mathbf{M}-\eta\mathbf{M}\|_{\rm
F}\leq\mathbb{E}\|\widehat{}\mathbf{M}\|_{\rm F}+\mathbb{E}\|\mathbf{M}\|_{\rm
F}\lesssim\frac{\lambda_{1}n^{1/4}}{\sqrt{d}}\left(\frac{1}{\sqrt{r}}\sum_{j=1}^{r}\lambda_{j}\right)+\sqrt{\sum_{j=1}^{r}\lambda_{j}^{2}}$

Then (39) implies that

$\displaystyle\mathbb{E}\min_{\eta\in\\{\pm
1\\}}\|\widehat{}\mathbf{M}-\eta\mathbf{M}\|_{\rm F}$
$\displaystyle\leq\sqrt{\frac{dr}{n}}+\frac{r}{\sum_{j=1}^{r}\lambda_{j}}\sqrt{\frac{d}{n}}+\left[\frac{\lambda_{1}n^{1/4}}{\sqrt{d}}\left(\frac{1}{\sqrt{r}}\sum_{j=1}^{r}\lambda_{j}\right)+\sqrt{\sum_{j=1}^{r}\lambda_{j}^{2}}\right]\exp(-c(r^{-1}d\wedge
n))$
$\displaystyle\lesssim\sqrt{\frac{dr}{n}}+\frac{r}{\sum_{j=1}^{r}\lambda_{j}}\sqrt{\frac{d}{n}}$

provided that $\lambda_{1}\asymp\lambda_{r}\asymp\lambda$ and
$\lambda\lesssim\exp(c(r^{-1}{d}\wedge n)-\log n)$.

###  A.3 Proof of Theorem 2

The main idea is to construct a set of sufficiently dissimilar hypotheses to
apply Fano’s method. To this end, we fix some
$\mathbf{U}_{0}\in\mathbb{O}_{d,r}$ and consider the ball centered at
$\mathbf{U}_{0}$ with radius of $\epsilon\in(0,\sqrt{2r}]$ under the chordal
Frobenius-norm metric
$\textsf{dist}(\mathbf{U}_{1},\mathbf{U}_{2}):=\min_{\mathbf{O}\in\mathbb{O}_{r}}\|\mathbf{U}_{1}-\mathbf{U}_{2}\mathbf{O}\|_{\rm
F}$ :

$B_{\epsilon}(\mathbf{U}_{0}):=\\{\mathbf{U}:\textsf{dist}(\mathbf{U},\mathbf{U}_{0})\leq\epsilon\\}$

By Lemma 1 in Cai et al. (2013) and the equivalence between
$\textsf{dist}(\cdot,\cdot)$ and $\|\sin\Theta(\cdot,\cdot)\|$, we have for
any $\alpha\in(0,1)$, there exists $\\{\mathbf{U}_{i}^{\prime}\\}_{i=1}^{m}$,
a packing of $B_{\epsilon}(\mathbf{U}_{0})$ such that for some absolute
constant $c_{0}>0$:

$m\geq\left(\frac{c_{0}}{\alpha}\right)^{r(d-r)},\quad\min_{i<j}\textsf{dist}(\mathbf{U}_{i}^{\prime},\mathbf{U}_{j}^{\prime})\geq\alpha\epsilon$

Denote
$\mathbf{O}_{i}=\arg\min_{\mathbf{O}\in\mathbb{O}_{r}}\|\mathbf{U}_{i}^{\prime}-\mathbf{U}_{0}\mathbf{O}\|_{\rm
F}$. Fix $\mathbf{\Sigma}=\text{diag}(\lambda_{1},\cdots,\lambda_{r})$ with
$\lambda_{1}=\cdots=\lambda_{r}=\lambda$ and $\mathbf{V}\in\mathbb{O}_{r}$, we
can construct
$\mathbf{M}_{i}=\mathbf{U}_{i}^{\prime}\mathbf{O}_{i}^{\top}\mathbf{\Sigma}\mathbf{V}^{\top}$
for $i=1,\cdots,m$. Notice that

$\displaystyle\min_{\eta\in\\{\pm
1\\}}\|\mathbf{M}_{i}-\eta\mathbf{M}_{j}\|_{\rm F}$
$\displaystyle=\min_{\eta\in\\{\pm
1\\}}\|\mathbf{U}_{i}^{\prime}\mathbf{O}_{i}^{\top}\mathbf{\Sigma}V^{\top}-\eta\mathbf{U}_{j}^{\prime}\mathbf{O}_{j}^{\top}\mathbf{\Sigma}\mathbf{V}^{\top}\|_{\rm
F}=\lambda\min_{\eta\in\\{\pm
1\\}}\|\mathbf{U}_{i}^{\prime}\mathbf{O}_{i}^{\top}-\eta\mathbf{U}_{j}^{\prime}\mathbf{O}_{j}^{\top}\|_{\rm
F}$
$\displaystyle\geq\lambda\cdot\textsf{dist}(\mathbf{U}_{i}^{\prime},\mathbf{U}_{j}^{\prime})\geq\lambda\alpha\epsilon$

Let $P_{\mathbf{M}}$ denote the distribution of
$\mathbf{X}=s\mathbf{M}+\mathbf{Z}$ and let ${P}_{j}^{1:n}$ denote the
distribution of
$\\{\mathbf{X}_{i}^{(j)}=s_{i}\mathbf{M}_{j}+\mathbf{Z}_{i},i=1,\cdots,n\\}$,
i.e, the $j$-th model parametrized by $\mathbf{M}_{j}$ for $j=1,\cdots,m$.
When $\|\mathbf{\Sigma}\|_{\rm F}=\sqrt{r}\lambda\geq 1$, since $s$ has a
Rademacher prior, using the log-sum inequality (see, e.g., Do (2003)) we have

$\displaystyle\text{D}_{\textsf{KL}}({P}_{j}^{1:n}||{P}_{k}^{1:n})$
$\displaystyle\leq\sum_{i=1}^{n}\frac{1}{2}\|\mathbf{M}_{j}-\mathbf{M}_{k}\|_{\rm
F}^{2}=\frac{1}{2}n\|\mathbf{U}_{j}^{\prime}\mathbf{O}_{j}^{\top}\mathbf{\Sigma}\mathbf{V}^{\top}-\mathbf{U}_{k}^{\prime}\mathbf{O}_{k}^{\top}\mathbf{\Sigma}\mathbf{V}^{\top}\|_{\rm
F}^{2}=\frac{1}{2}n\lambda^{2}\|\mathbf{U}_{j}^{\prime}\mathbf{O}_{j}^{\top}-\mathbf{U}_{k}^{\prime}\mathbf{O}_{k}^{\top}\|_{\rm
F}^{2}$ $\displaystyle\leq
n\lambda^{2}\left(\textsf{dist}^{2}(\mathbf{U}_{j}^{\prime},\mathbf{U}_{0})+\textsf{dist}^{2}(\mathbf{U}_{k}^{\prime},\mathbf{U}_{0})\right)\leq
2n\lambda^{2}\epsilon^{2}$

When $\|\mathbf{\Sigma}\|_{\rm F}=\sqrt{r}\lambda\leq 1$, by Lemma 27 in Wu
and Zhou (2019), there exists a universal constant $C>0$, such that for any
$\mathbf{U},\tilde{}\mathbf{U}\in\mathbb{O}_{d,r}$:

$\displaystyle\text{D}_{\textsf{KL}}(P_{\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}}||P_{\tilde{}\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}})$
$\displaystyle\leq C\min_{\eta\in\\{\pm
1\\}}\|\text{vec}(\|\mathbf{\Sigma}\|_{\rm
F}^{-1}\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top})-\eta\text{vec}(\|\mathbf{\Sigma}\|_{\rm
F}^{-1}\tilde{}\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top})\|_{\rm
F}^{2}\|\mathbf{\Sigma}\|_{\rm F}^{4}$ $\displaystyle\leq
C\|\mathbf{\Sigma}\|_{\rm
F}^{4}\frac{\|\mathbf{\Sigma}\|^{2}}{\|\mathbf{\Sigma}\|_{\rm
F}^{2}}\min_{\eta\in\\{\pm 1\\}}\|\mathbf{U}-\eta\tilde{}\mathbf{U}\|_{\rm
F}^{2}=Cr\lambda^{4}\min_{\eta\in\\{\pm
1\\}}\|\mathbf{U}-\eta\tilde{}\mathbf{U}\|_{\rm F}^{2}$

which implies that

$\displaystyle\text{D}_{\textsf{KL}}({P}_{j}^{1:n}||{P}_{k}^{1:n})$
$\displaystyle=n\text{D}_{\textsf{KL}}(P_{\mathbf{U}_{j}^{\prime}\mathbf{O}_{j}^{\top}\mathbf{\Sigma}\mathbf{V}^{\top}}||P_{\mathbf{U}_{k}^{\prime}\mathbf{O}_{k}^{\top}\mathbf{\Sigma}\mathbf{V}^{\top}})\leq
Cnr\lambda^{4}\min_{\eta\in\\{\pm
1\\}}\|\mathbf{U}_{j}^{\prime}\mathbf{O}_{j}^{\top}-\eta\mathbf{U}_{k}^{\prime}\mathbf{O}_{k}^{\top}\|_{\rm
F}^{2}$ $\displaystyle\leq
Cnr\lambda^{4}\left(\|\mathbf{U}_{j}^{\prime}\mathbf{O}_{j}^{\top}-\mathbf{U}_{0}\|_{\rm
F}^{2}+\min_{\eta\in\\{\pm
1\\}}\|\mathbf{U}_{0}-\eta\mathbf{U}_{k}^{\prime}\mathbf{O}_{k}^{\top}\|_{\rm
F}^{2}\right)$
$\displaystyle=Cnr\lambda^{4}\left(\textsf{dist}^{2}(\mathbf{U}_{j}^{\prime},\mathbf{U}_{0})+\textsf{dist}^{2}(\mathbf{U}_{k}^{\prime},\mathbf{U}_{0})\right)\leq
Cnr\lambda^{4}\epsilon^{2}$

Hence we have

$\text{D}_{\textsf{KL}}({P}_{j}^{1:n}||{P}_{k}^{1:n})\leq
Cn\lambda^{2}(r\lambda^{2}\wedge 1)\epsilon^{2}$

By Fano’s lower bound on minimax risk (see, e.g., Proposition 15.12 in
Wainwright (2019)), we have

$\inf_{\widehat{}\mathbf{M}}\sup_{\mathbf{M}\in\mathcal{M}_{d_{1},d_{2}}(r,\lambda)}\mathbb{E}\min_{\eta\in\\{\pm
1\\}}\|\widehat{}\mathbf{M}-\eta\mathbf{M}\|_{\rm
F}\geq\lambda\alpha\epsilon\left(1-\frac{Cn\lambda^{2}(r\lambda^{2}\wedge
1)\epsilon^{2}+\log 2}{r(d-r)\log(c_{0}/\alpha)}\right)$

By choosing $\epsilon=\sqrt{\frac{r(d-r)}{C_{0}n\lambda^{2}(r\lambda^{2}\wedge
1)}}\wedge\sqrt{2r}$ for some large absolute constant $C_{0}>0$ and
$\alpha=(c_{0}\wedge 1)/8$, we can guarantee that
$\left(1-\frac{Cn\lambda^{2}(r\lambda^{2}\wedge 1)\epsilon^{2}+\log
2}{r(d-r)\log(c_{0}/\alpha)}\right)\geq\frac{1}{2}$. Hence

$\displaystyle\inf_{\widehat{}\mathbf{M}}\sup_{\mathbf{M}\in\mathcal{M}_{d_{1},d_{2}}(r,\lambda)}\mathbb{E}\min_{\eta\in\\{\pm
1\\}}\|\widehat{}\mathbf{M}-\eta\mathbf{M}\|_{\rm
F}\gtrsim\lambda\left(\sqrt{\frac{dr/n}{\lambda^{2}(r\lambda^{2}\wedge
1)}}\wedge\sqrt{r}\right)\gtrsim\left(\frac{1}{\lambda}\sqrt{\frac{d}{n}}+\sqrt{\frac{dr}{n}}\right)\wedge\lambda\sqrt{r}$

∎

###  A.4 Proof of Theorem 4

Denote the prior distribution for $(\mathbf{M},\mathbf{s})$ defined in (12) as
$\Pi$, where $\mathbf{s}=(s_{1},\cdots,s_{n})$ is the latent label vector. Let
$(\mathbf{M}^{(1)},\mathbf{s}^{(1)}),(\mathbf{M}^{(2)},\mathbf{s}^{(2)})$ be
two independent copies from prior distribution $\Pi$. By Theorem 2.6 in
Kunisky et al. (2019), we have the following formula for $\|L_{n}^{\leq D}\|$
under the additive Gaussian noise model:

$\|L_{n}^{\leq
D}\|^{2}={\mathbb{E}_{\Pi}}\sum_{k=1}^{D}\frac{1}{k!}\langle\mathbf{s}^{(1)},\mathbf{s}^{(2)}\rangle^{k}\langle\mathbf{M}^{(1)},\mathbf{M}^{(2)}\rangle^{k}=1+\mathbb{E}_{\Pi}\sum_{k=1}^{\lfloor
D/2\rfloor}\frac{1}{(2k)!}\langle\mathbf{s}^{(1)},\mathbf{s}^{(2)}\rangle^{2k}\langle\mathbf{M}^{(1)},\mathbf{M}^{(2)}\rangle^{2k}$

The last inequality is due to the fact that
$\langle\mathbf{s}^{(1)},\mathbf{s}^{(2)}\rangle$ in distribution equals to
the sum of $n$ i.i.d. Rademacher random variables, denoted by
$\sum_{i=1}^{n}U_{i}$, and hence
$\mathbb{E}\langle\mathbf{s}^{(1)},\mathbf{s}^{(2)}\rangle^{k}=0$ for odd $k$.
Hence we have

$\mathbb{E}\langle\mathbf{s}^{(1)},\mathbf{s}^{(2)}\rangle^{2k}=\mathbb{E}\left(\sum_{i=1}^{n}U_{i}\right)^{2k}=\mathbb{E}\sum_{2k_{1}+\cdots+2k_{n}=2k}U_{1}^{2k_{1}}\cdots