# A Random Matrix Approach to Low-Multilinear-Rank Tensor Approximation

Hugo Lebeau<EMAIL_ADDRESS>
Université Grenoble Alpes, CNRS, Inria, Grenoble INP, LIG
Grenoble, 38000, France Florent Chatelain<EMAIL_ADDRESS>
Université Grenoble Alpes, CNRS, Grenoble INP, GIPSA-lab
Grenoble, 38000, France Romain Couillet<EMAIL_ADDRESS>
Université Grenoble Alpes, CNRS, Inria, Grenoble INP, LIG
Grenoble, 38000, France

###### Abstract

This work presents a comprehensive understanding of the estimation of a
planted low-rank signal from a general spiked tensor model near the
computational threshold. Relying on standard tools from the theory of large
random matrices, we characterize the large-dimensional spectral behavior of
the unfoldings of the data tensor and exhibit relevant signal-to-noise ratios
governing the detectability of the principal directions of the signal. These
results allow to accurately predict the reconstruction performance of
truncated multilinear SVD (MLSVD) in the non-trivial regime. This is
particularly important since it serves as an initialization of the higher-
order orthogonal iteration (HOOI) scheme, whose convergence to the best low-
multilinear-rank approximation depends entirely on its initialization. We give
a sufficient condition for the convergence of HOOI and show that the number of
iterations before convergence tends to $1$ in the large-dimensional limit.

Keywords: random tensors, random matrix theory, spiked tensor model, truncated
MLSVD, higher-order orthogonal iteration

##  1 Introduction

Information retrieval from large amounts of data has become a common task of
signal processing and machine learning in the past decades. Often, these data
have several modes as they may come from various sources, modalities, domains,
and so on. Tensors (multi-way arrays) are therefore a natural representation
for such datasets — they appear in multiple areas such as brain imaging (Zhou
et al., 2013), neurophysiological measurements (Seely et al., 2016), community
detection (Anandkumar et al., 2013), compression of hyperspectral images (Li
and Li, 2010), spatio-temporal gene expression (Liu et al., 2017), recommender
systems (Karatzoglou et al., 2010; Rendle and Schmidt-Thieme, 2010; Frolov and
Oseledets, 2017) and topic modelling (Anandkumar et al., 2014). Indeed,
tensors as multi-way arrays provide a more detailed representation of data
than mere matrices (two-way arrays) as they convey a structural information.
For instance, the modes of a data tensor can represent pixel $\times$ pixel
$\times$ wavelength $\times$ sample in hyperspectral imaging (Zhang et al.,
2013; Kanatsoulis et al., 2018), time $\times$ spatial scale $\times$
electrode in the EEG analysis by Acar et al. (2007) or neuron $\times$ time
$\times$ stimuli in the study of the visual cortex by Rabinowitz et al.
(2015).

In an information retrieval context, it is common to make use of tensor
decompositions in order to estimate a sought signal. In their fMRI study,
Hunyadi et al. (2017) perform a blind source separation via a joint tensor
decomposition on a channel $\times$ time $\times$ patient array, whereas
Williams et al. (2018) use a low-rank tensor approximation on a neuron
$\times$ time $\times$ trial array as a dimensionality reduction technique to
study neural dynamics. In fact, supposing that the signal has a low-rank
structure is a natural sparsity assumption (Kadmon and Ganguli, 2018;
Anandkumar et al., 2014), and low-rank tensor approximations are key tools to
extract information from multi-way data.

In the present work, we propose a random matrix analysis of a general low-rank
information $+$ noise tensor model and precisely quantify the amount of
information which can be recovered with a low-rank tensor approximation
depending on the signal-to-noise ratio (SNR). For a general introduction to
tensors, we refer the reader to Comon (2014, 2009); Landsberg (2011);
Hackbusch (2012) and, for an emphasis on statistical learning applications,
Sun et al. (2021). In the remainder of the introduction, the main concepts and
challenges behind low-rank tensor estimation are presented in Section 1.1.
Then, Section 1.2 introduces some important related works. Our main results
are finally summarized in Section 1.3. All the notations are properly defined
in Section 2.

###  1.1 Low-Rank Tensor Estimation

What is meant by a low-rank approximation of a tensor? And how is the _rank_
of a tensor actually defined? Let us start with a familiar matrix case: a
matrix ${\bm{M}}\in{\mathbb{R}}^{n_{1}\times n_{2}}$ is a two-way array (or
order-$2$ tensor). A singular value decomposition (SVD) allows us to write
${\bm{M}}$ in a compact way as the sum of $R$ rank-$1$ terms,
${\bm{M}}=\sum_{i=1}^{R}\sigma_{i}{\bm{u}}_{i}{\bm{v}}_{i}^{\top}={\bm{U}}{\bm{\Sigma}}{\bm{V}}^{\top}$
where ${\bm{U}}$, ${\bm{V}}$ are respectively $n_{1}\times R$ and $n_{2}\times
R$ matrices with _orthonormal_ columns and ${\bm{\Sigma}}$ is the $R\times R$
_diagonal_ matrix of singular values. The rank of ${\bm{M}}$ is here the
minimal number of rank-$1$ terms in which it can be exactly decomposed.
Extending this notion to tensors therefore seems straightforward: a
tensor111It is chosen of order $3$ for simplicity of exposure.
${\bm{\mathscr{T}}}\in{\mathbb{R}}^{n_{1}\times n_{2}\times n_{3}}$ has rank
$R$ if it is the minimal number of rank-$1$ terms in which it can be exactly
decomposed,
${\bm{\mathscr{T}}}=\sum_{i=1}^{R}\sigma_{i}{\bm{a}}_{i}\otimes{\bm{b}}_{i}\otimes{\bm{c}}_{i}$.
What we have just described is the canonical polyadic decomposition (CPD) of
${\bm{\mathscr{T}}}$, it dates back to Hitchcock (1927) and is unique under
very mild conditions (Kolda and Bader, 2009). However, we have lost an
important property in this process: the unit vectors ${\bm{a}}_{i}$ (resp.
${\bm{b}}_{i}$, ${\bm{c}}_{i}$) are, in general, no longer orthonormal.
Conversely, retaining the orthonormality property inevitably results in the
loss of the diagonality property,
${\bm{\mathscr{T}}}=\sum_{i=1}^{r_{1}}\sum_{j=1}^{r_{2}}\sum_{k=1}^{r_{3}}{\mathscr{G}}_{i,j,k}{\bm{u}}_{i}\otimes{\bm{v}}_{j}\otimes{\bm{w}}_{k}$.
This latter decomposition is called a Tucker decomposition and dates back to
Tucker (1966). In fact, the best way to represent ${\bm{\mathscr{T}}}$ with a
Tucker decomposition is to choose the ${\bm{u}}_{i}$ (resp. ${\bm{v}}_{i}$,
${\bm{w}}_{i}$) as the left singular vectors of the unfolding of
${\bm{\mathscr{T}}}$ along mode $1$ (resp. $2$, $3$)222This is properly
defined in Section 2.2.. This is called the multilinear SVD (MLSVD, De
Lathauwer et al., 2000b) and gives rise to a new definition of rank: the
_multilinear_ -rank $(r_{1},r_{2},r_{3})$. Note that, in the matrix case,
$r_{1}=r_{2}=R$ since both the diagonality and orthonormality properties are
verified. However, $r_{1},r_{2},r_{3}$ are, in general, not equal in the
tensor case, but $\max(r_{1},r_{2},r_{3})\leqslant
R\leqslant\min(r_{1}r_{2},r_{2}r_{3},r_{1}r_{3})$. See, e.g., Sidiropoulos et
al. (2017) for details. Other relevant references for the reader interested in
tensor decompositions are Kolda and Bader (2009); Cichocki et al. (2015);
Rabanser et al. (2017).

Given an order-$d$ tensor
${\bm{\mathscr{T}}}\in{\mathbb{R}}^{n_{1}\times\ldots\times n_{d}}$ of
possibly very high rank, we are interested in finding a low-rank
approximation, i.e., an $n_{1}\times\ldots\times n_{d}$ tensor
${\bm{\mathscr{X}}}$ which minimizes the distance
$\lVert{\bm{\mathscr{T}}}-{\bm{\mathscr{X}}}\rVert_{\mathrm{F}}$ on a set of
low-rank tensors. Yet, the problem of the best rank-$R$ approximation of a
tensor is ill-posed as soon as $R>1$ because the set of rank-$R$ tensors is
not closed (Kolda and Bader, 2009). Instead, we shall consider the best low-
_multilinear_ -rank problem, which is always well-posed,

$\min_{\operatorname{rank}({\bm{\mathscr{X}}})\leqslant(r_{1},\ldots,r_{d})}\left\lVert{\bm{\mathscr{T}}}-{\bm{\mathscr{X}}}\right\rVert_{\mathrm{F}}^{2}.$
(1)

It is well known in the matrix case that the best rank-$R$ approximation can
be easily computed by truncating the SVD to its $R$ most energetic terms
(Eckart and Young, 1936; Mirsky, 1960). Could this also be true for the MLSVD?
Unfortunately, counter-examples exist (Kolda, 2003), showing that there is no
tensor equivalent of the Eckart-Young-Mirsky theorem. Worse still, Problem (1)
is in fact NP-hard (Hillar and Lim, 2013). Nevertheless, despite not being the
best low-multilinear-rank approximation, the truncated MLSVD
$\hat{{\bm{\mathscr{T}}}}$ remains a very good “first guess” as it verifies
$\lVert{\bm{\mathscr{T}}}-{\bm{\mathscr{T}}}_{\star}\rVert_{\mathrm{F}}\leqslant\lVert{\bm{\mathscr{T}}}-\hat{{\bm{\mathscr{T}}}}\rVert_{\mathrm{F}}\leqslant\sqrt{d}\lVert{\bm{\mathscr{T}}}-{\bm{\mathscr{T}}}_{\star}\rVert_{\mathrm{F}}$
where ${\bm{\mathscr{T}}}_{\star}$ denotes a solution to Problem (1) and $d$
is the order of the tensor (Grasedyck et al., 2013; Hackbusch, 2012). It is a
cheap (it consists only in $d$ standard matrix SVDs) and quasi-optimal low-
multilinear-rank approximation of ${\bm{\mathscr{T}}}$. Moreover, it is often
used as an initialization of numerical methods which estimate a solution to
Problem (1), among which the most common is the higher-order orthogonal
iteration (HOOI) algorithm (Kroonenberg and de Leeuw, 1980; Kapteyn et al.,
1986; De Lathauwer et al., 2000a).

Another motivation for the analysis of the low-multilinear-rank approximation
problem is that it has also a practical interest for the numerical computation
of the canonical polyadic decomposition (CPD). Indeed, when dealing with large
tensors, it is computationally more efficient to first compress the tensor
with a low-multilinear-rank approximation and then compute the CPD on the
smaller core tensor rather than computing the CPD of the large tensor directly
(Bro and Andersson, 1998). This is done, e.g., by the cpd function of the
MATLAB toolbox Tensorlab (Vervliet et al., 2016).

###  1.2 Related Work

Multilinear SVD (MLSVD) has a wide range of applications and is often used to
extract relevant information from multi-way arrays. For instance, it has been
used in human motion recognition (Vasilescu, 2002), face recognition
(Vasilescu and Terzopoulos, 2003), handwritten digit classification (Savas and
Eldén, 2007) but also genomics (Omberg et al., 2007, 2009; Muralidhara et al.,
2011) and syndromic surveillance (Fanaee-T and Gama, 2015).

The analysis of _spiked_ tensor models — i.e., low-rank perturbations of large
random tensors — has started with the introduction by Montanari and Richard
(2014) of the rank-$1$ symmetric spiked tensor model,
${\bm{\mathscr{T}}}=\beta{\bm{x}}^{\otimes d}+{\bm{\mathscr{N}}}$ with
$\lVert{\bm{x}}\rVert=1$, ${\bm{\mathscr{N}}}$ Gaussian noise and $\beta$ a
parameter controlling the signal-to-noise ratio (SNR). They showed that
estimation of ${\bm{x}}$ from ${\bm{\mathscr{T}}}$ is _theoretically_ possible
as soon as $\beta$ is above a certain threshold $\beta_{c}$ behaving like
$\sqrt{d\log d}$, which is reminiscent of the now well-known spiked _matrix_
model where signal reconstruction is only possible above a critical threshold
(Péché, 2006) — a phenomenon called the BBP phase transition (Baik et al.,
2005). The behavior of singular values and singular vectors of spiked matrix
models was comprehensively studied by Benaych-Georges and Nadakuditi (2012).
Contrary to the matrix case however, Montanari and Richard (2014) made the
disturbing observation that none of the polynomial-time estimation algorithms
among tensor unfolding, power iteration and approximate message passing (AMP)
succeed unless $\beta$ diverges as the dimensions of the tensor grow large.
While Perry et al. (2017) showed that the information-theoretic threshold is
of order $1$, this suggests the existence of a computational-to-statistical
gap in spiked tensor estimation, as in a myriad of other problems (Bandeira et
al., 2018; Zdeborová and Krzakala, 2016).

The landscape of the rank-$1$ symmetric spiked tensor model was studied by Ben
Arous et al. (2019b), who showed that the number of local optima to Problem
(1) grows exponentially with the size of the tensor, but all lie close to a
subspace orthogonal to the sought solution, except for one if $\beta$ exceeds
a critical threshold $\beta_{c}$. Completing this analysis, Jagannath et al.
(2019) showed the existence of two close but different thresholds
$\beta_{s}<\beta_{c}$ such that the solution aligned with the underlying
signal is a local minimum of Problem (1) as soon as $\beta>\beta_{s}$ but
becomes a _global_ one only if $\beta>\beta_{c}$.

So far, we have only referred to works dealing with the rank-$1$ symmetric
case, but there has also been some studies about higher-(low-)rank spiked
models. Chevreuil and Loubaton (2018) give a sufficient (but not necessary)
condition for the _non-detectability_ of a rank-$R$ asymmetric signal
perturbed by an additive Gaussian noise. Chen et al. (2021) also discuss
signal detectability in the rank-$R$ symmetric case. Zhang and Xia (2020)
considered a general $\textit{low-multilinear-rank signal}\leavevmode\nobreak\
{\bm{\mathscr{P}}}+\textit{Gaussian noise}\leavevmode\nobreak\
{\bm{\mathscr{N}}}$ model and brought to light the same statistical-to-
computational gap: if $\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}$ is above a
statistical threshold of order $1$ then Problem (1) has a solution which is
aligned with the signal but is computationally intractable unless
$\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}$ is above a computational
threshold of order $N^{\frac{d-2}{4}}$, where $N$ scales as the dimensions of
the data tensor. In this strong SNR regime, the higher-order orthogonal
iteration (HOOI) algorithm (De Lathauwer et al., 2000a) is minimax-optimal. In
fact, it was already showed by Ben Arous et al. (2019a) that the computational
threshold should behave like $N^{\alpha}$ with $\alpha>0$. AMP and tensor
power iteration algorithms achieve $N^{\frac{d-1}{2}}$ (Lesieur et al., 2017;
Huang et al., 2020) while tensor unfolding methods (truncated MLSVD and HOOI
algorithm) achieve $N^{\frac{d-2}{4}}$ as already conjectured by Montanari and
Richard (2014) and later proven by Ben Arous et al. (2021). The convergence of
the HOOI algorithm towards a local maximum for a _sufficiently close_
initialization was proven by Xu (2018).

Recently, a new approach relying on tools from random matrix theory has
broaden the understanding of spiked tensor models. In particular, Goulart et
al. (2022) studied the rank-$1$ symmetric case and were able to recover
explicitly the same $\beta_{s}$ threshold as Jagannath et al. (2019) as well
as to precisely quantify the alignment between a solution to Problem (1) and
the signal. A similar analysis was carried out by Seddik et al. (2022) for the
more general asymmetric case, relying solely on classical techniques from
random matrix theory. Such tools show promise for the theoretical
understanding of learning from tensor data (Seddik et al., 2023).

###  1.3 Summary of Contributions

Consider the general spiked tensor model,

${\bm{\mathscr{T}}}={\bm{\mathscr{P}}}+\frac{1}{\sqrt{N}}{\bm{\mathscr{N}}}\quad\in{\mathbb{R}}^{n_{1}\times\ldots\times
n_{d}},\qquad{\mathscr{N}}_{i_{1},\ldots,i_{d}}\overset{\text{i.i.d.}}{\sim}{\mathcal{N}}(0,1),$
(2)

where ${\bm{\mathscr{N}}}$ is an additive Gaussian noise, $N$ is a parameter
controlling the size of the tensor such that
$n_{\ell}/N=\operatorname*{\Theta}(1)$333That is, $n_{\ell}/N$ neither
vanishes nor diverges as $N\to+\infty$. This ensures that the spectral norm of
$\frac{1}{\sqrt{N}}{\bm{\mathscr{N}}}$ is of order $1$ (Tomioka and Suzuki,
2014). for all $\ell\in\\{1,\ldots,d\\}$ (for instance, $N=n_{1}$ or
$N=\sum_{\ell=1}^{d}n_{\ell}$) and ${\bm{\mathscr{P}}}$ is a low-multilinear-
rank deterministic tensor, i.e., which can be decomposed as

${\bm{\mathscr{P}}}=\sum_{q_{1}=1}^{r_{1}}\ldots\sum_{q_{d}=1}^{r_{d}}{\mathscr{H}}_{q_{1},\ldots,q_{d}}[{\bm{x}}^{(1)}_{q_{1}}\otimes\ldots\otimes{\bm{x}}^{(d)}_{q_{d}}]\leavevmode\nobreak\
\overset{\text{def}}{=}\leavevmode\nobreak\
\left\llbracket{\bm{\mathscr{H}}};{\bm{X}}^{(1)},\ldots,{\bm{X}}^{(d)}\right\rrbracket,$
(3)

with ${\bm{\mathscr{H}}}\in{\mathbb{R}}^{r_{1}\times\ldots\times r_{d}}$ and
${\bm{X}}^{(\ell)}$ an $n_{\ell}\times r_{\ell}$ matrix with orthonormal
columns ${\bm{x}}^{(\ell)}_{q_{\ell}}$. The range of ${\bm{X}}^{(\ell)}$ is
the $\ell$-th singular subspace of ${\bm{\mathscr{P}}}$. This decomposition is
illustrated for the case $d=3$ in Figure 1. Model (2) with
${\bm{\mathscr{P}}}$ as in (3) is the most general spiked tensor model — i.e.,
low-rank perturbation of a large random tensor. Indeed, any of the models
referred to in the previous Section 1.2 fall into this definition since
decomposition (3) always exists and low CPD-rank is equivalent to low
multilinear-rank thanks to the inequality $\max_{\ell}\\{r_{\ell}\\}\leqslant
R\leqslant\min_{\ell\neq\ell^{\prime}}\\{r_{\ell}r_{\ell^{\prime}}\\}$
(Sidiropoulos et al., 2017).

$n_{1}$$n_{2}$$n_{3}$${\bm{\mathscr{P}}}$$=$$r_{1}$$r_{2}$$r_{3}$${\bm{\mathscr{H}}}$$r_{1}$$n_{1}$${\bm{X}}^{(1)}$$r_{2}$$n_{2}$${\bm{X}}^{(2)}$$r_{3}$$n_{3}$${\bm{X}}^{(3)}$
Figure 1: Illustration of the Tucker decomposition (3) of an $n_{1}\times
n_{2}\times n_{3}$ tensor ${\bm{\mathscr{P}}}$ with multilinear-rank
$(r_{1},r_{2},r_{3})$. ${\bm{\mathscr{H}}}$ is the $r_{1}\times r_{2}\times
r_{3}$ core tensor and ${\bm{X}}^{(1)},{\bm{X}}^{(2)},{\bm{X}}^{(3)}$ are
matrices with orthonormal columns spanning the singular subspaces of
${\bm{\mathscr{P}}}$.

In the regime where $N\to+\infty$ — representing the fact that, in practice,
the dimensions of the tensor are large compared to its rank —, we study the
estimation of ${\bm{\mathscr{P}}}$ from ${\bm{\mathscr{T}}}$ with a truncated
MLSVD, which serves as initialization of the HOOI algorithm. In particular, we
reveal that the interesting _non-trivial_ regime is characterized by the
$\operatorname*{\Theta}(N^{\frac{d-2}{2}})$ quantity
$\sigma_{N}=\frac{1}{N}\prod_{\ell=1}^{d}\sqrt{n_{\ell}}$.

  * •

If
$\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}^{2}/\sigma_{N}\xrightarrow[N\to+\infty]{}0$,
then the noise completely masks the signal, and truncated MLSVD fails to
recover ${\bm{\mathscr{P}}}$.

  * •

If
$\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}^{2}/\sigma_{N}\xrightarrow[N\to+\infty]{}+\infty$,
then the signal clearly stands out from the noise, and reconstruction of
${\bm{\mathscr{P}}}$ with a truncated MLSVD is easy.

  * •

If
$\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}^{2}/\sigma_{N}=\operatorname*{\Theta}(1)$
as $N\to+\infty$, then we are precisely in the _non-trivial_ regime between
the two previous situations, and truncated MLSVD may partially recover
${\bm{\mathscr{P}}}$.

It is the analysis of this last regime which is of practical interest. Given
the low-multilinear-rank approximation
$\hat{{\bm{\mathscr{T}}}}=\llbracket\hat{{\bm{\mathscr{G}}}};\hat{{\bm{U}}}^{(1)},\ldots,\hat{{\bm{U}}}^{(d)}\rrbracket$
obtained with a truncated MLSVD of ${\bm{\mathscr{T}}}$, we quantify how well
$\hat{{\bm{\mathscr{T}}}}$ reconstructs ${\bm{\mathscr{P}}}$ in this non-
trivial regime. To do so, we study the spectral properties of the unfoldings
(i.e., matricizations) of the tensor ${\bm{\mathscr{T}}}$, i.e., the
$n_{\ell}\times\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}$ matrices
${\bm{T}}^{(\ell)}$ whose columns are mode-$\ell$ fibers of
${\bm{\mathscr{T}}}$ and the columns of $\hat{{\bm{U}}}^{(\ell)}$ are its
dominant left singular vectors. Such _long_ matrices (the second dimension
grows faster than the first one) have already been studied by Ben Arous et al.
(2021) in order to analyze the properties of tensor-unfolding methods in the
particular setting of a rank-$1$ spike. Here, we tackle this problem with a
highly different approach relying solely on classical tools from random matrix
theory (Couillet and Liao, 2022) and give very general results on the spiked
tensor model that go beyond the specific rank-$1$ one case. Moreover, we
justify the practical use of truncated MLSVD as an initialization of the HOOI
algorithm by showing its optimality in the large $N$ regime.

Although the spectrum of ${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}$,
$\ell\in\\{1,\ldots,d\\}$, diverges as $N\to+\infty$, we show that its
eigenvalues (i.e., the squared singular values of ${\bm{T}}^{(\ell)}$) gather
in an interval $[\mu^{(\ell)}_{N}\pm 2\sigma_{N}]$ with
$\mu^{(\ell)}_{N}=\frac{1}{N}\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}=\operatorname*{\Theta}(N^{d-2})$.
More precisely, the empirical spectral distribution of the rescaled matrix
$\frac{1}{\sigma_{N}}[{\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}-\mu^{(\ell)}_{N}{\bm{I}}_{n_{\ell}}]$
converges weakly to the semicircle distribution on $[-2,+2]$ (see Figure 3,
Theorem 5 and Corollary 6). Furthermore, we show that a BBP phase transition
phenomenon occur: for each singular value of ${\bm{P}}^{(\ell)}$ (the
unfolding of ${\bm{\mathscr{P}}}$ along mode $\ell$) which is above the
threshold $\sqrt{\sigma_{N}}$, an eigenvalue of
${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}$ isolates itself on the right side of
the _bulk_ (see Figure 3) and its corresponding eigenvector (i.e., left
singular vector of ${\bm{T}}^{(\ell)}$) is aligned with the corresponding
singular subspace of ${\bm{\mathscr{P}}}$. The position of the isolated
eigenvalue and this alignment are efficiently predicted by Theorem 8 (see also
Figure 3).

$0$$2$$4$$6$$8$$10$$12$$14$$16$$18$$20$$0$$0.2$$0.4$$0.6$$0.8$$1$$\omega$
(SNR)Subspace alignmentMode $\ell$ $\rightarrow$$1$$2$$3$T-MLSVD (th)T-MLSVD
(sim)HOOI (sim) Figure 2: Alignments between singular subspaces (see Section
3.2) of the observation
${\bm{\mathscr{T}}}=\sqrt{\omega}{\bm{\mathscr{P}}}_{\circ}+\frac{1}{\sqrt{N}}{\bm{\mathscr{N}}}$
and of the signal ${\bm{\mathscr{P}}}_{\circ}$, with
$\lVert{\bm{\mathscr{P}}}_{\circ}\rVert_{\mathrm{F}}^{2}=\frac{\sqrt{n_{1}n_{2}n_{3}}}{N}$,
as a function of the signal-to-noise ratio $\omega$. Theoretical alignments
(Theorem 8) achieved with truncated MLSVD are compared with simulations and
those achieved with the HOOI algorithm. Empirical results are averaged over
$10$ trials, with error bars representing standard deviation. Experimental
setting: $d=3$, $(n_{1},n_{2},n_{3})=(100,200,300)$, $N=n_{1}+n_{2}+n_{3}$ and
$(r_{1},r_{2},r_{3})=(3,4,5)$.

As a result, Figure 2 plots, for an order-$3$ tensor, as a function of the
signal-to-noise ratio (SNR)
$\omega=\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}^{2}/\sigma_{N}$, the
alignments between the singular subspace of the signal ${\bm{\mathscr{P}}}$
spanned by ${\bm{X}}^{(\ell)}$ and the dominant singular subspace of the
observation ${\bm{\mathscr{T}}}$ spanned by $\hat{{\bm{U}}}^{(\ell)}$. Solid
curves are the alignments predicted by Theorem 8 while dotted curves are
empirical alignments computed on a $100\times 200\times 300$ tensor with
signal-rank $(3,4,5)$. If the SNR $\omega$ is too small, there is no
alignment, meaning that truncated MLSVD fails to recover ${\bm{\mathscr{P}}}$
— the signal is masked by the noise. When it exceeds a critical value (see
Theorem 8 and Section 3.2 for details), a phase transition phenomenon
occurs444In fact, we will see in Section 3.2 that there is one phase
transition for each principal direction of the singular subspaces of
${\bm{\mathscr{P}}}$, resulting in $\sum_{\ell=1}^{d}r_{\ell}$ phase
transitions. Their positions corresponds to sudden changes of slope in the
solid curves of Figure 2.: the alignment starts to grow — i.e., truncated
MLSVD now partially recovers ${\bm{\mathscr{P}}}$ — and converges to $1$ as
$\omega\to+\infty$.

Besides, Figure 2 also plots the empirical alignments between the singular
subspaces of ${\bm{\mathscr{P}}}$ and those estimated with the HOOI algorithm
(De Lathauwer et al., 2000a) given in Algorithm 1, whose truncated MLSVD
serves as initialization. This yields much better alignments, especially close
to the phase transition. In fact, we show in Theorem 13 that the HOOI
algorithm converges to a solution to Problem (2) as soon as its initialization
sufficiently preserves the underlying signal. This provides new insight into
the computational barrier: initialization is the limiting factor here. Had one
prior information on the solution, one could initialize the HOOI algorithm in
the right basin of attraction and still be able to perfectly (i.e., with
alignment $1$) reconstruct the signal in the regime
$1\ll\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}\ll N^{\frac{d-2}{4}}$, which
is computationally hard but statistically easy (see details in Section 4.1 and
discussion in Section 4.2).

In a nutshell, our contributions could be summarized as follows.

  * •

We characterize, in the large $N$ limit, the behavior of the singular values
of the unfoldings of the tensor ${\bm{\mathscr{T}}}$ — denoted
${\bm{T}}^{(\ell)}$, $\ell\in\\{1,\ldots,d\\}$ — when it follows the general
spiked tensor model (2) (Theorem 5 and Corollary 6). This is performed through
the analysis of the limiting spectral distribution of the symmetric matrix
${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}$ using standard tools from the theory
of large random matrices.

  * •

We give a precise condition, depending on a
$\operatorname*{\Theta}(N^{\frac{d-2}{4}})$ threshold on the signal, for the
detectability of a principal direction of the $\ell$-th singular subspace of
${\bm{\mathscr{P}}}$ from the unfolding ${\bm{T}}^{(\ell)}$. This corresponds
to the presence of an isolated eigenvalue in the spectrum of
${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}$ with associated eigenvector aligned
with the sought singular subspace. We specify the asymptotic position of this
isolated eigenvalue, as well as the quality of the alignment (Theorem 8).

  * •

Relying on our random matrix analysis, we characterize the performance of
truncated MLSVD in the reconstruction of the signal ${\bm{\mathscr{P}}}$ from
the observation ${\bm{\mathscr{T}}}$.

  * •

We show that exact reconstruction of ${\bm{\mathscr{P}}}$ from
${\bm{\mathscr{T}}}$ is possible with the HOOI algorithm (De Lathauwer et al.,
2000a) as long as $\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}\gg 1$ and it is
initialized in the right bassin of attraction (Theorem 13). Without prior
information this depends on the detectability of ${\bm{\mathscr{P}}}$ in the
truncated MLSVD of ${\bm{\mathscr{T}}}$, which is only possible above the
$\operatorname*{\Theta}(N^{\frac{d-2}{4}})$ computational threshold. Moreover,
as $N\to+\infty$, the number of iterations needed for the convergence of the
algorithm converges to $1$.

In section 2, we introduce our notations, tensor-related operations and
decompositions and useful tools from random matrix theory. Section 3 presents
the random matrix analysis of long matrices emerging from the unfoldings of
tensors following the general spiked tensor model (2). These results are
presented in the context of truncated MLSVD and exploited to quantitatively
explain its reconstruction performances. Then, relying on these results,
Section 4 deals with the numerical estimation of a solution to Problem (1)
with the HOOI algorithm. We show its asymptotic optimality and provide insight
into the limiting factors for numerical estimation below the computational
threshold. We conclude and discuss our results in Section 5. Most proofs are
deferred to the appendix.

##  2 Preliminaries on Tensors and Random Matrix Theory

We start by introducing some notation and the main tools which are needed to
expose the results of the next sections.

###  2.1 General Notations

$a$, , ${\bm{A}}$ and ${\bm{\mathscr{A}}}$ respectively denote a scalar, a
vector, a matrix and a tensor. $a_{i}$, $A_{i,j}$ and
${\mathscr{A}}_{i_{1},\ldots,i_{d}}$ are their entries. For
$x\in{\mathbb{R}}$, $[x]^{+}=\max(0,x)$. The imaginary part of
$z\in{\mathbb{C}}$ is $\Im[z]$. The set $\\{1,\ldots,n\\}$ of positive
integers smaller or equal to $n$ is denoted $[n]$. $\delta_{x}$ is the Dirac
measure at point $x$. The support of a probability measure $\mu$ is denoted
$\operatorname{Supp}\mu$. The notation $X\sim{\mathcal{L}}$ means that the
random variable $X$ is distributed according to the law ${\mathcal{L}}$. Given
a sequence of random variables $(X_{n})_{n\geqslant 0}$, its convergence in
distribution to ${\mathcal{L}}$ is denoted
$X_{n}\xrightarrow[n\to+\infty]{{\mathcal{D}}}{\mathcal{L}}$ and its almost
sure convergence to $L$ is denoted
$X_{n}\xrightarrow[n\to+\infty]{\text{a.s.}}L$. The normal distribution with
mean $\mu$ and variance $\sigma^{2}$ is denoted
${\mathcal{N}}(\mu,\sigma^{2})$. The span of an $n_{1}\times n_{2}$ matrix
${\bm{A}}$ is
$\operatorname{Span}{\bm{A}}=\\{{\bm{A}}{\bm{x}}\mid{\bm{x}}\in{\mathbb{R}}^{n_{2}}\\}\subset{\mathbb{R}}^{n_{1}}$.
The singular values of ${\bm{A}}$ in _non-increasing_ order are denoted
$s_{1}({\bm{A}})\geqslant s_{2}({\bm{A}})\geqslant\ldots\geqslant 0$. Given an
$n\times n$ matrix ${\bm{B}}$, its trace is
$\operatorname{Tr}{\bm{B}}=\sum_{i=1}^{n}B_{i,i}$ and its spectrum,
$\operatorname{Sp}{\bm{B}}$, is the set of all its eigenvalues.
$\lVert\cdot\rVert$ denotes the standard Euclidean norm for vectors and the
corresponding operator norm (spectral norm) for matrices.

Given two real sequences $(u_{n})_{n\geqslant 0}$ and $(v_{n})_{n\geqslant
0}$, we write $u_{n}=\operatorname*{\mathcal{O}}_{n\to+\infty}(v_{n})$ if
there exists a constant $C>0$ and an integer $n_{0}$ such that $\lvert
u_{n}\rvert\leqslant C\lvert v_{n}\rvert$ as soon as $n\geqslant n_{0}$. If
$u_{n}=\operatorname*{\mathcal{O}}_{n\to+\infty}(v_{n})$ and
$v_{n}=\operatorname*{\mathcal{O}}_{n\to+\infty}(u_{n})$555I.e., if there
exist two constants $c,C>0$ and an integer $n_{0}$ such that $c\lvert
v_{n}\rvert\leqslant\lvert u_{n}\rvert\leqslant C\lvert v_{n}\rvert$ as soon
as $n\geqslant n_{0}$., then we write
$u_{n}=\operatorname*{\Theta}_{n\to+\infty}(v_{n})$. We also write
$u_{n}\ll_{n\to+\infty}v_{n}$ (or $v_{n}\gg_{n\to+\infty}u_{n}$) if, for all
$\varepsilon>0$, there exists an integer $n_{0}$ such that $\lvert
u_{n}\rvert\leqslant\varepsilon\lvert v_{n}\rvert$ as soon as $n\geqslant
n_{0}$. When it is clear from context that $n\to+\infty$, we simply write
$u_{n}=\operatorname*{\Theta}(v_{n})$,
$u_{n}=\operatorname*{\mathcal{O}}(v_{n})$, $u_{n}\ll v_{n}$ or $u_{n}\gg
v_{n}$.

Unless stated otherwise, $d$ represents the order of a tensor; $\ell$ is an
index ranging from $1$ to $d$ and $i_{\ell}$, $q_{\ell}$ range from $1$ to
$n_{\ell}$, $r_{\ell}$ respectively.

###  2.2 Tensors, Related Operations and Decompositions

For our purposes, tensors are considered as multi-way arrays:
${\bm{\mathscr{T}}}\in{\mathbb{R}}^{n_{1}\times\ldots\times n_{d}}$ is a
collection of elements ${\mathscr{T}}_{i_{1},\ldots,i_{d}}\in{\mathbb{R}}$
with $i_{\ell}\in[n_{\ell}]$, $\ell\in[d]$. An _$\ell$ -fiber_ of
${\bm{\mathscr{T}}}$ is the vector of ${\mathbb{R}}^{n_{\ell}}$ obtained by
fixing all the indices of ${\bm{\mathscr{T}}}$ but the $\ell$-th. This is the
generalization of rows and columns of matrices, which are respectively $1$\-
and $2$-fibers. _Unfolding_ (or matricization) is the process by which a
matrix is built from a tensor — ${\bm{T}}^{(\ell)}$ is the unfolding of
${\bm{\mathscr{T}}}$ along mode $\ell$, i.e., the
$n_{\ell}\times\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}$ matrix whose
columns are $\ell$-fibers of ${\bm{\mathscr{T}}}$666The order of the columns
does not matter, as long as other operations, such as the Kronecker product,
are defined in a consistent manner..

Given two $n_{1}\times\ldots\times n_{d}$ tensors
${\bm{\mathscr{A}}},{\bm{\mathscr{B}}}$, their scalar product is

$\left\langle{\bm{\mathscr{A}}},{\bm{\mathscr{B}}}\right\rangle_{\mathrm{F}}\overset{\text{def}}{=}\sum_{i_{1}=1}^{n_{1}}\dots\sum_{i_{d}=1}^{n_{d}}{\mathscr{A}}_{i_{1},\ldots,i_{d}}{\mathscr{B}}_{i_{1},\ldots,i_{d}}.$

The Frobenius norm of ${\bm{\mathscr{A}}}$ is
$\lVert{\bm{\mathscr{A}}}\rVert_{\mathrm{F}}\overset{\text{def}}{=}\sqrt{\langle{\bm{\mathscr{A}}},{\bm{\mathscr{A}}}\rangle_{\mathrm{F}}}$.

Just as a matrix is said to have rank $1$ if it can be expressed as the _outer
product_ of two vectors,
${\bm{x}}{\bm{y}}^{\top}\overset{\text{def}}{=}{\bm{x}}\otimes{\bm{y}}$, an
order-$d$ tensor is said to have rank $1$ if it can be expressed as the outer
product of $d$ vectors, i.e.,
${\mathscr{T}}_{i_{1},\ldots,i_{d}}=x^{(1)}_{i_{1}}\ldots
x^{(d)}_{i_{d}}\iff{\bm{\mathscr{T}}}\overset{\text{def}}{=}\bigotimes_{\ell=1}^{d}{\bm{x}}^{(\ell)}$.

Given two matrices ${\bm{A}}$ and ${\bm{B}}$ of respective sizes $n_{1}\times
n_{2}$ and $p_{1}\times p_{2}$, their _Kronecker product_ , denoted
${\bm{A}}\boxtimes{\bm{B}}$, is the $n_{1}p_{1}\times n_{2}p_{2}$ matrix such
that
$[{\bm{A}}\boxtimes{\bm{B}}]_{p_{1}(i_{1}-1)+j_{1},p_{2}(i_{2}-1)+j_{2}}=A_{i_{1},i_{2}}B_{j_{1},j_{2}}$
for all
$(i_{1},i_{2})\in[n_{1}]\times[n_{2}],(j_{1},j_{2})\in[p_{1}]\times[p_{2}]$.
It is also defined for $n$-dimensional vectors, seen as $n\times 1$ matrices.
This product is useful to express unfoldings of tensors defined as outer
products. For instance,

$[{\bm{x}}\otimes{\bm{y}}\otimes{\bm{z}}]^{(1)}={\bm{x}}({\bm{y}}\boxtimes{\bm{z}})^{\top},\quad[{\bm{x}}\otimes{\bm{y}}\otimes{\bm{z}}]^{(2)}={\bm{y}}({\bm{x}}\boxtimes{\bm{z}})^{\top},\quad[{\bm{x}}\otimes{\bm{y}}\otimes{\bm{z}}]^{(3)}={\bm{z}}({\bm{x}}\boxtimes{\bm{y}})^{\top}.$
(4)

Among the various properties of the Kronecker product (Abadir and Magnus,
2005, Chapter 10), we highlight that it is bilinear, associative, _non_
-commutative,
$({\bm{A}}\boxtimes{\bm{B}})^{\top}={\bm{A}}^{\top}\boxtimes{\bm{B}}^{\top}$
and
$({\bm{A}}\boxtimes{\bm{B}})({\bm{C}}\boxtimes{\bm{D}})=({\bm{A}}{\bm{C}})\boxtimes({\bm{B}}{\bm{D}})$
when the matrix products ${\bm{A}}{\bm{C}}$ and ${\bm{B}}{\bm{D}}$ are
defined.

Because of their multi-dimensional nature, tensors can be very large and cause
computational and storage difficulties. Hence, tensor decompositions are an
appropriate way to provide a sparse representation of these objects and reveal
relevant knowledge from multi-way data arrays. We briefly present two of the
most common decompositions, namely CPD and MLSVD, although many others exist
(see, e.g., Kolda and Bader (2009); Vervliet et al. (2014)).

#### Canonical Polyadic Decomposition (CPD).

${\bm{\mathscr{T}}}$ is decomposed as a sum of $R$ rank-$1$ terms,
${\bm{\mathscr{T}}}=\sum_{i=1}^{R}\sigma_{i}\bigotimes_{\ell=1}^{d}{}^{(\ell)}_{i}$,
with ${}^{(\ell)}_{i}$ of unit norm, $\sigma_{i}>0$ and minimal $R$.
Introduced by Hitchcock (1927), and also called CANDECOMP or PARAFAC, this
decomposition is, up to permutation, essentially unique under mild conditions
(Kolda and Bader, 2009).

#### Multilinear Singular Value Decomposition (MLSVD).

${\bm{\mathscr{T}}}$ is written as the multilinear tensor-matrix product of an
$r_{1}\times\ldots\times r_{d}$ core tensor ${\bm{\mathscr{G}}}$ with $d$
factors ${\bm{U}}^{(\ell)}$, $\ell\in[d]$, of respective sizes $n_{\ell}\times
r_{\ell}$ and orthonormal columns ${\bm{u}}^{(\ell)}_{q_{\ell}}$,

${\bm{\mathscr{T}}}=\sum_{q_{1}=1}^{r_{1}}\dots\sum_{q_{d}=1}^{r_{d}}{\mathscr{G}}_{q_{1},\ldots,q_{d}}\bigotimes_{\ell=1}^{d}{\bm{u}}^{(\ell)}_{q_{\ell}}\leavevmode\nobreak\
\overset{\text{def}}{=}\leavevmode\nobreak\
\left\llbracket{\bm{\mathscr{G}}};{\bm{U}}^{(1)},\ldots,{\bm{U}}^{(d)}\right\rrbracket.$

$r_{\ell}$ is the dimension of the subspace spanned by the $\ell$-fibers of
${\bm{\mathscr{T}}}$, whose an orthonormal basis is formed by the columns of
${\bm{U}}^{(\ell)}$. This decomposition, introduced by Tucker (1966), is also
called higher-order SVD (HOSVD, De Lathauwer et al., 2000b). From this
decomposition, the unfoldings of ${\bm{\mathscr{T}}}$ can be expressed with
matrix and Kronecker products between the factor matrices and the unfoldings
of the core tensor ${\bm{\mathscr{G}}}$. For instance, if
${\bm{\mathscr{T}}}=\llbracket{\bm{\mathscr{G}}};{\bm{U}},{\bm{V}},{\bm{W}}\rrbracket$,

${\bm{T}}^{(1)}={\bm{U}}{\bm{G}}^{(1)}({\bm{V}}\boxtimes{\bm{W}})^{\top},\quad{\bm{T}}^{(2)}={\bm{V}}{\bm{G}}^{(2)}({\bm{U}}\boxtimes{\bm{W}})^{\top},\quad{\bm{T}}^{(3)}={\bm{W}}{\bm{G}}^{(3)}({\bm{U}}\boxtimes{\bm{V}})^{\top}.$

Note how these expressions generalize (4).

####

The _contraction_ of ${\bm{\mathscr{T}}}$ on
${\bm{A}}^{(1)},\ldots,{\bm{A}}^{(d)}$ with
${\bm{A}}^{(\ell)}\in{\mathbb{R}}^{n_{\ell}\times p_{\ell}}$, $\ell\in[d]$, is
the $p_{1}\times\ldots\times p_{d}$ tensor
${\bm{\mathscr{T}}}({\bm{A}}^{(1)},\ldots,{\bm{A}}^{(d)})$ whose
$(j_{1},\ldots,j_{d})$-entry is
$\sum_{i_{1}=1}^{n_{1}}\dots\sum_{i_{d}=1}^{n_{d}}{\mathscr{T}}_{i_{1},\ldots,i_{d}}\prod_{\ell=1}^{d}A^{(\ell)}_{i_{\ell},j_{\ell}}$.
If
${\bm{\mathscr{T}}}=\llbracket{\bm{\mathscr{G}}};{\bm{U}}^{(1)},\ldots,{\bm{U}}^{(d)}\rrbracket$
then
${\bm{\mathscr{G}}}=\llbracket{\bm{\mathscr{T}}};{\bm{U}}^{(1)\top},\ldots,{\bm{U}}^{(d)\top}\rrbracket={\bm{\mathscr{T}}}({\bm{U}}^{(1)},\ldots,{\bm{U}}^{(d)})$.
The tensor contraction generalizes the matrix operation
${\bm{A}}^{(1)\top}{\bm{M}}{\bm{A}}^{(2)}={\bm{M}}({\bm{A}}^{(1)},{\bm{A}}^{(2)})$.

######  Remark 1 (Non-unicity of the Tucker decomposition)

For all $r_{\ell}\times r_{\ell}$ orthogonal matrices ${\bm{O}}^{(\ell)}$,
$\ell\in[d]$, we have the equivalent decomposition

$\left\llbracket{\bm{\mathscr{G}}};{\bm{U}}^{(1)},\ldots,{\bm{U}}^{(d)}\right\rrbracket=\left\llbracket\breve{{\bm{\mathscr{G}}}};{\bm{U}}^{(1)}{\bm{O}}^{(1)},\ldots,{\bm{U}}^{(d)}{\bm{O}}^{(d)}\right\rrbracket\quad\text{with}\quad\breve{{\bm{\mathscr{G}}}}$

where
$\breve{{\bm{\mathscr{G}}}}={\bm{\mathscr{G}}}({\bm{O}}^{(1)},\ldots,{\bm{O}}^{(d)})$
is the contraction of ${\bm{\mathscr{G}}}$ on
${\bm{O}}^{(1)},\ldots,{\bm{O}}^{(d)}$.

The actual rank of a tensor is usually understood as its CPD-rank $R$ whereas
its MLSVD-rank $(r_{1},\ldots,r_{d})$ is called _multilinear-rank_. These two
ranks are constrained by the following inequalities,

$\max_{1\leqslant\ell\leqslant d}\\{r_{\ell}\\}\leqslant
R\leqslant\min_{1\leqslant\ell\neq\ell^{\prime}\leqslant
d}\\{r_{\ell}r_{\ell^{\prime}}\\},$ (5)

which implies one can generally speak of a “low-rank tensor” without having to
mention the rank to which it refers.

###  2.3 Tools from Random Matrix Theory

The results presented in this work rely on tools from the theory of large
random matrices. Given a certain $n\times n$ matrix ${\bm{M}}$ whose entries
are random variables, one is interested in the behavior its eigenvalues and
eigenvectors as $n\to+\infty$. Common questions are, in this regime, how does
the _empirical spectral distribution_ (ESD) of ${\bm{M}}$, i.e.,
$\frac{1}{n}\sum_{\lambda\in\operatorname{Sp}{\bm{M}}}\delta_{\lambda}$,
behave? Does it (weakly) converge to a _limiting spectral distribution_ (LSD)?
Are there any isolated eigenvalues? If so, are the corresponding eigenvectors
aligned with a relevant subspace?

In case ${\bm{M}}$ is symmetric, its eigenvalues and eigenvectors are real and
a key tool to answer the preceding questions is its resolvent,
${\bm{Q}}_{\bm{M}}(z)=\left({\bm{M}}-z{\bm{I}}_{n}\right)^{-1}$, defined for
all $z\in{\mathbb{C}}\setminus\operatorname{Sp}{\bm{M}}$. Indeed,
$\frac{1}{n}\operatorname{Tr}{\bm{Q}}_{\bm{M}}(z)=\frac{1}{n}\sum_{\lambda\in\operatorname{Sp}{\bm{M}}}\frac{1}{\lambda-z}$
is the _Stieltjes transform_ of its empirical spectral distribution.

######  Definition 2 (Stieltjes transform)

The Stieltjes transform of a real probability measure $\mu$ is defined for all
$z\in{\mathbb{C}}\setminus\operatorname{Supp}\mu$ as
$m_{\mu}(z)=\int_{\mathbb{R}}\frac{1}{t-z}{\mathrm{d}}\mu(t)$.

The knowledge of $m_{\mu}$ is equivalent to the knowledge of $\mu$ thanks to
the inversion formula $\mu([a,b])=\frac{1}{\pi}\lim_{\eta\downarrow
0}\int_{a}^{b}\Im[m_{\mu}(x+{\mathrm{i}}\eta)]\leavevmode\nobreak\
{\mathrm{d}}x$ for any continuity points $a,b$ of $\mu$ (Couillet and Liao,
2022, Theorem 2.1). Hence, following the behavior of
$\frac{1}{n}\operatorname{Tr}{\bm{Q}}_{\bm{M}}(z)$ as $n\to+\infty$ gives
insight into the limiting spectral distribution of ${\bm{M}}$ via its
Stieltjes transform. To do so, we seek first a _deterministic equivalent_ of
${\bm{Q}}_{\bm{M}}$.

######  Definition 3 (Matrix equivalent)

For two (random) matrices ${\bm{X}},{\bm{Y}}\in{\mathbb{R}}^{n\times n}$, we
write ${\bm{X}}\leftrightarrow{\bm{Y}}$ if, for any deterministic matrix
${\bm{A}}\in{\mathbb{R}}^{n\times n}$ and vectors $\in{\mathbb{R}}^{n}$ of
bounded norms (spectral and Euclidean norms respectively)777Strictly speaking,
we refer to _sequences_ of matrices and vectors whose size grow with the index
$n$. Yet, this mathematical aspect is made implicit in order to simplify
notation and be closer to practical considerations where these sequences, in
fact, do not exist — the assumption $n\to+\infty$ models the fact that $n$ is
large but remains finite.,

$\frac{1}{n}\operatorname{Tr}{\bm{A}}\left({\bm{X}}-{\bm{Y}}\right)\xrightarrow[n\to+\infty]{\text{a.s.}}0,\quad^{\top}\left({\bm{X}}-{\bm{Y}}\right)\xrightarrow[n\to+\infty]{\text{a.s.}}0\quad\text{and}\quad\left\lVert\mathbb{E}\left[{\bm{X}}-{\bm{Y}}\right]\right\rVert\xrightarrow[n\to+\infty]{\text{a.s.}}0.$

If ${\bm{Y}}$ is deterministic, then it is said to be a _deterministic
equivalent_ of ${\bm{X}}$.

The following lemma will be extensively used to derive such equivalents.

######  Lemma 4 (Stein, 1981)

Let $Z\sim{\mathcal{N}}(0,1)$ and $f:{\mathbb{R}}\to{\mathbb{R}}$ be a
continuously differentiable function. When the following expectations exist,
$\mathbb{E}[Zf(Z)]=\mathbb{E}[f^{\prime}(Z)]$.

Once a deterministic equivalent $\bar{{\bm{Q}}}_{\bm{M}}$ of
${\bm{Q}}_{\bm{M}}$ is found, the limiting spectral distribution of ${\bm{M}}$
is accessible through the limit of
$\frac{1}{n}\operatorname{Tr}\bar{{\bm{Q}}}_{\bm{M}}$ as $n\to+\infty$. In the
study of spiked models, eigenvalues are usually gathered in a _bulk_ ,
described by the limiting spectral distribution (LSD), with a finite number of
isolated eigenvalues, which do not appear in the LSD (see Figure 3 for
example). The position of these isolated eigenvalues are singular points of
the resolvent ${\bm{Q}}_{\bm{M}}$. Hence, our deterministic equivalent
provides equations defining their almost sure asymptotic position. Moreover,
thanks to Cauchy’s integral formula, for all $\in{\mathbb{R}}^{n}$,
$\lvert^{\top}{\bm{u}}\rvert^{2}=\frac{-1}{2{\mathrm{i}}\pi}\oint_{\gamma}^{\top}{\bm{Q}}_{\bm{M}}(z)\leavevmode\nobreak\
{\mathrm{d}}z$ where $\gamma$ is a positively-oriented complex contour
circling around an isolated eigenvalue (assuming it has multiplicity $1$) with
all other eigenvalues left outside and ${\bm{u}}$ is the eigenvector
associated to the corresponding eigenvalue. The almost sure asymptotic value
of this contour integral can be computed with our deterministic equivalent,
thus yielding formulae for the alignments of the _spike_ eigenvectors with
relevant directions of the model.

This briefly summarizes the main techniques we use in our proofs. A complete
presentation of these tools can be found in Couillet and Liao (2022). Other
valuable references for the reader interested in random matrix theory are
Potters and Bouchaud (2020); Bai and Silverstein (2010); Pastur and Shcherbina
(2011); Tao (2012).

##  3 Analysis of Truncated MLSVD under the General Spiked Tensor Model

This section presents a random matrix analysis of the general spiked tensor
model introduced in (2) using tools presented in Section 2.3. We give precise
results on the spectral behavior of the unfoldings of the observed tensor
${\bm{\mathscr{T}}}$, and specify the achievable performance in the estimation
of the underlying signal ${\bm{\mathscr{P}}}$ with a truncated MLSVD.
Although, as explained in Section 1.1, this approach is only _quasi-optimal_ ,
it is very easy to implement and represents an excellent “first guess” to
initialize a numerical scheme converging to a solution to Problem (1), which
is discussed in Section 4.

###  3.1 Random Matrix Results on the Model

Under the general spiked tensor model (2), we consider an
$n_{1}\times\ldots\times n_{d}$ tensor
${\bm{\mathscr{T}}}={\bm{\mathscr{P}}}+\frac{1}{\sqrt{N}}{\bm{\mathscr{N}}}$
of order $d\geqslant 3$, modeling a low-rank deterministic signal
${\bm{\mathscr{P}}}$ corrupted by an additive Gaussian noise tensor
${\bm{\mathscr{N}}}$ whose entries are independent ${\mathcal{N}}(0,1)$ random
variables888We highlight the fact that the _Gaussian_ noise assumption is not
restrictive at all and greatly simplifies the presentation. Firstly, the
universality result of Gurau (2014) shows that, as $N\to+\infty$, any random
tensor with i.i.d. entries has the same limit (in distribution) than that of a
tensor with i.i.d.  _Gaussian_ entries. Moreover, in order to extend our
results to non-Gaussian noise, the “interpolation trick” of Lytova and Pastur
(2009, Corollary 3.1) can be viewed as a generalized Stein’s lemma (Lemma 4)
for non-Gaussian distributions. All Gaussian expectations of our proofs can
then be expressed as their non-Gaussian form plus a residual term controlled
by the moments of the distribution.. We denote by $(r_{1},\ldots,r_{d})$ the
multilinear-rank of ${\bm{\mathscr{P}}}$ and study this model in the
asymptotic regime $N\to+\infty$ with $n_{\ell}=\operatorname*{\Theta}(N)$ and
$r_{\ell}=\operatorname*{\Theta}(1)$, $\ell\in[d]$.

The estimation of ${\bm{\mathscr{P}}}$ with a truncated MLSVD on
${\bm{\mathscr{T}}}$ is simply the computation of the dominant singular
subspaces of ${\bm{\mathscr{T}}}$. Specifically,
$\hat{{\bm{U}}}^{(\ell)}\in{\mathbb{R}}^{n_{\ell}\times r_{\ell}}$ gathers the
$r_{\ell}$ dominant left singular vectors of ${\bm{T}}^{(\ell)}$ — and thus,
$\hat{{\bm{U}}}^{(\ell)\top}\hat{{\bm{U}}}^{(\ell)}={\bm{I}}_{r_{\ell}}$.
Then, a low-multilinear-rank approximation of ${\bm{\mathscr{T}}}$ is
$\hat{{\bm{\mathscr{T}}}}=\llbracket\hat{{\bm{\mathscr{G}}}};\hat{{\bm{U}}}^{(1)},\ldots,\hat{{\bm{U}}}^{(d)}\rrbracket$
with an $r_{1}\times\ldots\times r_{d}$ core tensor
$\hat{{\bm{\mathscr{G}}}}={\bm{\mathscr{T}}}(\hat{{\bm{U}}}^{(1)},\ldots,\hat{{\bm{U}}}^{(d)})$.
An equivalent expression is
$\hat{{\bm{\mathscr{T}}}}=\llbracket{\bm{\mathscr{T}}};\hat{{\bm{U}}}^{(1)}\hat{{\bm{U}}}^{(1)\top},\ldots,\hat{{\bm{U}}}^{(d)}\hat{{\bm{U}}}^{(d)\top}\rrbracket$,
which explicitly shows that $\hat{{\bm{\mathscr{T}}}}$ is the projection of
${\bm{\mathscr{T}}}$ on its dominant singular subspaces999In the more familiar
matrix case ($d=2$), the expression
$\llbracket{\bm{T}};\hat{{\bm{U}}}^{(1)}\hat{{\bm{U}}}^{(1)\top},\hat{{\bm{U}}}^{(2)}\hat{{\bm{U}}}^{(2)\top}\rrbracket$
is equivalent to
$\hat{{\bm{U}}}^{(1)}\hat{{\bm{U}}}^{(1)\top}{\bm{T}}\hat{{\bm{U}}}^{(2)}\hat{{\bm{U}}}^{(2)\top}$..
Thus, the quality of this estimation hinges upon the alignments between the
singular subspaces of ${\bm{\mathscr{P}}}$ and the dominant singular subspaces
of ${\bm{\mathscr{T}}}$. Namely, denoting
$\hat{{\bm{u}}}^{(\ell)}_{q_{\ell}}$, for $q_{\ell}\in[r_{\ell}]$, the columns
of
$\hat{{\bm{U}}}^{(\ell)}=\begin{bmatrix}\hat{{\bm{u}}}^{(\ell)}_{1}&\dots&\hat{{\bm{u}}}^{(\ell)}_{r_{\ell}}\end{bmatrix}$
and given that
${\bm{\mathscr{P}}}=\llbracket{\bm{\mathscr{H}}};{\bm{X}}^{(1)},\ldots,{\bm{X}}^{(d)}\rrbracket$,
the quantities of interest are
$\lVert{\bm{X}}^{(\ell)\top}\hat{{\bm{u}}}^{(\ell)}_{q_{\ell}}\rVert^{2}$
since they represent how much of $\hat{{\bm{u}}}^{(\ell)}_{q_{\ell}}$ is in
the $\ell$-th singular subspace of the signal ${\bm{\mathscr{P}}}$.

In order to understand how the singular subspaces of ${\bm{\mathscr{P}}}$ are
perturbed by the addition of noise, we study the spectral properties of the
unfoldings
${\bm{T}}^{(\ell)}={\bm{P}}^{(\ell)}+\frac{1}{\sqrt{N}}{\bm{N}}^{(\ell)}$. In
fact, since we are only interested in the _left_ singular vectors of
${\bm{T}}^{(\ell)}$, it is more convenient to consider the $n_{\ell}\times
n_{\ell}$ _symmetric_ matrix ${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}$. Note
that this is different from standard spiked matrix models (Benaych-Georges and
Nadakuditi, 2011) because the second dimension of ${\bm{T}}^{(\ell)}$ grows at
a faster polynomial rate than the first one ($\operatorname*{\Theta}(N^{d-1})$
versus $\operatorname*{\Theta}(N)$). Hence, it is easy to see that the
spectrum of ${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}$ should diverge as
$N\to+\infty$: set ${\bm{\mathscr{P}}}=\mathbf{0}$ for simplicity and consider
the expected mean of the eigenvalues,

$\frac{1}{n_{\ell}}\mathbb{E}\left[\sum_{\lambda\in\operatorname{Sp}\left({\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}\right)}\lambda\right]=\frac{1}{n_{\ell}}\mathbb{E}\left[\operatorname{Tr}\left(\frac{1}{N}{\bm{N}}^{(\ell)}{\bm{N}}^{(\ell)\top}\right)\right]=\frac{1}{N}\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}\xrightarrow[N\to+\infty]{}+\infty.$

Hence, we need to consider instead a _rescaled_ version of our random matrix
${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}$ to properly study the behavior of its
spectrum. The quantities $\mu^{(\ell)}_{N}$ and $\sigma_{N}$ introduced in
Theorem 5 below are such that the eigenvalues of
$\frac{1}{\sigma_{N}}\left[{\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}-\mu^{(\ell)}_{N}{\bm{I}}_{n_{\ell}}\right]$
neither diverge nor vanish but stay at a $\operatorname*{\Theta}(1)$ scale as
$N\to+\infty$.

######  Theorem 5 (Deterministic equivalent)

For $\ell\in[d]$, define the following quantities,

$\mu^{(\ell)}_{N}=\frac{1}{N}\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}},\qquad\sigma_{N}=\frac{1}{N}\sqrt{\prod_{\ell\in[d]}n_{\ell}}.$

As $N\to+\infty$, if the ratio
$\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}^{2}/\sigma_{N}$ is bounded, then
the resolvent of the rescaled matrix
$\frac{1}{\sigma_{N}}\left[{\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}-\mu^{(\ell)}_{N}{\bm{I}}_{n_{\ell}}\right]$
has the following deterministic equivalent, defined for all
$\tilde{z}\in{\mathbb{C}}\setminus{\mathbb{R}}$,

$\widetilde{{\bm{Q}}}^{(\ell)}(\tilde{z})\overset{\text{def}}{=}\left(\frac{1}{\sigma_{N}}\left[{\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}-\mu^{(\ell)}_{N}{\bm{I}}_{n_{\ell}}\right]-\tilde{z}{\bm{I}}_{n_{\ell}}\right)^{-1}\\\
\longleftrightarrow\leavevmode\nobreak\
\bar{{\bm{Q}}}^{(\ell)}(\tilde{z})\overset{\text{def}}{=}\left(\frac{1}{\sigma_{N}}{\bm{P}}^{(\ell)}{\bm{P}}^{(\ell)\top}+\frac{1}{\tilde{m}(\tilde{z})}{\bm{I}}_{n_{\ell}}\right)^{-1}$

where
$\displaystyle\tilde{m}(\tilde{z})\overset{\text{def}}{=}\lim_{N\to+\infty}\frac{1}{n_{\ell}}\operatorname{Tr}\widetilde{{\bm{Q}}}^{(\ell)}(\tilde{z})$
does not depend on $\ell\in[d]$ and satisfies the following equation,

$\tilde{m}^{2}(\tilde{z})+\tilde{z}\tilde{m}(\tilde{z})+1=0.$ (6)

Proof  See Appendix A.

This theorem is fundamental. It gives a deterministic equivalent
$\bar{{\bm{Q}}}^{(\ell)}(\tilde{z})$ of the resolvent of the rescaled matrix
$\frac{1}{\sigma_{N}}\left[{\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}-\mu^{(\ell)}_{N}{\bm{I}}_{n_{\ell}}\right]$,
which is our entry point into the precise characterization of its spectral
behavior, following the approach presented in Section 2.3. First of all,
notice that the “scaling parameters” $\mu^{(\ell)}_{N}$ and $\sigma_{N}$ are
respectively $\operatorname*{\Theta}(N^{d-2})$ and
$\operatorname*{\Theta}(N^{\frac{d-2}{2}})$, meaning that the eigenvalues of
${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}$ grow at a speed $N^{d-2}$ and spread
over an interval whose length grows as $\sqrt{N^{d-2}}$ and does not depend on
the mode $\ell$.

Moreover, the relation given in (6) characterizes $\tilde{m}$, the Stieltjes
transform of the limiting spectral distribution (LSD) of
$\frac{1}{\sigma_{N}}\left[{\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}-\mu^{(\ell)}_{N}{\bm{I}}_{n_{\ell}}\right]$.
We see that this LSD is the same regardless of the low-rank perturbation
${\bm{P}}^{(\ell)}{\bm{P}}^{(\ell)\top}$, as it is expected that the
perturbation should only cause the presence of a _finite_ number of isolated
eigenvalues in the empirical spectral distribution. Notice that, if
${\bm{P}}^{(\ell)}=\mathbf{0}$, we simply have
$\bar{{\bm{Q}}}^{(\ell)}(\tilde{z})=\tilde{m}(\tilde{z}){\bm{I}}_{n_{\ell}}$.
We are now one step away from a precise description of the limiting spectral
distribution.

######  Corollary 6 (Limiting spectral distribution)

As $N\to+\infty$, the empirical spectral distribution of the rescaled matrix
$\frac{1}{\sigma_{N}}\left[{\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}-\mu^{(\ell)}_{N}{\bm{I}}_{n_{\ell}}\right]$
weakly converges to $\mu_{\text{SC}}$, the semicircle distribution on
$[-2,+2]$,

${\mathrm{d}}\mu_{\text{SC}}(x)\overset{\text{def}}{=}\frac{\sqrt{\left[4-x^{2}\right]^{+}}}{2\pi}{\mathrm{d}}x.$

Proof  Following (6),
$\tilde{m}(\tilde{z})=\frac{1}{2}\left[-z\pm\sqrt{z^{2}-4}\right]$ where the
$\pm$ sign is chosen so that $\tilde{m}$ satisfies the properties of a
Stieltjes transform, in particular $\Im[\tilde{z}]\Im[\tilde{m}(\tilde{z})]>0$
for all $\tilde{z}\in{\mathbb{C}}\setminus{\mathbb{R}}$. Then, the result
follows from the Stieltjes transform inversion formula introduced in Section
2.3 (see also Couillet and Liao (2022, Section 2.1.2)).

$0$$2$$4$$6$$\cdot 10^{-2}$Density
$\rightarrow$$\ell=1$$\ell=2$$\ell=3$ESDEigenvalueSemicircleSpike$220$$240$$260$$280$0.20.40.60.81.0$\leftarrow$
AlignmentObservedAsymptotic$140$$160$$180$$100$$120$$140$ Figure 3: Top:
empirical spectral distribution (ESD) of
${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}$. The orange curve is the density of
the stretched semicircle on $[\mu^{(\ell)}_{N}\pm 2\sigma_{N}]$ (Corollary 6).
Green dashed lines represent asymptotic positions of spikes
$\mu^{(\ell)}_{N}+\sigma_{N}\tilde{\xi}^{(\ell)}_{q_{\ell}}$ (Theorem 8).
Bottom: Observed alignments between the dominant eigenvectors of
${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}$ and
${\bm{P}}^{(\ell)}{\bm{P}}^{(\ell)\top}$ (purple bars) with their predicted
asymptotic values $[\zeta^{(\ell)}_{q_{\ell}}]^{+}$ (red curve, Theorem 8).
Experimental setting: $d=3$, $(n_{1},n_{2},n_{3})=(300,500,700)$,
$N=n_{1}+n_{2}+n_{3}$, $(r_{1},r_{2},r_{3})=(3,4,5)$ and
$\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}^{2}/\sigma_{N}=15$.

This result states that the limiting spectral distribution of
$\frac{1}{\sigma_{N}}\left[{\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}-\mu^{(\ell)}_{N}{\bm{I}}_{n_{\ell}}\right]$
is, in fact, the very-well-known semicircle distribution, first observed by
Wigner (1955, 1958) in the study of certain special classes of random matrices
arising in quantum mechanics. It indicates that, as $N\to+\infty$, the density
of eigenvalues of ${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}$ is a stretched
semicircle on $[\mu^{(\ell)}_{N}\pm 2\sigma_{N}]$. This phenomenon is
illustrated in the first row of Figure 3, where the empirical spectral
distribution (ESD) of ${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}$ is represented
with the corresponding stretched semicircle for every mode $\ell$ of an
order-$3$ tensor of size $300\times 500\times 700$ following the general
spiked tensor model (2).

######  Remark 7 (From Marčenko-Pastur to Wigner)

Given a random matrix ${\bm{X}}\in{\mathbb{R}}^{p_{1}\times p_{2}}$ with
i.i.d. ${\mathcal{N}}(0,\frac{1}{p_{2}})$ entries, it is well known that the
empirical spectral distribution of ${\bm{X}}{\bm{X}}^{\top}$ converges weakly
to the Marčenko-Pastur distribution as $p_{1},p_{2}\to+\infty$ with
$p_{1}/p_{2}\to c>0$ (Marčenko and Pastur, 1967; Pastur and Shcherbina, 2011,
Corollary 7.2.5; Couillet and Liao, 2022, Theorem 2.4; Potters and Bouchaud,
2020, Chapter 4). On the other hand, the standard semicircle distribution
$\mu_{\text{SC}}$ is known to be the limiting spectral distribution of
symmetric $p\times p$ random matrices with i.i.d. (up to symmetry)
${\mathcal{N}}(0,\frac{1}{p})$ entries (Pastur and Shcherbina, 2011, Corollary
2.2.8; Couillet and Liao, 2022, Theorem 2.5; Potters and Bouchaud, 2020,
Chapter 2). Here, our results show that if $p_{2}$ grows at a faster
polynomial rate than $p_{1}$, the matrix ${\bm{X}}{\bm{X}}^{\top}$ behaves
asymptotically (up to rescaling) like a Wigner matrix, even if its entries are
not independent. Experimentally, we observe that, if $n_{2}$ and $n_{3}$ are
chosen small compared to $n_{1}$ (in contradiction with our assumption
$n_{1},n_{2},n_{3}=\operatorname*{\Theta}(N)$), e.g.,
$(n_{1},n_{2},n_{3})=(1000,40,40)$, then the empirical spectral distribution
of ${\bm{T}}^{(1)}{\bm{T}}^{(1)\top}$ is, after rescaling, better modeled by a
Marčenko-Pastur distribution than by a Wigner semicircle.

The empirical spectral distributions of Figure 3 also show isolated
eigenvalues on the right side of each semicircle. They are caused by the low-
rank perturbation ${\bm{\mathscr{P}}}$ which, in this setting, has
multilinear-rank $(3,4,5)$. The estimate of ${\bm{\mathscr{P}}}$ given by a
truncated MLSVD on ${\bm{\mathscr{T}}}$ has its singular subspaces spanned by
the dominant eigenvectors of ${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}$, i.e.,
precisely those associated with these isolated eigenvalues. Hence, a precise
characterization of the behavior of these _spikes_ is needed to plainly
understand the recovery performance of this estimate. As explained in Section
2.3, this can be achieved with the deterministic equivalent given in Theorem
5.

######  Theorem 8 (Spike behavior)

For $\ell\in[d]$ and $q_{\ell}\in[r_{\ell}]$, define the following
quantities101010We recall that the notation $s_{i}({\bm{A}})$ denotes the
$i$-th singular value of ${\bm{A}}$ in non-increasing order.,

$\rho^{(\ell)}_{q_{\ell}}=\frac{s_{q_{\ell}}^{2}({\bm{P}}^{(\ell)})}{\sigma_{N}},\qquad\tilde{\xi}^{(\ell)}_{q_{\ell}}=\rho^{(\ell)}_{q_{\ell}}+\frac{1}{\rho^{(\ell)}_{q_{\ell}}}\qquad\text{and}\qquad\zeta^{(\ell)}_{q_{\ell}}=1-\frac{1}{\left[\rho^{(\ell)}_{q_{\ell}}\right]^{2}}.$

As $N\to+\infty$, if the ratio
$\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}^{2}/\sigma_{N}$ is bounded and
$\rho^{(\ell)}_{q_{\ell}}>1$, then

$s_{q_{\ell}}^{2}({\bm{T}}^{(\ell)})=\mu^{(\ell)}_{N}+\sigma_{N}\tilde{\xi}^{(\ell)}_{q_{\ell}}+\operatorname*{o}(1)\qquad\text{and}\qquad\left\lVert{\bm{X}}^{(\ell)\top}\hat{{\bm{u}}}^{(\ell)}_{q_{\ell}}\right\rVert^{2}=\zeta^{(\ell)}_{q_{\ell}}+\operatorname*{o}(1)$

where $\hat{{\bm{u}}}^{(\ell)}_{q_{\ell}}$ is the $q_{\ell}$-th dominant left
singular vector of ${\bm{T}}^{(\ell)}$ and both $\operatorname*{o}(1)$ are
quantities vanishing almost surely as $N\to+\infty$.

Proof  See Appendix B.

The first quantity defined in this theorem, $\rho^{(\ell)}_{q_{\ell}}$, should
be understood as a signal-to-noise ratio (SNR). Indeed, the squared
$q_{\ell}$-th singular value of ${\bm{P}}^{(\ell)}$ (i.e., the $q_{\ell}$-th
eigenvalue of ${\bm{P}}^{(\ell)}{\bm{P}}^{(\ell)\top}$),
$s_{q_{\ell}}^{2}({\bm{P}}^{(\ell)})$, measures the “strength” of the signal
in its $q_{\ell}$-th principal direction, whereas $\sigma_{N}$ measures the
spread of the noise, as seen in Theorem 5. The two quantities
$\tilde{\xi}^{(\ell)}_{q_{\ell}}$ and $\zeta^{(\ell)}_{q_{\ell}}$ depend only
on the value of this SNR and indicate respectively the position of an isolated
eigenvalue in the sprectrum of ${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}$ and
the alignment of the corresponding eigenvector with the sought signal. In
fact, we observe a _phase transition_ phenomenon: if the SNR is large enough,
i.e., if $\rho^{(\ell)}_{q_{\ell}}>1$, an eigenvalue of
${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}$ isolates itself from the
semicircle111111Indeed, note that
$\rho^{(\ell)}_{q_{\ell}}>1\implies\tilde{\xi}^{(\ell)}_{q_{\ell}}>2$. around
$\mu^{(\ell)}_{N}+\sigma_{N}\tilde{\xi}^{(\ell)}_{q_{\ell}}$. Moreover,
recalling that
${\bm{\mathscr{P}}}=\llbracket{\bm{\mathscr{H}}};{\bm{X}}^{(1)},\ldots,{\bm{X}}^{(d)}\rrbracket$,
the eigenvector associated with this isolated eigenvalue is aligned with the
subspace spanned by ${\bm{X}}^{(\ell)}$, which is the $\ell$-th singular
subspace of ${\bm{\mathscr{P}}}$. The quality of this alignment is given by
$0<\zeta^{(\ell)}_{q_{\ell}}\leqslant 1$.

Most importantly, this result reveals the _non-trivial regime_ for the
estimation of ${\bm{\mathscr{P}}}$ with a truncated MLSVD. Since
$\sigma_{N}=\operatorname*{\Theta}(N^{\frac{d-2}{2}})$, it shows that
$\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}^{2}=\sum_{q_{\ell}=1}^{r_{\ell}}s_{q_{\ell}}^{2}({\bm{P}}^{(\ell)})$
must also be of the same order. Indeed, if
$\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}^{2}\ll N^{\frac{d-2}{2}}$, then
$\rho^{(\ell)}_{q_{\ell}}\to 0$, the SNR is too small and no signal can be
recovered, whereas if $\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}^{2}\gg
N^{\frac{d-2}{2}}$, then $\rho^{(\ell)}_{q_{\ell}}\to+\infty$, the SNR is very
high and recovery of ${\bm{\mathscr{P}}}$ is easy. It is precisely between
these two regimes, i.e.,
$\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}=\operatorname*{\Theta}(N^{\frac{d-2}{4}})$,
that the recovery is non-trivial. Note that this observation is in line with
the results of Ben Arous et al. (2021) and Zhang and Xia (2020). In this non-
trivial regime, the quantities $\zeta^{(\ell)}_{q_{\ell}}$ given in Theorem 8
precisely quantify how well the dominant eigenvectors of
${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}$ are aligned with the sought signal,
i.e., the singular subspaces of ${\bm{\mathscr{P}}}$. In section 3.2 below,
this result is used to study the reconstruction performance of truncated
MLSVD.

######  Remark 9 (Link with the spiked Wigner model)

The reader familiar with the spiked Wigner model may have recognized the
expressions of $\tilde{\xi}^{(\ell)}_{q_{\ell}}$ and
$\zeta^{(\ell)}_{q_{\ell}}$ given in Theorem 8. Indeed, given a symmetric
$p\times p$ random matrix ${\bm{W}}$ with i.i.d. (up to symmetry)
${\mathcal{N}}(0,\frac{1}{p})$ entries, the spectrum of
$[\rho{\bm{x}}{\bm{x}}^{\top}+{\bm{W}}]$ with $\lVert{\bm{x}}\rVert=1$ follows
a semicircle distribution as $p\to+\infty$ with an isolated eigenvalue at
$\rho+\frac{1}{\rho}$ if, and only if, $\rho>1$ (Féral and Péché, 2007;
Edwards and Jones, 1976; Füredi and Komlós, 1981). Moreover, the corresponding
eigenvector ${\bm{u}}$ is such that
$\lvert{\bm{x}}^{\top}{\bm{u}}\rvert^{2}\to 1-\frac{1}{\rho^{2}}$ almost
surely as $p\to+\infty$ (Benaych-Georges and Nadakuditi, 2011). As discussed
in Remark 7, up to rescaling, ${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}$
asymptotically behaves like a spiked Wigner matrix.

Theorem 8 is illustrated in Figure 3. In the first row, asymptotic positions
of isolated eigenvalues
$\mu^{(\ell)}_{N}+\sigma_{N}\tilde{\xi}^{(\ell)}_{q_{\ell}}$ are represented
by the green dashed lines. In our experiment, ${\bm{\mathscr{P}}}$ has
multilinear-rank $(3,4,5)$. Hence $3$, $4$ and $5$ isolated eigenvalues are
expected in the spectrum of ${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}$ for
$\ell=1$, $2$ and $3$ respectively. This is indeed the case for $\ell=1$ and
$\ell=2$ but not $\ell=3$ where there are only $4$ spike eigenvalues. In fact,
$s_{5}({\bm{P}}^{(3)})$ is not “energetic enough” to extricate itself from the
bulk of eigenvalues, i.e, the SNR $\rho^{(3)}_{5}$ is below the phase
transition threshold. Hence the $5$-th dominant left singular vector of
${\bm{T}}^{(\ell)}$ is not informative as it is not aligned with the $3$-rd
singular subspace of ${\bm{\mathscr{P}}}$, spanned by ${\bm{X}}^{(3)}$. The
second row of Figure 3 depicts the alignments of the spiked eigenvectors
$\hat{{\bm{u}}}^{(\ell)}_{q_{\ell}}$ with the corresponding singular subspaces
of ${\bm{\mathscr{P}}}$ as well as the asymptotic alignment given by Theorem 8
as a function of the position of the associated eigenvalue. It appears that
the higher is the SNR $\rho^{(\ell)}_{q_{\ell}}$, the farther is the isolated
eigenvalue from the bulk and the more is the corresponding eigenvector aligned
with the span of ${\bm{X}}^{(\ell)}$. This assertion can be intuitively
understood in terms of “energy” $s_{q_{\ell}}^{2}({\bm{P}}^{(\ell)})$ of a
given principal direction. More energy pushes the eigenvalue farther from the
bulk and aligns the corresponding eigenvector with the corresponding singular
subspace of ${\bm{\mathscr{P}}}$.

###  3.2 Reconstruction Performance of Truncated MLSVD

Our random matrix results allow to accurately study the reconstruction
performance of truncated MLSVD. Given a data tensor ${\bm{\mathscr{T}}}$
following the general spiked tensor model (2), we consider its low-rank
approximation
$\hat{{\bm{\mathscr{T}}}}=\llbracket\hat{{\bm{\mathscr{G}}}};\hat{{\bm{U}}}^{(1)},\ldots,\hat{{\bm{U}}}^{(d)}\rrbracket$
where $\hat{{\bm{U}}}^{(\ell)}$ is the $n_{\ell}\times r_{\ell}$ matrix whose
columns are the $r_{\ell}$ dominant singular vectors of ${\bm{\mathscr{T}}}$
and
$\hat{{\bm{\mathscr{G}}}}={\bm{\mathscr{T}}}(\hat{{\bm{U}}}^{(1)},\ldots,\hat{{\bm{U}}}^{(d)})$.
$\hat{{\bm{\mathscr{T}}}}$ is the projection of ${\bm{\mathscr{T}}}$ on its
dominant singular subspaces, hence the name _truncated MLSVD_ as it
generalizes the truncated SVD of matrices. Given model (2), the underlying
signal estimated by $\hat{{\bm{\mathscr{T}}}}$ is
${\bm{\mathscr{P}}}=\llbracket{\bm{\mathscr{H}}};{\bm{X}}^{(1)},\ldots,{\bm{X}}^{(d)}\rrbracket$.
The reconstruction performance of $\hat{{\bm{\mathscr{T}}}}$ hence depends on
how well the _subspace_ 121212The object of importance is indeed the subspace
and not the matrix $\hat{{\bm{U}}}^{(\ell)}$ since any other matrix
$\hat{{\bm{U}}}^{(\ell)}{\bm{O}}^{(\ell)}$, with ${\bm{O}}^{(\ell)}$ an
$r_{\ell}\times r_{\ell}$ orthogonal matrix, would span the same subspace and
therefore give the same approximation. spanned by $\hat{{\bm{U}}}^{(\ell)}$
estimates the one spanned by ${\bm{X}}^{(\ell)}$.

Metrics between singular subspaces are often expressed in terms of _principal
angles_ (Björck and Golub, 1973; Stewart and Sun, 1990, II.4), which
generalize the concept of angle between lines. Given two subspaces (here,
$\operatorname{Span}{\bm{X}}^{(\ell)}$ and
$\operatorname{Span}\hat{{\bm{U}}}^{(\ell)}$), one can define a set of mutual
angles which are invariant under isometric transformation.

######  Definition 10 (Principal angles)

The principal angles $\theta^{(\ell)}_{q_{\ell}}\in[0,\frac{\pi}{2}]$ between
the subspaces spanned by ${\bm{X}}^{(\ell)}$ and $\hat{{\bm{U}}}^{(\ell)}$ are
recursively defined for $q_{\ell}=1,\ldots,r_{\ell}$ by

$\cos\theta^{(\ell)}_{q_{\ell}}={\bm{x}}_{q_{\ell}}^{\top}{\bm{u}}_{q_{\ell}}\quad\text{with}\quad({\bm{x}}_{q_{\ell}},{\bm{u}}_{q_{\ell}})\in\operatorname*{\arg\max}_{\begin{subarray}{c}({\bm{x}},{\bm{u}})\in\operatorname{Span}{\bm{X}}^{(\ell)}\times\operatorname{Span}\hat{{\bm{U}}}^{(\ell)}\\\
{\bm{x}}^{\top}{\bm{x}}_{q_{\ell}^{\prime}}=0,\leavevmode\nobreak\
{\bm{u}}^{\top}{\bm{u}}_{q_{\ell}^{\prime}}=0,\leavevmode\nobreak\ 1\leqslant
q_{\ell}^{\prime}<q_{\ell}\end{subarray}}{\bm{x}}^{\top}{\bm{u}}.$

Moreover, we have the following useful property.

######  Proposition 11 (Björck and Golub, 1973)

The $q_{\ell}$-th singular value of
${\bm{X}}^{(\ell)\top}\hat{{\bm{U}}}^{(\ell)}$ in non-increasing order equals
the cosine of the $q_{\ell}$-th principal angle,

$s_{q_{\ell}}({\bm{X}}^{(\ell)\top}\hat{{\bm{U}}}^{(\ell)})=\cos\theta^{(\ell)}_{q_{\ell}},\qquad\ell\in[d],\quad
q_{\ell}\in[r_{\ell}].$

In fact, information about the alignment between the subspaces induced by
${\bm{X}}^{(\ell)}$ and $\hat{{\bm{U}}}^{(\ell)}$ are contained entirely in
the $r_{\ell}\times r_{\ell}$ matrix
${\bm{X}}^{(\ell)\top}\hat{{\bm{U}}}^{(\ell)}$. Following Definition 10, we
know from Theorem 8 that, as $N\to+\infty$,
$\cos^{2}\theta^{(\ell)}_{q_{\ell}}\to[\zeta^{(\ell)}_{q_{\ell}}]^{+}$ almost
surely131313Indeed, since
$\zeta^{(\ell)}_{1}\geqslant\ldots\geqslant\zeta^{(\ell)}_{r_{\ell}}$, observe
that, in Definition 10,
$\lvert{\bm{x}}_{q_{\ell}}^{\top}{\bm{u}}_{q_{\ell}}\rvert^{2}$ is
asymptotically bounded by $[\zeta^{(\ell)}_{q_{\ell}}]^{+}$., where we have
used the handy notation $[\cdot]^{+}=\max(\cdot,0)$ since
$\rho^{(\ell)}_{q_{\ell}}>1\iff\zeta^{(\ell)}_{q_{\ell}}>0$. Hence, using
Proposition 11,

$\frac{1}{r_{\ell}}\left\lVert{\bm{X}}^{(\ell)\top}\hat{{\bm{U}}}^{(\ell)}\right\rVert_{\mathrm{F}}^{2}=\frac{1}{r_{\ell}}\sum_{q_{\ell}=1}^{r_{\ell}}\cos^{2}\theta^{(\ell)}_{q_{\ell}}\xrightarrow[N\to+\infty]{\text{a.s.}}\frac{1}{r_{\ell}}\sum_{q_{\ell}=1}^{r_{\ell}}\left[\zeta^{(\ell)}_{q_{\ell}}\right]^{+}.$
(7)

Therefore, the quantity
$\frac{1}{r_{\ell}}\lVert{\bm{X}}^{(\ell)\top}\hat{{\bm{U}}}^{(\ell)}\rVert_{\mathrm{F}}^{2}\in[0,1]$
appears as a relevant measure of alignment between the singular subspaces of
${\bm{\mathscr{P}}}$ and $\hat{{\bm{\mathscr{T}}}}$ and does not depend on the
chosen orthonormal bases ${\bm{X}}^{(\ell)}$ and $\hat{{\bm{U}}}^{(\ell)}$.
More details on metrics between subspaces can be found in Stewart and Sun
(1990, II.4).

In Figure 2 are represented the alignments between the singular subspaces of
$\hat{{\bm{\mathscr{T}}}}$ and
${\bm{\mathscr{P}}}=\sqrt{\omega}{\bm{\mathscr{P}}}_{\circ}$ with
$\lVert{\bm{\mathscr{P}}}_{\circ}\rVert^{2}=\sigma_{N}$ as a function of the
signal-to-noise ratio $\omega$. The fact that
$\lVert{\bm{\mathscr{P}}}_{\circ}\rVert^{2}=\sigma_{N}$ ensures that the
estimation problem is _non-trivial_ (neither too easy nor too hard) as
$\rho^{(\ell)}_{q_{\ell}}=\operatorname*{\Theta}(1)$. Plain curves are the
alignments given by Theorem 8 as $N\to+\infty$ (right-hand side of equation
(7)) whereas dotted curves are simulation results at finite $N$ (left-hand
side of equation (7)). In this setting, ${\bm{\mathscr{P}}}_{\circ}$ has
multilinear-rank $(3,4,5)$. Hence, its “energy”
$\lVert{\bm{\mathscr{P}}}_{\circ}\rVert^{2}$ is spread among $3$, $4$ and $5$
directions along modes $1$, $2$ and $3$ respectively. Each break in the plain
curves correspond to a value of $\omega$ such that
$\rho^{(\ell)}_{q_{\ell}}=1$, that is,
$\omega=\sigma_{N}/s_{q_{\ell}}^{2}({\bm{P}}_{\circ}^{(\ell)})$. In other
words, there are $r_{\ell}$ phase transitions along mode $\ell$ whose
positions depend on the singular values of ${\bm{P}}_{\circ}^{(\ell)}$. If
$\omega$ is too small (here, $\omega\lesssim 2$), truncated MLSVD fails to
recover any direction of the singular subspaces of
${\bm{\mathscr{P}}}_{\circ}$. As $\omega$ passes the first phase transition
(here, at $\omega=\sigma_{N}/s_{1}^{2}({\bm{P}}_{\circ}^{(1)})\approx 2$), a
first principal direction is partially reconstructed. Then, more and more
phase transitions occur, corresponding to more and more principal directions
being recovered as $\omega$ grows. Simultaneously, the reconstruction of
previous directions keeps improving. Eventually, as $\omega$ is large,
subspace alignments approach $1$, indicating that truncated MLSVD accurately
recovers the singular subspaces of ${\bm{\mathscr{P}}}_{\circ}$.

The reader has not missed the dashed lines in Figure 2 showing much better
subspace alignments than truncated MLSVD. They result from the numerical
estimation of the _best_ rank-$(3,4,5)$ approximation of ${\bm{\mathscr{T}}}$
with the HOOI algorithm, which is discussed below.

##  4 Numerical Estimation of the Best Low-Multilinear-Rank Approximation

In search of an efficient estimator of the planted signal
${\bm{\mathscr{P}}}$, one naturally considers the best low-multilinear-rank
approximation of ${\bm{\mathscr{T}}}$, that is, a solution to Problem (1). As
explained in the introduction, this is NP-hard in general but numerical
schemes can compute it in polynomial time above the computational threshold
(Montanari and Richard, 2014; Zhang and Xia, 2020). In this section, we
examine the most standard of these numerical schemes, namely the Higher Order
Orthogonal Iteration (HOOI) algorithm (De Lathauwer et al., 2000a; Kroonenberg
and de Leeuw, 1980; Kapteyn et al., 1986), and discuss the numerical
difficulties faced in the computation of a solution to Problem (1).

###  4.1 Higher-Order Orthogonal Iteration

Following De Lathauwer et al. (2000a, Theorem 4.2), the maximum likelihood
estimation formulated in Problem (1) is equivalent to

$\left({\bm{U}}^{(1)}_{\star},\ldots,{\bm{U}}^{(d)}_{\star}\right)\in\operatorname*{\arg\max}_{{\bm{U}}^{(\ell)}\in
V_{r_{\ell}}({\mathbb{R}}^{n_{\ell}}),\leavevmode\nobreak\
\ell\in[d]}\frac{1}{2}\left\lVert{\bm{\mathscr{T}}}({\bm{U}}^{(1)},\ldots,{\bm{U}}^{(d)})\right\rVert_{\mathrm{F}}^{2}$
(8)

where
$V_{r_{\ell}}({\mathbb{R}}^{n_{\ell}})=\\{{\bm{U}}^{(\ell)}\in{\mathbb{R}}^{n_{\ell}\times
r_{\ell}}\mid{\bm{U}}^{(\ell)\top}{\bm{U}}^{(\ell)}={\bm{I}}_{r_{\ell}}\\}$ is
the set of $r_{\ell}\times n_{\ell}$ matrices with orthonormal columns, known
as the _Stiefel manifold_ (Chikuse, 2003; Absil et al., 2009). Then, since the
Frobenius norm of ${\bm{\mathscr{T}}}({\bm{U}}^{(1)},\ldots,{\bm{U}}^{(d)})$
is equal to the Frobenius norm of any of its unfoldings,

$\left\lVert{\bm{\mathscr{T}}}({\bm{U}}^{(1)},\ldots,{\bm{U}}^{(d)})\right\rVert_{\mathrm{F}}=\left\lVert{\bm{U}}^{(\ell)\top}{\bm{T}}^{(\ell)}\operatorname*{\scalerel*{\boxtimes}{\bigotimes}}_{\ell^{\prime}\neq\ell}{\bm{U}}^{(\ell^{\prime})}\right\rVert_{\mathrm{F}},\qquad\ell\in[d].$

And we see from Problem (8) that ${\bm{U}}^{(\ell)}_{\star}$ is the matrix
gathering the $r_{\ell}$ dominant left singular vectors of
${\bm{T}}^{(\ell)}\operatorname*{\scalerel*{\boxtimes}{\bigotimes}}_{\ell^{\prime}\neq\ell}{\bm{U}}^{(\ell^{\prime})}_{\star}$.
This is precisely what motivates the HOOI algorithm presented in Algorithm 1.
It performs fixed-point iterations to compute a solution
${\bm{U}}^{(1)}_{\star},\ldots,{\bm{U}}^{(d)}_{\star}$ satisfying the previous
property141414In fact, this property corresponds to first-order optimality
conditions of Problem (8), with the squared singular values as Lagrange
multipliers..

for  _$\ell=1,\ldots,d$_ do ${\bm{U}}^{(\ell)}_{0}\leftarrow$ $r_{\ell}$
dominant left singular vectors of ${\bm{T}}^{(\ell)}$

repeat

for  _$\ell=1,\ldots,d$_ do ${\bm{U}}^{(\ell)}_{t+1}\leftarrow$ $r_{\ell}$
dominant left singular vectors of
${\bm{T}}^{(\ell)}\operatorname*{\scalerel*{\boxtimes}{\bigotimes}}_{\ell^{\prime}\neq\ell}{\bm{U}}^{(\ell^{\prime})}_{t}$

until  _convergence at $t=T$_

${\bm{\mathscr{G}}}_{\text{HOOI}}\leftarrow{\bm{\mathscr{T}}}({\bm{U}}^{(1)}_{T},\ldots,{\bm{U}}^{(d)}_{T})$

Algorithm 1 Higher-Order Orthogonal Iteration (De Lathauwer et al., 2000a)

The HOOI algorithm is initialized with
${\bm{U}}^{(1)}_{0},\ldots,{\bm{U}}^{(d)}_{0}$, the truncated
MLSVD151515Consistently with the notations of Section 3, this is
$\hat{{\bm{U}}}^{(1)},\ldots,\hat{{\bm{U}}}^{(d)}$. of ${\bm{\mathscr{T}}}$.
Given the results of Section 3, this is indeed a very good and easily
computable first guess. Then, fixed-point iterations are repeated in order to
find a solution ${\bm{U}}^{(1)}_{\star},\ldots,{\bm{U}}^{(d)}_{\star}$ such
that ${\bm{U}}^{(\ell)}_{\star}$ spans the left $r_{\ell}$-dimensional
dominant singular subspace of
${\bm{T}}^{(\ell)}\operatorname*{\scalerel*{\boxtimes}{\bigotimes}}_{\ell^{\prime}\neq\ell}{\bm{U}}^{(\ell^{\prime})}_{\star}$
for all $\ell\in[d]$, which corresponds to a solution to Problem (8). In
practice, the stopping criterion can be chosen as a negligible change in the
norm of the estimated core tensor,
$\lVert{\bm{\mathscr{T}}}({\bm{U}}^{(1)}_{t},\ldots,{\bm{U}}^{(d)}_{t})\rVert_{\mathrm{F}}$.

Xu (2018) showed that the convergence of this algorithm towards a local
minimum of Problem (8) is guaranteed as long as its initialization is
_sufficiently close_ to this local minimum. In light of our previous results,
we can provide further insight into this “sufficiently close” property.
Indeed, in Theorem 13 below, we show that an initialization with non-zero
alignment with the signal ${\bm{\mathscr{P}}}$ is sufficient to ensure that
the HOOI algorithm perfectly reconstructs it _after a single iteration_.

Before introducing Theorem 13, we formulate an important preliminary result
which essentially states that, given ${\bm{A}}^{(\ell)}\in
V_{r_{\ell}}({\mathbb{R}}^{n_{\ell}})$, $\ell\in[d]$, the quantity
$\lVert{\bm{\mathscr{N}}}({\bm{A}}^{(1)},\ldots,{\bm{A}}^{(d)})\rVert_{\mathrm{F}}/\sqrt{N}$
is almost surely bounded as $N\to+\infty$.

###### Lemma 12

With probability at least $1-\delta$,

$\sup_{{\bm{A}}^{(\ell)}\in
V_{r_{\ell}}({\mathbb{R}}^{n_{\ell}}),\leavevmode\nobreak\
\ell\in[d]}\left\lVert{\bm{\mathscr{N}}}({\bm{A}}^{(1)},\ldots,{\bm{A}}^{(d)})\right\rVert_{\mathrm{F}}^{2}\\\
\leqslant
16\left[\left(\sum_{\ell=1}^{d}r_{\ell}\left(n_{\ell}-\frac{r_{\ell}+1}{2}\right)\right)\log\frac{Cd}{\log\frac{3}{2}}+\log\left(\frac{1}{\delta}\max\left(1,e^{\frac{1}{2}\prod_{\ell=1}^{d}r_{\ell}-1}\right)\right)\right]$

where $C>0$ is a universal constant.

Proof  See Appendix C.
This result is crucial to handle the behavior of the noise in our analysis of
Algorithm 1 (Appendix D), which leads to the following result on the alignment
between the singular subspaces of the signal ${\bm{\mathscr{P}}}$ (spanned by
${\bm{X}}^{(\ell)}$) and those estimated from the observation
${\bm{\mathscr{T}}}$ after the first iteration of HOOI (spanned by
${\bm{U}}^{(\ell)}_{1}$).

######  Theorem 13 (Asymptotic optimality of HOOI)

As $N\to+\infty$, if $\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}\gg 1$ and

$\min_{\ell\in[d],\leavevmode\nobreak\
q_{\ell}\in[r_{\ell}]}\left\lVert{\bm{\mathscr{P}}}({\bm{U}}^{(1)}_{0},\ldots,{\bm{x}}^{(\ell)}_{q_{\ell}},\ldots,{\bm{U}}^{(d)}_{0})\right\rVert_{\mathrm{F}}\overset{\text{def}}{=}L_{N}\gg\left\lVert{\bm{\mathscr{P}}}\right\rVert_{\mathrm{F}}^{1/2}$

where ${\bm{x}}^{(\ell)}_{q_{\ell}}$ is the $q_{\ell}$-th column of
${\bm{X}}^{(\ell)}$, then,

$\frac{1}{r_{\ell}}\left\lVert{\bm{X}}^{(\ell)\top}{\bm{U}}^{(\ell)}_{1}\right\rVert_{\mathrm{F}}^{2}=1+\operatorname*{\mathcal{O}}\left(\frac{\left\lVert{\bm{\mathscr{P}}}\right\rVert_{\mathrm{F}}}{L_{N}^{2}}\right)\qquad\text{almost
surely.}$

Proof  See Appendix D.

It is important to carefully understand the assumptions of Theorem 13.
Firstly, it assumes that $\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}\gg 1$ as
$N\to+\infty$, that is, the signal is not necessarily in the non-trivial
regime $\operatorname*{\Theta}(N^{\frac{d-2}{4}})$ but can be smaller or
bigger as long as $\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}\to+\infty$,
regardless its speed. Then, the second assumption
$\lVert{\bm{\mathscr{P}}}({\bm{U}}^{(1)}_{0},\ldots,{\bm{x}}^{(\ell)}_{q_{\ell}},\ldots,{\bm{U}}^{(d)}_{0})\rVert_{\mathrm{F}}\gg\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}^{1/2}$
means that each principal directions of the $\ell$-th singular subspace are
sufficiently preserved after contraction on
$\\{{\bm{U}}^{(\ell^{\prime})}_{0}\\}_{\ell^{\prime}\neq\ell}$. When these
assumptions are verified, Theorem 13 states that the matrices
${\bm{U}}^{(1)}_{1},\ldots,{\bm{U}}^{(d)}_{1}$ computed after the _first_
iteration of HOOI _perfectly_ reconstruct the singular subspaces of the sought
signal ${\bm{\mathscr{P}}}$ as the dimensions of the tensor,
$n_{1},\ldots,n_{d}$, grow large. More formally, as $N\to+\infty$, the
alignment
$\frac{1}{r_{\ell}}\lVert{\bm{X}}^{(\ell)\top}{\bm{U}}^{(\ell)}_{1}\rVert_{\mathrm{F}}^{2}$
approaches $1$ almost surely. Furthermore, the speed of this convergence
behaves like $\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}/L_{N}^{2}$.

Theorem 13 does not assume a particular choice of initialization
${\bm{U}}^{(1)}_{0},\ldots,{\bm{U}}^{(d)}_{0}$ and gives a sufficient
condition for it to ensure the convergence of the algorithm. Nevertheless, as
it is presented in Algorithm 1, truncated MLSVD is a standard choice of
initialization. In this case, the assumption
$L_{N}=\operatorname*{\Theta}(\lVert{\bm{\mathscr{P}}}\rVert)$ is verified as
soon as enough principal directions are recovered. According to Theorem 8,
this is only possible if
$\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}\geqslant\sqrt{\sigma_{N}}=\operatorname*{\Theta}(N^{\frac{d-2}{4}})$
since a _necessary_ condition is $\rho^{(\ell)}_{1}>1$ for all $\ell\in[d]$,
while a _sufficient_ condition is $\rho^{(\ell)}_{q_{\ell}}>1$ for all
$\ell\in[d]$ and $q_{\ell}\in[d]$. In other words, convergence at speed
$\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}^{-1}$ as assured above a critical
signal-to-noise ratio lying between the first and the last phase transition of
each mode, and which depends on the particular structure of the core tensor
${\bm{\mathscr{H}}}$. Yet, in most cases, this happens quite early, right
after the first phase transitions, see for example Figure 2.

We emphasize the fact that the assumption of Theorem 13 can already be
verified as soon as $\min_{\ell}\rho^{(\ell)}_{1}>1$. That is, there is no
need for all the principal directions to be reconstructed at initialization.
In fact, it could very well be that $\max_{\ell}\rho^{(\ell)}_{2}<1$. If the
singular subspaces of ${\bm{\mathscr{P}}}$ are sufficiently preserved with the
initialization ${\bm{U}}^{(1)}_{0},\ldots,{\bm{U}}^{(d)}_{0}$ — i.e., if
$L_{N}\gg\lVert{\bm{\mathscr{P}}}\rVert^{1/2}$ —, then the other principal
directions still emerge after the first iteration.

These technical considerations vanish in the simpler rank-$1$ case: if
${\bm{\mathscr{P}}}=\beta_{N}\bigotimes_{\ell\in[d]}{\bm{x}}^{(\ell)}$ then a
_necessary and sufficient_ condition for
$L_{N}=\operatorname*{\Theta}(\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}})$ is
simply $\beta_{N}^{2}>\sigma_{N}$. Indeed,
$\rho^{(\ell)}_{1}=\beta_{N}^{2}/\sigma_{N}$ for all $\ell\in[d]$. Therefore,
when this assumption if verified, Theorem 13 ensures the asymptotic exact
reconstruction of ${\bm{x}}^{(1)},\ldots,{\bm{x}}^{(d)}$ in a single iteration
with a $\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}^{-1}=\beta_{N}^{-1}$ speed
of convergence.

######  Remark 14 (Practical implications)

In practice, one should still run several iterations of Algorithm 1 until a
certain stopping criterion is verified as this effectively improves the final
estimate and converges to a solution to the maximum likelihood estimation (1)
(Xu, 2018). Theorem 13 states that the reconstruction performance of HOOI
after the first iteration increases as we consider larger tensors, until it
reaches perfect recovery in the large $N$ limit. In other words, the number of
iterations required to achieve a specific level of accuracy in maximum
likelihood estimation tends to $1$ as $N\to+\infty$.

$0$$500$$1000$$1500$$0$$0.2$$0.4$$0.6$$0.8$$1$$N$InitializationMode $\ell$
$\rightarrow$$1$$2$$3$ObservedPredicted$0$$500$$1000$$1500$$0$$0.2$$0.4$$0.6$$0.8$$1$$N$First
iteration$0$$500$$1000$$1500$$1\cdot 10^{-1}$$0.2$$0.3$$0.4$$N$First iteration
(rescaled) Figure 4: Alignments between singular subspaces of the observation
${\bm{\mathscr{T}}}={\bm{\mathscr{P}}}+\frac{1}{\sqrt{N}}{\bm{\mathscr{N}}}$
and of the signal ${\bm{\mathscr{P}}}$, with
$\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}^{2}/\sigma_{N}=10$, at
initialization of Algorithm 1 (i.e., truncated MLSVD) and after the first
iteration, as a function of the size of the tensor given by the parameter $N$.
Left:
$\frac{1}{r_{\ell}}\lVert{\bm{X}}^{(\ell)\top}{\bm{U}}^{(\ell)}_{0}\rVert_{\mathrm{F}}^{2}$.
Middle:
$\frac{1}{r_{\ell}}\lVert{\bm{X}}^{(\ell)\top}{\bm{U}}^{(\ell)}_{1}\rVert_{\mathrm{F}}^{2}$.
Right:
$(1-\frac{1}{r_{\ell}}\lVert{\bm{X}}^{(\ell)\top}{\bm{U}}^{(\ell)}_{1}\rVert_{\mathrm{F}}^{2})\times\sqrt{\sigma_{N}}$.
Experimental setting: $d=3$,
$(\frac{n_{1}}{N},\frac{n_{2}}{N},\frac{n_{3}}{N})=(\frac{1}{6},\frac{2}{6},\frac{3}{6})$,
$N=n_{1}+n_{2}+n_{3}$ and $(r_{1},r_{2},r_{3})=(3,4,5)$.

Theorem 13 is illustrated in Figure 4. As a function of $N$ — the size of the
tensor — we represent the subspace alignments observed at initialization and
after the first iteration for a fixed signal-to-noise ratio
$\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}^{2}/\sigma_{N}=10$. The left
panel compares the observed alignments achieved with truncated MLSVD
(initialization of Algorithm 1) with the asymptotic alignments predicted by
Theorem 8. As $N$ grows, the observed alignments remain around their
asymptotic values, with only a decrease in variance. The middle panel presents
the alignments after the first iteration. Here, as $N$ increases, we observe
an increase in the values of the alignments, which approach $1$, consistently
with Theorem 13. This is specified in the right panel where the value
$(1-\frac{1}{r_{\ell}}\lVert{\bm{X}}^{(\ell)\top}{\bm{U}}^{(\ell)}_{1}\rVert_{\mathrm{F}}^{2})\times\sqrt{\sigma_{N}}$
is plotted. According to Theorem 13, this value should be
$\operatorname*{\mathcal{O}}(1)$ since
$L_{N}=\operatorname*{\Theta}(\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}})=\operatorname*{\Theta}(\sqrt{\sigma_{N}})$
here. The observed behavior confirms the
$\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}^{-1}$ speed of convergence
asserted in Theorem 13.

###  4.2 Discussion on Signal Reconstructibility

Our results on truncated MLSVD (Section 3) and HOOI (Section 4.1) bring
insight into the computational-to-statistical gap observed in the low-
multilinear-rank approximation problem. Truncated MLSVD can only work
efficiently if $\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}$ is at least
$\operatorname*{\Theta}(N^{\frac{d-2}{4}})$ and its reconstruction performance
has been discussed in Section 3.2. However, Theorem 13 suggests that it is
possible to perfectly reconstruct the signal ${\bm{\mathscr{P}}}$ from the
observation ${\bm{\mathscr{T}}}$ as long as
$\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}\gg 1$ and HOOI is accurately
initialized. Yet, it is known that, without prior information on
${\bm{\mathscr{P}}}$, maximum likelihood estimation (1) is NP-hard below the
$\operatorname*{\Theta}(N^{\frac{d-2}{4}})$ computational threshold (Zhang and
Xia, 2020), which lies precisely in the non-trivial regime of truncated MLSVD.

In fact, what can be understood from Theorem 13 is that it suffices to have an
initialization ${\bm{U}}^{(\ell)}_{0}$, $\ell\in[d]$, slightly aligned with
the underlying signal ${\bm{\mathscr{P}}}$ to be in the right basin of
attraction and allow the convergence of Algorithm 1 towards a solution to
Problem (1). This complements the results of Xu (2018) on the conditions of
convergence of HOOI. Furthermore, a solution to Problem (1) is aligned with
${\bm{\mathscr{P}}}$ as soon as
$\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}=\operatorname*{\Theta}(1)$ (Ben
Arous et al., 2019b; Jagannath et al., 2019; Zhang and Xia, 2020). Hence, as
the HOOI algorithm is meant to compute a maximum likelihood estimator, with
the assumption $\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}\gg 1$ made in
Theorem 13, it is expected that these iterations allow to perfectly recover
the signal asymptotically. Maximum likelihood estimation is indeed
(theoretically) trivial if
$\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}\to+\infty$. It is more surprising
however that this already happens at the first iteration.

As said previously, the choice of initialization ${\bm{U}}^{(\ell)}_{0}$,
$\ell\in[d]$, does not matter in Theorem 13. In fact, without prior
information, a truncated MLSVD is the best choice as it allows to partially
reconstruct the signal at the $\operatorname*{\Theta}(N^{\frac{d-2}{4}})$
computational threshold. Nevertheless, had one prior information allowing such
a reconstruction in the regime $1\ll{\bm{\mathscr{P}}}\ll N^{\frac{d-2}{4}}$ —
where truncated MLSVD would not be fruitful —, HOOI would still be able to
perfectly reconstruct the signal ${\bm{\mathscr{P}}}$ given this
initialization.

Hence, HOOI initialized with a truncated MLSVD, as it is presented in
Algorithm 1, allows to numerically compute a maximum likelihood estimator
(solution to Problem (1)) but only above the phase transition of truncated
MLSVD. Indeed, its initialization plays a determining role: it must place
${\bm{U}}^{(\ell)}_{0}$, $\ell\in[d]$, in the right basin of attraction,
which, without prior information, is only possible above the
$\operatorname*{\Theta}(N^{\frac{d-2}{4}})$ computational threshold.

Finally, we highlight the fact that these results concern the large $N$ limit.
In practice, it makes no sense to talk about $\operatorname*{\Theta}(1)$ or
$\operatorname*{\Theta}(N^{\frac{d-2}{4}})$ regimes at finite $N$. Figure 2
also depicts the subspace alignments achieved with the maximum likelihood
estimator computed with Algorithm 1 on
${\bm{\mathscr{T}}}=\sqrt{\omega}{\bm{\mathscr{P}}}_{\circ}+\frac{1}{\sqrt{N}}{\bm{\mathscr{N}}}$.
Although
$\lVert{\bm{\mathscr{P}}}_{\circ}\rVert_{\mathrm{F}}=\sqrt{\sigma_{N}}=\operatorname*{\Theta}(N^{\frac{d-2}{4}})$,
HOOI does not achieve perfect recovery of ${\bm{\mathscr{P}}}$ as one might
expect from Theorem 13 (even if several iterations were performed here). In
fact, at finite $N$, perfect recovery is not feasible. But, as $N$ grows, the
dashed line would approach $1$ above the (computational) phase transition
determined by the truncated MLSVD and stay close to $0$ below.

##  5 Concluding Remarks

The analysis presented in this work yield theoretical and practical insight
into the estimation of a low-rank signal from an observation
${\bm{\mathscr{T}}}={\bm{\mathscr{P}}}+\frac{1}{\sqrt{N}}{\bm{\mathscr{N}}}$
following the most general spiked tensor model. While Zhang and Xia (2020)
gave a general overview of the different regimes governing the estimation of
${\bm{\mathscr{P}}}$ with a low-multilinear-rank approximation of
${\bm{\mathscr{T}}}$ — thereby confirming the existence of a computational-to-
statistical gap —, our results shed light on the non-trivial aspects at stake
around the $\operatorname*{\Theta}(N^{\frac{d-2}{4}})$ computational
threshold. This is of particular importance as practical applications lie in
this non-trivial regime where signal and noise have the same magnitude and
must be decoupled. In particular, truncated MLSVD and HOOI are very standard
and efficient algorithms to compute low-multilinear-rank approximations.
Performances of the latter rely strongly on the quality of its initialization,
which is usually performed with a truncated MLSVD in the absence of prior
information. This approach allows the detection of the underlying signal as
early as the computational threshold contrary to other methods such as AMP or
tensor power iteration, which are efficient above a
$\operatorname*{\Theta}(N^{\frac{d-1}{2}})$ algorithmic threshold (Montanari
and Richard, 2014).

Relying on standard tools and methods from the theory of large random
matrices, we have characterized the spectral behavior of the unfoldings of
${\bm{\mathscr{T}}}$ in the large $N$ limit. Specifically, our first main
result shows that, when properly rescaled, the eigenvalues of
${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}$ asymptotically follow a semicircle
distribution. The rescaling exhibits their mean
$\mu^{(\ell)}_{N}=\operatorname*{\Theta}(N^{d-2})$ and a quantity
$\sigma_{N}=\operatorname*{\Theta}(N^{\frac{d-2}{2}})$ governing their spread.
From our denoising perspective, $\sigma_{N}$ indicates the _strength_ of the
noise. Indeed, while the global behavior of the eigenvalues is controlled by
the noise, the addition of a low-rank signal causes the presence of a finite
number of eigenvalues outside the limiting semicircle distribution with
corresponding eigenvectors aligned with the singular subspaces of the sought
signal ${\bm{\mathscr{P}}}$. Yet, the existence of these outlier eigenvalues
hinges on the values of the signal-to-noise ratios
$\rho^{(\ell)}_{q_{\ell}}=s_{q_{\ell}}^{2}({\bm{P}}^{(\ell)})/\sigma_{N}$,
manifesting a BBP phase transition phenomenon. When they exist, the positions
of these isolated eigenvalues and the quality of the corresponding alignments
are completely determined by $\rho^{(\ell)}_{q_{\ell}}$. These results justify
the use of a truncated MLSVD to estimate ${\bm{\mathscr{P}}}$ from the
observation ${\bm{\mathscr{T}}}$ and allow the precise characterization of the
achievable reconstruction performances in the non-trivial regime, i.e., close
to the computational threshold. In particular, we have seen that each singular
value of ${\bm{P}}^{(\ell)}$ determines the position of a phase transition
corresponding to the detectability of the corresponding principal direction.

Although truncated MLSVD does not yield the best low-multilinear-rank
approximation — i.e., a maximum likelihood solution —, it serves as an
excellent initialization for the HOOI algorithm, which converges to such an
estimator if it is initialized _sufficiently close_ to it (Xu, 2018). In fact,
we precise this last assertion by showing that, as long as the initialization
preserves the singular subspaces of ${\bm{\mathscr{P}}}$ in a sense precised
in Theorem 13, HOOI converges to a maximum likelihood solution in a number of
iterations which tends to $1$ as $N\to+\infty$. Hence, when it is initialized
with a truncated MLSVD, it shares the same phase transition, whose position
depends on the singular values $s^{(\ell)}_{q_{\ell}}({\bm{P}}^{(\ell)})$ of
the unfoldings of ${\bm{\mathscr{P}}}$. Yet, given prior information, HOOI can
still reconstruct the maximum likelihood solution below the computational
threshold
$\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}=\operatorname*{\Theta}(N^{\frac{d-2}{4}})$,
where its success depends entirely on the quality of its initialization.

This work gives a comprehensive understanding of the low-multilinear-rank
approximation problem near the computational threshold, which has both
practical and theoretical implications. Besides, from a theoretical
perspective, the behavior of the maximum likelihood estimator is still unclear
near the statistical threshold — that is, in the regime where
$\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}=\operatorname*{\Theta}(1)$.
Several works have studied the rank-$1$ symmetric case (Ben Arous et al.,
2019b; Jagannath et al., 2019) and the approach developed by Seddik et al.
(2022) in their analysis of the rank-$1$ _asymmetric_ case may be an
attractive direction to consider. Relying solely on classical tools from
random matrix theory, they bring the study of the best rank-$1$ tensor
approximation down to that of a structured matrix defined from the
contractions of the data tensor on its dominant singular vectors. Extending
this procedure to our general spiked tensor model (2) presents no conceptual
difficulty, despite being computationally cumbersome due to the multiple
dimensions of the singular subspaces. It is an interesting line of
investigation to refine our understanding of the statistical limits of spiked
tensor models.

Acknowledgments and Disclosure of Funding

This work is supported by the MIAI LargeDATA Chair at Université Grenoble
Alpes.

##  Appendix A Proof of Theorem 5

The resolvent of ${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}$ is defined for all
$z\in{\mathbb{C}}\setminus\operatorname{Sp}{\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}$
as

${\bm{Q}}^{(\ell)}(z)\overset{\text{def}}{=}\left({\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}-z{\bm{I}}_{n_{\ell}}\right)^{-1}.$

We will often drop the dependence in $z$ to simplify notations.

Since ${\bm{Q}}^{(\ell)-1}{\bm{Q}}^{(\ell)}={\bm{I}}_{n_{\ell}}$ and
${\bm{T}}^{(\ell)}={\bm{P}}^{(\ell)}+\frac{1}{\sqrt{N}}{\bm{N}}^{(\ell)}$, we
have,

${\bm{P}}^{(\ell)}{\bm{T}}^{(\ell)\top}{\bm{Q}}^{(\ell)}+\frac{1}{\sqrt{N}}{\bm{N}}^{(\ell)}{\bm{T}}^{(\ell)\top}{\bm{Q}}^{(\ell)}-z{\bm{Q}}^{(\ell)}={\bm{I}}_{n_{\ell}}.$
(9)

###  A.1 Expressions with Stein’s lemma

Using Stein’s lemma (Lemma 4), we find the following expressions.

$\displaystyle\mathbb{E}\left[{\bm{P}}^{(\ell)}{\bm{T}}^{(\ell)\top}{\bm{Q}}^{(\ell)}\right]=\mathbb{E}\left[{\bm{P}}^{(\ell)}{\bm{P}}^{(\ell)\top}{\bm{Q}}^{(\ell)}\right]-\frac{1}{N}\mathbb{E}\left[{\bm{P}}^{(\ell)}{\bm{T}}^{(\ell)\top}{\bm{Q}}^{(\ell)}\operatorname{Tr}{\bm{Q}}^{(\ell)}+{\bm{P}}^{(\ell)}{\bm{T}}^{(\ell)\top}{\bm{Q}}^{(\ell)2}\right],$
(10)
$\displaystyle\mathbb{E}\left[{\bm{N}}^{(\ell)}{\bm{T}}^{(\ell)\top}{\bm{Q}}^{(\ell)}\right]=\frac{\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}}{\sqrt{N}}\mathbb{E}\left[{\bm{Q}}^{(\ell)}\right]-\frac{1}{\sqrt{N}}\mathbb{E}\left[\left(n_{\ell}+1\right){\bm{Q}}^{(\ell)}-z\left({\bm{Q}}^{(\ell)2}+{\bm{Q}}^{(\ell)}\operatorname{Tr}{\bm{Q}}^{(\ell)}\right)\right].$
(11)

#### Derivative of ${\bm{Q}}^{(\ell)}$

Firstly, we need to show that

$\frac{\partial Q^{(\ell)}_{a,b}}{\partial
N^{(\ell)}_{c,d}}=-\frac{1}{\sqrt{N}}\left(Q^{(\ell)}_{a,c}\left[{\bm{T}}^{(\ell)\top}{\bm{Q}}^{(\ell)}\right]_{d,b}+\left[{\bm{Q}}^{(\ell)}{\bm{T}}^{(\ell)}\right]_{a,d}Q^{(\ell)}_{c,b}\right).$
(12)

Indeed, using the fact that
$\partial{\bm{Q}}^{(\ell)}=-{\bm{Q}}^{(\ell)}\partial({\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}){\bm{Q}}^{(\ell)}$,
we have,

$\frac{\partial Q^{(\ell)}_{a,b}}{\partial
N^{(\ell)}_{c,d}}=-\left[{\bm{Q}}^{(\ell)}\frac{\partial{\bm{T}}^{(\ell)}}{\partial
N^{(\ell)}_{c,d}}{\bm{T}}^{(\ell)\top}{\bm{Q}}^{(\ell)}\right]_{a,b}-\left[{\bm{Q}}^{(\ell)}{\bm{T}}^{(\ell)}\frac{\partial{\bm{T}}^{(\ell)\top}}{\partial
N^{(\ell)}_{c,d}}{\bm{Q}}^{(\ell)}\right]_{a,b}$

and, since
${\bm{T}}^{(\ell)}={\bm{P}}^{(\ell)}+\frac{1}{\sqrt{N}}{\bm{N}}^{(\ell)}$,

$\displaystyle\left[{\bm{Q}}^{(\ell)}\frac{\partial{\bm{T}}^{(\ell)}}{\partial
N^{(\ell)}_{c,d}}{\bm{T}}^{(\ell)\top}{\bm{Q}}^{(\ell)}\right]_{a,b}=\frac{1}{\sqrt{N}}Q^{(\ell)}_{a,c}\left[{\bm{T}}^{(\ell)\top}{\bm{Q}}^{(\ell)}\right]_{d,b},$
$\displaystyle\left[{\bm{Q}}^{(\ell)}{\bm{T}}^{(\ell)}\frac{\partial{\bm{T}}^{(\ell)\top}}{\partial
N^{(\ell)}_{c,d}}{\bm{Q}}^{(\ell)}\right]_{a,b}=\frac{1}{\sqrt{N}}\left[{\bm{Q}}^{(\ell)}{\bm{T}}^{(\ell)}\right]_{a,d}Q^{(\ell)}_{c,b}.$

#### Proof of equation (10)

Since
${\bm{T}}^{(\ell)}={\bm{P}}^{(\ell)}+\frac{1}{\sqrt{N}}{\bm{N}}^{(\ell)}$, we
have,

$\mathbb{E}[{\bm{P}}^{(\ell)}{\bm{T}}^{(\ell)\top}{\bm{Q}}^{(\ell)}]=\mathbb{E}[{\bm{P}}^{(\ell)}{\bm{P}}^{(\ell)\top}{\bm{Q}}^{(\ell)}]+\frac{1}{\sqrt{N}}\mathbb{E}[{\bm{P}}^{(\ell)}{\bm{N}}^{(\ell)\top}{\bm{Q}}^{(\ell)}].$

To deal with the rightmost term, we successively use Stein’s lemma (Lemma 4)
and equation (12).

$\displaystyle\mathbb{E}\left[{\bm{P}}^{(\ell)}{\bm{N}}^{(\ell)\top}{\bm{Q}}^{(\ell)}\right]_{i,j}$
$\displaystyle=\sum_{k=1}^{\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}}\sum_{l=1}^{n_{\ell}}\mathbb{E}\left[P^{(\ell)}_{i,k}N^{(\ell)}_{l,k}Q^{(\ell)}_{l,j}\right]$
$\displaystyle=\sum_{k=1}^{\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}}\sum_{l=1}^{n_{\ell}}\mathbb{E}\left[P^{(\ell)}_{i,k}\frac{\partial
Q^{(\ell)}_{l,j}}{\partial N^{(\ell)}_{l,k}}\right]$
$\displaystyle=\begin{multlined}-\frac{1}{\sqrt{N}}\sum_{k=1}^{\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}}\sum_{l=1}^{n_{\ell}}\mathbb{E}\left[P^{(\ell)}_{i,k}Q^{(\ell)}_{l,l}\left[{\bm{T}}^{(\ell)\top}{\bm{Q}}^{(\ell)}\right]_{k,j}\right]\\\
-\frac{1}{\sqrt{N}}\sum_{k=1}^{\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}}\sum_{l=1}^{n_{\ell}}\mathbb{E}\left[P^{(\ell)}_{i,k}\left[{\bm{Q}}^{(\ell)}{\bm{T}}^{(\ell)}\right]_{l,k}Q^{(\ell)}_{l,j}\right]\end{multlined}-\frac{1}{\sqrt{N}}\sum_{k=1}^{\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}}\sum_{l=1}^{n_{\ell}}\mathbb{E}\left[P^{(\ell)}_{i,k}Q^{(\ell)}_{l,l}\left[{\bm{T}}^{(\ell)\top}{\bm{Q}}^{(\ell)}\right]_{k,j}\right]\\\
-\frac{1}{\sqrt{N}}\sum_{k=1}^{\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}}\sum_{l=1}^{n_{\ell}}\mathbb{E}\left[P^{(\ell)}_{i,k}\left[{\bm{Q}}^{(\ell)}{\bm{T}}^{(\ell)}\right]_{l,k}Q^{(\ell)}_{l,j}\right]$
$\displaystyle=-\frac{1}{\sqrt{N}}\mathbb{E}\left[{\bm{P}}^{(\ell)}{\bm{T}}^{(\ell)\top}{\bm{Q}}^{(\ell)}\operatorname{Tr}{\bm{Q}}^{(\ell)}+{\bm{P}}^{(\ell)}{\bm{T}}^{(\ell)\top}{\bm{Q}}^{(\ell)2}\right]_{i,j}.$

#### Proof of equation (11)

We proceed similarly with Stein’s lemma (Lemma 4) and equation (12).

$\displaystyle\mathbb{E}\left[{\bm{N}}^{(\ell)}{\bm{T}}^{(\ell)\top}{\bm{Q}}^{(\ell)}\right]_{i,j}$
$\displaystyle\qquad=\sum_{k=1}^{\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}}\sum_{l=1}^{n_{\ell}}\mathbb{E}\left[N^{(\ell)}_{i,k}T^{(\ell)}_{l,k}Q^{(\ell)}_{l,j}\right]$
$\displaystyle\qquad=\sum_{k=1}^{\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}}\sum_{l=1}^{n_{\ell}}\mathbb{E}\left[\frac{\partial
T^{(\ell)}_{l,k}}{\partial
N^{(\ell)}_{i,k}}Q^{(\ell)}_{l,j}+T^{(\ell)}_{l,k}\frac{\partial
Q^{(\ell)}_{l,j}}{\partial N^{(\ell)}_{i,k}}\right]$
$\displaystyle\qquad=\begin{multlined}\frac{\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}}{\sqrt{N}}\mathbb{E}\left[{\bm{Q}}^{(\ell)}\right]_{i,j}-\frac{1}{\sqrt{N}}\sum_{k=1}^{\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}}\sum_{l=1}^{n_{\ell}}\mathbb{E}\left[T^{(\ell)}_{l,k}Q^{(\ell)}_{l,i}\left[{\bm{T}}^{(\ell)\top}{\bm{Q}}^{(\ell)}\right]_{k,j}\right]\\\
-\frac{1}{\sqrt{N}}\sum_{k=1}^{\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}}\sum_{l=1}^{n_{\ell}}\mathbb{E}\left[T^{(\ell)}_{l,k}\left[{\bm{Q}}^{(\ell)}{\bm{T}}^{(\ell)}\right]_{l,k}Q^{(\ell)}_{i,j}\right]\end{multlined}\frac{\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}}{\sqrt{N}}\mathbb{E}\left[{\bm{Q}}^{(\ell)}\right]_{i,j}-\frac{1}{\sqrt{N}}\sum_{k=1}^{\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}}\sum_{l=1}^{n_{\ell}}\mathbb{E}\left[T^{(\ell)}_{l,k}Q^{(\ell)}_{l,i}\left[{\bm{T}}^{(\ell)\top}{\bm{Q}}^{(\ell)}\right]_{k,j}\right]\\\
-\frac{1}{\sqrt{N}}\sum_{k=1}^{\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}}\sum_{l=1}^{n_{\ell}}\mathbb{E}\left[T^{(\ell)}_{l,k}\left[{\bm{Q}}^{(\ell)}{\bm{T}}^{(\ell)}\right]_{l,k}Q^{(\ell)}_{i,j}\right]$
$\displaystyle\qquad=\frac{\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}}{\sqrt{N}}\mathbb{E}\left[{\bm{Q}}^{(\ell)}\right]_{i,j}-\frac{1}{\sqrt{N}}\mathbb{E}\left[{\bm{Q}}^{(\ell)}{\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}{\bm{Q}}^{(\ell)}+{\bm{Q}}^{(\ell)}\operatorname{Tr}{\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}{\bm{Q}}^{(\ell)}\right]_{i,j}.$

Since
${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}{\bm{Q}}^{(\ell)}-z{\bm{Q}}^{(\ell)}={\bm{I}}_{n_{\ell}}$,
we find the result stated in equation (11).

###  A.2 Deterministic Equivalent

Taking the expectation of equation (9) and injecting equation (11) yields

$\mathbb{E}\left[{\bm{P}}^{(\ell)}{\bm{T}}^{(\ell)\top}{\bm{Q}}^{(\ell)}\right]+\frac{\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}}{N}\mathbb{E}\left[{\bm{Q}}^{(\ell)}\right]\\\
-\frac{1}{N}\mathbb{E}\left[\left(n_{\ell}+1\right){\bm{Q}}^{(\ell)}-z\left({\bm{Q}}^{(\ell)2}+{\bm{Q}}^{(\ell)}\operatorname{Tr}{\bm{Q}}^{(\ell)}\right)\right]-z\mathbb{E}\left[{\bm{Q}}^{(\ell)}\right]={\bm{I}}_{n_{\ell}}.$

We rearrange this expression into the more convenient following form,

$z\frac{n_{\ell}}{N}\mathbb{E}\left[{\bm{Q}}^{(\ell)}\frac{\operatorname{Tr}{\bm{Q}}^{(\ell)}}{n_{\ell}}\right]+\left(z+\frac{n_{\ell}-\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}}{N}\right)\mathbb{E}\left[{\bm{Q}}^{(\ell)}\right]+{\bm{I}}_{n_{\ell}}\\\
=-\frac{1}{N}\mathbb{E}\left[{\bm{Q}}^{(\ell)}-z{\bm{Q}}^{(\ell)2}\right]+\mathbb{E}\left[{\bm{P}}^{(\ell)}{\bm{T}}^{(\ell)\top}{\bm{Q}}^{(\ell)}\right].$

Here, the divergence of the spectrum of
${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}$ becomes problematic: its resolvent
${\bm{Q}}^{(\ell)}$ vanish asymptotically, allowing the presence of the
diverging coefficient
$\frac{1}{N}\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}$ in the previous
equation. To bypass this difficulty, we proceed to a rescaling,

$\displaystyle\tilde{z}\overset{\text{def}}{=}\frac{z-\mu^{(\ell)}_{N}}{\sigma_{N}},\qquad\widetilde{{\bm{Q}}}^{(\ell)}(\tilde{z})\overset{\text{def}}{=}\left(\frac{{\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}-\mu^{(\ell)}_{N}{\bm{I}}_{n_{\ell}}}{\sigma_{N}}-\tilde{z}{\bm{I}}_{n_{\ell}}\right)^{-1}=\sigma_{N}{\bm{Q}}^{(\ell)}(z),$
$\displaystyle\text{with}\qquad\mu^{(\ell)}_{N}=\frac{1}{N}\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}},\qquad\sigma_{N}=\frac{1}{N}\sqrt{\prod_{\ell\in[d]}n_{\ell}}.$

This changes our equation into

$\frac{\mu^{(\ell)}_{N}+\sigma_{N}\tilde{z}}{\sigma_{N}^{2}}\frac{n_{\ell}}{N}\mathbb{E}\left[\widetilde{{\bm{Q}}}^{(\ell)}\frac{\operatorname{Tr}\widetilde{{\bm{Q}}}^{(\ell)}}{n_{\ell}}\right]+\left(\frac{\mu^{(\ell)}_{N}}{\sigma_{N}}+\tilde{z}+\frac{n_{\ell}-\prod_{\ell^{\prime}\neq\ell}n_{\ell^{\prime}}}{\sigma_{N}N}\right)\mathbb{E}\left[\widetilde{{\bm{Q}}}^{(\ell)}\right]+{\bm{I}}_{n_{\ell}}\\\
=-\frac{1}{\sigma_{N}N}\mathbb{E}\left[\widetilde{{\bm{Q}}}^{(\ell)}-\left(\frac{\mu^{(\ell)}_{N}}{\sigma_{N}}+\tilde{z}\right)\widetilde{{\bm{Q}}}^{(\ell)2}\right]+\frac{1}{\sigma_{N}}\mathbb{E}\left[{\bm{P}}^{(\ell)}{\bm{T}}^{(\ell)\top}\widetilde{{\bm{Q}}}^{(\ell)}\right].$

Replacing the vanishing terms by their asymptotic order of magnitude, we have,

$\left(\frac{\mu^{(\ell)}_{N}}{\sigma_{N}^{2}}+\operatorname*{\Theta}(N^{-\frac{d-2}{2}})\right)\frac{n_{\ell}}{N}\mathbb{E}\left[\widetilde{{\bm{Q}}}^{(\ell)}\frac{\operatorname{Tr}\widetilde{{\bm{Q}}}^{(\ell)}}{n_{\ell}}\right]+\left(\tilde{z}+\operatorname*{\Theta}(N^{-\frac{d-2}{2}})\right)\mathbb{E}\left[\widetilde{{\bm{Q}}}^{(\ell)}\right]+{\bm{I}}_{n_{\ell}}\\\
=-\operatorname*{\Theta}(N^{-\frac{d}{2}})\times\mathbb{E}\left[\widetilde{{\bm{Q}}}^{(\ell)}\right]+\operatorname*{\Theta}(N^{-1})\times\mathbb{E}\left[\widetilde{{\bm{Q}}}^{(\ell)2}\right]+\frac{1}{\sigma_{N}}\mathbb{E}\left[{\bm{P}}^{(\ell)}{\bm{T}}^{(\ell)\top}\widetilde{{\bm{Q}}}^{(\ell)}\right].$

We see that there is no longer any diverging coefficient, which justifies our
choice of rescaling. Thus, keeping only the non-vanishing terms, we can now
proceed to the following matrix equivalent formula161616This is justified by
the use of standard concentration arguments such as Poincaré-Nash inequality
(Chen, 1982; Ledoux, 2001) and Borel-Cantelli lemma (Billingsley, 2012). See,
e.g., Couillet and Liao (2022, §2.2.2). This reasoning is applied whenever we
move from an equality with expectations to a matrix equivalent formula.

$\widetilde{{\bm{Q}}}^{(\ell)}\frac{\operatorname{Tr}\widetilde{{\bm{Q}}}^{(\ell)}}{n_{\ell}}+\tilde{z}\widetilde{{\bm{Q}}}^{(\ell)}+{\bm{I}}_{n_{\ell}}\longleftrightarrow\frac{1}{\sigma_{N}}{\bm{P}}^{(\ell)}{\bm{T}}^{(\ell)\top}\widetilde{{\bm{Q}}}^{(\ell)}$
(13)

where we have used the fact that
$\frac{\mu^{(\ell)}_{N}}{\sigma_{N}^{2}}\frac{n_{\ell}}{N}\xrightarrow[N\to+\infty]{}1$.

At this point, we still need to handle the term on the right-hand-side, which
has been purposely ignored until now. From equation (10), we have,

$\mathbb{E}\left[{\bm{P}}^{(\ell)}{\bm{T}}^{(\ell)\top}\left({\bm{Q}}^{(\ell)}+\frac{n_{\ell}}{N}{\bm{Q}}^{(\ell)}\frac{\operatorname{Tr}{\bm{Q}}^{(\ell)}}{n_{\ell}}+\frac{1}{N}{\bm{Q}}^{(\ell)2}\right)\right]=\mathbb{E}\left[{\bm{P}}^{(\ell)}{\bm{P}}^{(\ell)\top}{\bm{Q}}^{(\ell)}\right].$

Applying the same rescaling, this becomes

$\mathbb{E}\left[{\bm{P}}^{(\ell)}{\bm{T}}^{(\ell)\top}\left(\widetilde{{\bm{Q}}}^{(\ell)}+\operatorname*{\Theta}(N^{-\frac{d-2}{2}})\times\widetilde{{\bm{Q}}}^{(\ell)}\frac{\operatorname{Tr}\widetilde{{\bm{Q}}}^{(\ell)}}{n_{\ell}}+\operatorname*{\Theta}(N^{-\frac{d}{2}})\times\widetilde{{\bm{Q}}}^{(\ell)2}\right)\right]\\\
=\mathbb{E}\left[{\bm{P}}^{(\ell)}{\bm{P}}^{(\ell)\top}\widetilde{{\bm{Q}}}^{(\ell)}\right].$

Therefore,

${\bm{P}}^{(\ell)}{\bm{T}}^{(\ell)\top}\widetilde{{\bm{Q}}}^{(\ell)}\longleftrightarrow{\bm{P}}^{(\ell)}{\bm{P}}^{(\ell)\top}\widetilde{{\bm{Q}}}^{(\ell)}$

and equation (13) becomes

$\widetilde{{\bm{Q}}}^{(\ell)}\frac{\operatorname{Tr}\widetilde{{\bm{Q}}}^{(\ell)}}{n_{\ell}}+\tilde{z}\widetilde{{\bm{Q}}}^{(\ell)}+{\bm{I}}_{n_{\ell}}\longleftrightarrow\frac{1}{\sigma_{N}}{\bm{P}}^{(\ell)}{\bm{P}}^{(\ell)\top}\widetilde{{\bm{Q}}}^{(\ell)}.$

Thus, under the assumption that
$\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}}^{2}/\sigma_{N}$ stays bounded
asymptotically, we have a relation for the Stieltjes transform of the limiting
spectral distribution,

$\frac{\operatorname{Tr}\widetilde{{\bm{Q}}}^{(\ell)}}{n_{\ell}}\xrightarrow[N\to+\infty]{\text{a.s.}}\tilde{m}(\tilde{z})\qquad\text{with}\qquad\tilde{m}^{2}(\tilde{z})+\tilde{z}\tilde{m}(\tilde{z})+1=0.$

Indeed,
$\frac{1}{n_{\ell}}\frac{1}{\sigma_{N}}\operatorname{Tr}{\bm{P}}^{(\ell)}{\bm{P}}^{(\ell)\top}\widetilde{{\bm{Q}}}^{(\ell)}\to
0$ almost surely as $N\to+\infty$ because ${\bm{P}}^{(\ell)}$ is a low-rank
matrix. The indeterminacy in the solution to the quadratic equation is raised
by choosing $\tilde{m}(\tilde{z})$ such that
$\tilde{z}\mapsto\tilde{m}(\tilde{z})$ follows the properties of a Stieltjes
transform (in particular $\Im[\tilde{z}]\Im[\tilde{m}(\tilde{z})>0$ for all
$\tilde{z}\in{\mathbb{C}}\setminus{\mathbb{R}}$, see Couillet and Liao (2022,
§2.1.2)).

Moreover, we can define $\bar{{\bm{Q}}}^{(\ell)}$, a _deterministic
equivalent_ of $\widetilde{{\bm{Q}}}^{(\ell)}$, such that

$\tilde{m}(\tilde{z})\bar{{\bm{Q}}}^{(\ell)}+\tilde{z}\bar{{\bm{Q}}}^{(\ell)}+{\bm{I}}_{n_{\ell}}=\frac{1}{\sigma_{N}}{\bm{P}}^{(\ell)}{\bm{P}}^{(\ell)\top}\bar{{\bm{Q}}}^{(\ell)}\\\
\iff\bar{{\bm{Q}}}^{(\ell)}=\left(\frac{1}{\sigma_{N}}{\bm{P}}^{(\ell)}{\bm{P}}^{(\ell)\top}-\left(\tilde{m}(\tilde{z})+\tilde{z}\right){\bm{I}}_{n_{\ell}}\right)^{-1}.$

Since $\tilde{m}^{2}(\tilde{z})+\tilde{z}\tilde{m}(\tilde{z})+1=0$, we have
$\tilde{m}(\tilde{z})+\tilde{z}=-\frac{1}{\tilde{m}(\tilde{z})}$ and we find
the definition of $\bar{{\bm{Q}}}^{(\ell)}$ given in Theorem 5.

##  Appendix B Proof of Theorem 8

###  B.1 Isolated Eigenvalues

The asymptotic positions of isolated eigenvalues are singular points of the
deterministic equivalent
$\bar{{\bm{Q}}}^{(\ell)}:\tilde{z}\mapsto\left(\frac{1}{\sigma_{N}}{\bm{P}}^{(\ell)}{\bm{P}}^{(\ell)\top}+\frac{1}{\tilde{m}(\tilde{z})}{\bm{I}}_{n_{\ell}}\right)^{-1}$.
Thus, we seek $\tilde{\xi}^{(\ell)}_{q_{\ell}}\in{\mathbb{R}}\setminus[-2,+2]$
such that

$\frac{s_{q_{\ell}}^{2}({\bm{P}}^{(\ell)})}{\sigma_{N}}+\frac{1}{\tilde{m}(\tilde{\xi}^{(\ell)}_{q_{\ell}})}=0,\qquad
q_{\ell}\in[r_{\ell}].$

Injecting the expression
$\tilde{m}(\tilde{\xi}^{(\ell)}_{q_{\ell}})=-\frac{\sigma_{N}}{s_{q_{\ell}}^{2}({\bm{P}}^{(\ell)})}$
into the equation
$\tilde{m}^{2}(\tilde{\xi}^{(\ell)}_{q_{\ell}})+\tilde{\xi}^{(\ell)}_{q_{\ell}}\tilde{m}(\tilde{\xi}^{(\ell)}_{q_{\ell}})+1=0$
yields

$\frac{\sigma_{N}^{2}}{s_{q_{\ell}}^{4}({\bm{P}}^{(\ell)})}-\tilde{\xi}^{(\ell)}_{q_{\ell}}\frac{\sigma_{N}}{s_{q_{\ell}}^{2}({\bm{P}}^{(\ell)})}+1=0\iff\tilde{\xi}^{(\ell)}_{q_{\ell}}=\frac{s_{q_{\ell}}^{2}({\bm{P}}^{(\ell)})}{\sigma_{N}}+\frac{\sigma_{N}}{s_{q_{\ell}}^{2}({\bm{P}}^{(\ell)})}.$

As $\tilde{\xi}^{(\ell)}_{q_{\ell}}>0$ by definition, to be a singular point
of $\bar{{\bm{Q}}}^{(\ell)}(\tilde{z})$, it must be strictly greater than $2$
(the right edge of the semicircle). This is true only if
$\rho^{(\ell)}_{q_{\ell}}\overset{\text{def}}{=}\frac{s_{q_{\ell}}^{2}({\bm{P}}^{(\ell)})}{\sigma_{N}}>1$.

###  B.2 Eigenvector Alignments

Let $\hat{{\bm{u}}}^{(\ell)}_{i_{\ell}}$, $i_{\ell}\in[n_{\ell}]$, denote the
$i_{\ell}$-th left singular vector of ${\bm{T}}^{(\ell)}$ (sorted in non-
increasing order of its corresponding singular value). From the spectral
decomposition
${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}=\sum_{i_{\ell}=1}^{n_{\ell}}s_{i_{\ell}}^{2}({\bm{T}}^{(\ell)})\hat{{\bm{u}}}^{(\ell)}_{i_{\ell}}\hat{{\bm{u}}}^{(\ell)\top}_{i_{\ell}}$,
we have,

$\widetilde{{\bm{Q}}}^{(\ell)}(\tilde{z})=\sum_{i_{\ell}=1}^{n_{\ell}}\frac{\hat{{\bm{u}}}^{(\ell)}_{q_{\ell}}\hat{{\bm{u}}}^{(\ell)\top}_{q_{\ell}}}{s_{i_{\ell}}^{2}({\bm{T}}^{(\ell)})-\tilde{z}}.$

If $\rho^{(\ell)}_{q_{\ell}}>1$, $q_{\ell}\in[r_{\ell}]$, then
$s_{q_{\ell}}^{2}({\bm{T}}^{(\ell)})$ is an isolated eigenvalue in the
spectrum of ${\bm{T}}^{(\ell)}{\bm{T}}^{(\ell)\top}$ and
$\frac{1}{\sigma_{N}}\left[s_{q_{\ell}}^{2}({\bm{T}}^{(\ell)})-\mu^{(\ell)}_{N}\right]\xrightarrow[N\to+\infty]{\text{a.s.}}\tilde{\xi}^{(\ell)}_{q_{\ell}}$.
Hence, for any positively-oriented complex contour
$\gamma^{(\ell)}_{q_{\ell}}$ circling around
$\tilde{\xi}^{(\ell)}_{q_{\ell}}$, leaving all the other
$\tilde{\xi}^{(\ell)}_{q_{\ell}^{\prime}}$, $q_{\ell}^{\prime}\neq q_{\ell}$,
outside and not crossing $[-2,+2]$, Cauchy’s integral formula yields, for $N$
large enough and any $\in{\mathbb{R}}^{n_{\ell}}$,

$\left\lvert{}^{\top}\hat{{\bm{u}}}^{(\ell)}_{q_{\ell}}\right\rvert^{2}=-\frac{1}{2{\mathrm{i}}\pi}\oint_{\gamma^{(\ell)}_{q_{\ell}}}^{\top}\widetilde{{\bm{Q}}}^{(\ell)}(\tilde{z})\leavevmode\nobreak\
{\mathrm{d}}\tilde{z}\xrightarrow[N\to+\infty]{\text{a.s.}}-\frac{1}{2{\mathrm{i}}\pi}\oint_{\gamma^{(\ell)}_{q_{\ell}}}^{\top}\bar{{\bm{Q}}}^{(\ell)}(\tilde{z})\leavevmode\nobreak\
{\mathrm{d}}\tilde{z}.$

Using residue calculus, we can compute,

$-\frac{1}{2{\mathrm{i}}\pi}\oint_{\gamma^{(\ell)}_{q_{\ell}}}^{\top}\bar{{\bm{Q}}}^{(\ell)}(\tilde{z})\leavevmode\nobreak\
{\mathrm{d}}\tilde{z}=-\lim_{\tilde{z}\to\tilde{\xi}^{(\ell)}_{q_{\ell}}}\left(\tilde{z}-\tilde{\xi}^{(\ell)}_{q_{\ell}}\right)^{\top}\left(\frac{1}{\sigma_{N}}{\bm{P}}^{(\ell)}{\bm{P}}^{(\ell)\top}+\frac{1}{\tilde{m}(\tilde{z})}{\bm{I}}_{n_{\ell}}\right)^{-1}$

Note that
${\bm{P}}^{(\ell)}{\bm{P}}^{(\ell)\top}={\bm{X}}^{(\ell)}{\bm{H}}^{(\ell)}{\bm{H}}^{(\ell)\top}{\bm{X}}^{(\ell)\top}$
and there exist an $r_{\ell}\times r_{\ell}$ orthogonal matrix
${\bm{O}}^{(\ell)}$ such that
${\bm{H}}^{(\ell)}{\bm{H}}^{(\ell)\top}={\bm{O}}^{(\ell)}{\bm{L}}^{(\ell)}{\bm{O}}^{(\ell)\top}$
with
${\bm{L}}^{(\ell)}=\operatorname{Diag}(s_{1}^{2}({\bm{P}}^{(\ell)}),\ldots,s_{r_{\ell}}^{2}({\bm{P}}^{(\ell)}))$.
Hence,

${}^{\top}\left(\frac{1}{\sigma_{N}}{\bm{P}}^{(\ell)}{\bm{P}}^{(\ell)\top}+\frac{1}{\tilde{m}(\tilde{z})}{\bm{I}}_{n_{\ell}}\right)^{-1}{}^{\top}{\bm{X}}^{(\ell)}{\bm{O}}^{(\ell)}\left(\frac{1}{\sigma_{N}}{\bm{L}}^{(\ell)}+\frac{1}{\tilde{m}(\tilde{z})}{\bm{I}}_{n_{\ell}}\right)^{-1}{\bm{O}}^{(\ell)\top}{\bm{X}}^{(\ell)\top}$

Let us therefore compute the following quantity,

$-\lim_{\tilde{z}\to\tilde{\xi}^{(\ell)}_{q_{\ell}}}\left(\tilde{z}-\tilde{\xi}^{(\ell)}_{q_{\ell}}\right)\left(\frac{s_{q_{\ell}^{\prime}}^{2}({\bm{P}}^{(\ell)})}{\sigma_{N}}+\frac{1}{\tilde{m}(\tilde{z})}\right)^{-1}=\left\\{\begin{array}[]{ll}0&\text{if}\leavevmode\nobreak\
q_{\ell}^{\prime}\neq q_{\ell}\\\
\zeta^{(\ell)}_{q_{\ell}}&\text{if}\leavevmode\nobreak\
q_{\ell}^{\prime}=q_{\ell}\end{array}\right.,\qquad
q_{\ell}^{\prime}\in[r_{\ell}],$

where we have used the fact that
$\tilde{m}(\tilde{\xi}^{(\ell)}_{q_{\ell}^{\prime}})=-\frac{\sigma_{N}}{s_{q_{\ell}^{\prime}}^{2}({\bm{P}}^{(\ell)})}$.
The case $q_{\ell}^{\prime}=q_{\ell}$ can be handled with L’Hôpital’s rule,

$\displaystyle\zeta^{(\ell)}_{q_{\ell}}$
$\displaystyle=-\left(\frac{{\mathrm{d}}}{{\mathrm{d}}\tilde{z}}\left[\frac{s_{q_{\ell}}^{2}({\bm{P}}^{(\ell)})}{\sigma_{N}}+\frac{1}{\tilde{m}(\tilde{z})}\right]_{\tilde{z}=\tilde{\xi}^{(\ell)}_{q_{\ell}}}\right)^{-1}$
$\displaystyle=\frac{\tilde{m}^{2}(\tilde{\xi}^{(\ell)}_{q_{\ell}})}{\tilde{m}^{\prime}(\tilde{\xi}^{(\ell)}_{q_{\ell}})}$
$\displaystyle=\frac{\sigma_{N}^{2}}{s_{q_{\ell}}^{4}({\bm{P}}^{(\ell)})\tilde{m}^{\prime}(\tilde{\xi}^{(\ell)}_{q_{\ell}})}.$

In order to compute $\tilde{m}^{\prime}(\tilde{\xi}^{(\ell)}_{q_{\ell}})$, let
us differentiate the relation
$\tilde{m}^{2}(\tilde{z})+\tilde{z}\tilde{m}(\tilde{z})+1=0$,

$\displaystyle
2\tilde{m}^{\prime}(\tilde{z})\tilde{m}(\tilde{z})+\tilde{m}(\tilde{z})+\tilde{z}\tilde{m}^{\prime}(\tilde{z})=0,$
$\displaystyle\tilde{m}^{\prime}(\tilde{z})=-\frac{\tilde{m}(\tilde{z})}{2\tilde{m}(\tilde{z})+\tilde{z}}.$

Hence,

$\displaystyle\tilde{m}^{\prime}(\tilde{\xi}^{(\ell)}_{q_{\ell}})$
$\displaystyle=-\frac{-\frac{\sigma_{N}}{s_{q_{\ell}}^{2}({\bm{P}}^{(\ell)})}}{-2\frac{\sigma_{N}}{s_{q_{\ell}}^{2}({\bm{P}}^{(\ell)})}+\frac{s_{q_{\ell}}^{2}({\bm{P}}^{(\ell)})}{\sigma_{N}}+\frac{\sigma_{N}}{s_{q_{\ell}}^{2}({\bm{P}}^{(\ell)})}}$
$\displaystyle=\frac{1}{\frac{s_{q_{\ell}}^{4}({\bm{P}}^{(\ell)})}{\sigma_{N}^{2}}-1}.$

Back to our previous expression of $\zeta^{(\ell)}_{q_{\ell}}$, we now have,

$\zeta^{(\ell)}_{q_{\ell}}=1-\frac{\sigma_{N}^{2}}{s_{q_{\ell}}^{4}({\bm{P}}^{(\ell)})}.$

Therefore, for all $\in{\mathbb{R}}^{n_{\ell}}$,

$-\frac{1}{2{\mathrm{i}}\pi}\oint_{\gamma^{(\ell)}_{q_{\ell}}}^{\top}\bar{{\bm{Q}}}^{(\ell)}(\tilde{z})\leavevmode\nobreak\
{\mathrm{d}}\tilde{z}=^{\top}{\bm{X}}^{(\ell)}{\bm{O}}^{(\ell)}{\bm{Z}}^{(\ell)}_{q_{\ell}}{\bm{O}}^{(\ell)\top}{\bm{X}}^{(\ell)\top}$

where ${\bm{Z}}^{(\ell)}_{q_{\ell}}$ is an $r_{\ell}\times r_{\ell}$ matrix
with all its entries equal to $0$ except
$[{\bm{Z}}^{(\ell)}_{q_{\ell}}]_{q_{\ell},q_{\ell}}=\zeta^{(\ell)}_{q_{\ell}}$.
Thus, summing the alignments of $\hat{{\bm{u}}}^{(\ell)}_{q_{\ell}}$ with each
column of ${\bm{X}}^{(\ell)}$ yields

$\left\lVert{\bm{X}}^{(\ell)\top}\hat{{\bm{u}}}^{(\ell)}_{q_{\ell}}\right\rVert^{2}\xrightarrow[N\to+\infty]{\text{a.s.}}\zeta^{(\ell)}_{q_{\ell}}\sum_{q_{\ell}^{\prime}=1}^{r_{\ell}}O^{(\ell)2}_{q_{\ell}^{\prime},q_{\ell}}=\zeta^{(\ell)}_{q_{\ell}}.$

##  Appendix C Proof of Lemma 12

Our proof of Lemma 12 uses the notion of $\varepsilon$-covering. An
$\varepsilon$-covering of a _compact_ set ${\mathcal{S}}$ for the norm
$\lVert\cdot\rVert$ is a _finite_ set ${\mathcal{C}}\subset{\mathcal{S}}$ such
that for all $x\in{\mathcal{S}}$, there exists $\bar{x}\in{\mathcal{C}}$ such
that $\lVert x-\bar{x}\rVert\leqslant\varepsilon$. We also define the covering
number $N(\varepsilon,{\mathcal{S}},\lVert\cdot\rVert)$ as the smallest
possible number of elements in ${\mathcal{C}}$.

Moreover, we recall the definition of the Gamma function
$\Gamma(s)=\int_{0}^{+\infty}t^{s-1}e^{-t}{\mathrm{d}}t$ and the (upper)
incomplete Gamma function
$\Gamma(s,x)=\int_{x}^{+\infty}t^{s-1}e^{-t}{\mathrm{d}}t$ for $s>0$ and
$x\geqslant 0$.

For our proof, we need to introduce a few preliminary results which are stated
and proven below (except Lemma 17 for which a reference is given).

###### Lemma 15

For $\ell\in[d]$ and $\varepsilon>0$, let
${\bm{\Delta}}^{(\ell)}\in{\mathbb{R}}^{n_{\ell}\times r_{\ell}}$ be such that
$\lVert{\bm{\Delta}}^{(\ell)}\rVert\leqslant\varepsilon$ and
${\bm{V}}^{(\ell)}\in V_{r_{\ell}}({\mathbb{R}}^{n_{\ell}})$ be the matrix of
its left singular vectors. For all
${\bm{A}}^{(\ell^{\prime})}\in{\mathbb{R}}^{n_{\ell^{\prime}}\times
r_{\ell^{\prime}}}$, $\ell^{\prime}\neq\ell$,

$\left\lVert{\bm{\mathscr{N}}}({\bm{A}}^{(1)},\ldots,{\bm{\Delta}}^{(\ell)},\ldots,{\bm{A}}^{(d)})\right\rVert_{\mathrm{F}}\leqslant\varepsilon\left\lVert{\bm{\mathscr{N}}}({\bm{A}}^{(1)},\ldots,{\bm{V}}^{(\ell)},\ldots,{\bm{A}}^{(d)})\right\rVert_{\mathrm{F}}.$

Proof  Let ${\bm{V}}^{(\ell)}{\bm{\Sigma}}^{(\ell)}{\bm{W}}^{(\ell)\top}$ be
the singular value decomposition of ${\bm{\Delta}}^{(\ell)}$. We have,

$\left\lVert{\bm{\mathscr{N}}}({\bm{A}}^{(1)},\ldots,{\bm{\Delta}}^{(\ell)},\ldots,{\bm{A}}^{(d)})\right\rVert_{\mathrm{F}}=\left\lVert{\bm{\Delta}}^{(\ell)\top}{\bm{N}}^{(\ell)}\operatorname*{\scalerel*{\boxtimes}{\bigotimes}}_{\ell^{\prime}\neq\ell}{\bm{A}}^{(\ell^{\prime})}\right\rVert_{\mathrm{F}}\\\
=\left\lVert{\bm{W}}^{(\ell)\top}{\bm{\Sigma}}^{(\ell)}{\bm{V}}^{(\ell)\top}{\bm{N}}^{(\ell)}\operatorname*{\scalerel*{\boxtimes}{\bigotimes}}_{\ell^{\prime}\neq\ell}{\bm{A}}^{(\ell^{\prime})}\right\rVert_{\mathrm{F}}\leqslant\underbrace{\left\lVert{\bm{W}}^{(\ell)\top}{\bm{\Sigma}}^{(\ell)}\right\rVert}_{=\varepsilon}\underbrace{\left\lVert{\bm{V}}^{(\ell)\top}{\bm{N}}^{(\ell)}\operatorname*{\scalerel*{\boxtimes}{\bigotimes}}_{\ell^{\prime}\neq\ell}{\bm{A}}^{(\ell^{\prime})}\right\rVert_{\mathrm{F}}}_{=\left\lVert{\bm{\mathscr{N}}}({\bm{A}}^{(1)},\ldots,{\bm{V}}^{(\ell)},\ldots,{\bm{A}}^{(d)})\right\rVert_{\mathrm{F}}}$

using the fact that
$\lVert{\bm{A}}{\bm{B}}\rVert_{\mathrm{F}}\leqslant\lVert{\bm{A}}\rVert\lVert{\bm{B}}\rVert_{\mathrm{F}}$.

###### Lemma 16

Given ${\bm{A}}{(\ell)}\in V_{r_{\ell}}({\mathbb{R}}^{n_{\ell}})$,
$\ell\in[d]$,

$\left\lVert{\bm{\mathscr{N}}}({\bm{A}}^{(1)},\ldots,{\bm{A}}^{(d)})\right\rVert_{\mathrm{F}}^{2}\sim\chi^{2}\left(\prod_{\ell\in[d]}r_{\ell}\right).$

Proof  Firstly, observe that, for all
$(q_{1},\ldots,q_{d})\in\bigtimes_{\ell\in[d]}[r_{\ell}]$,

$[{\bm{\mathscr{N}}}({\bm{A}}^{(1)},\ldots,{\bm{A}}^{(d)})]_{q_{1},\ldots,q_{d}}=\sum_{i_{1},\ldots,i_{d}=1}^{n_{1},\ldots,n_{d}}{\mathscr{N}}_{i_{1},\ldots,i_{d}}A^{(1)}_{i_{1},q_{1}}\ldots
A^{(d)}_{i_{d},q_{d}}\sim{\mathcal{N}}(0,1)$

since $\sum_{i_{\ell}=1}^{n_{\ell}}A^{(\ell)2}_{i_{\ell},q_{\ell}}=1$ for all
$\ell\in[d]$. Then, we show that all the entries of
${\bm{\mathscr{N}}}({\bm{A}}^{(1)},\ldots,{\bm{A}}^{(d)})$ are independent
because their covariance is identity,

$\displaystyle\mathbb{E}\left[[{\bm{\mathscr{N}}}({\bm{A}}^{(1)},\ldots,{\bm{A}}^{(d)})]_{q_{1},\ldots,q_{d}}[{\bm{\mathscr{N}}}({\bm{A}}^{(1)},\ldots,{\bm{A}}^{(d)})]_{q_{1}^{\prime},\ldots,q_{d}^{\prime}}\right]$
$\displaystyle=\sum_{i_{1},\ldots,i_{d}=1}^{n_{1},\ldots,n_{d}}\sum_{i_{1}^{\prime},\ldots,i_{d}^{\prime}=1}^{n_{1},\ldots,n_{d}}\mathbb{E}\left[{\mathscr{N}}_{i_{1},\ldots,i_{d}}{\mathscr{N}}_{i_{1}^{\prime},\ldots,i_{d}^{\prime}}\right]A^{(1)}_{i_{1},q_{1}}A^{(1)}_{i_{1}^{\prime},q_{1}^{\prime}}\ldots
A^{(d)}_{i_{d},q_{d}}A^{(d)}_{i_{d}^{\prime},q_{d}^{\prime}}$
$\displaystyle=\sum_{i_{1},\ldots,i_{d}=1}^{n_{1},\ldots,n_{d}}A^{(1)}_{i_{1},q_{1}}A^{(1)}_{i_{1},q_{1}^{\prime}}\ldots
A^{(d)}_{i_{d},q_{d}}A^{(d)}_{i_{d},q_{d}^{\prime}}$
$\displaystyle=\left\\{\begin{array}[]{ll}1&\text{if}\leavevmode\nobreak\
(q_{1},\ldots,q_{d})=(q_{1}^{\prime},\ldots,q_{d}^{\prime})\\\
0&\text{otherwise}\end{array}\right..$

Hence, the result follows from the fact that
$\lVert{\bm{\mathscr{N}}}({\bm{A}}^{(1)},\ldots,{\bm{A}}^{(d)})\rVert_{\mathrm{F}}^{2}$
is the sum of $\prod_{\ell\in[d]}r_{\ell}$ squared independent
${\mathcal{N}}(0,1)$ variables.

######  Lemma 17 (Hinrichs et al., 2016, Lemma 4.1)

For $0<\varepsilon<1$, we have the following upper bound on the
$\varepsilon$-covering number of the Stiefel manifold
$V_{r}({\mathbb{R}}^{n})$ for the spectral norm $\lVert\cdot\rVert$,

$N(\varepsilon,V_{r}({\mathbb{R}}^{n}),\left\lVert\cdot\right\rVert)\leqslant\left[\frac{C}{\varepsilon}\right]^{r\left(n-\frac{r+1}{2}\right)}.$

where $C>0$ is a universal constant.

###### Lemma 18

$\Gamma(s,x)\leqslant\max(1,e^{s-1})\Gamma(s)e^{-x/2}$ for all $x\geqslant 0$
and $s>0$.

Proof  Given $s>0$, consider the function
$f:x\in[0,+\infty[\mapsto\frac{\Gamma(s,x)}{Ce^{-x/2}}$ with $C>0$. Our goal
is to show that $0<f\leqslant 1$ when $C$ is well chosen. $f$ is continuously
differentiable on $[0,+\infty[$ and

$f^{\prime}(x)=\frac{1}{Ce^{-x/2}}\left(-x^{s-1}e^{-x}+\frac{1}{2}\Gamma(s,x)\right)\geqslant
0\iff\Gamma(s,x)-2x^{s-1}e^{-x}\geqslant 0.$

Consider the function $g:x\in[0,+\infty[\mapsto\Gamma(s,x)-2x^{s-1}e^{-x}$.
$g$ is also continuously differentiable on $[0,+\infty[$ and

$g^{\prime}(x)=x^{s-1}e^{-x}-2\left(s-1\right)x^{s-2}e^{-x}\geqslant 0\iff
x\geqslant 2\left(s-1\right).$

We distinguish two cases.

  1. 1.

If $0<s\leqslant 1$, then $g$ increases monotonically on $[0,+\infty[$. Since
$\lim_{x\to+\infty}g(x)=0$, we necessarily have $g(x)\leqslant 0$ for all
$x\in[0,+\infty[$. Hence, $f(x)\leqslant f(0)=\frac{\Gamma(s)}{C}$ and we can
choose $C=\Gamma(s)$.

  2. 2.

If $s>1$, our conclusion stems from the following table.

$x$$g^{\prime}(x)$$g(x)$$0$$2\left(s-1\right)$$+\infty$$\Gamma(s)$$g(2\left(s-1\right))$

Since $g$ is strictly increasing on $\left[2\left(s-1\right),+\infty\right[$
and $\lim_{x\to+\infty}g(x)=0$, we necessarily have $g(2\left(s-1\right))<0$.
Hence, since $\Gamma(s)>0$, the equation $g(x)=0$ has a unique solution on
$[0,+\infty[$ and it lies between $0$ and $2\left(s-1\right)$. Let
$x_{0}(s)\in\left[0,2\left(s-1\right)\right]$ be this unique solution. Then,
$\sup_{[0,+\infty[}f=f(x_{0}(s))=\frac{\Gamma(s,x_{0}(s))}{Ce^{-x_{0}(s)/2}}$
and we can choose $C=\Gamma(s,x_{0}(s))e^{x_{0}(s)/2}$. The final result
follows from $\Gamma(s,x_{0}(s))\leqslant\Gamma(s)$ and
$e^{x_{0}(s)/2}\leqslant e^{s-1}$.

We are now ready to prove Lemma 12.

Let $\varepsilon>0$ and ${\mathcal{C}}_{\ell}$ be an $\varepsilon$-covering of
$V_{r_{\ell}}({\mathbb{R}}^{n_{\ell}})$ for the spectral norm
$\lVert\cdot\rVert$, $\ell\in[d]$. Since
$\bigtimes_{\ell\in[d]}V_{r_{\ell}}({\mathbb{R}}^{n_{\ell}})$ is compact, it
contains an element $({\bm{A}}^{(1)}_{\star},\ldots,{\bm{A}}^{(d)}_{\star})$
such that

$\sup_{{\bm{A}}^{(\ell)}\in
V_{r_{\ell}}({\mathbb{R}}^{n_{\ell}}),\leavevmode\nobreak\
\ell\in[d]}\left\lVert{\bm{\mathscr{N}}}({\bm{A}}^{(1)},\ldots,{\bm{A}}^{(d)})\right\rVert_{\mathrm{F}}=\left\lVert{\bm{\mathscr{N}}}({\bm{A}}^{(1)}_{\star},\ldots,{\bm{A}}^{(d)}_{\star})\right\rVert_{\mathrm{F}}.$

Let $\bar{{\bm{A}}}^{(\ell)}\in{\mathcal{C}}_{\ell}$ be such that
${\bm{A}}^{(\ell)}_{\star}=\bar{{\bm{A}}}^{(\ell)}+{\bm{\Delta}}^{(\ell)}$
with $\lVert{\bm{\Delta}}^{(\ell)}\rVert\leqslant\varepsilon$. Then, using the
triangle inequality, Lemma 15 and the optimality of
$({\bm{A}}^{(1)}_{\star},\ldots,{\bm{A}}^{(d)}_{\star})$, we have,

$\left\lVert{\bm{\mathscr{N}}}(\bar{{\bm{A}}}^{(1)}+{\bm{\Delta}}^{(1)},\ldots,\bar{{\bm{A}}}^{(d)}+{\bm{\Delta}}^{(d)})\right\rVert_{\mathrm{F}}\leqslant\left\lVert{\bm{\mathscr{N}}}(\bar{{\bm{A}}}^{(1)},\ldots,\bar{{\bm{A}}}^{(d)})\right\rVert_{\mathrm{F}}+S\left\lVert{\bm{\mathscr{N}}}({\bm{A}}^{(1)}_{\star},\ldots,{\bm{A}}^{(d)}_{\star})\right\rVert_{\mathrm{F}}$

with
$S\overset{\text{def}}{=}\sum_{k=1}^{d}\binom{d}{k}\varepsilon^{k}\leqslant\sum_{k=1}^{d}\varepsilon^{k}\frac{d^{k}}{k!}\leqslant
e^{\varepsilon d}-1$. Hence, choosing
$\varepsilon=\frac{1}{d}\log\frac{3}{2}$, we get,

$\left\lVert{\bm{\mathscr{N}}}({\bm{A}}^{(1)}_{\star},\ldots,{\bm{A}}^{(d)}_{\star})\right\rVert_{\mathrm{F}}\leqslant
2\left\lVert{\bm{\mathscr{N}}}(\bar{{\bm{A}}}^{(1)},\ldots,\bar{{\bm{A}}}^{(d)})\right\rVert_{\mathrm{F}}$

and, from the union bound, for any $t\geqslant 0$,

$\displaystyle{\mathbb{P}}\left(\left\lVert{\bm{\mathscr{N}}}({\bm{A}}^{(1)}_{\star},\ldots,{\bm{A}}^{(d)}_{\star})\right\rVert_{\mathrm{F}}\geqslant
t\right)$
$\displaystyle\leqslant{\mathbb{P}}\left(\bigcup_{{\bm{A}}^{(\ell)}\in{\mathcal{C}}_{\ell},\leavevmode\nobreak\
\ell\in[d]}\left\\{\left\lVert{\bm{\mathscr{N}}}({\bm{A}}^{(1)},\ldots,{\bm{A}}^{(d)})\right\rVert_{\mathrm{F}}\geqslant\frac{t}{2}\right\\}\right)$
$\displaystyle\leqslant\sum_{{\bm{A}}^{(\ell)}\in{\mathcal{C}}_{\ell},\leavevmode\nobreak\
\ell\in[d]}{\mathbb{P}}\left(\left\lVert{\bm{\mathscr{N}}}({\bm{A}}^{(1)},\ldots,{\bm{A}}^{(d)})\right\rVert_{\mathrm{F}}\geqslant\frac{t}{2}\right).$

Thus, combining Lemma 16 and 17, we have,

${\mathbb{P}}\left(\left\lVert{\bm{\mathscr{N}}}({\bm{A}}^{(1)}_{\star},\ldots,{\bm{A}}^{(d)}_{\star})\right\rVert_{\mathrm{F}}\geqslant
t\right)\leqslant\left[\frac{Cd}{\log\frac{3}{2}}\right]^{\sum_{\ell=1}^{d}r_{\ell}\left(n_{\ell}-\frac{r_{\ell}+1}{2}\right)}{\mathbb{P}}\left(X\geqslant\frac{t^{2}}{4}\right)$

where $X$ is a random variable following a
$\chi^{2}(\prod_{\ell\in[d]}r_{\ell})$ distribution. Eventually, the
probability on the right-hand side can be bounded using Lemma 18,

${\mathbb{P}}\left(X\geqslant\frac{t^{2}}{4}\right)=\frac{\Gamma\left(\frac{1}{2}\prod_{\ell\in[d]}r_{\ell},\frac{t^{2}}{8}\right)}{\Gamma\left(\frac{1}{2}\prod_{\ell\in[d]}r_{\ell}\right)}\leqslant\max(1,e^{\frac{1}{2}\prod_{\ell\in[d]}r_{\ell}-1})e^{-t^{2}/16}.$

We get the result stated in Lemma 12 with

$t^{2}=16\left[\left(\sum_{\ell=1}^{d}r_{\ell}\left(n_{\ell}-\frac{r_{\ell}+1}{2}\right)\right)\log\frac{Cd}{\log\frac{3}{2}}+\log\left(\frac{1}{\delta}\max\left(1,e^{\frac{1}{2}\prod_{\ell=1}^{d}r_{\ell}-1}\right)\right)\right].$

##  Appendix D Proof of Theorem 13

Recall that
${\bm{\mathscr{P}}}=\llbracket{\bm{\mathscr{H}}};{\bm{X}}^{(1)},\ldots,{\bm{X}}^{(d)}\rrbracket$.
We use the following lemma whenever we state that
$\lVert{\bm{\mathscr{P}}}({\bm{U}}^{(1)},\ldots,{\bm{U}}^{(d)})\rVert_{\mathrm{F}}=\operatorname*{\mathcal{O}}_{N\to+\infty}(\lVert{\bm{\mathscr{P}}}\rVert_{\mathrm{F}})$.

###### Lemma 19

For all
$({\bm{A}}^{(1)},\ldots,{\bm{A}}^{(d)})\in\bigtimes_{\ell\in[d]}V_{r_{\ell}}({\mathbb{R}}^{n_{\ell}})$,